Fixing orphaned holes in R

Question

I'm trying to perform a union on a common field after merging two adjacent shapefiles. The shapefiles end up with at least one thin sliver of space between them. When I attempt a union I get the following orphaned hole error:

Error in createPolygonsComment(p) : rgeos_PolyCreateComment: orphaned hole, cannot find containing polygon for hole at index 17

I've uploaded a reproducible example to Dropbox at this link.

Here is the code to recreate the problem:

#loading required packages
require(sp)    
require(rgdal)
require(maptools)
require(rgeos)

#load example data, set "dsn=" to your working directory or specify the path
example <- readOGR(dsn=".",layer="ReproducibleExample")

#Attempting a UnionSpatialPolygons based on the COUNTY field
example.df <- as(example, "data.frame")
countycol <- example.df$COUNTY
example.diss <- unionSpatialPolygons(example, countycol)

Returns:

Error in createPolygonsComment(p) : rgeos_PolyCreateComment: orphaned hole, cannot find containing polygon for hole at index 17

Trying the fix proposed here and here:

slot(example, "polygons") <- lapply(slot(example, "polygons"), checkPolygonsHoles)

This returns the same error that comes from the union attempt but with different index number:

rgeos_PolyCreateComment: orphaned hole, cannot find containing polygon for hole at index 30

Trying the fix proposed in Roger Bivand's helpful tutorial

fix <- slot(example, "polygons")
fixa <- lapply(fix, checkPolygonsHoles)

Returns same error at index 30 as above.

Others have raised this problem here and here, and while the solutions posited above appear to work for some cases, other cases are not resolved. One user used QGIS to address the problem, and the other had 2 of 3 items fixed, but no resolution for the final one.

It appears that people continue to have problems despite this code working from time to time. Has anybody found a solution within R?

I've performed "repair geometry" tool in ArcGIS, and it corrected the problem, but it seems like there should be a fix in R.

Answer 1

I've analysed the geometry issues in the attached data, and it seems it does not ONLY have orphaned holes but also geometry validity issues. It's true that an orphaned hole is somehow a geometry validity issue, but rgeos does not handle it in the same way, as for orphaned holes, an error is raised, instead of a simple warning. As you indicate, they are hints to check polygon holes, but it is not always successfull when applied in order to fix orphaned holes.

So, let's:

1. clean your data (which is required if you wish to do geoprocessing like union)

2. use the cleaned data with your union process

1. Cleaning geometry Fixing geometries in R can be sometimes challenging, so i've tried to built an experimental R package (see https://github.com/eblondel/cleangeo) that intends to facilitate cleaning of sp objects (at now limited on polygonal shapes). You can install the package with:

require(devtools)
install_github("eblondel/cleangeo")
require(cleangeo)

To start, it's good that you see what are the geometry issues with your source data. For this, you can run the following (your data is large so it can take some time):

#get a report of geometry validity & issues for a sp spatial object
report <- clgeo_CollectionReport(sp)
summary <- clgeo_SummaryReport(report)
issues <- report[report$valid == FALSE,]

With this, you will see that your data has 2 kinds of issues: orphaned holes and geometry validity issues. Both (and not only the orphaned holes) are likely to make the union process failing, so the data should be cleaned before, in an automate way when possible. For a fast reproduction, the first sample code below only takes the subset of features that are tagged as suspicious (except the latest one, with index = 9002 in the original data - see my note below on this)

#get suspicious features (indexes)
nv <- clgeo_SuspiciousFeatures(report)
mysp <- sp[nv[-14],]

#try to clean data
mysp.clean <- clgeo_Clean(mysp, print.log = TRUE)

#check if they are still errors
report.clean <- clgeo_CollectionReport(mysp.clean)
summary.clean <- clgeo_SummaryReport(report.clean)

If clgeo_Clean does well the job, you should get all geometries valid now. You can apply this to the complete dataset (except feature index = 9002)

#try to clean data
mysp <- sp[-9002,]
mysp.clean <- clgeo_Clean(mysp, print.log = TRUE)

#check if they are still errors
report.clean <- clgeo_CollectionReport(mysp.clean)
summary.clean <- clgeo_SummaryReport(report.clean)

2. Union process Now, let's see if the union works on this dataset:

#Attempting a UnionSpatialPolygons based on the COUNTY field
mysp.df <- as(mysp, "data.frame")
countycol <- mysp.df$COUNTY
mysp.diss <- unionSpatialPolygons(mysp.clean, countycol)

Note: as said before, i've remove one feature (index = 9002).By plotting it: plot(sp[9002,]), you will see that this feature is very (very) complex. I've excluded it from the sample only because checking holes was taking too much time. Let's see now if the same problem occurs using readShapePoly (from maptools) for reading the data...

3. Switch to readShapePoly vs. readOGR for reading data (UPDATE)

readOGR is not the only function available to read shapefiles. You can also use readShapePoly from maptools package, generally more performant than the first one:

require(maptools)
mysp <- readShapePoly("ReproducibleExample.shp")

Apart from running faster:

• if you use the above code based on clgeo_CollectionReport, there is no problem of orphaned holes, but still problems of geometry.

• Cleaning the geometry with clgeo_Clean also runs well, and now it doesn't get stuck with the feature index 9002

• And... the union process works.

See below the plot result:

#plot the result
plot(mysp, border= "lightgray")
plot(mysp.diss, border="red", add = TRUE)

Conclusion: prefer maptools to read your shapefile data, and consider using cleangeo to clean your data before any geoprocessing.

Answer 2

A convenient solution that keeps working for me in R is to apply a zero-width buffer:

#loading required packages
require(sp)    
require(rgdal)
require(maptools)
require(rgeos)

#load example data, set "dsn=" to your working directory or specify the path
example <- readOGR(dsn=".",layer="ReproducibleExample")

#project your data (I'm using California Albers here) and apply a zero-width buffer
example <- spTransform(example, CRS("+init=epsg:3310"))
example <- gBuffer(example, byid = T, width = 0)

#Attempting a UnionSpatialPolygons based on the COUNTY field
example.df <- as(example, "data.frame")
countycol <- example.df$COUNTY
example.diss <- unionSpatialPolygons(example, countycol)

unionSpatialPolygons takes a while with this data set, but seems to work just fine.