Repository: shapely/shapely

Issue:
Expose TPVWSimplifier
Does shapely already expose the GEOS implementation of [Topology-Preserving Visvalingam-Whyatt Line Simplification](https://libgeos.org/doxygen/classgeos_1_1coverage_1_1TPVWSimplifier.html) through its Python API? (I couldn’t find it in the docs, but I might have looked at the wrong place). If it’s not exposed yet, would a Python wrapper be in scope for shapely?

Discussion:
jorisvandenbossche: @martinfleis I finally took the time to properly re-read our discussion earlier to understand again what we decided, and your last comment and updated implementation. And I think your implementation is actually just fine, and it's only a matter of updating the tests to follow this new behaviour (4th option).

> Right now, it works in principle but does not work like a ufunc, so casting it on a GeoSeries does not preserve GeoSeries output etc... And I think that this logic should happen in C but I think I've hit the limits of my C and wasn't able to make it work.


I think this is a bit the consequence of choosing the 4th option, because that option does not fit in the numpy ufunc model (it is not a simple element-wise ufunc, but also not a generalized ufunc with eg a reducing or reused axis). It is a function that preserves the shape, but looks at all the elements together, and so not a ufunc. 

I suppose we could still make it work that it behaves like a ufunc in the sense of preserving the original object (I think that would then rather use the `__array_function__` protocol instead of `__array_ufunc__`, see https://numpy.org/doc/stable/user/basics.interoperability.html#operating-on-foreign-objects-without-converting) 
But this seems like less critical, and something I could take a look at.


martinfleis: I am getting a strange bare `GEOSException` when passing a scalar Polygon. Everything else is fine. I thought that it may be some GEOS issue but using command line `geosop` is fine. No clue where that comes from.

Patch:
diff --git a/CHANGES.txt b/CHANGES.txt
index b6f8dca..3220959 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -33,6 +33,8 @@ Improvements:
 - Add function ``minimum_clearance_line`` (#2106)
 - Add function ``maximum_inscribed_circle`` (#1307)
 - Add function ``orient_polygons`` (#2147)
+- Add function ``coverage_simplify`` to allow topological simplification of polygonal
+  coverages (#1969)
 
 Breaking changes in GEOS 3.12:
 
diff --git a/docs/release/1.x.rst b/docs/release/1.x.rst
index 6f673bb..356f462 100644
--- a/docs/release/1.x.rst
+++ b/docs/release/1.x.rst
@@ -361,7 +361,7 @@ Bug fixes:
 1.6.1 (2017-09-01)
 ------------------
 
-- Avoid ``STRTree`` crashes due to dangling references (#505) by maintaining
+- Avoid ``STRtree`` crashes due to dangling references (#505) by maintaining
   references to added geometries.
 - Reduce log level to debug when reporting on calls to ctypes ``CDLL()`` that
   don't succeed and are retried (#515).
diff --git a/docs/strtree.rst b/docs/strtree.rst
index 26679a8..d4803ae 100644
--- a/docs/strtree.rst
+++ b/docs/strtree.rst
@@ -1,4 +1,4 @@
-STRTree
+STRtree
 =======
 
 .. autoclass:: shapely.STRtree
diff --git a/shapely/constructive.py b/shapely/constructive.py
index 1bbbbfb..915e2ed 100644
--- a/shapely/constructive.py
+++ b/shapely/constructive.py
@@ -2,8 +2,9 @@
 
 import numpy as np
 
-from shapely import lib
+from shapely import Geometry, GeometryType, lib
 from shapely._enum import ParamEnum
+from shapely._geometry import get_parts
 from shapely.algorithms._oriented_envelope import _oriented_envelope_min_area_vectorized
 from shapely.decorators import multithreading_enabled, requires_geos
 from shapely.errors import UnsupportedGEOSVersionError
@@ -18,6 +19,7 @@ __all__ = [
     "clip_by_rect",
     "concave_hull",
     "convex_hull",
+    "coverage_simplify",
     "delaunay_triangles",
     "segmentize",
     "envelope",
@@ -1031,6 +1033,71 @@ def simplify(geometry, tolerance, preserve_topology=True, **kwargs):
         return lib.simplify(geometry, tolerance, **kwargs)
 
 
+@requires_geos("3.12.0")
+@multithreading_enabled
+def coverage_simplify(geometry, tolerance, simplify_boundary=True):
+    """Return a simplified version of an input geometry using coverage simplification.
+
+    Assumes that the geometry forms a polygonal coverage. Under this assumption, the
+    function simplifies the edges using the Visvalingam-Whyatt algorithm, while
+    preserving a valid coverage. In the most simplified case, polygons are reduced to
+    triangles.
+
+    A collection of valid polygons is considered a coverage if the polygons are:
+
+    * **Non-overlapping** - polygons do not overlap (their interiors do not intersect)
+    * **Edge-Matched** - vertices along shared edges are identical
+
+    The function allows simplification of all edges including the outer boundaries of
+    the coverage or simplification of only the inner (shared) edges.
+
+    If there are other geometry types than Polygons or MultiPolygons present, the
+    array will not undergo simplification and geometries are returned unchanged.
+
+    If the geometry is polygonal but does not form a valid coverage due to overlaps,
+    it will be simplified but it may result in invalid topology.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    tolerance : float or array_like
+        The degree of simplification roughly equal to the square root of the area
+        of triangles that will be removed.
+    simplify_boundary : bool, optional
+        By default (True), simplifies both internal edges of the coverage as well
+        as its boundary. If set to False, only simplifies internal edges.
+
+    Returns
+    -------
+    numpy.ndarray | shapely.Geometry
+
+    Examples
+    --------
+    >>> from shapely import Polygon, MultiPolygon
+    >>> poly = Polygon([(0, 0), (20, 0), (20, 10), (10, 5), (0, 10), (0, 0)])
+    >>> coverage_simplify(poly, tolerance=2)
+    <POLYGON ((0 0, 20 0, 20 10, 10 5, 0 10, 0 0))>
+    """
+    scalar = False
+    if isinstance(geometry, Geometry):
+        scalar = True
+
+    geometries = np.asarray(geometry)
+    shape = geometries.shape
+    geometries = geometries.ravel()
+
+    # create_collection acts on the inner axis
+    collections = lib.create_collection(
+        geometries, np.intc(GeometryType.GEOMETRYCOLLECTION)
+    )
+
+    simplified = lib.coverage_simplify(collections, tolerance, simplify_boundary)
+    parts = get_parts(simplified).reshape(shape)
+    if scalar:
+        return parts.item()
+    return parts
+
+
 @multithreading_enabled
 def snap(geometry, reference, tolerance, **kwargs):
     """Snap an input geometry to reference geometry's vertices.
diff --git a/shapely/geometry/polygon.py b/shapely/geometry/polygon.py
index 9347a18..489b2c0 100644
--- a/shapely/geometry/polygon.py
+++ b/shapely/geometry/polygon.py
@@ -321,6 +321,9 @@ shapely.lib.registry[3] = Polygon
 
 def orient(polygon, sign=1.0):
     """Return an oriented polygon."""
+    if polygon.is_empty:
+        return polygon
+
     s = float(sign)
     rings = []
     ring = polygon.exterior
diff --git a/src/ufuncs.c b/src/ufuncs.c
index e9c7b8d..3530180 100644
--- a/src/ufuncs.c
+++ b/src/ufuncs.c
@@ -3612,6 +3612,70 @@ static PyUFuncGenericFunction to_geojson_funcs[1] = {&to_geojson_func};
 
 #endif  // GEOS_SINCE_3_10_0
 
+#if GEOS_SINCE_3_12_0
+static char coverage_simplify_dtypes[4] = {NPY_OBJECT, NPY_DOUBLE, NPY_BOOL, NPY_OBJECT};
+static void coverage_simplify_func(char** args, const npy_intp* dimensions, const npy_intp* steps,
+                                    void* data) {
+  GEOSGeometry* in1 = NULL;
+  GEOSGeometry** geom_arr;
+
+  CHECK_NO_INPLACE_OUTPUT(3);
+
+  // allocate a temporary array to store output GEOSGeometry objects
+  geom_arr = malloc(sizeof(void*) * dimensions[0]);
+  CHECK_ALLOC(geom_arr);
+
+  GEOS_INIT_THREADS;
+
+  TERNARY_LOOP {
+    CHECK_SIGNALS_THREADS(i);
+    if (errstate == PGERR_PYSIGNAL) {
+      destroy_geom_arr(ctx, geom_arr, i - 1);
+      break;
+    }
+    // get the geometry: return on error
+    if (!get_geom(*(GeometryObject**)ip1, &in1)) {
+      errstate = PGERR_NOT_A_GEOMETRY;
+      destroy_geom_arr(ctx, geom_arr, i - 1);
+      break;
+    }
+    double in2 = *(double*)ip2;
+    npy_bool in3 = !(*(npy_bool*)ip3);
+
+    int isValid = 1;
+    int numGeoms = GEOSGetNumGeometries_r(ctx, in1);
+    for (int j = 0; j < numGeoms; j++) {
+      GEOSGeometry* geom = GEOSGetGeometryN_r(ctx, in1, j);
+      if (GEOSGeomTypeId_r(ctx, geom) != GEOS_POLYGON && GEOSGeomTypeId_r(ctx, geom) != GEOS_MULTIPOLYGON) {
+          isValid = 0;
+          break;
+      }
+    }
+    if (isValid) {
+      geom_arr[i] = GEOSCoverageSimplifyVW_r(ctx, in1, in2, (int)in3);
+      if (geom_arr[i] == NULL) {
+        errstate = PGERR_GEOS_EXCEPTION;
+        destroy_geom_arr(ctx, geom_arr, i - 1);
+        break;
+      }
+    } else {
+      geom_arr[i] = GEOSGeom_clone_r(ctx, in1);
+    }
+
+  }
+
+
+  GEOS_FINISH_THREADS;
+
+  // fill the numpy array with PyObjects while holding the GIL
+  if (errstate == PGERR_SUCCESS) {
+    geom_arr_to_npy(geom_arr, args[3], steps[3], dimensions[0]);
+  }
+  free(geom_arr);
+}
+static PyUFuncGenericFunction coverage_simplify_funcs[1] = {&coverage_simplify_func};
+
+#endif  // GEOS_SINCE_3_12_0
 /*
 TODO polygonizer functions
 TODO prepared geometry predicate functions
@@ -3874,6 +3938,7 @@ int init_ufuncs(PyObject* m, PyObject* d) {
 #endif
 
 #if GEOS_SINCE_3_12_0
+  DEFINE_CUSTOM(coverage_simplify, 3);
   DEFINE_Y_Y(disjoint_subset_union);
   DEFINE_Y_b(has_m);
   DEFINE_Y_d(get_m);


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index ffd454d..67affa6 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -100,7 +100,7 @@ jobs:
         if: ${{ matrix.geos == 'main' }}
 
       - name: Set up Python ${{ matrix.python }}
-        uses: Quansight-Labs/setup-python@b9ab292c751a42bcd2bb465b7fa202ea2c3f5796  # v5.3.1
+        uses: Quansight-Labs/setup-python@869aeafb7eeb9dc48ba68acc0479e6fc3fd7ce5e  # v5.4.0
         with:
           python-version: ${{ matrix.python }}
           architecture: ${{ matrix.architecture }}
diff --git a/shapely/tests/legacy/test_orient.py b/shapely/tests/legacy/test_orient.py
index 71dce26..e9df865 100644
--- a/shapely/tests/legacy/test_orient.py
+++ b/shapely/tests/legacy/test_orient.py
@@ -39,6 +39,10 @@ class OrientTestCase(unittest.TestCase):
         assert orient(linearring, 1) == linearring
         assert orient(linearring, -1) == linearring
 
+    def test_empty_polygon(self):
+        polygon = Polygon()
+        assert orient(polygon) == polygon
+
     def test_polygon(self):
         polygon = Polygon([(0, 0), (0, 1), (1, 0)])
         polygon_reversed = Polygon(polygon.exterior.coords[::-1])
diff --git a/shapely/tests/test_constructive.py b/shapely/tests/test_constructive.py
index 352f6ef..0c6b5e2 100644
--- a/shapely/tests/test_constructive.py
+++ b/shapely/tests/test_constructive.py
@@ -1097,6 +1097,86 @@ def test_concave_hull_kwargs():
     assert shapely.get_num_coordinates(result4) < shapely.get_num_coordinates(result3)
 
 
+@pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12")
+@pytest.mark.parametrize("geometry", all_types)
+def test_coverage_simplify_scalars(geometry):
+    actual = shapely.coverage_simplify(geometry, 0.0)
+    assert isinstance(actual, Geometry)
+    assert shapely.get_type_id(actual) == shapely.get_type_id(geometry)
+    # Anything other than MultiPolygon or a GeometryCollection is returned as-is
+    if shapely.get_type_id(geometry) not in (3, 6):
+        assert actual.equals(geometry)
+
+
+@pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12")
+@pytest.mark.parametrize("geometry", all_types)
+def test_coverage_simplify_geom_types(geometry):
+    actual = shapely.coverage_simplify([geometry, geometry], 0.0)
+    assert isinstance(actual, np.ndarray)
+    assert actual.shape == (2,)
+    assert (shapely.get_type_id(actual) == shapely.get_type_id(geometry)).all()
+
+
+@pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12")
+def test_coverage_simplify_multipolygon():
+    mp = MultiPolygon(
+        [
+            Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)]),
+            Polygon([(2, 2), (2, 3), (3, 3), (3, 2), (2, 2)]),
+        ]
+    )
+    actual = shapely.coverage_simplify(mp, 1)
+    assert actual.equals(
+        shapely.from_wkt(
+            "MULTIPOLYGON (((0 1, 1 1, 1 0, 0 1)), ((2 3, 3 3, 3 2, 2 3)))"
+        )
+    )
+
+
+@pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12")
+def test_coverage_simplify_array():
+    polygons = np.array(
+        [
+            shapely.Polygon([(0, 0), (20, 0), (20, 10), (10, 5), (0, 10), (0, 0)]),
+            shapely.Polygon([(0, 10), (10, 5), (20, 10), (20, 20), (0, 20), (0, 10)]),
+        ]
+    )
+    low_tolerance = shapely.coverage_simplify(polygons, 1)
+    mid_tolerance = shapely.coverage_simplify(polygons, 8)
+    high_tolerance = shapely.coverage_simplify(polygons, 10)
+
+    assert shapely.equals(low_tolerance, shapely.normalize(polygons)).all()
+    assert shapely.equals(
+        mid_tolerance,
+        shapely.from_wkt(
+            [
+                "POLYGON ((20 10, 0 10, 0 0, 20 0, 20 10))",
+                "POLYGON ((20 10, 0 10, 0 20, 20 20, 20 10))",
+            ]
+        ),
+    ).all()
+    assert shapely.equals(
+        high_tolerance,
+        shapely.from_wkt(
+            [
+                "POLYGON ((20 10, 0 10, 20 0, 20 10))",
+                "POLYGON ((20 10, 0 10, 0 20, 20 10))",
+            ]
+        ),
+    ).all()
+
+    no_boundary = shapely.coverage_simplify(polygons, 10, simplify_boundary=False)
+    assert shapely.equals(
+        no_boundary,
+        shapely.from_wkt(
+            [
+                "POLYGON ((20 10, 0 10, 0 0, 20 0, 20 10))",
+                "POLYGON ((20 10, 0 10, 0 20, 20 20, 20 10))",
+            ]
+        ),
+    ).all()
+
+
 @pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12")
 def test_voronoi_polygons_ordered():
     mp = MultiPoint([(3.0, 1.0), (3.0, 2.0), (1.0, 2.0), (1.0, 1.0)])