Merge 8de9d6e into 38a9021

shapely/constructive.py

@@ -34,6 +34,7 @@
     "oriented_envelope",
     "minimum_rotated_rectangle",
     "minimum_bounding_circle",
+    "coverage_simplify",
 ]
 
 
@@ -91,7 +92,7 @@ def buffer(
     join_style="round",
     mitre_limit=5.0,
     single_sided=False,
-    **kwargs
+    **kwargs,
 ):
     """
     Computes the buffer of a geometry for positive and negative buffer distance.
@@ -185,7 +186,7 @@ def buffer(
         np.intc(join_style),
         mitre_limit,
         np.bool_(single_sided),
-        **kwargs
+        **kwargs,
     )
 
 
@@ -251,7 +252,7 @@ def offset_curve(
         np.intc(quad_segs),
         np.intc(join_style),
         np.double(mitre_limit),
-        **kwargs
+        **kwargs,
     )
 
 
@@ -330,7 +331,7 @@ def clip_by_rect(geometry, xmin, ymin, xmax, ymax, **kwargs):
         np.double(ymin),
         np.double(xmax),
         np.double(ymax),
-        **kwargs
+        **kwargs,
     )
 
 
@@ -861,6 +862,43 @@ 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, **kwargs):
+    """Returns a simplified version of an input geometry using the coverage simplification
+
+    Assumes that the geometry forms a polygonal coverage, i.e. that the boundaries of
+    the polygons are shared with neighboring polygons. 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.
+
+    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 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.
+    **kwargs
+        See :ref:`NumPy ufunc docs <ufuncs.kwargs>` for other keyword arguments.
+
+    """
+    return lib.coverage_simplify(geometry, tolerance, simplify_boundary, **kwargs)
+
+
 @multithreading_enabled
 def snap(geometry, reference, tolerance, **kwargs):
     """Snaps an input geometry to reference geometry's vertices.

shapely/tests/test_constructive.py

@@ -982,3 +982,66 @@ def test_concave_hull_kwargs():
     result3 = shapely.concave_hull(mp, ratio=0)
     result4 = shapely.concave_hull(mp, ratio=1)
     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_geom_types(geometry):
+    actual = shapely.coverage_simplify([geometry, geometry], 0.0)
+    assert actual.shape == (2,)
+    assert isinstance(actual[0], Geometry)
+    # 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)
+
+
+@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(
+            "GEOMETRYCOLLECTION (POLYGON ((0 1, 1 1, 1 0, 0 1)), "
+            "POLYGON ((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(
+        [
+            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)
+
+    assert low_tolerance.equals(mp.normalize())
+    assert mid_tolerance.equals(
+        shapely.from_wkt(
+            "GEOMETRYCOLLECTION (POLYGON ((20 10, 0 10, 0 0, 20 0, 20 10)), "
+            "POLYGON ((20 10, 0 10, 0 20, 20 20, 20 10)))"
+        )
+    )
+    assert high_tolerance.equals(
+        shapely.from_wkt(
+            "GEOMETRYCOLLECTION (POLYGON ((20 10, 0 10, 20 0, 20 10)), "
+            "POLYGON ((20 10, 0 10, 0 20, 20 10)))"
+        )
+    )
+
+    no_boundary = shapely.coverage_simplify(mp, 10, simplify_boundary=False)
+    assert no_boundary.equals(
+        shapely.from_wkt(
+            "GEOMETRYCOLLECTION (POLYGON ((20 10, 0 10, 0 0, 20 0, 20 10)), "
+            "POLYGON ((20 10, 0 10, 0 20, 20 20, 20 10)))"
+        )
+    )

src/ufuncs.c

@@ -2457,7 +2457,7 @@ static void points_func(char** args, const npy_intp* dimensions, const npy_intp*
 
   GEOS_INIT_THREADS;
 
-  char *ip1 = args[0];               
+  char *ip1 = args[0];
   npy_intp is1 = steps[0], cs1 = steps[3];
   npy_intp n = dimensions[0], n_c1 = dimensions[1];
   npy_intp i;
@@ -3474,6 +3474,80 @@ 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);
+    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);
+          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
@@ -3733,6 +3807,10 @@ int init_ufuncs(PyObject* m, PyObject* d) {
   DEFINE_CUSTOM(concave_hull, 3);
 #endif
 
+#if GEOS_SINCE_3_12_0
+  DEFINE_CUSTOM(coverage_simplify, 3);
+#endif
+
   Py_DECREF(ufunc);
   return 0;
 }