Skip to content

Handling invalid geometries from the ‘mapdata’ package

7 messages · Karl Ove Hufthammer, Roger Bivand

Karl Ove Hufthammer
# 
I have discovered what looks like a geometry error in data from ?worldHires?
data set. Here?s an illustration:

library(rgeos)
library(rgdal)o 
library(mapdata)
library(maptools)

x.map = map("worldHires", "Norway", fill=TRUE, col="transparent", plot=FALSE, xlim=c(0,40), ylim=c(50,73))
x.sp = map2SpatialPolygons(x.map, x.map$names, proj4string=CRS("+init=epsg:4326"))

plot(x.sp)            # *Looks* OK
gPointOnSurface(x.sp) # Error: TopologyException: side location conflict at 3.97911 0.0769463
gIsValid(x.sp)        # Error: Self-intersection at or near point 27.9364 70.0864
points(27.9364, 70.0864, col="red", pch=19)

I guess this is to be be expected, as the ?worldHires? data isn?t of very
high quality, and is not guaranteed to be consistent. My questions are 
1) is it possible to automatically ?fix? this problem inside R, and 
2) what other problems may I expect from continuing to use the invalid data?
Roger Bivand
#
On Mon, 28 Mar 2011, Karl Ove Hufthammer wrote:

            

      
      
      
    

        
To start with, it is always advisible to promote the sp SpatialPolygons 
object Polygons components to GEOS compatibility:

pls <- slot(x.sp, "polygons")
pls1 <- lapply(pls, checkPolygonsHoles)
slot(x.sp, "polygons") <- pls1

which adds a comment attribute to each Polygons object assigning hole 
interior rings to their containing exterior rings.

Your problem cannot be solved automatically, but here from zooming in on 
the problem area:

gIsValid(x.sp)
plot(x.sp, xlim=c(27.9, 28), ylim=c(69.9, 70.1))

and eyeballing the bounding boxes of the Polygons objects:

lapply(pls1, bbox)

it looks as though 26 is the problem:

plot(x.sp[26], add=TRUE, border="orange")

x.sp1 <- x.sp[-26]
gIsValid(x.sp1)
gPointOnSurface(x.sp1)

Fresh releases of rgeos and maptools are on their way to CRAN, but source 
may be checked out from R-Forge now (the package builds there are of 
yesterday's version).

Hope this helps,

Roger

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Karl Ove Hufthammer
    

    

      
      #
    

  


  

      
Roger Bivand wrote:

        

            

      
      
      
    

        
Thanks! That helped. ?Norway:Tana? is really a very strange polygon ... :)

I also tried running your procedure on the whole ?worldHires?, to see if 
this is a frequent problem. The idea was to run

  status=gIsValid(x.sp, byid=TRUE)
  status[!status]

to get a list of invalid polygons (though perhaps the SpatialPolygons object 
may be invalid even if all the sub-polygons are valid?)

First I had to remove the proj4string, as some of the polygons are outside 
the [-180, 180] range. And it turned out that I got an error message already 
at the

  pls <- slot(x.sp, "polygons")
  pls1 <- lapply(pls, checkPolygonsHoles)
  slot(x.sp, "polygons") <- pls1

step.

So I guess it?s better to stay away from the ?worldHires? data set,
and use GSHHS instead ? :-/

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Mon, 28 Mar 2011, Karl Ove Hufthammer wrote:

        

            

      
      
      
    

        
If the error was:

            

      
      
      
    

        
Error in RGEOSBinPredFunc(spgeom1, spgeom2, byid, "rgeos_equals") :
   TopologyException: side location conflict at 78.9364 33.4081

then I can reproduce it, and will examine whether it can be handled - the 
checkPolygonsHoles() function also uses GEOS, but probably should warn 
after failing try(gEquals(...)).

Was this your error?

Throwing bad data at maptools/rgeos is a nice way of stress-testing it.

Roger

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Karl Ove Hufthammer
    

    

      
      #
    

  


  

      
Roger Bivand wrote:

        

            

      
      
      
    

        
Yes.

I also eventually managed to get a segfault when trying various stuff like 
this, but it may have been unrelated. (At least, I can?t reproduce it.)

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Mon, 28 Mar 2011, Karl Ove Hufthammer wrote:

        

            

      
      
      
    

        
Thanks. Updating to maptools 0_8-5 and rgeos 0_1-3 (sources on CRAN now) 
should provide better protection against non-conformant data causing a 
seg-fault in GEOS - GEOS really likes its polygons *very* tidy.

Roger

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Mon, 28 Mar 2011, Roger Bivand wrote:

        

            

      
      
      
    

        
The attached plot shows this issue, which is in the area between northern 
India, Kashmir, and China, so perhaps reality is more complicated than 
GEOS prefers? I've only plotted the border given in the data set for 
China, which includes the odd reversed straight line. If I include other 
borders, I see a line starting at right angles to the northern end of this 
line and heading WSW, which may be another demarcation line. Apart from 
this case, the rest of worldHires got through OK - the maptools code on 
SVN now reports:

            

      
      
      
    

        
Warning message:
In checkPolygonsGEOS(x, properly = properly) :
   Polygons object China, Polygon 1: Error in RGEOSBinPredFunc(spgeom1, 
spgeom2, byid, "rgeos_equals") :
   TopologyException: side location conflict at 78.9364 33.4081

instead of failing.

Roger