Splitting polygons at narrowest part 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.

This would take you rather far afield, but one way to define the most narrow waist of a shape is to identify the event at which the straight skeleton partitions the shrinking shape into two pieces.

---

         
          ^{Image from here (and probably copied from elsewhere).}

---

You can view the polygon boundary edges moving inward at a fixed rate—a type of inward offset. A narrow neck will pinch off the shape into two pieces. In the above image, this happens a bit beyond t=4.