gIntersection fails
One has to find the appropriate scale manually in such cases, curiously
I have no idea to generate the scale programmatically.
Well, the EuroStat dataset goes up to French Guiana to the west, then to
Reunion Island to the South etc. It is deprived of any index to identify
those overseas regions.
Either we create a new index of we clip the area to be plotted. In both
cases we have to manually retrieve the extreme points of this area. In the
present case, one can start from "Cabo de Roca" in Portugal etc. A much
simple solution might be to exclude France, Spain and Portugal from the EU
and plot the remaining countries...
I slightly changed the script so that the outcome sp object is now of class
SpatialPolygonsDataFrame.
I couldn't find out Andorra in the Eurostat dataset. No NUTS_ID starting
with AD.
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("BE","CH","DE","DK","ES","IE","FR","IT","LU","NL","PT","UK"),]
# 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))
clip.extent <- as(extent(-10, 19, 35, 59), "SpatialPolygons")
proj4string(clip.extent) <- CRS(proj4string(s1))
clip.extent <- spTransform(clip.extent,proj.laea)
summary(clip.extent)
EuropeW <- gIntersection(s1.laea, clip.extent, byid=c(TRUE, FALSE))
# side effect: the data frale is removed
# You may need a SpatialPolygonDataFrame for future use
EuropeW <- s1.laea[sapply(slot(s1.laea,"polygons"),slot,"ID") %in%
sapply(slot(EuropeW,"polygons"),slot,"ID"),]
EuropeW$NUTS_ID <- EuropeW$NUTS_ID[drop=T]
plot(EuropeW,col=brewer.pal(10,"Set3"))
plot(EuropeW,col=brewer.pal(12,"Set3")[as.factor(substring(EuropeWdf$NUTS_ID
,1,2))])
-----Message d'origine-----
De?: Roger Bivand [mailto:Roger.Bivand at nhh.no]
Envoy??: mercredi 17 d?cembre 2014 13:19
??: Philippe Liege
Cc?: r-sig-geo at r-project.org
Objet?: RE: [R-sig-Geo] gIntersection fails
On Wed, 17 Dec 2014, Philippe Liege wrote:
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
OK, thanks for this. The resolution is either as you chose, or: ce_ll <- spTransform(clip.extent, CRS(proj4string(s1))) oS <- getScale() oS # [1] 1e+08 EuropeW <- gIntersection(s1, ce_ll, byid=c(TRUE, FALSE)) # Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, # "rgeos_intersection") : # TopologyException: Input geom 0 is invalid: Self-intersection at or # near point 10.648729489383225 57.7454124997018 at 10.648729489383225 # 57.7454124997018 setScale(1e+9) EuropeW <- gIntersection(s1, ce_ll, byid=c(TRUE, FALSE)) # OK setScale(oS) One has to find the appropriate scale manually in such cases, curiously setScale(1e+7) also works, but for geographical coordinates one probably wants to increase rather than decrease. Roger (I'm assuming that NL hasn't been flooded?)
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
-- 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