Repository: geopandas/geopandas
Pull request #3552 (2025-04-23 20:08:38)

ENH: add constrained_delaunay_triangles
Add a `constrained_delaunay_triangles` method on GeoSeries/GeoDataFrame as this function became available in shapely 2.1.

xref #2010


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index a797e06d..2c2361cd 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -15,12 +15,15 @@ New features and improvements:
 - Improve performance of `overlay` with `how=identity` (#3504).
 - Fix ambiguous error when GeoDataFrame is initialised with a column called "crs" (#3502).
 - Added `disjoint_subset` union algorithm for `union_all` and `dissolve` (#3534).
+- Added `constrained_delaunay_triangles` method to GeoSeries/GeoDataframe (#3552).
 - Added `to_pandas_kwargs` argument to `from_arrow`, `read_parquet` and `read_feather`
   to allow better control of conversion of non-geometric Arrow data to DataFrames (#3466).
 - Added `is_valid_coverage` and `invalid_coverage_edges` to GeoSeries/GeoDataFrame to allow validation of polygonal coverage (#3545).
 - Added `maximum_inscribed_circle` method from shapely to GeoSeries/GeoDataframe (#3544).
 - Added `minimum_clearance_line` method from shapely to GeoSeries/GeoDataframe (#3543).
 - Added ``method`` and ``keep_collapsed`` argument to ``make_valid`` (#3548).
+- Added `simplify_coverage` method for topological simplification of polygonal coverages
+  to GeoSeries/GeoDataframe (#3541).
 
 Bug fixes:
 
diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index 48ebea22..5020e51f 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -123,6 +123,7 @@ Constructive methods and attributes
    GeoSeries.segmentize
    GeoSeries.shortest_line
    GeoSeries.simplify
+   GeoSeries.simplify_coverage
    GeoSeries.snap
    GeoSeries.transform
 
@@ -156,6 +157,7 @@ Aggregating and exploding
    :toctree: api/
 
    GeoSeries.build_area
+   GeoSeries.constrained_delaunay_triangles
    GeoSeries.delaunay_triangles
    GeoSeries.explode
    GeoSeries.intersection_all
diff --git a/doc/source/docs/user_guide/geometric_manipulations.rst b/doc/source/docs/user_guide/geometric_manipulations.rst
index b9a121eb..4e2a8498 100644
--- a/doc/source/docs/user_guide/geometric_manipulations.rst
+++ b/doc/source/docs/user_guide/geometric_manipulations.rst
@@ -38,6 +38,14 @@ Constructive methods
   number of points in the object is less than three. For two points,
   the convex hull collapses to a `LineString`; for 1, a `Point`.
 
+.. method:: GeoSeries.constrained_delaunay_triangles
+
+  Returns a :class:`~geopandas.GeoSeries` with 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.
+
 .. method:: GeoSeries.delaunay_triangles(tolerance, preserve_topology=True)
 
   Returns a :class:`~geopandas.GeoSeries` consisting of polygons (default) or linestrings
diff --git a/geopandas/array.py b/geopandas/array.py
index 44011c88..51df9dc4 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -597,6 +597,14 @@ class GeometryArray(ExtensionArray):
     def concave_hull(self, ratio, allow_holes):
         return shapely.concave_hull(self._data, ratio=ratio, allow_holes=allow_holes)
 
+    def constrained_delaunay_triangles(self):
+        if not SHAPELY_GE_21:
+            raise ImportError("'constrained_delaunay_triangles' requires shapely>=2.1.")
+
+        return GeometryArray(
+            shapely.constrained_delaunay_triangles(self._data), crs=self.crs
+        )
+
     @property
     def convex_hull(self):
         return GeometryArray(shapely.convex_hull(self._data), crs=self.crs)
@@ -871,6 +879,18 @@ class GeometryArray(ExtensionArray):
             crs=self.crs,
         )
 
+    def simplify_coverage(self, tolerance, simplify_boundary=True):
+        if not (SHAPELY_GE_21 and GEOS_GE_312):
+            raise ImportError(
+                "'simplify_coverage' requires shapely>=2.1 and GEOS>=3.12."
+            )
+        return GeometryArray(
+            shapely.coverage_simplify(
+                self._data, tolerance, simplify_boundary=simplify_boundary
+            ),
+            crs=self.crs,
+        )
+
     def project(self, other, normalized=False):
         if isinstance(other, GeometryArray):
             other = other._data
diff --git a/geopandas/base.py b/geopandas/base.py
index 554e0027..010c2fca 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -1094,6 +1094,35 @@ GeometryCollection
             "concave_hull", self, ratio=ratio, allow_holes=allow_holes
         )
 
+    def constrained_delaunay_triangles(self):
+        """Returns a :class:`~geopandas.GeoSeries` with 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.
+
+        Examples
+        --------
+
+        >>> from shapely.geometry import Polygon
+        >>> s = geopandas.GeoSeries([Polygon([(0, 0), (1, 1), (0, 1)])])
+        >>> s
+        0                       POLYGON ((0 0, 1 1, 0 1, 0 0))
+        dtype: geometry
+
+        >>> s.constrained_delaunay_triangles()
+        0         GEOMETRYCOLLECTION (POLYGON ((0 0, 0 1, 1 1, 0...
+        dtype: geometry
+
+        See also
+        --------
+        GeoSeries.delaunay_triangles : Delaunay triangulation
+
+        """
+        return _delegate_geo_method("constrained_delaunay_triangles", self)
+
     @property
     def convex_hull(self):
         """Returns a ``GeoSeries`` of geometries representing the convex hull
@@ -1229,6 +1258,7 @@ GeometryCollection
         See also
         --------
         GeoSeries.voronoi_polygons : Voronoi diagram around vertices
+        GeoSeries.constrained_delaunay_triangles : constrained Delaunay triangulation
         """
         from .geoseries import GeoSeries
 
@@ -5289,9 +5319,14 @@ GeometryCollection
         to the straight line is smaller than `tolerance`. It does not
         move any points and it always preserves endpoints of
         the original line or polygon.
-        See http://shapely.readthedocs.io/en/latest/manual.html#object.simplify
+        See https://shapely.readthedocs.io/en/latest/manual.html#object.simplify
         for details
 
+        Simplifies individual geometries independently, without considering
+        the topology of a potential polygonal coverage. If you would like to treat
+        the ``GeoSeries`` as a coverage and simplify its edges, while preserving the
+        coverage topology, see :meth:`simplify_coverage`.
+
         Parameters
         ----------
         tolerance : float
@@ -5311,6 +5346,10 @@ GeometryCollection
         coordinates: two geometries differing only in order of coordinates may be
         simplified differently.
 
+        See also
+        --------
+        simplify_coverage : simplify geometries using coverage simplification
+
         Examples
         --------
         >>> from shapely.geometry import Point, LineString
@@ -5331,6 +5370,80 @@ GeometryCollection
             "simplify", self, tolerance=tolerance, preserve_topology=preserve_topology
         )
 
+    def simplify_coverage(self, tolerance, simplify_boundary=True):
+        """Returns a ``GeoSeries`` containing a simplified representation of
+        polygonal coverage.
+
+        Assumes that the ``GeoSeries`` forms a polygonal coverage. Under this
+        assumption, the method simplifies the edges using the Visvalingam-Whyatt
+        algorithm, while preserving a valid coverage. In the most simplified case,
+        polygons are reduced to triangles.
+
+        A ``GeoSeries`` of valid polygons is considered a coverage if the polygons are:
+
+        * **Non-overlapping** - polygons do not overlap (their interiors do not
+          intersect)
+        * **Edge-Matched** - vertices along shared edges are identical
+
+        The method allows simplification of all edges including the outer boundaries of
+        the coverage or simplification of only the inner (shared) edges.
+
+        If there are other geometry types than Polygons or MultiPolygons present, the
+        method will raise an error.
+
+        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 coverage topology.
+
+        .. versionadded:: 1.1.0
+
+        Parameters
+        ----------
+        tolerance : float
+            The degree of simplification roughly equal to the square root of the area
+            of triangles that will be removed. It has the same units
+            as the coordinate reference system of the GeoSeries.
+            For example, using `tolerance=100` in a projected CRS with meters
+            as units means a distance of 100 meters in reality.
+        simplify_boundary: bool (default True)
+            By default (True), simplifies both internal edges of the coverage as well
+            as its boundary. If set to False, only simplifies internal edges.
+
+
+        See also
+        --------
+        simplify : simplification of individual geometries
+
+        Examples
+        --------
+        >>> from shapely.geometry import Polygon
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...         Polygon([(0, 0), (1, 1.1), (2, 0), (1.5, 1), (2, 2), (0, 2)]),
+        ...         Polygon([(2, 0), (4, 0), (4, 2), (2, 2), (1.5, 1)]),
+        ...     ]
+        ... )
+        >>> s
+        0    POLYGON ((0 0, 1 1.1, 2 0, 1.5 1, 2 2, 0 2, 0 0))
+        1           POLYGON ((2 0, 4 0, 4 2, 2 2, 1.5 1, 2 0))
+        dtype: geometry
+
+        >>> s.simplify_coverage(1)
+        0         POLYGON ((2 0, 2 2, 0 2, 2 0))
+        1    POLYGON ((2 0, 4 0, 4 2, 2 2, 2 0))
+        dtype: geometry
+
+        >>> s.simplify_coverage(1, simplify_boundary=False)
+        0    POLYGON ((2 0, 2 2, 0 2, 0 0, 1 1.1, 2 0))
+        1           POLYGON ((2 0, 4 0, 4 2, 2 2, 2 0))
+        dtype: geometry
+        """
+        return _delegate_geo_method(
+            "simplify_coverage",
+            self,
+            tolerance=tolerance,
+            simplify_boundary=simplify_boundary,
+        )
+
     def relate(self, other, align=None):
         """
         Returns the DE-9IM intersection matrices for the geometries


Test patch:
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index 772ad6ab..42d39af9 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -1198,6 +1198,16 @@ class TestGeomMethods:
                 ratio=0.1, allow_holes=Series([True, False], index=[99, 98])
             )
 
+    @pytest.mark.skipif(not SHAPELY_GE_21, reason="requires shapely 2.1")
+    def test_constrained_delaunay_triangles(self):
+        input = GeoSeries([Polygon([(0, 0), (1, 1), (0, 1)])])
+        expected = GeoSeries(
+            [GeometryCollection([Polygon([(0, 0), (0, 1), (1, 1), (0, 0)])])]
+        )
+        assert_geoseries_equal(
+            input.constrained_delaunay_triangles(), expected, check_geom_type=True
+        )
+
     def test_convex_hull(self):
         # the convex hull of a square should be the same as the square
         assert_geoseries_equal(self.squares, self.squares.convex_hull)
@@ -1483,6 +1493,38 @@ class TestGeomMethods:
         ):
             self.g1.simplify(Series([0.1], index=[99]))
 
+    @pytest.mark.skipif(
+        not (GEOS_GE_312 and SHAPELY_GE_21), reason="GEOS 3.12 and shapely 2.1 needed."
+    )
+    def test_simplify_coverage(self):
+        s = GeoSeries(
+            [
+                shapely.Polygon(
+                    [(0, 0), (10, 1), (20, 0), (20, 10), (10, 5), (0, 10), (0, 0)]
+                ),
+                shapely.Polygon(
+                    [(0, 10), (10, 5), (20, 10), (20, 20), (0, 20), (0, 10)]
+                ),
+            ]
+        )
+        e = GeoSeries(
+            [
+                shapely.Polygon([(0, 0), (20, 0), (20, 10), (0, 10)]),
+                shapely.Polygon([(0, 10), (20, 10), (20, 20), (0, 20)]),
+            ]
+        )
+        assert_geoseries_equal(s.simplify_coverage(8), e.normalize())
+
+        e_boundary = GeoSeries(
+            [
+                shapely.Polygon([(0, 0), (10, 1), (20, 0), (20, 10), (0, 10)]),
+                shapely.Polygon([(0, 10), (20, 10), (20, 20), (0, 20)]),
+            ]
+        )
+        assert_geoseries_equal(
+            s.simplify_coverage(8, simplify_boundary=False), e_boundary.normalize()
+        )
+
     def test_envelope(self):
         e = self.g3.envelope
         assert np.all(e.geom_equals(self.sq))