I'm wondering if this is the right way to do this, or if there's an easier way.
For illustration sake, take this
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).