# Eliminating island sub-polygons

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.