Repository: shapely/shapely
Pull request #2248 (2025-03-03 08:42:11)

ENH - coverage validation: add coverage_invalid_edges and expose gap_width
Follow-up on https://github.com/shapely/shapely/pull/2156, adding the `gap_width` keyword to `coverage_is_valid()`, and also exposing the ability to return the invalid edges as a separate `coverage_invalid_edges()` (both use the same GEOS function under the hood, `GEOSCoverageIsValid_r`, which does both at the same time).


Still need to clean-up the implementation, add more tests, and add/update the docstring.


coveralls: ## Pull Request Test Coverage Report for [Build 13726397921](https://coveralls.io/builds/72601096)


### Details

* **9** of **9**   **(100.0%)**  changed or added relevant lines in **1** file are covered.
* No unchanged relevant lines lost coverage.
* Overall coverage increased (+**0.3%**) to **88.138%**

---



|  Totals | [![Coverage Status](https://coveralls.io/builds/72601096/badge)](https://coveralls.io/builds/72601096) |
| :-- | --: |
| Change from base [Build 13719146731](https://coveralls.io/builds/72594364): |  0.3% |
| Covered Lines: | 2712 |
| Relevant Lines: | 3077 |

---
##### 💛  - [Coveralls](https://coveralls.io)


Patch:
diff --git a/shapely/_coverage.py b/shapely/_coverage.py
index 3f9e60c..70a20d6 100644
--- a/shapely/_coverage.py
+++ b/shapely/_coverage.py
@@ -4,22 +4,32 @@ from shapely import Geometry, GeometryType, lib
 from shapely._geometry import get_parts
 from shapely.decorators import multithreading_enabled, requires_geos
 
-__all__ = ["coverage_is_valid", "coverage_simplify"]
+__all__ = ["coverage_invalid_edges", "coverage_is_valid", "coverage_simplify"]
 
 
 @requires_geos("3.12.0")
-def coverage_is_valid(geometry, **kwargs):
+@multithreading_enabled
+def coverage_is_valid(geometry, gap_width=0.0, **kwargs):
     """Verify if a coverage is valid.
 
     The coverage is represented by an array of polygonal geometries with
     exactly matching edges and no overlap.
 
+    A valid coverage may contain holes (regions of no coverage). However,
+    sometimes it might be desirable to detect narrow gaps as invalidities in
+    the coverage. The `gap_width` parameter allows to specify the maximum
+    width of gaps to detect. When gaps are detected, this function will
+    return False and the `coverage_invalid_edges` function can be used to
+    find the edges of those gaps.
+
     Geometries that are not Polygon or MultiPolygon are ignored.
 
     Parameters
     ----------
     geometry : array_like
         Array of geometries to verify.
+    gap_width : float, default 0.0
+        The maximum width of gaps to detect.
     **kwargs
         See :ref:`NumPy ufunc docs <ufuncs.kwargs>` for other keyword arguments.
 
@@ -27,11 +37,61 @@ def coverage_is_valid(geometry, **kwargs):
     -------
     bool
 
+    See Also
+    --------
+    coverage_invalid_edges, coverage_simplify
+
+    """
+    geometries = np.asarray(geometry)
+    # we always consider the full array as a single coverage -> ravel the input
+    # to pass a 1D array
+    return lib.coverage_is_valid(geometries.ravel(order="K"), gap_width, **kwargs)
+
+
+@requires_geos("3.12.0")
+@multithreading_enabled
+def coverage_invalid_edges(geometry, gap_width=0.0, **kwargs):
+    """Verify if a coverage is valid and return invalid edges.
+
+    This functions returns linear indicators showing the location of invalid
+    edges (if any) in each polygon in the input array.
+
+    The coverage is represented by an array of polygonal geometries with
+    exactly matching edges and no overlap.
+
+    A valid coverage may contain holes (regions of no coverage). However,
+    sometimes it might be desirable to detect narrow gaps as invalidities in
+    the coverage. The `gap_width` parameter allows to specify the maximum
+    width of gaps to detect. When gaps are detected, the `coverage_is_valid`
+    function will return False and this function can be used to find the
+    edges of those gaps.
+
+    Geometries that are not Polygon or MultiPolygon are ignored.
+
+    .. versionadded:: 2.1.0
+
+    Parameters
+    ----------
+    geometry : array_like
+        Array of geometries to verify.
+    gap_width : float, default 0.0
+        The maximum width of gaps to detect.
+    **kwargs
+        See :ref:`NumPy ufunc docs <ufuncs.kwargs>` for other keyword arguments.
+
+    Returns
+    -------
+    numpy.ndarray | shapely.Geometry
+
+    See Also
+    --------
+    coverage_is_valid, coverage_simplify
+
     """
     geometries = np.asarray(geometry)
     # we always consider the full array as a single coverage -> ravel the input
     # to pass a 1D array
-    return lib.coverage_is_valid(geometries.ravel(order="K"), **kwargs)
+    return lib.coverage_invalid_edges(geometries.ravel(order="K"), gap_width, **kwargs)
 
 
 @requires_geos("3.12.0")
@@ -72,6 +132,10 @@ def coverage_simplify(geometry, tolerance, *, simplify_boundary=True):
     -------
     numpy.ndarray | shapely.Geometry
 
+    See Also
+    --------
+    coverage_is_valid, coverage_invalid_edges
+
     Examples
     --------
     >>> import shapely
diff --git a/shapely/_geometry_helpers.pyx b/shapely/_geometry_helpers.pyx
index 5e2d7d1..69b9cbb 100644
--- a/shapely/_geometry_helpers.pyx
+++ b/shapely/_geometry_helpers.pyx
@@ -11,6 +11,7 @@ import shapely
 
 from shapely._geos cimport (
     GEOSContextHandle_t,
+    GEOSCoordSeq_clone_r,
     GEOSCoordSeq_getSize_r,
     GEOSCoordSequence,
     GEOSGeom_clone_r,
@@ -21,6 +22,7 @@ from shapely._geos cimport (
     GEOSGeom_createPoint_r,
     GEOSGeom_createPolygon_r,
     GEOSGeom_destroy_r,
+    GEOSGeom_getCoordSeq_r,
     GEOSGeometry,
     GEOSGeomTypeId_r,
     GEOSGetExteriorRing_r,
@@ -686,3 +688,43 @@ def _geom_factory(uintptr_t g):
         geom = PyGEOS_CreateGeometry(<GEOSGeometry *>g, geos_handle)
 
     return geom
+
+
+def linestring_to_linearring(object line):
+    cdef GEOSGeometry *geom = NULL
+    cdef const GEOSCoordSequence *seq = NULL
+    cdef GEOSCoordSequence *seq_cloned = NULL
+    cdef GEOSGeometry *ring = NULL
+    cdef int geom_type
+
+    with get_geos_handle() as geos_handle:
+        if PyGEOS_GetGEOSGeometry(<PyObject*> line, &geom) == 0:
+            raise TypeError(
+                "The argument is of incorrect type. Please provide a Geometry object."
+            )
+        if geom == NULL:
+            raise TypeError(
+                "The argument is of incorrect type. Please provide a Geometry object."
+            )
+
+        geom_type = GEOSGeomTypeId_r(geos_handle, geom)
+        if geom_type != 1:
+            raise TypeError(
+                "The argument is of incorrect type. Please provide a LineString object."
+            )
+
+        seq = GEOSGeom_getCoordSeq_r(geos_handle, geom)
+        if seq == NULL:
+            return
+
+        seq_cloned = GEOSCoordSeq_clone_r(geos_handle, seq)
+        if seq_cloned == NULL:
+            return
+
+        ring = GEOSGeom_createLinearRing_r(geos_handle, seq_cloned)
+        if ring == NULL:
+            return
+
+        result = PyGEOS_CreateGeometry(ring, geos_handle)
+
+    return result
diff --git a/shapely/_geos.pxd b/shapely/_geos.pxd
index 30ad384..6a7eb68 100644
--- a/shapely/_geos.pxd
+++ b/shapely/_geos.pxd
@@ -40,6 +40,8 @@ cdef extern from "geos_c.h":
     void GEOSGeom_destroy_r(GEOSContextHandle_t handle, GEOSGeometry* g) nogil
 
     # Coordinate sequences
+    const GEOSCoordSequence* GEOSGeom_getCoordSeq_r(GEOSContextHandle_t handle, const GEOSGeometry* g)
+    GEOSCoordSequence* GEOSCoordSeq_clone_r(GEOSContextHandle_t handle, const GEOSCoordSequence* s)
     GEOSCoordSequence* GEOSCoordSeq_create_r(GEOSContextHandle_t handle, unsigned int size, unsigned int dims) nogil
     void GEOSCoordSeq_destroy_r(GEOSContextHandle_t handle, GEOSCoordSequence* s) nogil
     int GEOSCoordSeq_setX_r(GEOSContextHandle_t handle, GEOSCoordSequence* s, unsigned int idx, double val) nogil
diff --git a/shapely/_ragged_array.py b/shapely/_ragged_array.py
index 73ff162..df157c0 100644
--- a/shapely/_ragged_array.py
+++ b/shapely/_ragged_array.py
@@ -347,6 +347,8 @@ def _point_from_flatcoords(coords):
 
 def _multipoint_from_flatcoords(coords, offsets):
     # recreate points
+    if len(offsets):
+        coords = coords[offsets[0] :]
     points = creation.points(coords)
 
     # recreate multipoints
@@ -364,6 +366,8 @@ def _multipoint_from_flatcoords(coords, offsets):
 
 def _linestring_from_flatcoords(coords, offsets):
     # recreate linestrings
+    if len(offsets):
+        coords = coords[offsets[0] :]
     linestring_n = np.diff(offsets)
     linestring_indices = np.repeat(np.arange(len(linestring_n)), linestring_n)
 
diff --git a/shapely/constructive.py b/shapely/constructive.py
index abce8fd..bf97bb6 100644
--- a/shapely/constructive.py
+++ b/shapely/constructive.py
@@ -22,6 +22,7 @@ __all__ = [
     "centroid",
     "clip_by_rect",
     "concave_hull",
+    "constrained_delaunay_triangles",
     "convex_hull",
     "delaunay_triangles",
     "envelope",
@@ -490,8 +491,8 @@ def delaunay_triangles(geometry, tolerance=0.0, only_edges=False, **kwargs):
     """Compute a Delaunay triangulation around the vertices of an input geometry.
 
     The output is a geometrycollection containing polygons (default)
-    or linestrings (see only_edges). Returns an None if an input geometry
-    contains less than 3 vertices.
+    or linestrings (see ``only_edges``). Returns an empty geometry for input
+    geometries that contain less than 3 vertices.
 
     Parameters
     ----------
@@ -505,6 +506,14 @@ def delaunay_triangles(geometry, tolerance=0.0, only_edges=False, **kwargs):
     **kwargs
         See :ref:`NumPy ufunc docs <ufuncs.kwargs>` for other keyword arguments.
 
+    Returns
+    -------
+    GeometryCollection or array of GeometryCollections
+
+    See Also
+    --------
+    constrained_delaunay_triangles
+
     Examples
     --------
     >>> import shapely
@@ -532,6 +541,54 @@ def delaunay_triangles(geometry, tolerance=0.0, only_edges=False, **kwargs):
     return lib.delaunay_triangles(geometry, tolerance, only_edges, **kwargs)
 
 
+@requires_geos("3.10.0")
+@multithreading_enabled
+def constrained_delaunay_triangles(geometry, **kwargs):
+    """Compute the constrained Delaunay triangulation of polygons.
+
+    A constrained Delaunay triangulation requires the edges of the input
+    polygon(s) to be in the set of resulting triangle edges. An unconstrained
+    delaunay triangulation only triangulates based on the vertices, hence
+    triangle edges could cross polygon boundaries.
+
+    .. versionadded:: 2.1.0
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Returns
+    -------
+    GeometryCollection or array of GeometryCollections
+        * GeometryCollection of polygons, given polygonal input
+        * Empty GeometryCollection, given non-polygonal input
+
+    See Also
+    --------
+    delaunay_triangles
+
+    Examples
+    --------
+    >>> import shapely
+    >>> from shapely import MultiPoint, MultiPolygon, Polygon
+    >>> shapely.constrained_delaunay_triangles(Polygon([(10, 10), (20, 40), (90, 90), (90, 10), (10, 10)]))
+    <GEOMETRYCOLLECTION (POLYGON ((90 10, 20 40, 90 90, 90 10)), POLYGON ((20 40...>
+    >>> shapely.constrained_delaunay_triangles(Polygon())
+    <GEOMETRYCOLLECTION EMPTY>
+    >>> shapely.constrained_delaunay_triangles(MultiPolygon([Polygon(((50, 30), (60, 30), (100, 100), (50, 30))), Polygon(((10, 10), (20, 40), (90, 90), (90, 10), (10, 10)))]))
+    <GEOMETRYCOLLECTION (POLYGON ((50 30, 100 100, 60 30, 50 30)), POLYGON ((90 ...>
+    >>> shapely.constrained_delaunay_triangles(MultiPolygon())
+    <GEOMETRYCOLLECTION EMPTY>
+    >>> shapely.constrained_delaunay_triangles(MultiPoint([(50, 30), (51, 30), (60, 30), (100, 100)]))
+    <GEOMETRYCOLLECTION EMPTY>
+
+    """  # noqa: E501
+    return lib.constrained_delaunay_triangles(geometry, **kwargs)
+
+
 @multithreading_enabled
 def envelope(geometry, **kwargs):
     """Compute the minimum bounding box that encloses an input geometry.
@@ -1124,13 +1181,13 @@ def simplify(geometry, tolerance, preserve_topology=True, **kwargs):
 
 @multithreading_enabled
 def snap(geometry, reference, tolerance, **kwargs):
-    """Snap an input geometry to reference geometry's vertices.
+    """Snap the vertices and segments of the geometry to vertices of the reference.
 
-    Vertices of the first geometry are snapped to vertices of the second.
-    geometry, returning a new geometry; the input geometries are not modified.
-    The result geometry is the input geometry with the vertices snapped.
-    If no snapping occurs then the input geometry is returned unchanged.
-    The tolerance is used to control where snapping is performed.
+    Vertices and segments of the input geometry are snapped to vertices of the
+    reference geometry, returning a new geometry; the input geometries are not
+    modified. The result geometry is the input geometry with the vertices and
+    segments snapped. If no snapping occurs then the input geometry is returned
+    unchanged. The tolerance is used to control where snapping is performed.
 
     Where possible, this operation tries to avoid creating invalid geometries;
     however, it does not guarantee that output geometries will be valid. It is
diff --git a/shapely/geometry/polygon.py b/shapely/geometry/polygon.py
index f62621a..c7793e7 100644
--- a/shapely/geometry/polygon.py
+++ b/shapely/geometry/polygon.py
@@ -3,6 +3,7 @@
 import numpy as np
 
 import shapely
+from shapely import _geometry_helpers
 from shapely.algorithms.cga import signed_area  # noqa
 from shapely.errors import TopologicalError
 from shapely.geometry.base import BaseGeometry
@@ -15,7 +16,7 @@ __all__ = ["LinearRing", "Polygon", "orient"]
 def _unpickle_linearring(wkb):
     linestring = shapely.from_wkb(wkb)
     srid = shapely.get_srid(linestring)
-    linearring = shapely.linearrings(shapely.get_coordinates(linestring))
+    linearring = _geometry_helpers.linestring_to_linearring(linestring)
     if srid:
         linearring = shapely.set_srid(linearring, srid)
     return linearring
diff --git a/src/ufuncs.c b/src/ufuncs.c
index 2f97cf2..e1f6ccc 100644
--- a/src/ufuncs.c
+++ b/src/ufuncs.c
@@ -737,6 +737,7 @@ static void* maximum_inscribed_circle_data[1] = {GEOSMaximumInscribedCircleWithD
 
 #if GEOS_SINCE_3_10_0
 static void* segmentize_data[1] = {GEOSDensify_r};
+static void* constrained_delaunay_triangles_data[1] = {GEOSConstrainedDelaunayTriangulation_r};
 #endif
 
 #if GEOS_SINCE_3_11_0
@@ -2445,7 +2446,7 @@ static PyUFuncGenericFunction polygonize_full_funcs[1] = {&polygonize_full_func}
 
 #if GEOS_SINCE_3_12_0
 
-static char coverage_is_valid_dtypes[2] = {NPY_OBJECT, NPY_BOOL};
+static char coverage_is_valid_dtypes[3] = {NPY_OBJECT, NPY_DOUBLE, NPY_BOOL};
 static void coverage_is_valid_func(char** args, const npy_intp* dimensions, const npy_intp* steps,
                                    void* data) {
   GEOSGeometry* geom = NULL;
@@ -2455,6 +2456,12 @@ static void coverage_is_valid_func(char** args, const npy_intp* dimensions, cons
   unsigned int n_geoms;
   int ret;
 
+  if (steps[1] != 0) {
+    PyErr_Format(PyExc_ValueError, "coverage_is_valid function called with non-scalar gap_width");
+    return;
+  }
+  double gap_width = *(double*)args[1];
+
   GEOS_INIT_THREADS;
 
   GEOSGeometry** geoms = malloc(sizeof(void*) * dimensions[1]);
@@ -2463,13 +2470,14 @@ static void coverage_is_valid_func(char** args, const npy_intp* dimensions, cons
     goto finish;
   }
 
-  SINGLE_COREDIM_LOOP_OUTER {
+  BINARY_SINGLE_COREDIM_LOOP_OUTER {
     CHECK_SIGNALS(i);
     if (errstate == PGERR_PYSIGNAL) {
       goto finish;
     }
+    cp1 = ip1;
     n_geoms = 0;
-    SINGLE_COREDIM_LOOP_INNER {
+    BINARY_SINGLE_COREDIM_LOOP_INNER {
       if (!get_geom(*(GeometryObject**)cp1, &geom)) {
         errstate = PGERR_NOT_A_GEOMETRY;
         goto finish;
@@ -2488,7 +2496,7 @@ static void coverage_is_valid_func(char** args, const npy_intp* dimensions, cons
       goto finish;
     }
 
-    ret = GEOSCoverageIsValid_r(ctx, collection, 0, NULL);
+    ret = GEOSCoverageIsValid_r(ctx, collection, gap_width, NULL);
     if (ret == 2) {
       errstate = PGERR_GEOS_EXCEPTION;
       goto finish;
@@ -2513,6 +2521,96 @@ finish:
 }
 static PyUFuncGenericFunction coverage_is_valid_funcs[1] = {&coverage_is_valid_func};
 
+
+static char coverage_invalid_edges_dtypes[3] = {NPY_OBJECT, NPY_DOUBLE, NPY_OBJECT};
+static void coverage_invalid_edges_func(char** args, const npy_intp* dimensions, const npy_intp* steps,
+                                   void* data) {
+  GEOSGeometry* geom = NULL;
+  GEOSGeometry* collection = NULL;
+  GEOSGeometry** collection_parts;
+  GEOSGeometry* result_collection = NULL;
+  GEOSGeometry** result_collection_parts;
+  unsigned int n_parts, n_parts_result;
+  unsigned int n_geoms;
+  int ret;
+
+  if (steps[1] != 0) {
+    PyErr_Format(PyExc_ValueError, "coverage_invalid_edges function called with non-scalar gap_width");
+    return;
+  }
+  double gap_width = *(double*)args[1];
+
+  // allocate a temporary array to store input GEOSGeometry objects
+  GEOSGeometry** geoms = malloc(sizeof(void*) * dimensions[1]);
+  CHECK_ALLOC(geoms);
+
+  GEOS_INIT;
+
+  npy_intp ocs1 = steps[4];
+  BINARY_SINGLE_COREDIM_LOOP_OUTER {
+    Py_BEGIN_ALLOW_THREADS;
+    CHECK_SIGNALS(i);
+    if (errstate == PGERR_PYSIGNAL) {
+      goto finish;
+    }
+    cp1 = ip1;
+    n_geoms = 0;
+    BINARY_SINGLE_COREDIM_LOOP_INNER {
+      if (!get_geom(*(GeometryObject**)cp1, &geom)) {
+        errstate = PGERR_NOT_A_GEOMETRY;
+        goto finish;
+      }
+      if (geom == NULL) {
+        continue;
+      }
+      // we do not clone the geometries, so have to release the collection later
+      geoms[n_geoms] = geom;
+      n_geoms++;
+    }
+    collection =
+        GEOSGeom_createCollection_r(ctx, GEOS_GEOMETRYCOLLECTION, geoms, n_geoms);
+    if (collection == NULL) {
+      errstate = PGERR_GEOS_EXCEPTION;
+      goto finish;
+    }
+
+    ret = GEOSCoverageIsValid_r(ctx, collection, gap_width, &result_collection);
+    if ((ret == 2) || (result_collection == NULL)) {
+      errstate = PGERR_GEOS_EXCEPTION;
+      goto finish;
+    }
+
+    Py_END_ALLOW_THREADS;
+    result_collection_parts = GEOSGeom_releaseCollection_r(ctx, result_collection, &n_parts_result);
+    geom_arr_to_npy(result_collection_parts, op1, ocs1, n_parts_result);
+    GEOSFree_r(ctx, result_collection_parts);
+    GEOSGeom_destroy_r(ctx, result_collection);
+    result_collection = NULL;
+
+    collection_parts = GEOSGeom_releaseCollection_r(ctx, collection, &n_parts);
+    GEOSFree_r(ctx, collection_parts);
+    GEOSGeom_destroy_r(ctx, collection);
+    collection = NULL;
+  }
+
+finish:
+  if (collection != NULL) {
+    collection_parts = GEOSGeom_releaseCollection_r(ctx, collection, &n_parts);
+    GEOSFree_r(ctx, collection_parts);
+    GEOSGeom_destroy_r(ctx, collection);
+  }
+  if (result_collection != NULL) {
+    result_collection_parts = GEOSGeom_releaseCollection_r(ctx, result_collection, &n_parts_result);
+    GEOSFree_r(ctx, result_collection_parts);
+    GEOSGeom_destroy_r(ctx, result_collection);
+  }
+  if (geoms != NULL) {
+    free(geoms);
+  }
+  GEOS_FINISH;
+}
+static PyUFuncGenericFunction coverage_invalid_edges_funcs[1] = {&coverage_invalid_edges_func};
+
 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) {
@@ -2546,7 +2644,7 @@ static void coverage_simplify_func(char** args, const npy_intp* dimensions, cons
     // Validate the geometries in the collection
     int num_geoms = GEOSGetNumGeometries_r(ctx, in1);
     for (int j = 0; j < num_geoms; j++) {
-      GEOSGeometry* geom = GEOSGetGeometryN_r(ctx, in1, j);
+      const GEOSGeometry* geom = GEOSGetGeometryN_r(ctx, in1, j);
       geom_type = GEOSGeomTypeId_r(ctx, geom);
       if (geom_type != GEOS_POLYGON && geom_type != GEOS_MULTIPOLYGON) {
         errstate = PGERR_GEOMETRY_TYPE;
@@ -2574,7 +2672,7 @@ finish:
 }
 static PyUFuncGenericFunction coverage_simplify_funcs[1] = {&coverage_simplify_func};
 
-#endif
+#endif  // GEOS_SINCE_3_12_0
 
 static char shortest_line_dtypes[3] = {NPY_OBJECT, NPY_OBJECT, NPY_OBJECT};
 static void shortest_line_func(char** args, const npy_intp* dimensions, const npy_intp* steps,
@@ -3992,6 +4090,7 @@ int init_ufuncs(PyObject* m, PyObject* d) {
   DEFINE_CUSTOM(dwithin, 3);
   DEFINE_CUSTOM(from_geojson, 2);
   DEFINE_CUSTOM(to_geojson, 2);
+  DEFINE_Y_Y(constrained_delaunay_triangles);
 #endif
 
 #if GEOS_SINCE_3_11_0
@@ -4001,7 +4100,8 @@ int init_ufuncs(PyObject* m, PyObject* d) {
 #endif
 
 #if GEOS_SINCE_3_12_0
-  DEFINE_GENERALIZED(coverage_is_valid, 1, "(d)->()");
+  DEFINE_GENERALIZED(coverage_is_valid, 2, "(d),()->()");
+  DEFINE_GENERALIZED(coverage_invalid_edges, 2, "(d),()->(d)");
   DEFINE_CUSTOM(coverage_simplify, 3);
   DEFINE_Y_Y(disjoint_subset_union);
   DEFINE_Y_b(has_m);


Test patch:
diff --git a/shapely/tests/test_constructive.py b/shapely/tests/test_constructive.py
index 2a7da71..1b3f495 100644
--- a/shapely/tests/test_constructive.py
+++ b/shapely/tests/test_constructive.py
@@ -51,6 +51,12 @@ CONSTRUCTIVE_NO_ARGS = (
     shapely.node,
     shapely.normalize,
     shapely.point_on_surface,
+    pytest.param(
+        shapely.constrained_delaunay_triangles,
+        marks=pytest.mark.skipif(
+            shapely.geos_version < (3, 10, 0), reason="GEOS < 3.10"
+        ),
+    ),
 )
 
 CONSTRUCTIVE_FLOAT_ARG = (
@@ -1109,6 +1115,44 @@ def test_concave_hull_kwargs():
     assert shapely.get_num_coordinates(result4) < shapely.get_num_coordinates(result3)
 
 
+@pytest.mark.skipif(shapely.geos_version < (3, 10, 0), reason="GEOS < 3.10")
+class TestConstrainedDelaunayTriangulation:
+    """
+    Only testing the number of triangles and their type here.
+    This doesn't actually test the points in the resulting geometries.
+
+    """
+
+    def test_poly(self):
+        polys = shapely.constrained_delaunay_triangles(
+            Polygon([(10, 10), (20, 40), (90, 90), (90, 10), (10, 10)])
+        )
+        assert len(polys.geoms) == 2
+        for p in polys.geoms:
+            assert isinstance(p, Polygon)
+
+    def test_multi_polygon(self):
+        multipoly = MultiPolygon(
+            [
+                Polygon(((50, 30), (60, 30), (100, 100), (50, 30))),
+                Polygon(((10, 10), (20, 40), (90, 90), (90, 10), (10, 10))),
+            ]
+        )
+        polys = shapely.constrained_delaunay_triangles(multipoly)
+        assert len(polys.geoms) == 3
+        for p in polys.geoms:
+            assert isinstance(p, Polygon)
+
+    def test_point(self):
+        p = Point(1, 1)
+        polys = shapely.constrained_delaunay_triangles(p)
+        assert len(polys.geoms) == 0
+
+    def test_empty_poly(self):
+        polys = shapely.constrained_delaunay_triangles(Polygon())
+        assert len(polys.geoms) == 0
+
+
 @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)])
diff --git a/shapely/tests/test_coverage.py b/shapely/tests/test_coverage.py
index 3ed2c4f..db6f370 100644
--- a/shapely/tests/test_coverage.py
+++ b/shapely/tests/test_coverage.py
@@ -9,9 +9,11 @@ from shapely import (
     Polygon,
 )
 from shapely.errors import UnsupportedGEOSVersionError
+from shapely.testing import assert_geometries_equal
 from shapely.tests.common import (
     all_types,
     all_types_z,
+    empty_line_string,
 )
 
 
@@ -23,6 +25,10 @@ def test_coverage_is_valid(geometry):
     assert actual.dtype == np.bool_
     assert actual.item() is True
 
+    actual = shapely.coverage_invalid_edges([geometry])
+    expected = np.array([empty_line_string], dtype=object)
+    assert_geometries_equal(actual, expected)
+
 
 @pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="requires >= 3.12")
 def test_coverage_is_valid_non_polygonal():
@@ -34,31 +40,111 @@ def test_coverage_is_valid_non_polygonal():
         Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)]),
     ]
     assert shapely.coverage_is_valid(geoms)
+    assert (shapely.coverage_invalid_edges(geoms) == empty_line_string).all()
 
 
 @pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="requires >= 3.12")
 def test_coverage_is_valid_polygonal():
-    geoms = [
-        Polygon([(0, 0), (1, 1), (1, 0), (0, 0)]),
-        Polygon([(0, 0), (1, 1), (0, 1), (0, 0)]),
-    ]
-    assert shapely.coverage_is_valid(geoms)
+    # adjacent triangles
+    poly1 = Polygon([(0, 0), (1, 1), (1, 0), (0, 0)])
+    poly2 = Polygon([(0, 0), (1, 1), (0, 1), (0, 0)])
+    assert shapely.coverage_is_valid([poly1, poly2])
+    assert shapely.is_empty(shapely.coverage_invalid_edges([poly1, poly2])).all()
+
+    # shared egde but without identical vertices
+    poly2b = Polygon([(0, 0), (0.5, 0.5), (1, 1), (0, 1), (0, 0)])
+    assert not shapely.coverage_is_valid([poly1, poly2b])
+    result = shapely.coverage_invalid_edges([poly1, poly2b])
+    expected = [LineString([(0, 0), (1, 1)]), LineString([(0, 0), (0.5, 0.5), (1, 1)])]
+    assert_geometries_equal(result, expected)
+
+    # overlap
+    poly3 = Polygon([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])
+    assert not shapely.coverage_is_valid([poly1, poly3])
+    result = shapely.coverage_invalid_edges([poly1, poly3])
+    assert not shapely.is_empty(result).any()
 
-    geoms = [
-        Polygon([(0, 0), (1, 1), (1, 0), (0, 0)]),
-        Polygon([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]),
-    ]
-    assert not shapely.coverage_is_valid(geoms)
 
+@pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="requires >= 3.12")
+def test_coverage_is_valid_gap_width():
+    # shared edge of boxes with multiple vertices
+    poly1 = shapely.from_wkt("POLYGON ((0 10, 10 10, 10 7, 10 3, 10 0, 0 0, 0 10))")
+    poly2 = shapely.from_wkt("POLYGON ((10 10, 20 10, 20 0, 10 0, 10 3, 10 7, 10 10))")
+
+    # extra vertex in middle of shared edge
+    poly2_extra = shapely.from_wkt(
+        "POLYGON ((10 10, 20 10, 20 0, 10 0, 10 3, 10 5, 10 7, 10 10))"
+    )
+    # extra vertex shifted -> gap of max 1 wide
+    poly2_shift = shapely.from_wkt(
+        "POLYGON ((10 10, 20 10, 20 0, 10 0, 10 3, 11 5, 10 7, 10 10))"
+    )
+
+    # valid coverage -> gap_width value does not matter
+    assert shapely.coverage_is_valid([poly1, poly2], gap_width=0.0)
+    assert shapely.coverage_is_valid([poly1, poly2], gap_width=2.0)
 
-@pytest.mark.skipif(shapely.geos_version >= (3, 12, 0), reason="requires >= 3.12")
-def test_coverage_is_valid_unsupported_geos():
-    geoms = [
-        Polygon([(0, 0), (1, 1), (1, 0), (0, 0)]),
-        Polygon([(0, 0), (1, 1), (0, 1), (0, 0)]),
-    ]
-    with pytest.raises(UnsupportedGEOSVersionError):
-        shapely.coverage_is_valid(geoms)
+    result = shapely.coverage_invalid_edges([poly1, poly2], gap_width=0.0)
+    assert_geometries_equal(result, [empty_line_string] * 2)
+    result = shapely.coverage_invalid_edges([poly1, poly2], gap_width=2.0)
+    assert_geometries_equal(result, [empty_line_string] * 2)
+
+    # invalid coverage -> gap_width value does not matter
+    assert not shapely.coverage_is_valid([poly1, poly2_extra], gap_width=0.0)
+    assert not shapely.coverage_is_valid([poly1, poly2_extra], gap_width=2.0)
+
+    expected = shapely.from_wkt(
+        ["LINESTRING (10 7, 10 3)", "LINESTRING (10 3, 10 5, 10 7)"]
+    )
+    result = shapely.coverage_invalid_edges([poly1, poly2_extra], gap_width=0.0)
+    assert_geometries_equal(result, expected)
+    result = shapely.coverage_invalid_edges([poly1, poly2_extra], gap_width=2.0)
+    assert_geometries_equal(result, expected)
+
+    # coverage with gap of 1 unit wide
+    assert shapely.coverage_is_valid([poly1, poly2_shift], gap_width=0.0)
+    assert shapely.coverage_is_valid([poly1, poly2_shift], gap_width=0.5)
+    assert not shapely.coverage_is_valid([poly1, poly2_shift], gap_width=1.0)
+    assert not shapely.coverage_is_valid([poly1, poly2_shift], gap_width=1.5)
+    # TODO why this behaviour?
+    assert shapely.coverage_is_valid([poly1, poly2_shift], gap_width=2.0)
+
+    assert_geometries_equal(
+        shapely.coverage_invalid_edges([poly1, poly2_shift], gap_width=0.0),
+        [empty_line_string] * 2,
+    )
+    assert_geometries_equal(
+        shapely.coverage_invalid_edges([poly1, poly2_shift], gap_width=1.0),
+        shapely.from_wkt(["LINESTRING (10 7, 10 3)", "LINESTRING (10 3, 11 5, 10 7)"]),
+    )
+
+
+@pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="requires >= 3.12")
+def test_coverage_invalid_edges_gufunc():
+    poly1 = shapely.from_wkt("POLYGON ((0 10, 10 10, 10 7, 10 3, 10 0, 0 0, 0 10))")
+    poly2 = shapely.from_wkt("POLYGON ((10 10, 20 10, 20 0, 10 0, 10 3, 10 7, 10 10))")
+    poly2_extra = shapely.from_wkt(
+        "POLYGON ((10 10, 20 10, 20 0, 10 0, 10 3, 10 5, 10 7, 10 10))"
+    )
+    poly3 = shapely.from_wkt("POLYGON ((20 10, 30 10, 30 7, 30 3, 30 0, 20 0, 20 10))")
+
+    arr = np.array([[poly1, poly2, poly3], [poly1, poly2_extra, poly3]])
+    result = shapely.lib.coverage_invalid_edges(arr, 0.0)
+    expected = shapely.from_wkt(
+        [
+            ["LINESTRING EMPTY"] * 3,
+            [
+                "LINESTRING (10 7, 10 3)",
+                "LINESTRING (10 3, 10 5, 10 7)",
+                "LINESTRING EMPTY",
+            ],
+        ]
+    )
+    assert_geometries_equal(result, expected)
+
+    arr2 = np.array(arr, order="F")
+    result = shapely.lib.coverage_invalid_edges(arr2, 0.0)
+    assert_geometries_equal(result, expected)
 
 
 @pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12")
@@ -145,3 +231,19 @@ def test_coverage_simplify_array():
             ]
         ),
     ).all()
+
+
+@pytest.mark.skipif(shapely.geos_version >= (3, 12, 0), reason="requires >= 3.12")
+def test_coverage_unsupported_geos():
+    geoms = [
+        Polygon([(0, 0), (1, 1), (1, 0), (0, 0)]),
+        Polygon([(0, 0), (1, 1), (0, 1), (0, 0)]),
+    ]
+    with pytest.raises(UnsupportedGEOSVersionError):
+        shapely.coverage_is_valid(geoms)
+
+    with pytest.raises(UnsupportedGEOSVersionError):
+        shapely.coverage_invalid_edges(geoms)
+
+    with pytest.raises(UnsupportedGEOSVersionError):
+        shapely.coverage_simplify(geoms, 1.0)
diff --git a/shapely/tests/test_io.py b/shapely/tests/test_io.py
index 1b403cb..f41d61e 100644
--- a/shapely/tests/test_io.py
+++ b/shapely/tests/test_io.py
@@ -1153,9 +1153,6 @@ def test_pickle(geom):
 
 @pytest.mark.parametrize("geom", all_types_z)
 def test_pickle_z(geom):
-    if shapely.get_type_id(geom) == shapely.GeometryType.LINEARRING:
-        # https://github.com/shapely/shapely/pull/1808#discussion_r1182174077
-        pytest.xfail("TODO: Z gets dropped for LinearRing")
     pickled = pickle.dumps(geom)
     actual = pickle.loads(pickled)
     assert_geometries_equal(actual, geom, tolerance=0)
@@ -1171,8 +1168,6 @@ def test_pickle_z(geom):
 )
 @pytest.mark.parametrize("geom", all_types_m)
 def test_pickle_m(geom):
-    if shapely.get_type_id(geom) == shapely.GeometryType.LINEARRING:
-        pytest.xfail("TODO: M gets dropped for LinearRing types; see GH-2005")
     pickled = pickle.dumps(geom)
     actual = pickle.loads(pickled)
     assert_geometries_equal(actual, geom, tolerance=0)
@@ -1187,8 +1182,6 @@ def test_pickle_m(geom):
 )
 @pytest.mark.parametrize("geom", all_types_zm)
 def test_pickle_zm(geom):
-    if shapely.get_type_id(geom) == shapely.GeometryType.LINEARRING:
-        pytest.xfail("TODO: ZM gets dropped for LinearRing types; see GH-2005")
     pickled = pickle.dumps(geom)
     actual = pickle.loads(pickled)
     assert_geometries_equal(actual, geom, tolerance=0)
diff --git a/shapely/tests/test_ragged_array.py b/shapely/tests/test_ragged_array.py
index 0e4d5df..32bc9e8 100644
--- a/shapely/tests/test_ragged_array.py
+++ b/shapely/tests/test_ragged_array.py
@@ -246,6 +246,15 @@ def test_linestrings():
     arr[-2] = shapely.from_wkt("LINESTRING EMPTY")
     assert_geometries_equal(result, arr)
 
+    # sliced
+    offsets_sliced = (offsets[0][1:],)
+    result = shapely.from_ragged_array(typ, coords, offsets_sliced)
+    assert_geometries_equal(result, arr[1:])
+
+    offsets_sliced = (offsets[0][:-1],)
+    result = shapely.from_ragged_array(typ, coords, offsets_sliced)
+    assert_geometries_equal(result, arr[:-1])
+
 
 def test_polygons():
     arr = shapely.from_wkt(
@@ -299,6 +308,17 @@ def test_polygons():
     arr[-2] = shapely.from_wkt("POLYGON EMPTY")
     assert_geometries_equal(result, arr)
 
+    # sliced:
+    # - indices into the coordinate array for the whole buffer
+    # - indices into the ring array for *just* the sliced part
+    offsets_sliced = (offsets[0], offsets[1][1:])
+    result = shapely.from_ragged_array(typ, coords, offsets_sliced)
+    assert_geometries_equal(result, arr[1:])
+
+    offsets_sliced = (offsets[0], offsets[1][:-1])
+    result = shapely.from_ragged_array(typ, coords, offsets_sliced)
+    assert_geometries_equal(result, arr[:-1])
+
 
 def test_multipoints():
     arr = shapely.from_wkt(
@@ -338,6 +358,15 @@ def test_multipoints():
     arr[-2] = shapely.from_wkt("MULTIPOINT EMPTY")
     assert_geometries_equal(result, arr)
 
+    # sliced:
+    offsets_sliced = (offsets[0][1:],)
+    result = shapely.from_ragged_array(typ, coords, offsets_sliced)
+    assert_geometries_equal(result, arr[1:])
+
+    offsets_sliced = (offsets[0][:-1],)
+    result = shapely.from_ragged_array(typ, coords, offsets_sliced)
+    assert_geometries_equal(result, arr[:-1])
+
 
 def test_multilinestrings():
     arr = shapely.from_wkt(
@@ -389,6 +418,17 @@ def test_multilinestrings():
     arr[-2] = shapely.from_wkt("MULTILINESTRING EMPTY")
     assert_geometries_equal(result, arr)
 
+    # sliced:
+    # - indices into the coordinate array for the whole buffer
+    # - indices into the line parts for *just* the sliced part
+    offsets_sliced = (offsets[0], offsets[1][1:])
+    result = shapely.from_ragged_array(typ, coords, offsets_sliced)
+    assert_geometries_equal(result, arr[1:])
+
+    offsets_sliced = (offsets[0], offsets[1][:-1])
+    result = shapely.from_ragged_array(typ, coords, offsets_sliced)
+    assert_geometries_equal(result, arr[:-1])
+
 
 def test_multipolygons():
     arr = shapely.from_wkt(
@@ -453,6 +493,16 @@ def test_multipolygons():
     arr[-2] = shapely.from_wkt("MULTIPOLYGON EMPTY")
     assert_geometries_equal(result, arr)
 
+    # sliced:
+    offsets_sliced = (offsets[0], offsets[1], offsets[2][1:])
+    result = shapely.from_ragged_array(typ, coords, offsets_sliced)
+    assert_geometries_equal(result, arr[1:])
+
+    offsets_sliced = (offsets[0], offsets[1], offsets[2][:-3])
+    result = shapely.from_ragged_array(typ, coords, offsets_sliced)
+    assert_geometries_equal(result, arr[:-3])
+    print(result)
+
 
 def test_mixture_point_multipoint():
     typ, coords, offsets = shapely.to_ragged_array([point, multi_point])