Repository: geopandas/geopandas
Pull request #3211 (2024-03-06 10:15:14)

ENH: add relate_pattern method
xref #2010


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2f725460..fdd3e6a0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,7 +12,8 @@ Notes on dependencies:
 
 API changes:
 
-- `unary_union` is now deprecated and replaced by the `union_all` method (#3007).
+- `unary_union` attribute is now deprecated and replaced by the `union_all()` method (#3007) allowing
+  opting for a faster union algorithm for coverages (#3151).
 - `align` keyword in binary methods now defaults to `None`, treated as True. Explicit True
   will silence the warning about mismachted indices. (#3212)
 - The `sjoin` method will now preserve the name of the index of the right
@@ -23,6 +24,7 @@ API changes:
 
 New methods:
 
+- Added `relate_pattern` method from shapely to GeoSeries/GeoDataframe (#3211).
 - Added `intersection_all` method from shapely to GeoSeries/GeoDataframe (#3228).
 - Added `line_merge` method from shapely to GeoSeries/GeoDataframe (#3214).
 - Added `set_precision` and `get_precision` methods from shapely to GeoSeries/GeoDataframe (#3175).
diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index 74d61f0a..d17ef7c7 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -74,6 +74,7 @@ Binary predicates
    GeoSeries.covers
    GeoSeries.covered_by
    GeoSeries.relate
+   GeoSeries.relate_pattern
 
 
 Set-theoretic methods
diff --git a/geopandas/array.py b/geopandas/array.py
index d1bcdb26..d745712f 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -829,6 +829,11 @@ class GeometryArray(ExtensionArray):
             other = other._data
         return shapely.relate(self._data, other)
 
+    def relate_pattern(self, other, pattern):
+        if isinstance(other, GeometryArray):
+            other = other._data
+        return shapely.relate_pattern(self._data, other, pattern)
+
     #
     # Reduction operations that return a Shapely geometry
     #
@@ -842,8 +847,15 @@ class GeometryArray(ExtensionArray):
         )
         return self.union_all()
 
-    def union_all(self):
-        return shapely.union_all(self._data)
+    def union_all(self, method="unary"):
+        if method == "coverage":
+            return shapely.coverage_union_all(self._data)
+        elif method == "unary":
+            return shapely.union_all(self._data)
+        else:
+            raise ValueError(
+                f"Method '{method}' not recognized. Use 'coverage' or 'unary'."
+            )
 
     def intersection_all(self):
         return shapely.intersection_all(self._data)
diff --git a/geopandas/base.py b/geopandas/base.py
index 8fa660db..40eab3c5 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -1749,10 +1749,26 @@ GeometryCollection
 
         return self.geometry.values.union_all()
 
-    def union_all(self):
+    def union_all(self, method="unary"):
         """Returns a geometry containing the union of all geometries in the
         ``GeoSeries``.
 
+        By default, the unary union algorithm is used. If the geometries are
+        non-overlapping (forming a coverage), GeoPandas can use a significantly faster
+        algorithm to perform the union using the ``method="coverage"`` option.
+
+        Parameters
+        ----------
+        method : str (default ``"unary"``)
+            The method to use for the union. Options are:
+
+            * ``"unary"``: use the unary union algorithm. This option is the most robust
+              but can be slow for large numbers of geometries (default).
+            * ``"coverage"``: use the coverage union algorithm. This option is optimized
+              for non-overlapping polygons and can be significantly faster than the
+              unary union algorithm. However, it can produce invalid geometries if the
+              polygons overlap.
+
         Examples
         --------
 
@@ -1766,7 +1782,7 @@ GeometryCollection
         >>> s.union_all()
         <POLYGON ((0 1, 0 2, 2 2, 2 0, 1 0, 0 0, 0 1))>
         """
-        return self.geometry.values.union_all()
+        return self.geometry.values.union_all(method=method)
 
     def intersection_all(self):
         """Returns a geometry containing the intersection of all geometries in
@@ -4732,6 +4748,118 @@ GeometryCollection
         """
         return _binary_op("relate", self, other, align)
 
+    def relate_pattern(self, other, pattern, align=None):
+        """
+        Returns True if the DE-9IM string code for the relationship between
+        the geometries satisfies the pattern, else False.
+
+        This function compares the DE-9IM code string for two geometries
+        against a specified pattern. If the string matches the pattern then
+        ``True`` is returned, otherwise ``False``. The pattern specified can
+        be an exact match (``0``, ``1`` or ``2``), a boolean match
+        (uppercase ``T`` or ``F``), or a wildcard (``*``). For example,
+        the pattern for the ``within`` predicate is ``'T*F**F***'``
+
+        The operation works on a 1-to-1 row-wise manner:
+
+        .. image:: ../../../_static/binary_op-01.svg
+           :align: center
+
+        Parameters
+        ----------
+        other : BaseGeometry or GeoSeries
+            The other geometry to be tested agains the pattern.
+        pattern : str
+            The DE-9IM pattern to test against.
+        align : bool | None (default None)
+            If True, automatically aligns GeoSeries based on their indices.
+            If False, the order of elements is preserved. None defaults to True.
+
+        Returns
+        -------
+        Series
+
+        Examples
+        --------
+        >>> from shapely.geometry import Polygon, LineString, Point
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...         Polygon([(0, 0), (2, 2), (0, 2)]),
+        ...         Polygon([(0, 0), (2, 2), (0, 2)]),
+        ...         LineString([(0, 0), (2, 2)]),
+        ...         LineString([(2, 0), (0, 2)]),
+        ...         Point(0, 1),
+        ...     ],
+        ... )
+        >>> s2 = geopandas.GeoSeries(
+        ...     [
+        ...         Polygon([(0, 0), (1, 1), (0, 1)]),
+        ...         LineString([(1, 0), (1, 3)]),
+        ...         LineString([(2, 0), (0, 2)]),
+        ...         Point(1, 1),
+        ...         Point(0, 1),
+        ...     ],
+        ...     index=range(1, 6),
+        ... )
+
+        >>> s
+        0    POLYGON ((0 0, 2 2, 0 2, 0 0))
+        1    POLYGON ((0 0, 2 2, 0 2, 0 0))
+        2             LINESTRING (0 0, 2 2)
+        3             LINESTRING (2 0, 0 2)
+        4                       POINT (0 1)
+        dtype: geometry
+
+        >>> s2
+        1    POLYGON ((0 0, 1 1, 0 1, 0 0))
+        2             LINESTRING (1 0, 1 3)
+        3             LINESTRING (2 0, 0 2)
+        4                       POINT (1 1)
+        5                       POINT (0 1)
+        dtype: geometry
+
+        We can check the relate pattern of each geometry and a single
+        shapely geometry:
+
+        .. image:: ../../../_static/binary_op-03.svg
+           :align: center
+
+        >>> s.relate_pattern(Polygon([(0, 0), (1, 1), (0, 1)]), "2*T***F**")
+        0     True
+        1     True
+        2    False
+        3    False
+        4    False
+        dtype: bool
+
+        We can also check two GeoSeries against each other, row by row.
+        The GeoSeries above have different indices. We can either align both GeoSeries
+        based on index values and compare elements with the same index using
+        ``align=True`` or ignore index and compare elements based on their matching
+        order using ``align=False``:
+
+        .. image:: ../../../_static/binary_op-02.svg
+
+        >>> s.relate_pattern(s2, "TF******T", align=True)
+        0    False
+        1    False
+        2     True
+        3     True
+        4    False
+        5    False
+        dtype: bool
+
+        >>> s.relate_pattern(s2, "TF******T", align=False)
+        0    False
+        1     True
+        2     True
+        3     True
+        4     True
+        dtype: bool
+
+        """
+        return _binary_op("relate_pattern", self, other, pattern=pattern, align=align)
+
     def project(self, other, normalized=False, align=None):
         """
         Return the distance along each geometry nearest to *other*
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 170dc56a..3a9791b8 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -1755,6 +1755,7 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
         sort=True,
         observed=False,
         dropna=True,
+        method="unary",
         **kwargs,
     ):
         """
@@ -1787,26 +1788,28 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
         level : int or str or sequence of int or sequence of str, default None
             If the axis is a MultiIndex (hierarchical), group by a
             particular level or levels.
-
-            .. versionadded:: 0.9.0
         sort : bool, default True
             Sort group keys. Get better performance by turning this off.
             Note this does not influence the order of observations within
             each group. Groupby preserves the order of rows within each group.
-
-            .. versionadded:: 0.9.0
         observed : bool, default False
             This only applies if any of the groupers are Categoricals.
             If True: only show observed values for categorical groupers.
             If False: show all values for categorical groupers.
-
-            .. versionadded:: 0.9.0
         dropna : bool, default True
             If True, and if group keys contain NA values, NA values
             together with row/column will be dropped. If False, NA
             values will also be treated as the key in groups.
+        method : str (default ``"unary"``)
+            The method to use for the union. Options are:
+
+            * ``"unary"``: use the unary union algorithm. This option is the most robust
+              but can be slow for large numbers of geometries (default).
+            * ``"coverage"``: use the coverage union algorithm. This option is optimized
+              for non-overlapping polygons and can be significantly faster than the
+              unary union algorithm. However, it can produce invalid geometries if the
+              polygons overlap.
 
-            .. versionadded:: 0.9.0
         **kwargs :
             Keyword arguments to be passed to the pandas `DataFrameGroupby.agg` method
             which is used by `dissolve`. In particular, `numeric_only` may be
@@ -1881,7 +1884,7 @@ properties': {'col1': 'name1'}, 'geometry': {'type': 'Point', 'coordinates': (1.
 
         # Process spatial component
         def merge_geometries(block):
-            merged_geom = block.union_all()
+            merged_geom = block.union_all(method=method)
             return merged_geom
 
         g = self.groupby(group_keys=False, **groupby_kwargs)[self.geometry.name].agg(


Test patch:
diff --git a/geopandas/tests/test_array.py b/geopandas/tests/test_array.py
index 9de3679b..d5b88bd7 100644
--- a/geopandas/tests/test_array.py
+++ b/geopandas/tests/test_array.py
@@ -685,6 +685,12 @@ def test_union_all():
     expected = shapely.geometry.Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
     assert u.equals(expected)
 
+    u_cov = G.union_all(method="coverage")
+    assert u_cov.equals(expected)
+
+    with pytest.raises(ValueError, match="Method 'invalid' not recognized."):
+        G.union_all(method="invalid")
+
 
 @pytest.mark.parametrize(
     "attr, arg",
diff --git a/geopandas/tests/test_dissolve.py b/geopandas/tests/test_dissolve.py
index f9f7cfce..3bcdff9e 100644
--- a/geopandas/tests/test_dissolve.py
+++ b/geopandas/tests/test_dissolve.py
@@ -349,3 +349,25 @@ def test_dissolve_multi_agg(nybb_polydf, merged_shapes):
         )
     assert_geodataframe_equal(test, merged_shapes)
     assert len(record) == 0
+
+
+def test_coverage_dissolve(nybb_polydf):
+    manhattan_bronx = nybb_polydf.loc[3:4]
+    others = nybb_polydf.loc[0:2]
+
+    collapsed = [
+        others.geometry.union_all(method="coverage"),
+        manhattan_bronx.geometry.union_all(method="coverage"),
+    ]
+    merged_shapes = GeoDataFrame(
+        {"myshapes": collapsed},
+        geometry="myshapes",
+        index=pd.Index([5, 6], name="manhattan_bronx"),
+        crs=nybb_polydf.crs,
+    )
+
+    merged_shapes["BoroName"] = ["Staten Island", "Manhattan"]
+    merged_shapes["BoroCode"] = [5, 1]
+
+    test = nybb_polydf.dissolve("manhattan_bronx", method="coverage")
+    assert_frame_equal(merged_shapes, test, check_column_type=False)
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index e03c1128..e8af60de 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -485,6 +485,8 @@ class TestGeomMethods:
         g3 = GeoSeries([None, None])
         assert g3.union_all().equals(shapely.GeometryCollection())
 
+        assert g.union_all(method="coverage").equals(expected)
+
     def test_unary_union_deprecated(self):
         p1 = self.t1
         p2 = Polygon([(2, 0), (3, 0), (3, 1)])
@@ -631,6 +633,31 @@ class TestGeomMethods:
         )
         assert_array_dtype_equal(expected, self.g0.relate(self.g9, align=False))
 
+    def test_relate_pattern(self):
+        expected = Series([True] * 4 + [False] * 3, index=self.g0.index, dtype=bool)
+        assert_array_dtype_equal(
+            expected, self.g0.relate_pattern(self.inner_sq, "2********")
+        )
+
+        expected = Series([True, False], index=self.g6.index, dtype=bool)
+        assert_array_dtype_equal(
+            expected, self.g6.relate_pattern(self.na_none, "FF0******")
+        )
+
+        expected = Series(
+            [False] + [True] * 5 + [False, False], index=range(8), dtype=bool
+        )
+        with pytest.warns(UserWarning, match="The indices of the left and right"):
+            assert_array_dtype_equal(
+                expected, self.g0.relate_pattern(self.g9, "T********", align=None)
+            )
+        expected = Series(
+            [False] + [True] * 2 + [False] * 4, index=self.g0.index, dtype=bool
+        )
+        assert_array_dtype_equal(
+            expected, self.g0.relate_pattern(self.g9, "T********", align=False)
+        )
+
     def test_distance(self):
         expected = Series(
             np.array([np.sqrt((5 - 1) ** 2 + (5 - 1) ** 2), np.nan]), self.na_none.index