Repository: geopandas/geopandas
Pull request #3151 (2024-01-18 21:37:17)

ENH: include coverage_union_all in union_all and dissolve
Xref #2010

I was surprised how faster this actually is. On geodatasets' `geoda.south` it goes from 121 ms to 4.43 ms. It will likely depend on the data but it is impressive nonetheless. 

The formulations in the docstring feel a bit sloppy so I'd appreciate some native speaker to check if it makes sense.

martinfleis: Change the keyword to `method` with string options to allow `GEOSDisjointSubsetUnion` to be exposed in future. 

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2f725460..0517b198 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
diff --git a/geopandas/array.py b/geopandas/array.py
index d1bcdb26..7a6d9e24 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -842,8 +842,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..b89d8f63 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
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..bdec0770 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)])