Repository: geopandas/geopandas
Pull request #3534 (2025-03-24 20:14:36)

ENH: expose disjoint_subset union in union_all and dissolve
Mimics the way we exposed `coverage_union`.


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index ca42f3b5..2972cbf3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -14,6 +14,7 @@ New features and improvements:
 - GeoDataFrame no longer hard-codes the class internally, allowing easier subclassing (#3505).
 - 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).
 
 Bug fixes:
 
diff --git a/geopandas/_compat.py b/geopandas/_compat.py
index 53ed07e5..fae2fbd9 100644
--- a/geopandas/_compat.py
+++ b/geopandas/_compat.py
@@ -20,9 +20,11 @@ PANDAS_GE_30 = Version(pd.__version__) >= Version("3.0.0.dev0")
 # -----------------------------------------------------------------------------
 
 SHAPELY_GE_204 = Version(shapely.__version__) >= Version("2.0.4")
+SHAPELY_GE_21 = Version(shapely.__version__) >= Version("2.1rc1")
 
 GEOS_GE_390 = shapely.geos_version >= (3, 9, 0)
 GEOS_GE_310 = shapely.geos_version >= (3, 10, 0)
+GEOS_GE_312 = shapely.geos_version >= (3, 12, 0)
 
 
 def import_optional_dependency(name: str, extra: str = ""):
diff --git a/geopandas/array.py b/geopandas/array.py
index 3af29449..51c4e7d5 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -19,7 +19,14 @@ import shapely.ops
 import shapely.wkt
 from shapely.geometry.base import BaseGeometry
 
-from ._compat import HAS_PYPROJ, PANDAS_GE_21, PANDAS_GE_22, requires_pyproj
+from ._compat import (
+    GEOS_GE_312,
+    HAS_PYPROJ,
+    PANDAS_GE_21,
+    PANDAS_GE_22,
+    SHAPELY_GE_21,
+    requires_pyproj,
+)
 from .sindex import SpatialIndex
 
 if HAS_PYPROJ:
@@ -849,15 +856,22 @@ class GeometryArray(ExtensionArray):
         return self.union_all()
 
     def union_all(self, method="unary", grid_size=None):
+        if method != "unary" and grid_size is not None:
+            raise ValueError(f"grid_size is not supported for method '{method}'.")
         if method == "coverage":
-            if grid_size is not None:
-                raise ValueError("grid_size is not supported for method 'coverage'.")
             return shapely.coverage_union_all(self._data)
         elif method == "unary":
             return shapely.union_all(self._data, grid_size=grid_size)
+        elif method == "disjoint_subset":
+            if not (SHAPELY_GE_21 and GEOS_GE_312):
+                raise ImportError(
+                    "Method 'disjoin_subset' requires shapely>=2.1 and GEOS>=3.12."
+                )
+            return shapely.disjoint_subset_union_all(self._data)
         else:
             raise ValueError(
-                f"Method '{method}' not recognized. Use 'coverage' or 'unary'."
+                f"Method '{method}' not recognized. Use 'coverage', 'unary' or "
+                "'disjoint_subset'."
             )
 
     def intersection_all(self):
diff --git a/geopandas/base.py b/geopandas/base.py
index e2996365..71c76178 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -2140,6 +2140,8 @@ GeometryCollection
         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.
+        Alternatively, for situations which can be divided into many disjoint subsets,
+        ``method="disjoint_subset"`` may be preferable.
 
         Parameters
         ----------
@@ -2152,6 +2154,10 @@ GeometryCollection
               for non-overlapping polygons and can be significantly faster than the
               unary union algorithm. However, it can produce invalid geometries if the
               polygons overlap.
+            * ``"disjoint_subset:``: use the disjoint subset union algorithm. This
+              option is optimized for inputs that can be divided into subsets that do
+              not intersect. If there is only one such subset, performance can be
+              expected to be worse than ``"unary"``.
 
         grid_size : float, default None
             When grid size is specified, a fixed-precision space is used to perform the
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 1c98bc99..93a4b3db 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -2130,7 +2130,7 @@ default 'snappy'
         sort: bool = True,
         observed: bool = False,
         dropna: bool = True,
-        method: Literal["unary", "coverage"] = "unary",
+        method: Literal["unary", "coverage", "disjoint_subset"] = "unary",
         grid_size: float | None = None,
         **kwargs,
     ) -> GeoDataFrame:
@@ -2185,6 +2185,11 @@ default 'snappy'
               for non-overlapping polygons and can be significantly faster than the
               unary union algorithm. However, it can produce invalid geometries if the
               polygons overlap.
+            * ``"disjoint_subset:``: use the disjoint subset union algorithm. This
+              option is optimized for inputs that can be divided into subsets that do
+              not intersect. If there is only one such subset, performance can be
+              expected to be worse than ``"unary"``.
+
 
         grid_size : float, default None
             When grid size is specified, a fixed-precision space is used to perform the


Test patch:
diff --git a/geopandas/tests/test_array.py b/geopandas/tests/test_array.py
index 66a5ea87..0fb707cb 100644
--- a/geopandas/tests/test_array.py
+++ b/geopandas/tests/test_array.py
@@ -13,7 +13,7 @@ from shapely import MultiPolygon, Polygon, geos_version
 from shapely.geometry.base import CAP_STYLE, JOIN_STYLE
 
 import geopandas
-from geopandas._compat import HAS_PYPROJ
+from geopandas._compat import GEOS_GE_312, HAS_PYPROJ, SHAPELY_GE_21
 from geopandas.array import (
     GeometryArray,
     _check_crs,
@@ -681,6 +681,10 @@ def test_union_all():
     u_cov = G.union_all(method="coverage")
     assert u_cov.equals(expected)
 
+    if GEOS_GE_312 and SHAPELY_GE_21:
+        u_disjoint = G.union_all(method="disjoint_subset")
+        assert u_disjoint.equals(expected)
+
     with pytest.raises(ValueError, match="Method 'invalid' not recognized."):
         G.union_all(method="invalid")
 
diff --git a/geopandas/tests/test_dissolve.py b/geopandas/tests/test_dissolve.py
index 0d8e5e62..382e54dd 100644
--- a/geopandas/tests/test_dissolve.py
+++ b/geopandas/tests/test_dissolve.py
@@ -7,7 +7,7 @@ from shapely import MultiPolygon, Polygon
 
 import geopandas
 from geopandas import GeoDataFrame, read_file
-from geopandas._compat import HAS_PYPROJ, PANDAS_GE_30
+from geopandas._compat import GEOS_GE_312, HAS_PYPROJ, PANDAS_GE_30, SHAPELY_GE_21
 
 import pytest
 from geopandas.testing import assert_geodataframe_equal, geom_almost_equals
@@ -353,13 +353,16 @@ def test_dissolve_multi_agg(nybb_polydf, merged_shapes):
     assert_geodataframe_equal(test, merged_shapes)
 
 
-def test_coverage_dissolve(nybb_polydf):
+@pytest.mark.parametrize("method", ["coverage", "disjoint_subset"])
+def test_dissolve_method(nybb_polydf, method):
+    if method == "disjoint_subset" and not (GEOS_GE_312 and SHAPELY_GE_21):
+        pytest.skip("Unsupported shapely/GEOS.")
     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"),
+        others.geometry.union_all(method=method),
+        manhattan_bronx.geometry.union_all(method=method),
     ]
     merged_shapes = GeoDataFrame(
         {"myshapes": collapsed},
@@ -371,5 +374,5 @@ def test_coverage_dissolve(nybb_polydf):
     merged_shapes["BoroName"] = ["Staten Island", "Manhattan"]
     merged_shapes["BoroCode"] = [5, 1]
 
-    test = nybb_polydf.dissolve("manhattan_bronx", method="coverage")
+    test = nybb_polydf.dissolve("manhattan_bronx", method=method)
     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 f64d8bbc..2ef4b70a 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -20,7 +20,7 @@ from shapely.geometry.collection import GeometryCollection
 from shapely.ops import unary_union
 
 from geopandas import GeoDataFrame, GeoSeries
-from geopandas._compat import HAS_PYPROJ
+from geopandas._compat import GEOS_GE_312, HAS_PYPROJ, SHAPELY_GE_21
 from geopandas.base import GeoPandasBase
 
 import pytest
@@ -500,6 +500,8 @@ class TestGeomMethods:
         assert g3.union_all().equals(shapely.GeometryCollection())
 
         assert g.union_all(method="coverage").equals(expected)
+        if GEOS_GE_312 and SHAPELY_GE_21:
+            assert g.union_all(method="disjoint_subset").equals(expected)
 
     def test_unary_union_deprecated(self):
         p1 = self.t1