[R-sig-Geo] gIntersection fails

[Roger Bivand]

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