Repository: shapely/shapely
Pull request #1982 (2024-02-06 15:57:24)

ENH: add disjoint_subset_union and disjoint_subset_union_all
Exposing `GEOSDisjointSubsetUnion` as `disjoint_subset_union` and `disjoint_subset_union_all`, matching the API of union and coverage_union (assuming coverage_unoin will be fixed as discussed in #1963).

Tests are currently only reusing the tooling for other set operations. Not sure if there's a need to go beyond that.

xref #1510



coveralls: ## Pull Request Test Coverage Report for [Build 7806019234](https://coveralls.io/builds/65540605)


* **-3** of **20**   **(85.0%)**  changed or added relevant lines in **1** file are covered.
* No unchanged relevant lines lost coverage.
* Overall coverage increased (+**0.3%**) to **87.995%**

---

|  Changes Missing Coverage | Covered Lines | Changed/Added Lines | % |
| :-----|--------------|--------|---: |
| [shapely/set_operations.py](https://coveralls.io/builds/65540605/source?filename=shapely%2Fset_operations.py#L549) | 17 | 20 | 85.0%
<!-- | **Total:** | **17** | **20** | **85.0%** | -->


|  Totals | [![Coverage Status](https://coveralls.io/builds/65540605/badge)](https://coveralls.io/builds/65540605) |
| :-- | --: |
| Change from base [Build 7796641302](https://coveralls.io/builds/65522606): |  0.3% |
| Covered Lines: | 2602 |
| Relevant Lines: | 2957 |

---
##### 💛  - [Coveralls](https://coveralls.io)


Patch:
diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml
index f2def91..0138cb3 100644
--- a/.github/workflows/lint.yml
+++ b/.github/workflows/lint.yml
@@ -11,7 +11,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v5
         with:
-          python-version: 3.8
+          python-version: 3.12
 
       - name: Run ruff
         uses: pre-commit/action@v3.0.1
diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
index d6151d0..79f9282 100644
--- a/.github/workflows/release.yml
+++ b/.github/workflows/release.yml
@@ -118,7 +118,7 @@ jobs:
         uses: pypa/cibuildwheel@v2.20.0
         env:
           CIBW_ARCHS: ${{ matrix.arch }}
-          CIBW_SKIP: cp36-* pp* *musllinux* *-manylinux_i686
+          CIBW_SKIP: pp* *musllinux* *-manylinux_i686
           CIBW_TEST_SKIP: "cp38-macosx_arm64"
           CIBW_ENVIRONMENT_LINUX:
             GEOS_VERSION=${{ env.GEOS_VERSION }}
diff --git a/CHANGES.txt b/CHANGES.txt
index c00cba9..d602d07 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -10,7 +10,7 @@ Bug fixes:
 
 Improvements:
 
-- Require GEOS >= 3.8, NumPy >= 1.16, and Python >= 3.8 (#1802, #1885)
+- Require GEOS >= 3.9, NumPy >= 1.20, and Python >= 3.9 (#1802, #1885, #2124)
 - Handle ``Feature`` type in ``shapely.geometry.shape`` (#1815)
 - Add a ``handle_nan`` parameter to ``shapely.linestrings()`` and ``shapely.linearrings()``
   to allow, skip, or error on nonfinite (NaN / Inf) coordinates. The default
@@ -21,6 +21,9 @@ Improvements:
   lets it automatically detect the dimensionality of each input geometry.
 - Add parameters ``method`` and ``keep_collapsed`` to ``shapely.make_valid()`` (#1941)
 - Upgraded the GEOS version in the binary wheel distributions to 3.12.1.
+- Add ``disjoint_subset_union`` and ``disjoint_subset_union_all`` as an optimized
+  version of union and union_all, assuming inputs can be divided into subsets that do
+  not intersect. Requires at least GEOS 3.12.
 - The ``voronoi_polygons`` now accepts the ``ordered`` keyword, optionally forcing the
   order of polygons within the GeometryCollection to follow the order of input
   coordinates. Requires at least GEOS 3.12. (#1968)
diff --git a/README.rst b/README.rst
index 79ea2de..f8a58d1 100644
--- a/README.rst
+++ b/README.rst
@@ -111,9 +111,9 @@ Requirements
 
 Shapely 2.1 requires
 
-* Python >=3.8
+* Python >=3.9
 * GEOS >=3.9
-* NumPy >=1.16
+* NumPy >=1.20
 
 Installing Shapely
 ==================
diff --git a/asv.conf.json b/asv.conf.json
index d6dd865..6597043 100644
--- a/asv.conf.json
+++ b/asv.conf.json
@@ -55,7 +55,7 @@
 
     // The Pythons you'd like to test against.  If not provided, defaults
     // to the current version of Python used to run `asv`.
-    "pythons": ["3.8"],
+    "pythons": ["3.12"],
 
     // The list of conda channel names to be searched for benchmark
     // dependency packages in the specified order
diff --git a/docs/release/2.x.rst b/docs/release/2.x.rst
index 051364a..37a9369 100644
--- a/docs/release/2.x.rst
+++ b/docs/release/2.x.rst
@@ -8,7 +8,7 @@ Version 2.1.0 (unreleased)
 
 Improvements:
 
-- Require GEOS >= 3.7, NumPy >= 1.16, and Python >= 3.8 (#1802)
+- Require GEOS >= 3.9, NumPy >= 1.20, and Python >= 3.9 (#1802, #1885, #2124)
 - Add a ``handle_nan`` parameter to ``shapely.linestrings()`` and ``shapely.linearrings()``
   to allow, skip, or error on nonfinite (NaN / Inf) coordinates. The default
   behaviour (allow) is backwards compatible (#1594).
diff --git a/pyproject.toml b/pyproject.toml
index bc3970a..5a2dec7 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -33,15 +33,16 @@ classifiers = [
     "Operating System :: MacOS",
     "Operating System :: Microsoft :: Windows",
     "Programming Language :: Python :: 3",
-    "Programming Language :: Python :: 3.8",
     "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
+    "Programming Language :: Python :: 3.12",
+    "Programming Language :: Python :: 3.13",
     "Topic :: Scientific/Engineering :: GIS",
 ]
-requires-python = ">=3.8"
+requires-python = ">=3.9"
 dependencies = [
-    "numpy>=1.16",
+    "numpy>=1.20",
 ]
 
 [project.optional-dependencies]
diff --git a/shapely/geometry/base.py b/shapely/geometry/base.py
index e40babd..aa7ca0e 100644
--- a/shapely/geometry/base.py
+++ b/shapely/geometry/base.py
@@ -214,11 +214,8 @@ class BaseGeometry(shapely.Geometry):
         """Return True if geometries are equal, else False."""
         if not isinstance(other, BaseGeometry):
             return NotImplemented
-        # equal_nan=False is the default, but not yet available for older numpy
-        # TODO updated once we require numpy >= 1.19
         return type(other) is type(self) and np.array_equal(
-            self.coords,
-            other.coords,  # , equal_nan=False
+            self.coords, other.coords, equal_nan=False
         )
 
     def __ne__(self, other):
diff --git a/shapely/geometry/polygon.py b/shapely/geometry/polygon.py
index 11dd68f..071225a 100644
--- a/shapely/geometry/polygon.py
+++ b/shapely/geometry/polygon.py
@@ -283,10 +283,9 @@ class Polygon(BaseGeometry):
         ]
         if not len(my_coords) == len(other_coords):
             return False
-        # equal_nan=False is the default, but not yet available for older numpy
         return np.all(
             [
-                np.array_equal(left, right)  # , equal_nan=False)
+                np.array_equal(left, right, equal_nan=False)
                 for left, right in zip(my_coords, other_coords)
             ]
         )
diff --git a/shapely/set_operations.py b/shapely/set_operations.py
index 95ba384..5920b7e 100644
--- a/shapely/set_operations.py
+++ b/shapely/set_operations.py
@@ -2,8 +2,8 @@
 
 import numpy as np
 
-from shapely import GeometryType, lib
-from shapely.decorators import multithreading_enabled
+from shapely import Geometry, GeometryType, lib
+from shapely.decorators import multithreading_enabled, requires_geos
 
 __all__ = [
     "difference",
@@ -16,6 +16,8 @@ __all__ = [
     "union_all",
     "coverage_union",
     "coverage_union_all",
+    "disjoint_subset_union",
+    "disjoint_subset_union_all",
 ]
 
 
@@ -490,3 +492,114 @@ def coverage_union_all(geometries, axis=None, **kwargs):
         geometries, np.intc(GeometryType.GEOMETRYCOLLECTION)
     )
     return lib.coverage_union(collections, **kwargs)
+
+
+@requires_geos("3.12.0")
+@multithreading_enabled
+def disjoint_subset_union(a, b, **kwargs):
+    """Merge multiple polygons into one using algorithm optimised for subsets.
+
+    This is an optimized version of union which assumes inputs can be
+    divided into subsets that do not intersect.
+
+    If there is only one such subset, performance can be expected to be worse than
+    :func:`union`. As such, it is recommeded to use ``disjoint_subset_union`` with
+    GeometryCollections rather than individual geometries.
+
+    Parameters
+    ----------
+    a, b : Geometry or array_like
+        Geometry or geometries to merge (union).
+    **kwargs
+        See :ref:`NumPy ufunc docs <ufuncs.kwargs>` for other keyword arguments.
+
+    See Also
+    --------
+    union
+    coverage_union
+    disjoint_subset_union_all
+
+    Examples
+    --------
+    >>> from shapely import normalize, Polygon
+    >>> polygon_1 = Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)])
+    >>> polygon_2 = Polygon([(1, 0), (1, 1), (2, 1), (2, 0), (1, 0)])
+    >>> normalize(disjoint_subset_union(polygon_1, polygon_2))
+    <POLYGON ((0 0, 0 1, 1 1, 2 1, 2 0, 1 0, 0 0))>
+
+    Union with None returns same polygon:
+
+    >>> normalize(disjoint_subset_union(polygon_1, None))
+    <POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))>
+    """
+    if (isinstance(a, Geometry) or a is None) and (
+        isinstance(b, Geometry) or b is None
+    ):
+        pass
+    elif isinstance(a, Geometry) or a is None:
+        a = np.full_like(b, a)
+    elif isinstance(b, Geometry) or b is None:
+        b = np.full_like(a, b)
+    elif len(a) != len(b):
+        raise ValueError("Arrays a and b must have the same length")
+    return disjoint_subset_union_all([a, b], axis=0, **kwargs)
+
+
+@requires_geos("3.12.0")
+@multithreading_enabled
+def disjoint_subset_union_all(geometries, axis=None, **kwargs):
+    """Return the union of multiple polygons.
+
+    This is an optimized version of union which assumes inputs can be divided into
+    subsets that do not intersect.
+
+    If there is only one such subset, performance can be expected to be worse than
+    :func:`union_all`.
+
+    This function ignores None values when other Geometry elements are present.
+    If all elements of the given axis are None, an empty GeometryCollection is
+    returned.
+
+    Parameters
+    ----------
+    geometries : array_like
+        Geometries to union.
+    axis : int, optional
+        Axis along which the operation is performed. The default (None)
+        performs the operation over all axes, returning a scalar value.
+        Axis may be negative, in which case it counts from the last to the
+        first axis.
+    **kwargs
+        See :ref:`NumPy ufunc docs <ufuncs.kwargs>` for other keyword arguments.
+
+    See Also
+    --------
+    coverage_union_all
+    union_all
+    disjoint_subset_union
+
+    Examples
+    --------
+    >>> from shapely import normalize, Polygon
+    >>> polygon_1 = Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)])
+    >>> polygon_2 = Polygon([(1, 0), (1, 1), (2, 1), (2, 0), (1, 0)])
+    >>> normalize(disjoint_subset_union_all([polygon_1, polygon_2]))
+    <POLYGON ((0 0, 0 1, 1 1, 2 1, 2 0, 1 0, 0 0))>
+    >>> normalize(disjoint_subset_union_all([polygon_1, None]))
+    <POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))>
+    >>> normalize(disjoint_subset_union_all([None, None]))
+    <GEOMETRYCOLLECTION EMPTY>
+    """
+    geometries = np.asarray(geometries)
+    if axis is None:
+        geometries = geometries.ravel()
+    else:
+        geometries = np.rollaxis(
+            np.asarray(geometries), axis=axis, start=geometries.ndim
+        )
+    # create_collection acts on the inner axis
+    collections = lib.create_collection(
+        geometries, np.intc(GeometryType.GEOMETRYCOLLECTION)
+    )
+
+    return lib.disjoint_subset_union(collections, **kwargs)
diff --git a/shapely/strtree.py b/shapely/strtree.py
index faaffb5..0d28467 100644
--- a/shapely/strtree.py
+++ b/shapely/strtree.py
@@ -1,6 +1,7 @@
 """STRtree spatial index for efficient spatial queries."""
 
-from typing import Any, Iterable, Union
+from collections.abc import Iterable
+from typing import Any, Union
 
 import numpy as np
 
diff --git a/src/ufuncs.c b/src/ufuncs.c
index 195a3be..6e31da5 100644
--- a/src/ufuncs.c
+++ b/src/ufuncs.c
@@ -557,6 +557,9 @@ static void* oriented_envelope_data[1] = {GEOSMinimumRotatedRectangle_r};
 #if GEOS_SINCE_3_11_0
 static void* line_merge_directed_data[1] = {GEOSLineMergeDirected_r};
 #endif
+#if GEOS_SINCE_3_12_0
+static void* disjoint_subset_union_data[1] = {GEOSDisjointSubsetUnion_r};
+#endif  // GEOS_SINCE_3_12_0
 typedef void* FuncGEOS_Y_Y(void* context, void* a);
 static char Y_Y_dtypes[2] = {NPY_OBJECT, NPY_OBJECT};
 static void Y_Y_func(char** args, const npy_intp* dimensions, const npy_intp* steps, void* data) {
@@ -1588,7 +1591,7 @@ static char make_valid_with_params_inner(void* ctx, GEOSMakeValidParams* params,
   return PGERR_SUCCESS;
 }
 
-static char make_valid_with_params_dtypes[4] = {NPY_OBJECT, NPY_INT, NPY_BOOL, 
+static char make_valid_with_params_dtypes[4] = {NPY_OBJECT, NPY_INT, NPY_BOOL,
                                                 NPY_OBJECT};
 static void make_valid_with_params_func(char** args, const npy_intp* dimensions,
                                         const npy_intp* steps, void* data) {
@@ -3807,6 +3810,7 @@ int init_ufuncs(PyObject* m, PyObject* d) {
 #endif
 
 #if GEOS_SINCE_3_12_0
+  DEFINE_Y_Y(disjoint_subset_union);
   DEFINE_Y_b(has_m);
   DEFINE_Y_d(get_m);
 #endif


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index c88fa7a..938649f 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -16,14 +16,10 @@ jobs:
         architecture: [x64]
         geos: [3.9.5, 3.10.6, 3.11.4, 3.12.2, 3.13.0beta1, main]
         include:
-          # 2019
-          - python: 3.8
-            geos: 3.9.5  # this is here temporarily to allow matrix to work
-            numpy: 1.17.3
           # 2020
           - python: 3.9
             geos: 3.9.5
-            numpy: 1.19.5
+            numpy: 1.20.3
           # 2021
           - python: "3.10"
             geos: 3.10.6
@@ -55,9 +51,9 @@ jobs:
           # enable two 32-bit windows builds:
           - os: windows-2019
             architecture: x86
-            python: 3.8
+            python: 3.9
             geos: 3.9.5
-            numpy: 1.16.2
+            numpy: 1.20.3
           - os: windows-2019
             architecture: x86
             python: "3.11"
@@ -65,7 +61,7 @@ jobs:
             numpy: 1.24.4
           # pypy (use explicit ubuntu version to not overwrite existing ubuntu-latest + geos 3.11.0 build)
           - os: ubuntu-22.04
-            python: "pypy3.8"
+            python: "pypy3.10"
             geos: 3.12.2
             numpy: 1.24.4
 
@@ -157,7 +153,7 @@ jobs:
 
       # Windows requires special treatment:
       # - geos-config does not exist, so we specify include and library paths
-      # - Python >=3.8 ignores the PATH for finding DLLs, so we copy them into the package
+      # - Python ignores the PATH for finding DLLs, so we copy them into the package
       - name: Set environment variables + copy DLLs (Windows)
         run: |
           cp geosinstall/geos-${{ matrix.geos }}/bin/*.dll shapely
diff --git a/shapely/tests/test_set_operations.py b/shapely/tests/test_set_operations.py
index 52e91c0..4ce0f35 100644
--- a/shapely/tests/test_set_operations.py
+++ b/shapely/tests/test_set_operations.py
@@ -31,6 +31,11 @@ REDUCE_SET_OPERATIONS = (
 # operations that support fixed precision
 REDUCE_SET_OPERATIONS_PREC = ((shapely.union_all, shapely.union),)
 
+if shapely.geos_version >= (3, 12, 0):
+    SET_OPERATIONS += (shapely.disjoint_subset_union,)
+    REDUCE_SET_OPERATIONS += (
+        (shapely.disjoint_subset_union_all, shapely.disjoint_subset_union),
+    )
 
 reduce_test_data = [
     shapely.box(0, 0, 5, 5),
@@ -64,6 +69,8 @@ def test_set_operation_array(a, func):
 
 @pytest.mark.parametrize("func", SET_OPERATIONS)
 def test_set_operation_prec_nonscalar_grid_size(func):
+    if func is shapely.disjoint_subset_union:
+        pytest.skip("disjoint_subset_union does not support grid_size")
     with pytest.raises(
         ValueError, match="grid_size parameter only accepts scalar values"
     ):
@@ -74,6 +81,8 @@ def test_set_operation_prec_nonscalar_grid_size(func):
 @pytest.mark.parametrize("func", SET_OPERATIONS)
 @pytest.mark.parametrize("grid_size", [0, 1, 2])
 def test_set_operation_prec_array(a, func, grid_size):
+    if func is shapely.disjoint_subset_union:
+        pytest.skip("disjoint_subset_union does not support grid_size")
     actual = func([a, a], point, grid_size=grid_size)
     assert actual.shape == (2,)
     assert isinstance(actual[0], Geometry)