Repository: shapely/shapely
Pull request #2229 (2025-02-18 09:13:00)

TST: fix testing with older GEOS versions (3.9 to 3.11)
This PR allows older point-releases of GEOS to be successfully tested with shapely, usually by skipping over older upstream bugs and expected behaviours. For example [GEOS 3.10.2 on Ubuntu jammy](https://packages.ubuntu.com/jammy/libgeos-c1v5) previously raise a segfault.

I've manually checked this PR with a few GEOS version ranges (linux x64 only):

- GEOS 3.9.0 to 3.9.5 ([NEWS](https://github.com/libgeos/geos/blob/3.9/NEWS))
- GEOS 3.10.1 to 3.10.6 (not including GEOS 3.10.0) ([NEWS](https://github.com/libgeos/geos/blob/3.10/NEWS))
- GEOS 3.11.0 to 3.11.4 ([NEWS.md](https://github.com/libgeos/geos/blob/3.11/NEWS.md))


Patch:
diff --git a/shapely/_ragged_array.py b/shapely/_ragged_array.py
index 20cfd8a..2546110 100644
--- a/shapely/_ragged_array.py
+++ b/shapely/_ragged_array.py
@@ -63,7 +63,15 @@ def _get_arrays_point(arr, include_z):
 
 
 def _indices_to_offsets(indices, n):
-    offsets = np.insert(np.bincount(indices).cumsum(), 0, 0)
+    # default to int32 offsets if possible (to prefer the non-large arrow list variants)
+    # n_coords is the length of the array the indices are poin
+    if len(indices) > 2147483647:
+        dtype = np.int64
+    else:
+        dtype = np.int32
+
+    offsets = np.insert(np.bincount(indices).cumsum(dtype=dtype), 0, 0)
+
     if len(offsets) != n + 1:
         # last geometries might be empty or missing
         offsets = np.pad(
@@ -103,7 +111,7 @@ def _get_arrays_multilinestring(arr, include_z):
 
     # the coords and offsets into the coordinates of the linestrings
     coords, indices = get_coordinates(arr_flat, return_index=True, include_z=include_z)
-    offsets1 = np.insert(np.bincount(indices).cumsum(), 0, 0)
+    offsets1 = _indices_to_offsets(indices, len(arr_flat))
 
     return coords, (offsets1, offsets2)
 
@@ -116,7 +124,7 @@ def _get_arrays_polygon(arr, include_z):
 
     # the coords and offsets into the coordinates of the rings
     coords, indices = get_coordinates(arr_flat, return_index=True, include_z=include_z)
-    offsets1 = np.insert(np.bincount(indices).cumsum(), 0, 0)
+    offsets1 = _indices_to_offsets(indices, len(arr_flat))
 
     return coords, (offsets1, offsets2)
 
@@ -130,11 +138,11 @@ def _get_arrays_multipolygon(arr, include_z):
     # explode/flatten the Polygons into Rings
     arr_flat2, ring_indices = get_rings(arr_flat, return_index=True)
     # the offsets into the exterior/interior rings of the multipolygon parts
-    offsets2 = np.insert(np.bincount(ring_indices).cumsum(), 0, 0)
+    offsets2 = _indices_to_offsets(ring_indices, len(arr_flat))
 
     # the coords and offsets into the coordinates of the rings
     coords, indices = get_coordinates(arr_flat2, return_index=True, include_z=include_z)
-    offsets1 = np.insert(np.bincount(indices).cumsum(), 0, 0)
+    offsets1 = _indices_to_offsets(indices, len(arr_flat2))
 
     return coords, (offsets1, offsets2, offsets3)
 
@@ -179,6 +187,8 @@ def to_ragged_array(geometries, include_z=None):
             offset arrays depends on the geometry type. See
             https://github.com/geoarrow/geoarrow/blob/main/format.md
             for details.
+            Uses int32 dtype offsets if possible, otherwise int64 for
+            large inputs (coordinates > 32GB).
 
     Notes
     -----
@@ -233,7 +243,7 @@ def to_ragged_array(geometries, include_z=None):
            [ 2.,  2.]])
 
     >>> offsets
-    (array([0, 5, 9]), array([0, 2]))
+    (array([0, 5, 9], dtype=int32), array([0, 2], dtype=int32))
 
     As an example how to interpret the offsets: the i-th ring in the
     coordinates is represented by ``offsets[0][i]`` to ``offsets[0][i+1]``:


Test patch:
diff --git a/shapely/tests/test_constructive.py b/shapely/tests/test_constructive.py
index 0c6b5e2..89992d3 100644
--- a/shapely/tests/test_constructive.py
+++ b/shapely/tests/test_constructive.py
@@ -13,6 +13,7 @@ from shapely import (
     MultiPolygon,
     Point,
     Polygon,
+    geos_version,
 )
 from shapely.errors import UnsupportedGEOSVersionError
 from shapely.testing import assert_geometries_equal
@@ -30,6 +31,10 @@ from shapely.tests.common import (
     point_z,
 )
 
+geos39 = geos_version[0:2] == (3, 9)
+geos310 = geos_version[0:2] == (3, 10)
+geos311 = geos_version[0:2] == (3, 11)
+
 CONSTRUCTIVE_NO_ARGS = (
     shapely.boundary,
     shapely.centroid,
@@ -64,9 +69,12 @@ def test_no_args_array(geometry, func):
         geometry.is_empty
         and shapely.get_num_geometries(geometry) > 0
         and func is shapely.node
-        and shapely.geos_version < (3, 9, 3)
-    ):
-        pytest.xfail("GEOS < 3.9.3 crashes with empty geometries")  # GEOS GH-601
+        and (
+            (geos39 and geos_version < (3, 9, 3))
+            or (geos310 and geos_version < (3, 10, 3))
+        )
+    ):  # GEOS GH-601
+        pytest.xfail("GEOS < 3.9.3 or GEOS < 3.10.3 crashes with empty geometries")
     actual = func([geometry, geometry])
     assert actual.shape == (2,)
     assert actual[0] is None or isinstance(actual[0], Geometry)
@@ -630,8 +638,12 @@ def test_clip_by_rect_polygon(geom, rect, expected):
 def test_clip_by_rect_array(geometry):
     if (
         geometry.is_empty
-        and shapely.get_type_id(geometry) == 0
-        and shapely.geos_version < (3, 9, 5)
+        and shapely.get_type_id(geometry) == shapely.GeometryType.POINT
+        and (
+            (geos39 and geos_version < (3, 9, 5))
+            or (geos310 and geos_version < (3, 10, 6))
+            or (geos311 and geos_version < (3, 11, 3))
+        )
     ):
         # GEOS GH-913
         with pytest.raises(GEOSException):
@@ -1207,6 +1219,7 @@ def test_maximum_inscribed_circle_all_types(geometry):
         with pytest.raises(
             GEOSException,
             match=(
+                "Argument must be Polygonal or LinearRing|"  # GEOS < 3.10.4
                 "Input geometry must be a Polygon or MultiPolygon|"
                 "Operation not supported by GeometryCollection"
             ),
@@ -1247,7 +1260,11 @@ def test_maximum_inscribed_circle(geometry, expected):
 def test_maximum_inscribed_circle_empty():
     geometry = shapely.from_wkt("POINT EMPTY")
     with pytest.raises(
-        GEOSException, match="Input geometry must be a Polygon or MultiPolygon"
+        GEOSException,
+        match=(
+            "Argument must be Polygonal or LinearRing|"  # GEOS < 3.10.4
+            "Input geometry must be a Polygon or MultiPolygon"
+        ),
     ):
         shapely.maximum_inscribed_circle(geometry)
 
diff --git a/shapely/tests/test_io.py b/shapely/tests/test_io.py
index da71119..1b403cb 100644
--- a/shapely/tests/test_io.py
+++ b/shapely/tests/test_io.py
@@ -1319,7 +1319,7 @@ def test_to_geojson_exceptions():
         shapely.to_geojson(1)
 
 
-@pytest.mark.skipif(shapely.geos_version < (3, 10, 0), reason="GEOS < 3.10")
+@pytest.mark.skipif(shapely.geos_version < (3, 10, 2), reason="GEOS < 3.10.2")
 @pytest.mark.parametrize(
     "geom",
     [
@@ -1341,6 +1341,12 @@ def test_geojson_all_types(geom):
     type_id = shapely.get_type_id(geom)
     if type_id == shapely.GeometryType.LINEARRING:
         pytest.skip("Linearrings are not preserved in GeoJSON")
+    elif (
+        geom.is_empty
+        and type_id == shapely.GeometryType.POINT
+        and shapely.geos_version < (3, 10, 2)
+    ):
+        pytest.skip("GEOS < 3.10.2 with POINT EMPTY")  # TRAC-1139
     geojson = shapely.to_geojson(geom)
     actual = shapely.from_geojson(geojson)
     assert not actual.has_z
diff --git a/shapely/tests/test_ragged_array.py b/shapely/tests/test_ragged_array.py
index 8507152..023b16c 100644
--- a/shapely/tests/test_ragged_array.py
+++ b/shapely/tests/test_ragged_array.py
@@ -164,10 +164,11 @@ def test_linestrings():
             [10.0, 40.0],
         ]
     )
-    expected_offsets = np.array([0, 3, 7, 7, 7, 10, 10, 10])
+    expected_offsets = np.array([0, 3, 7, 7, 7, 10, 10, 10], dtype="int32")
     assert typ == shapely.GeometryType.LINESTRING
     assert_allclose(coords, expected)
     assert len(offsets) == 1
+    assert offsets[0].dtype == np.int32
     assert_allclose(offsets[0], expected_offsets)
 
     result = shapely.from_ragged_array(typ, coords, offsets)
@@ -218,6 +219,8 @@ def test_polygons():
     assert typ == shapely.GeometryType.POLYGON
     assert_allclose(coords, expected)
     assert len(offsets) == 2
+    assert offsets[0].dtype == np.int32
+    assert offsets[1].dtype == np.int32
     assert_allclose(offsets[0], expected_offsets1)
     assert_allclose(offsets[1], expected_offsets2)
 
@@ -257,6 +260,7 @@ def test_multipoints():
     assert typ == shapely.GeometryType.MULTIPOINT
     assert_allclose(coords, expected)
     assert len(offsets) == 1
+    assert offsets[0].dtype == np.int32
     assert_allclose(offsets[0], expected_offsets)
 
     result = shapely.from_ragged_array(typ, coords, offsets)
@@ -305,6 +309,8 @@ def test_multilinestrings():
     assert typ == shapely.GeometryType.MULTILINESTRING
     assert_allclose(coords, expected)
     assert len(offsets) == 2
+    assert offsets[0].dtype == np.int32
+    assert offsets[1].dtype == np.int32
     assert_allclose(offsets[0], expected_offsets1)
     assert_allclose(offsets[1], expected_offsets2)
 
@@ -365,6 +371,9 @@ def test_multipolygons():
     assert typ == shapely.GeometryType.MULTIPOLYGON
     assert_allclose(coords, expected)
     assert len(offsets) == 3
+    assert offsets[0].dtype == np.int32
+    assert offsets[1].dtype == np.int32
+    assert offsets[2].dtype == np.int32
     assert_allclose(offsets[0], expected_offsets1)
     assert_allclose(offsets[1], expected_offsets2)
     assert_allclose(offsets[2], expected_offsets3)