[R-sig-Geo] gIntersection fails
Philippe Liege
pliege at anisteme.net
Wed Dec 17 16:31:24 CET 2014
> 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
More information about the R-sig-Geo
mailing list