gIntersection fails
Please never post HTML, only plain text. Always include a link to data, but not this data, as this requires acceptance of conditions to download, so is not automatic. Note that the SO posting you quote was from April 2013 - it is always important to give your sessionInfo() output, as how things work changes over time, leaving legacy advice stale. Your results might match the SO ones if your platform is the same as that was then. In particular, the: http://cran.r-project.org/web/packages/rgeos/ChangeLog entry for 2014-02-08 indicates that non-polygon output from polygon intersections is now handled more gracefully. The underlying problem may be that these coordinates are geographical, not projected, and the integer discretisation used by GEOS is failing for some input - have you tried projecting to planar coordinates in a different metric and doing the same operation? Or using setScale()? For my 20M 2006: library(rgdal) s <- readOGR('.', layer="NUTS_RG_20M_2006") 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"),] library(raster) clip.extent <- as(extent(-10, 19, 36, 51), "SpatialPolygons") proj4string(clip.extent) <- CRS(proj4string(s)) EuropeW <- gIntersection(s1, clip.extent, byid=c(TRUE, FALSE)) plot(EuropeW) just works. Hope this clarifies, Roger PS:
sessionInfo()
R version 3.1.2 (2014-10-31) Platform: x86_64-unknown-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_US.utf8 LC_NUMERIC=C [3] LC_TIME=en_US.utf8 LC_COLLATE=en_US.utf8 [5] LC_MONETARY=en_US.utf8 LC_MESSAGES=en_US.utf8 [7] LC_PAPER=en_US.utf8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.utf8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] rgeos_0.3-8 raster_2.3-12 rgdal_0.9-1 sp_1.0-16 loaded via a namespace (and not attached): [1] grid_3.1.2 lattice_0.20-29
version_GEOS()
[1] "3.4.2-CAPI-1.8.2 r3921"
On Tue, 16 Dec 2014, Philippe Liege wrote:
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]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 91 00 e-mail: Roger.Bivand at nhh.no