Repository: shapely/shapely (python)

Issue:
Triangulation not working for non-convex polygon
## Expected behavior and actual behavior.

Triangulation of a *Polygon* should not produce triangles outside the Polygon. While they have to be present for a *MultiPoint* with the same coords, it is not useful to have them for the Polygon.

## Steps to reproduce the problem.

Repro data in https://github.com/jpmaterial/test_data/tree/master/Shapely/772 

## Operating system

Windows 10

## Shapely version and provenance

installed dependencies with conda, installed shapely 1.7a2 by git cloning from github and setup.py develop


Pull request description:
ENH: Add constrained_delaunay_triangles function
One of the items in https://github.com/shapely/shapely/issues/1570

Resolves #772

Patch:
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 46462d6..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.
diff --git a/src/ufuncs.c b/src/ufuncs.c
index 2f97cf2..48ebb4b 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
@@ -3992,6 +3993,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


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_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])