[R-sig-Geo] gIntersection fails

[Philippe Liege]

Hello,

 

I'm looking to clip the Europe map after reading this post:

http://stackoverflow.com/questions/15881455/how-to-clip-worldmap-with-polygo
n-in-r

 

The procedure in running fine when using the shape file NUTS_RG_60M_2010.shp

 

library(rgdal)

s <- readOGR('NUTS_2010_60M_SH/data', layer="NUTS_RG_60M_2010")

s1 <- s[s$STAT_LEVL_==2 & substring(s$NUTS_ID,1,2) %in%
c("FR","BE","ES","IT","PT","DE","CH","UK","IE","DK","AT"),]

# plot(s1)

clip.extent <- as(extent(-10, 19, 36, 51), "SpatialPolygons")

# clip.extent <- as(extent(-10, 40, 60, 72), "SpatialPolygons")

proj4string(clip.extent) <- CRS(proj4string(s))

gI <- gIntersects( s1 , clip.extent , byid = TRUE )

out <- lapply( which(gI) , function(x){ gIntersection( s1[x,] , clip.extent)
} )

 

# But let's look at what is returned

table( sapply( out , class ) )

 

# We want to keep only objects of class SpatialPolygons                 

keep <- sapply(out, class)

out <- out[keep == "SpatialPolygons"]

 

# Coerce list back to SpatialPolygons object

EuropeW <- SpatialPolygons( lapply( 1:length( out ) , function(i) { Pol <-
slot(out[[i]], "polygons")[[1]];   slot(Pol, "ID") <- as.character(i)

   Pol

}))

plot( EuropeW )

 

 

 

However I'm getting an error when using the shape files NUTS_RG_10M_2010.shp
or NUTS_RG_03M_2010.shp

 

Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td,
"rgeos_intersection"): TopologyException: Input geom 0 is invalid:
Self-intersection at or near point 0.68803650000000005 51.7317125 at
0.68803650000000005 51.7317125

 

 

The conflicts originate from the UK map, somewhere above lat 50. Steps to
reproduce:

 

library(rgdal)

  s <- readOGR('EUmaps/NUTS_2010_10M_SH/data',layer="NUTS_RG_10M_2010")

s1 <- s[s$STAT_LEVL_==2 & substring(s$NUTS_ID,1,2) %in% c("UK"),]

# plot(s1)

clip.extent <- as(extent(-10, 19, 36, 53), "SpatialPolygons")

proj4string(clip.extent) <- CRS(proj4string(s))

gI <- gIntersects( s1 , clip.extent , byid = TRUE )

out <- lapply( which(gI) , function(x){ gIntersection( s1[x,] , clip.extent)
} )

 

# But let's look at what is returned

table( sapply( out , class ) )

 

# We want to keep only objects of class SpatialPolygons                 

keep <- sapply(out, class)

out <- out[keep == "SpatialPolygons"]

 

# Coerce list back to SpatialPolygons object

UK <- SpatialPolygons( lapply( 1:length( out ) , function(i) { Pol <-
slot(out[[i]], "polygons")[[1]];   slot(Pol, "ID") <- as.character(i)

   Pol

}))

plot( UK )

 

 

What does this self-intersection means?

 

 

 

Regards,

 

 

Philippe

 


	[[alternative HTML version deleted]]