I'm wondering if this is the right way to do this, or if there's an easier way.

For illustration sake, take this SpatialPolygons object:

main = cbind(c(0, 0, 1, 1),
             c(0, 1, 1, 0))
hole = main/3 + 1/3
island = cbind(c(1.05, 1.05, 1.55, 1.55),
               c(0, .5, .5, 0))
P = Polygons(list(Polygon(main),
                  Polygon(hole, hole = TRUE),
                  Polygon(island)), ID = 1)
SP = SpatialPolygons(list(P))
plot(SP, col = 'red')


I'd like to clip the "island" polygon while retaining the topology of the "main" polygon -- namely, that is has a hole (in my case, it has several holes, but this should do).

Here's what I came up with:

#extract sub-polygons of this polygon
subP = SP@polygons[[1L]]@Polygons
#use the area to identify the largest
areas = sapply(subP, slot, 'area')
#keep holes or the largest actual polygon
keep = sapply(subP, function(p) p@hole | p@area == max(areas))
#construct a new SpatialPolygons object from scratch using index
rem = SpatialPolygons(list(Polygons(subP[keep], ID = 1)))
plot(rem, col = 'red')


That strikes me as being a bit circuitous -- am I missing a method that would handle this more directly? Other related questions I've seen speak more about a SpatialPolygons object that's made up of many polygons instead of one polygon and many Polygons (alas rendering this nuanced search all but impossible).


I spent way too much time thinking about this before remembering my unproblematic fave, the rmapshaper package. I think you should try using ms_explode to get to singlepart objects (or sp::disaggregate, per that tool's own recommendation, for SPDFs), and then experiment with the tolerances on ms_filter_islands.

The brute force method would involve looking for non-hole polygon parts whose coordinates aren't in use elsewhere in the object, but all that recursive searching seems like it would be pretty inefficient, even on a fortified object.

| improve this answer | |

When I saw your question I was thinking about this question here ("Transform raster donuts to circles") and the corresponding answer from @Kersten. The solution is based on a Mathematical Morphology application. Python is used in this answer.

There is also a CRAN - Package called mmand (https://cran.r-project.org/web/packages/mmand/index.html). Hopefully doing the same stuff.

Depending on your data you could convert you polygons to a raster, doing some 'Mathematical Morphology' and convert the raster back to polygons.

| improve this answer | |

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.