collect all geoms prior simplification

shapely/constructive.py

@@ -1,6 +1,6 @@
 import numpy as np
 
-from shapely import lib
+from shapely import GeometryType, lib
 from shapely._enum import ParamEnum
 from shapely.algorithms._oriented_envelope import _oriented_envelope_min_area_vectorized
 from shapely.decorators import multithreading_enabled, requires_geos
@@ -925,7 +925,7 @@ def simplify(geometry, tolerance, preserve_topology=True, **kwargs):
 
 @requires_geos("3.12.0")
 @multithreading_enabled
-def coverage_simplify(geometry, tolerance, simplify_boundary=True, **kwargs):
+def coverage_simplify(geometry, tolerance, simplify_boundary=True, axis=None, **kwargs):
     """Returns a simplified version of an input geometry using the coverage simplification
 
     Assumes that the geometry forms a polygonal coverage. Under this assumption, the
@@ -941,13 +941,13 @@ def coverage_simplify(geometry, tolerance, simplify_boundary=True, **kwargs):
     The function allows simplification of all edges including the outer boundaries of the
     coverage or simplification of only the inner (shared) edges.
 
-    Invalid polygons within the collections and geometry types other than MultiPolygons
-    or GeometryCollections of Polygons and MultiPolygons are returned unchanged.
+    If there are other geometry types than Polygons or MultiPolygons present, the
+    resulting GeometryCollection will not undergo simplification and geometries are
+    returned collected but 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
@@ -957,11 +957,25 @@ def coverage_simplify(geometry, tolerance, simplify_boundary=True, **kwargs):
     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.
+    axis : int, optional
+        Axis along which the operation is performed. The default (None)
+        performs the operation over all axes, returning a scalar value.
+        Axis may be negative, in which case it counts from the last to the
+        first axis.
     **kwargs
         See :ref:`NumPy ufunc docs <ufuncs.kwargs>` for other keyword arguments.
 
     """
-    return lib.coverage_simplify(geometry, tolerance, simplify_boundary, **kwargs)
+    geometries = np.asarray(geometry)
+    if axis is None:
+        geometries = geometries.ravel()
+    else:
+        geometries = np.rollaxis(geometries, axis=axis, start=geometries.ndim)
+
+    # create_collection acts on the inner axis
+    collections = lib.create_collection(geometries, GeometryType.GEOMETRYCOLLECTION)
+
+    return lib.coverage_simplify(collections, tolerance, simplify_boundary, **kwargs)
 
 
 @multithreading_enabled

shapely/tests/test_constructive.py

@@ -1101,11 +1101,11 @@ def test_concave_hull_kwargs():
 @pytest.mark.parametrize("geometry", all_types)
 def test_coverage_simplify_geom_types(geometry):
     actual = shapely.coverage_simplify([geometry, geometry], 0.0)
-    assert actual.shape == (2,)
-    assert isinstance(actual[0], Geometry)
+    assert isinstance(actual, Geometry)
+    assert shapely.get_type_id(actual) == 7
     # Anything other than MultiPolygon or a GeometryCollection is returned as-is
-    if shapely.get_type_id(geometry) not in (6, 7):
-        assert actual[0].equals(geometry)
+    if shapely.get_type_id(geometry) not in (3, 6):
+        assert actual.geoms[0].equals(geometry)
 
 
 @pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12")
@@ -1119,25 +1119,25 @@ def test_coverage_simplify_multipolygon():
     actual = shapely.coverage_simplify(mp, 1)
     assert actual.equals(
         shapely.from_wkt(
-            "GEOMETRYCOLLECTION (POLYGON ((0 1, 1 1, 1 0, 0 1)), "
-            "POLYGON ((2 3, 3 3, 3 2, 2 3)))"
+            "GEOMETRYCOLLECTION (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_collection():
-    mp = shapely.GeometryCollection(
+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(mp, 1)
-    mid_tolerance = shapely.coverage_simplify(mp, 8)
-    high_tolerance = shapely.coverage_simplify(mp, 10)
+    low_tolerance = shapely.coverage_simplify(polygons, 1)
+    mid_tolerance = shapely.coverage_simplify(polygons, 8)
+    high_tolerance = shapely.coverage_simplify(polygons, 10)
 
-    assert low_tolerance.equals(mp.normalize())
+    assert low_tolerance.equals(shapely.GeometryCollection(list(polygons)).normalize())
     assert mid_tolerance.equals(
         shapely.from_wkt(
             "GEOMETRYCOLLECTION (POLYGON ((20 10, 0 10, 0 0, 20 0, 20 10)), "
@@ -1151,7 +1151,7 @@ def test_coverage_simplify_collection():
         )
     )
 
-    no_boundary = shapely.coverage_simplify(mp, 10, simplify_boundary=False)
+    no_boundary = shapely.coverage_simplify(polygons, 10, simplify_boundary=False)
     assert no_boundary.equals(
         shapely.from_wkt(
             "GEOMETRYCOLLECTION (POLYGON ((20 10, 0 10, 0 0, 20 0, 20 10)), "

src/ufuncs.c

@@ -3604,37 +3604,27 @@ static void coverage_simplify_func(char** args, const npy_intp* dimensions, cons
     }
     double in2 = *(double*)ip2;
     npy_bool in3 = !(*(npy_bool*)ip3);
-    if ((in1 == NULL) || npy_isnan(in2)) {
-      // in case of a missing value: return NULL (None)
-      geom_arr[i] = NULL;
-    }
-    else if
-       (GEOSGeomTypeId_r(ctx, in1) != GEOS_MULTIPOLYGON &&
-        GEOSGeomTypeId_r(ctx, in1) != GEOS_GEOMETRYCOLLECTION) {
-        geom_arr[i] = GEOSGeom_clone_r(ctx, in1);
-    } else {
-      int isValid = 1;
-      if (GEOSGeomTypeId_r(ctx, in1) == GEOS_GEOMETRYCOLLECTION) {
-        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);
+
+    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;
-        }
-      } else {
-        geom_arr[i] = GEOSGeom_clone_r(ctx, in1);
       }
     }
+    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);
+    }
+
   }