ENH: Add parameters method and keep_collapsed to make_valid + improve doc

CHANGES.txt

@@ -19,6 +19,7 @@ Improvements:
   signature in the ``transformation`` function.
 - The ``include_z`` in ``shapely.transform()`` now also allows ``None``, which
   lets it automatically detect the dimensionality of each input geometry.
+- Add parameters ``method`` and ``keep_collapsed`` to ``shapely.make_valid()`` (#1941)
 - Upgraded the GEOS version in the binary wheel distributions to 3.12.0.
 
 Breaking changes in GEOS 3.12.0:

shapely/constructive.py

@@ -91,7 +91,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 +185,7 @@ def buffer(
         np.intc(join_style),
         mitre_limit,
         np.bool_(single_sided),
-        **kwargs
+        **kwargs,
     )
 
 
@@ -251,7 +251,7 @@ def offset_curve(
         np.intc(quad_segs),
         np.intc(join_style),
         np.double(mitre_limit),
-        **kwargs
+        **kwargs,
     )
 
 
@@ -330,7 +330,7 @@ def clip_by_rect(geometry, xmin, ymin, xmax, ymax, **kwargs):
         np.double(ymin),
         np.double(xmax),
         np.double(ymax),
-        **kwargs
+        **kwargs,
     )
 
 
@@ -508,12 +508,22 @@ def build_area(geometry, **kwargs):
 
 
 @multithreading_enabled
-def make_valid(geometry, **kwargs):
+def make_valid(geometry, method="linework", keep_collapsed=True, **kwargs):
     """Repairs invalid geometries.
 
     Parameters
     ----------
     geometry : Geometry or array_like
+    method:  {'linework', 'structure'}, default 'linework'
+        Algorithm to use when repairing geometry. The linework algorithm combines all
+        rings into a set of noded lines and then extracts valid polygons from that
+        linework. The structure algorithm first makes all rings valid, then merges
+        shells and subtracts holes from shells to generate valid result. It assumes that
+        holes and shells are correctly categorized. 'structure' requires GEOS >= 3.10.
+    keep_collapsed: bool, default True
+        For the 'structure' method, True will keep components that have collapsed into a
+        lower dimensionality. For example, a ring collapsing to a line, or a line
+        collapsing to a point. Must be True for the 'linework' method.
     **kwargs
         See :ref:`NumPy ufunc docs <ufuncs.kwargs>` for other keyword arguments.
 
@@ -525,8 +535,40 @@ def make_valid(geometry, **kwargs):
     False
     >>> make_valid(polygon)
     <MULTILINESTRING ((0 0, 1 1), (1 1, 1 2))>
+    >>> make_valid(polygon, method="structure", keep_collapsed=True)
+    <LINESTRING (0 0, 1 1, 1 2, 1 1, 0 0)>
+    >>> make_valid(polygon, method="structure", keep_collapsed=False)
+    <POLYGON EMPTY>
     """
-    return lib.make_valid(geometry, **kwargs)
+    if not np.isscalar(method):
+        raise TypeError("method only accepts scalar values")
+    if not np.isscalar(keep_collapsed):
+        raise TypeError("keep_collapsed only accepts scalar values")
+
+    if method == "linework":
+        if keep_collapsed is False:
+            raise ValueError(
+                "The 'linework' method does not support 'keep_collapsed=False'"
+            )
+
+        # The make_valid code can be removed once support for GEOS < 3.10 is dropped.
+        # In GEOS >= 3.10, make_valid just calls make_valid_with_params with
+        # method="linework" and keep_collapsed=True, so there is no advantage to keep
+        # both code paths in shapely on long term.
+        return lib.make_valid(geometry, **kwargs)
+
+    elif method == "structure":
+        if lib.geos_version < (3, 10, 0):
+            raise ValueError(
+                "The 'structure' method is only available in GEOS >= 3.10.0"
+            )
+
+        return lib.make_valid_with_params(
+            geometry, np.intc(1), np.bool_(keep_collapsed), **kwargs
+        )
+
+    else:
+        raise ValueError(f"Unknown method: {method}")
 
 
 @multithreading_enabled

shapely/tests/test_constructive.py

@@ -244,6 +244,106 @@ def test_make_valid_1d(geom, expected):
     assert np.all(shapely.normalize(actual) == shapely.normalize(expected))
 
 
+@pytest.mark.skipif(shapely.geos_version < (3, 10, 0), reason="GEOS < 3.10")
+@pytest.mark.parametrize(
+    "geom,expected",
+    [
+        (point, point),  # a valid geometry stays the same (but is copied)
+        # an L shaped polygon without area is converted to a linestring
+        (
+            Polygon([(0, 0), (1, 1), (1, 2), (1, 1), (0, 0)]),
+            LineString([(0, 0), (1, 1), (1, 2), (1, 1), (0, 0)]),
+        ),
+        # a polygon with self-intersection (bowtie) is converted into polygons
+        (
+            Polygon([(0, 0), (2, 2), (2, 0), (0, 2), (0, 0)]),
+            MultiPolygon(
+                [
+                    Polygon([(1, 1), (2, 2), (2, 0), (1, 1)]),
+                    Polygon([(0, 0), (0, 2), (1, 1), (0, 0)]),
+                ]
+            ),
+        ),
+        (empty, empty),
+        ([empty], [empty]),
+    ],
+)
+def test_make_valid_structure(geom, expected):
+    actual = shapely.make_valid(geom, method="structure")
+    assert actual is not expected
+    # normalize needed to handle variation in output across GEOS versions
+    assert shapely.normalize(actual) == expected
+
+
+@pytest.mark.skipif(shapely.geos_version < (3, 10, 0), reason="GEOS < 3.10")
+@pytest.mark.parametrize(
+    "geom,expected",
+    [
+        (point, point),  # a valid geometry stays the same (but is copied)
+        # an L shaped polygon without area is converted to Empty Polygon
+        (
+            Polygon([(0, 0), (1, 1), (1, 2), (1, 1), (0, 0)]),
+            Polygon(),
+        ),
+        # a polygon with self-intersection (bowtie) is converted into polygons
+        (
+            Polygon([(0, 0), (2, 2), (2, 0), (0, 2), (0, 0)]),
+            MultiPolygon(
+                [
+                    Polygon([(1, 1), (2, 2), (2, 0), (1, 1)]),
+                    Polygon([(0, 0), (0, 2), (1, 1), (0, 0)]),
+                ]
+            ),
+        ),
+        (empty, empty),
+        ([empty], [empty]),
+    ],
+)
+def test_make_valid_structure_keep_collapsed_false(geom, expected):
+    actual = shapely.make_valid(geom, method="structure", keep_collapsed=False)
+    assert actual is not expected
+    # normalize needed to handle variation in output across GEOS versions
+    assert shapely.normalize(actual) == expected
+
+
+@pytest.mark.skipif(shapely.geos_version >= (3, 10, 0), reason="GEOS >= 3.10")
+def test_make_valid_structure_unsupported_geos():
+    with pytest.raises(
+        ValueError, match="The 'structure' method is only available in GEOS >= 3.10.0"
+    ):
+        _ = shapely.make_valid(Point(), method="structure")
+
+
+@pytest.mark.skipif(shapely.geos_version < (3, 10, 0), reason="GEOS < 3.10")
+@pytest.mark.parametrize(
+    "method, keep_collapsed, error_type, error",
+    [
+        (
+            np.array(["linework", "structure"]),
+            True,
+            TypeError,
+            "method only accepts scalar values",
+        ),
+        (
+            "linework",
+            [True, False],
+            TypeError,
+            "keep_collapsed only accepts scalar values",
+        ),
+        ("unknown", True, ValueError, "Unknown method: unknown"),
+        (
+            "linework",
+            False,
+            ValueError,
+            "The 'linework' method does not support 'keep_collapsed=False'",
+        ),
+    ],
+)
+def test_make_valid_invalid_params(method, keep_collapsed, error_type, error):
+    with pytest.raises(error_type, match=error):
+        _ = shapely.make_valid(Point(), method=method, keep_collapsed=keep_collapsed)
+
+
 @pytest.mark.parametrize(
     "geom,expected",
     [

src/ufuncs.c

@@ -1506,7 +1506,7 @@ static void buffer_func(char** args, const npy_intp* dimensions, const npy_intp*
   GEOS_INIT_THREADS;
 
   GEOSBufferParams* params = GEOSBufferParams_create_r(ctx);
-  if (params != 0) {
+  if (params != NULL) {
     if (!GEOSBufferParams_setQuadrantSegments_r(ctx, params, *(int*)ip3)) {
       errstate = PGERR_GEOS_EXCEPTION;
     }
@@ -1541,7 +1541,7 @@ static void buffer_func(char** args, const npy_intp* dimensions, const npy_intp*
     }
   }
 
-  if (params != 0) {
+  if (params != NULL) {
     GEOSBufferParams_destroy_r(ctx, params);
   }
 
@@ -1555,6 +1555,96 @@ static void buffer_func(char** args, const npy_intp* dimensions, const npy_intp*
 }
 static PyUFuncGenericFunction buffer_funcs[1] = {&buffer_func};
 
+#if GEOS_SINCE_3_10_0
+
+static char make_valid_with_params_inner(void* ctx, GEOSMakeValidParams* params,
+                                         void* ip1, GEOSGeometry** geom_arr,
+                                         npy_intp i) {
+  GEOSGeometry* in1 = NULL;
+
+  /* get the geometry: return on error */
+  if (!get_geom(*(GeometryObject**)ip1, &in1)) {
+    return PGERR_NOT_A_GEOMETRY;
+  }
+  /* handle NULL geometries */
+  if (in1 == NULL) {
+    geom_arr[i] = NULL;
+  } else {
+    geom_arr[i] = GEOSMakeValidWithParams_r(ctx, in1, params);
+    if (geom_arr[i] == NULL) {
+      return PGERR_GEOS_EXCEPTION;
+    }
+  }
+  return PGERR_SUCCESS;
+}
+
+static char make_valid_with_params_dtypes[4] = {NPY_OBJECT, NPY_INT, NPY_BOOL, 
+                                                NPY_OBJECT};
+static void make_valid_with_params_func(char** args, const npy_intp* dimensions,
+                                        const npy_intp* steps, void* data) {
+  char *ip1 = args[0], *ip2 = args[1], *ip3 = args[2];
+  npy_intp is1 = steps[0], is2 = steps[1], is3 = steps[2];
+  npy_intp n = dimensions[0];
+  npy_intp i;
+  GEOSGeometry** geom_arr;
+
+  CHECK_NO_INPLACE_OUTPUT(3);
+
+  if ((is2 != 0) || (is3 != 0)) {
+    PyErr_Format(PyExc_ValueError,
+                 "make_valid_with_params function called with non-scalar parameters");
+    return;
+  }
+
+  // allocate a temporary array to store output GEOSGeometry objects
+  geom_arr = malloc(sizeof(void*) * n);
+  CHECK_ALLOC(geom_arr);
+
+  GEOS_INIT_THREADS;
+
+  GEOSMakeValidParams* params = GEOSMakeValidParams_create_r(ctx);
+  if (params != NULL) {
+    if (!GEOSMakeValidParams_setMethod_r(ctx, params, *(int*)ip2)) {
+      errstate = PGERR_GEOS_EXCEPTION;
+    }
+    if (!GEOSMakeValidParams_setKeepCollapsed_r(ctx, params, *(npy_bool*)ip3)) {
+      errstate = PGERR_GEOS_EXCEPTION;
+    }
+  } else {
+    errstate = PGERR_GEOS_EXCEPTION;
+  }
+
+  if (errstate == PGERR_SUCCESS) {
+    for (i = 0; i < n; i++, ip1 += is1) {
+      CHECK_SIGNALS_THREADS(i);
+      if (errstate == PGERR_PYSIGNAL) {
+        destroy_geom_arr(ctx, geom_arr, i - 1);
+        break;
+      }
+      errstate = make_valid_with_params_inner(ctx, params, ip1, geom_arr, i);
+      if (errstate != PGERR_SUCCESS) {
+        destroy_geom_arr(ctx, geom_arr, i - 1);
+        break;
+      }
+    }
+  }
+
+  if (params != NULL) {
+    GEOSMakeValidParams_destroy_r(ctx, params);
+  }
+
+  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 make_valid_with_params_funcs[1] = {&make_valid_with_params_func};
+
+#endif  // GEOS_SINCE_3_10_0
+
 static char offset_curve_dtypes[6] = {NPY_OBJECT, NPY_DOUBLE, NPY_INT,
                                       NPY_INT,    NPY_DOUBLE, NPY_OBJECT};
 static void offset_curve_func(char** args, const npy_intp* dimensions, const npy_intp* steps,
@@ -3721,6 +3811,7 @@ int init_ufuncs(PyObject* m, PyObject* d) {
 #endif
 
 #if GEOS_SINCE_3_10_0
+  DEFINE_CUSTOM(make_valid_with_params, 3);
   DEFINE_Yd_Y(segmentize);
   DEFINE_CUSTOM(dwithin, 3);
   DEFINE_CUSTOM(from_geojson, 2);