[R-sig-Geo] gIntersection fails

Wed Dec 17 13:18:44 CET 2014

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