Keep largest ring of multi-ring polygon

Question

I have a SpatialPolygons object. This object has a single feature, a Polygons object, that is itself made up of multiple Polygon objects.

I'd like to subset the SpatialPolygons object so that the Polygons object only has a single Polygon object, which is the Polygon with the largest area.

I've tried many different approaches and can't figure out how to do the subsetting at a sub-Polygons level.

Here is an example:

#create polygon
Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))
Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))
Srs = Polygons(list(Sr1,Sr2), "s1/2")

SpP = SpatialPolygons(list(Srs))

In this example, SpP has a single Polygons object. The first sub-polygon of the Polygons object has area 5.5, and the second has area 1.5. I'd like to subset the SpatialPolygons object such that the Polygons object only has one Polygon, the one with area 5.5.

Is this possible?

EDIT

Calum You thank you for your answer. I've never used the sf package, but it looks cool. Especially since it plays so nicely with dplyr!

I ended up finding a different solution though. My main issue was that I didn't grasp how the SpatialPolygons object structure worked. A SpatialPolygons object contains 1..* Polygons objects, so the constructor takes a list of Polygons objects. A Polygons object contains 1..* Polygon objects, so the constructor takes a list of Polygon objects. With this in mind, I used the following solution.

plst <- SpP@polygons[[1]]@Polygons #get the list of Polygon objects
plst <- plst[which.max(sapply(plst,function(p) return(p@area)))] #filter to the max area
spoly2 <- SpatialPolygons(list(Polygons(plst,'id')))

Plot shows this will filter the rings by the largest area.

plot(SpP)
plot(spoly2,col=alpha('red',0.1),add=T)

Is there just the one object? Are you creating all of them manually like this or is this just an example? Trying to figure out what the best way is — Calum You

Calum You Calum You · Accepted Answer · 2018-05-09T22:20:10

Here's one approach that uses the sf package. Unfortunately I am not too familiar with sp, but perhaps this will be an incentive to switch!

The problem is that as given, the geometry is a single polygon with an external outer ring. To work on the geometry separately, we need to split each ring into its own polygon. In sf, a POLYGON is a list of matrices, each matrix corresponding to the points in a ring; a MULTIPOLYGON is a list of theses POLYGON lists of matrices. So to convert we do the following:

(convert the SpatialPolygons into an sfc with st_as_sfc)
map(list) across the POLYGON wrapping each ring in a list, and then convert the list of lists with st_multipolygon. Note I also map across the sfc column in case there are more geometries, though here you only have one.
Convert this list back into an sfc object and then an sf object. The sfc is just the list of geometries with additional properties like a coordinate system, the sf makes the sfc just one column in a data frame. This is needed because we want to know the area of each geometry.
Use the wonderful function st_cast to split the MULTIPOLYGON into individual polygons. We added a rowid so that we know which new POLYGONs belonged to the original MULTIPOLYGON.
Add an area column with mutate and st_area, group_by it, and then use top_n to filter down to only the largest polygon per group. Note how nicely sf objects work with dplyr verbs!

The plot shows the result. The original polygon is in red, the largest ring of it is in blue.

library(tidyverse)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3
library(sp)
Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))
Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))
Srs = Polygons(list(Sr1,Sr2), "s1/2")

SpP = SpatialPolygons(list(Srs))
sfc <- st_as_sfc(SpP)
sfc
#> Geometry set for 1 feature 
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 1 ymin: 2 xmax: 5 ymax: 5
#> epsg (SRID):    NA
#> proj4string:    NA
#> POLYGON ((2 2, 1 4, 4 5, 4 3, 2 2), (5 2, 4 3, ...

sf_obj <- sfc %>%
  map(
    .f = function(polygon) polygon %>% map(list) %>% st_multipolygon()
  ) %>%
  st_sfc() %>%
  st_sf() %>%
  rowid_to_column() %>%
  st_cast("POLYGON") %>% 
  mutate(area = st_area(geometry)) %>%
  group_by(rowid) %>%
  top_n(1, area)
#> Warning in st_cast.sf(., "POLYGON"): repeating attributes for all sub-
#> geometries for which they may not be constant

plot(sfc, col = "red")
plot(sf_obj$geometry, col = "blue", add = TRUE)

Created on 2018-05-09 by the reprex package (v0.2.0).

Keep largest ring of multi-ring polygon

1 Answers