18

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.

  • Without your data, it is hard to say where is the problem. – user32309 Sep 16 '14 at 9:16
  • @Pascal, I just uploaded a dropbox link with a slimmed down shapefile of 10mb zipped and 16mb unzipped that will reproduce the problem. I wasn't sure how to provide the data as the original was 1.5 gb, but managed to use ArcGIS to narrow the issue to a smaller file. Is there a good protocol for creating and sharing manageable sized reproducible examples? – Luke Macaulay Sep 16 '14 at 19:22
  • Trying different approaches with R didn't work. And Qgis is freezing when checking geometries. – user32309 Sep 18 '14 at 9:28
25
+50

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)

Union result

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

  • Thanks eblondel! I will try this out. Thanks for the package development! – Luke Macaulay Sep 29 '14 at 20:37
  • Hi eblondel, This worked well, but I wanted to let you know that in correcting the geometry, it would oftentimes create a very large polygon when dealing with intricate and complex features. For example a road network was corrected to a large polygon that was basically the extent of the network. I'm not sure how easy that is to correct, but wanted to let you know. – Luke Macaulay Feb 6 '15 at 9:14
  • Wow. Very impressive package. That must have been a lot of work. – nograpes Jul 7 '15 at 17:25
  • 3
    Thanks @nograpes for your feedback. I've built this package from scratch when this issue was posted, also because cleaning geometries is not always an easy task. If you are on Github, i would welcome your 'star' :-), I would like to further improve the package in the future, and possibly release it on CRAN. – eblondel Jul 7 '15 at 18:55
  • 7
    Just to let you know that cleangeo has been published in CRAN (cran.r-project.org/package=cleangeo), to all people that use it, feel free to report enhancement requests or bugs in Github. – eblondel Sep 27 '15 at 9:43
0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.