[R-sig-Geo] gIntersection fails
Philippe Liege
pliege at anisteme.net
Wed Dec 17 11:58:58 CET 2014
Roger,
Thanks for your detailed answer. Sure it helps!
I have projected the SpatialPolygonsDataFrame object to planar coordinates
using etrs89/etrs-laea, and it works fine now.
The Eurostat website has been redesigned and links to publically available
reference data have not yet been updated in Google. The datasets for EU
Administrative Units/Statistical Units (the so-called NUTS files) can be
found here:
http://ec.europa.eu/eurostat/c/portal/layout?p_l_id=6033084&p_v_l_s_g_id=0
And the file I was referring to in my previous e-mail can be downloaded
from:
http://ec.europa.eu/eurostat/cache/GISCO/geodatafiles/NUTS_2010_10M_SH.zip
library(rgdal)
library(RColorBrewer)
library(raster)
library(rgeos)
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("FR","BE","ES","IT","PT","DE","CH","UK","IE","DK"),]
# ETRS89 / ETRS-LAEA
# as in: http://spatialreference.org/ref/epsg/3035/proj4/
proj.laea <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000
+ellps=GRS80 +units=m +no_defs")
s1.laea <-spTransform(s1,proj.laea)
summary(s1.laea)
clip.extent <- as(extent(2.5e6, +5e6, 1.5e6, 4.1e6), "SpatialPolygons")
proj4string(clip.extent) <- CRS(proj4string(s1.laea))
EuropeW <- gIntersection(s1.laea, clip.extent, byid=c(TRUE, FALSE))
plot(EuropeW,col=brewer.pal(10,"Set3"))
Regards,
Philippe
-----Message d'origine-----
De : Roger Bivand [mailto:Roger.Bivand at nhh.no]
Envoyé : mardi 16 décembre 2014 19:36
À : Philippe Liege
Cc : r-sig-geo at r-project.org
Objet : Re: [R-sig-Geo] 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
More information about the R-sig-Geo
mailing list