R Polygon Splitting – How to Split Polygons at Narrowest Part Using R

polygonrsplitting

I have a dataset of polygons like this one:

enter image description here

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?

Best Answer

Given a polygon pol, like this:

enter image description here

then:

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

produces this:

enter image description here

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.

Related Question