Add function to split geometries

[dbaston]

Many applications require the ability to split a geometry by a line. This has been implemented in PostGIS, QGIS, and when using the dateline as a splitting geometry, GDAL. Implementing this functionality in GEOS would allow applications to avoid maintaining their own implementations, and could likely provide better performance by avoiding the limitations of the C API.

[dr-jts]

This has been on the ToDo list for a while. The machinery in OverlayNG should accomodate this, with some reworking.

[jorisvandenbossche]

Shapely also has a custom split() function (https://shapely.readthedocs.io/en/stable/manual.html#shapely.ops.split, https://github.com/shapely/shapely/blob/70fcae0561264f8cc2175e27cf5c0c287af8fcf0/shapely/ops.py#L361-L547), so adding this to GEOS would be welcome!

[dbaston]

I'm looking into how we might accomplish this using OverlayNG, using a simple case of a square split with a line. It is straightfoward to mark a set of OverlayEdges to include in the result:

But we run into trouble when assembling a geometry from these edge. The PolygonBuilder logic isn't designed for a case like this -- its job is to build valid geometries, and one cannot be constructed from these edges. In order to produce the desired result (two polygons), it seems like the MaximalEdgeRing logic would need to link each edge not to the next edge with isInResultArea() (as it does now), but instead give a preference to edges on the boundary of the splitting geometry.

@dr-jts, I'm curious if this is something you've given thought to, and how you might approach it.

[dr-jts]

I'm curious if this is something you've given thought to, and how you might approach it.

I've thought about this at a high level. The problem of determining the set of faces in the result is isomorphic to polygonization, so that code provides a pattern to follow. Either the OverlayNG polygon building logic needs to be generalized, or else made polymorphic, with an alternative strategy used for splitting as opposed to overlay.

This probably requires some generalization of the edge labelling as well. I assume the desired result is a pair of valid polygonal geometries, one for either side of the splitting line (can be produced as a GeometryCollection with 2 elements). This will require the labelling to include information about the line side, so the polygonized faces can be assigned to the appropriate result geometry. (Some thought should be given to how to handled self-intersecting splitting lines. The polygonization is straightforward, but the labelling and assignment requires more logic).

There may be other impacts on the OverlayNG code as well. It contains some finicky logic to deal with mixed-type results, which might need to be bypassed?