Xjyuk

I have a dataset of polygons like this one:

I would like to split these polygons in separate parts at their most narrow location (if they have one). For example the two small polygons should not be split. So I would need to identify narrow locations and then split the polygon there.

How could this be done using R?

Given a polygon pol, like this:

then:

> library(sf)
> sdist = -0.055168
> ppol = splitnarrow(pol, sdist, 1e-3)
> plot(ppol, col=1:2)

produces this:

Here's the source code for splitnarrow. There's a zillion places where this can go wrong, and first you have to determine sdist and eps for your polygons.

splitnarrow <- function(pol, sdist, eps)
###
### split a polygon at its narrowest point.
###

### sdist is the smallest value for internal buffering that splits the
### polygon into a MULTIPOLYGON and needs computing before running this.

### eps is another tolerance that is needed to get the points at which the
### narrowest point is to be cut.

 ## split the polygon into two separate polygons
 bparts = st_buffer(pol, sdist)
 features = st_cast(st_sfc(bparts), "POLYGON")

 ## find where the two separate polygons are closest, this is where
 ## the internal buffering pinched off into two polygons.

 pinch = st_nearest_points(features[1],features[2])

 ## buffering the pinch point by a slightly larger buffer length should intersect with
 ## the polygon at the narrow point. 
 inter = st_intersection(
 st_cast(pol,"MULTILINESTRING"),
 st_buffer(pinch,-(sdist-(eps))
 )
 )
 join = st_cast(st_sfc(inter), "LINESTRING")

 ## join is now two small line segments of the polygon across the "waist".
 ## find the line of closest approach of them:
 splitline = st_nearest_points(join[1], join[2])

 ## that's our cut line. Now put that with the polygon and make new polygons:
 mm = st_union(splitline, st_cast(pol, "LINESTRING"))
 parts = st_collection_extract(st_polygonize(mm))
 parts

sdist is the smallest value that splits the single polygon into a multipolygon, and eps is the smallest value that touches both sides when buffered from the waist intersection point. Finding these could be automated.