Repository: Toblerity/Shapely (python)

Issue:
NumPy 1.21 __array_interface__ behavior
NumPy 1.21 has a new behavior described in [this release note](https://github.com/numpy/numpy/releases/tag/v1.21.0):
> Exceptions will be raised during array-like creation: When an object raised an exception during access of the special
attributes `__array__` or `__array_interface__`, this exception was usually ignored. A warning is now given when the exception is anything but AttributeError. To silence the warning, the type raising the exception has to be adapted to raise an AttributeError.

Shapely<2.0 versions still have a `__array_interface__` property for most (all?) geometry types, and some of them, e.g. Polygon raise `NotImplementedError("A polygon does not itself provide the array interface. Its rings do.")`. This is problematic with this example:
```
import numpy as np
import warnings
from random import random
from shapely.geometry import Polygon
warnings.filterwarnings('always')

geoms = [Polygon((random(), random()) for _ in range(3)) for _ in range(3)]
arr = np.empty(len(geoms), dtype=object)
arr[:] = geoms
```
which will show this message three times (one for each geometry)
> `<stdin>`:1: DeprecationWarning: An exception was ignored while fetching the attribute `__array_interface__` from an object of type 'Polygon'.  With the exception of `AttributeError` NumPy will always raise this exception in the future.  Raise this deprecation warning to see the original exception. (Warning added NumPy 1.21)

This is also an issue for [GeoPandas](https://github.com/geopandas/geopandas/issues/2013), which will flood a terminal with warnings for large GeoSeries.

One very quick fix is to change the exception raised by `Polygon.__array_interface__` to `AttributeError`. Other geometries may also need to be adjusted, but the BaseGeometry class should still raise NotImplementedError. If this seems reasonable, I can prepare a PR for 1.8.

Pull request:
Silence NumPy 1.21 DeprecationWarning by raising AttributeError
Possible fix for NumPy 1.21 DeprecationWarning for exceptions raised in `__array_interface__`.  Only `BaseGeometry.__array_interface__` still raise NotImplementedError, as it is an abstract class.

Tests for each geometry type (except GeometryCollection) are added, which creates/assigns an object array.

Comments/reviews are welcome!

Closes #1173

Patch:
diff --git a/shapely/coords.py b/shapely/coords.py
index df673b0..ce526ef 100644
--- a/shapely/coords.py
+++ b/shapely/coords.py
@@ -116,7 +116,10 @@ class CoordinateSequence:
     def _ctypes(self):
         self._update()
         has_z = self._ndim == 3
-        n = self._ndim
+        n = self._ndim or 0
+        if n == 0:
+            # ignore with NumPy 1.21 __array_interface__
+            raise AttributeError("empty geometry sequence")
         m = self.__len__()
         array_type = c_double * (m * n)
         data = array_type()
diff --git a/shapely/geometry/base.py b/shapely/geometry/base.py
index 1ad0d07..2a5eacb 100644
--- a/shapely/geometry/base.py
+++ b/shapely/geometry/base.py
@@ -849,8 +849,8 @@ class BaseMultipartGeometry(BaseGeometry):
     @property
     def __array_interface__(self):
         """Provide the Numpy array protocol."""
-        raise NotImplementedError("Multi-part geometries do not themselves "
-                                  "provide the array interface")
+        raise AttributeError("Multi-part geometries do not themselves "
+                             "provide the array interface")
 
     def _get_coords(self):
         raise NotImplementedError("Sub-geometries may have coordinate "
diff --git a/shapely/geometry/polygon.py b/shapely/geometry/polygon.py
index 4386eb4..fabb6f2 100644
--- a/shapely/geometry/polygon.py
+++ b/shapely/geometry/polygon.py
@@ -317,7 +317,7 @@ class Polygon(BaseGeometry):
 
     @property
     def __array_interface__(self):
-        raise NotImplementedError(
+        raise AttributeError(
         "A polygon does not itself provide the array interface. Its rings do.")
 
     def _get_coords(self):
diff --git a/shapely/speedups/_speedups.pyx b/shapely/speedups/_speedups.pyx
index f5b6dc2..0b49959 100644
--- a/shapely/speedups/_speedups.pyx
+++ b/shapely/speedups/_speedups.pyx
@@ -421,7 +421,10 @@ def coordseq_ctypes(self):
     cdef GEOSCoordSequence *cs
     cdef double *data_p
     self._update()
-    n = self._ndim
+    n = self._ndim or 0
+    if n == 0:
+        # ignore with NumPy 1.21 __array_interface__
+        raise AttributeError("empty geometry sequence")
     m = self.__len__()
     array_type = ctypes.c_double * (m * n)
     data = array_type()
diff --git a/shapely/strtree.py b/shapely/strtree.py
index ac1008f..30c79a5 100644
--- a/shapely/strtree.py
+++ b/shapely/strtree.py
@@ -20,7 +20,7 @@ References
 
 import ctypes
 import logging
-from typing import Any, ItemsView, Iterable, Iterator, Sequence, Tuple, Union
+from typing import Any, ItemsView, Iterable, Iterator, Optional, Sequence, Tuple, Union
 import sys
 from warnings import warn
 
@@ -108,27 +108,24 @@ class STRtree:
         self.node_capacity = node_capacity
         self._rev = {
             item: geom
-            for geom, item in self._iterinitdata(
-                zip(geoms, items) if items is not None else geoms
-            )
+            for geom, item in self._iterinitdata(geoms, items)
             if not geom.is_empty
         }
-
         if self._rev:
             self._init_tree(self._rev.items())
 
     def _iterinitdata(
         self,
-        initdata: Union[Iterable[Tuple[BaseGeometry, Any]], Iterable[BaseGeometry]],
+        geoms: Iterable[BaseGeometry], items: Optional[Iterable[BaseGeometry]],
     ) -> Iterator[Tuple[BaseGeometry, Any]]:
-        if not initdata:
-            return
-
-        for enum_idx, item in enumerate(initdata):
-            if isinstance(item, tuple):
-                yield item[0], item[1]
-            elif isinstance(item, BaseGeometry):
-                yield (item, enum_idx)
+        if items is not None:
+            for geom, item in zip(geoms, items):
+                if isinstance(geom, BaseGeometry):
+                    yield (geom, item)
+        else:
+            for enum_idx, geom in enumerate(geoms):
+                if isinstance(geom, BaseGeometry):
+                    yield (geom, enum_idx)
 
     def _init_tree(self, rev_initdata: ItemsView[Any, BaseGeometry]):
         if rev_initdata:


Test patch:
diff --git a/tests/test_collection.py b/tests/test_collection.py
index 09b5c11..5ba176b 100644
--- a/tests/test_collection.py
+++ b/tests/test_collection.py
@@ -114,3 +114,14 @@ def test_len_deprecated(geometrycollection_geojson):
     geom = shape(geometrycollection_geojson)
     with pytest.warns(ShapelyDeprecationWarning, match="__len__"):
         assert len(geom) == 2
+
+
+@shapely20_deprecated
+@pytest.mark.filterwarnings("error:An exception was ignored")  # NumPy 1.21
+def test_numpy_object_array():
+    np = pytest.importorskip("numpy")
+
+    geom = GeometryCollection([LineString([(0, 0), (1, 1)])])
+    ar = np.empty(1, object)
+    ar[:] = [geom]
+    assert ar[0] == geom
diff --git a/tests/test_linestring.py b/tests/test_linestring.py
index 2b6e7f9..a452ff7 100644
--- a/tests/test_linestring.py
+++ b/tests/test_linestring.py
@@ -87,6 +87,7 @@ def test_from_numpy():
     assert line.coords[:] == [(1.0, 2.0), (3.0, 4.0)]
 
 
+@pytest.mark.filterwarnings("error:An exception was ignored")  # NumPy 1.21
 def test_numpy_empty_linestring_coords():
     np = pytest.importorskip("numpy")
 
@@ -97,6 +98,17 @@ def test_numpy_empty_linestring_coords():
     assert la.shape == (0,)
 
 
+@shapely20_deprecated
+@pytest.mark.filterwarnings("error:An exception was ignored")  # NumPy 1.21
+def test_numpy_object_array():
+    np = pytest.importorskip("numpy")
+
+    geom = LineString([(0.0, 0.0), (0.0, 1.0)])
+    ar = np.empty(1, object)
+    ar[:] = [geom]
+    assert ar[0] == geom
+
+
 def test_from_invalid_dim():
     with pytest.raises(ValueError, match="at least 2 coordinate tuples"):
         LineString([(1, 2)])
diff --git a/tests/test_multilinestring.py b/tests/test_multilinestring.py
index bf79ddf..c6475dc 100644
--- a/tests/test_multilinestring.py
+++ b/tests/test_multilinestring.py
@@ -117,3 +117,14 @@ def test_getitem_deprecated():
     geom = MultiLineString([[[5.0, 6.0], [7.0, 8.0]]])
     with pytest.warns(ShapelyDeprecationWarning, match="__getitem__"):
         part = geom[0]
+
+
+@shapely20_deprecated
+@pytest.mark.filterwarnings("error:An exception was ignored")  # NumPy 1.21
+def test_numpy_object_array():
+    np = pytest.importorskip("numpy")
+
+    geom = MultiLineString([[[5.0, 6.0], [7.0, 8.0]]])
+    ar = np.empty(1, object)
+    ar[:] = [geom]
+    assert ar[0] == geom
diff --git a/tests/test_multipoint.py b/tests/test_multipoint.py
index afbffaf..a3e2914 100644
--- a/tests/test_multipoint.py
+++ b/tests/test_multipoint.py
@@ -129,6 +129,17 @@ def test_multipoint_array_interface_numpy_deprecated():
         np.array(geom)
 
 
+@shapely20_deprecated
+@pytest.mark.filterwarnings("error:An exception was ignored")  # NumPy 1.21
+def test_numpy_object_array():
+    np = pytest.importorskip("numpy")
+
+    geom = MultiPoint(((1.0, 2.0), (3.0, 4.0)))
+    ar = np.empty(1, object)
+    ar[:] = [geom]
+    assert ar[0] == geom
+
+
 def test_iteration_deprecated():
     geom = MultiPoint([[5.0, 6.0], [7.0, 8.0]])
     with pytest.warns(ShapelyDeprecationWarning, match="Iteration"):
diff --git a/tests/test_multipolygon.py b/tests/test_multipolygon.py
index c3982a8..878211a 100644
--- a/tests/test_multipolygon.py
+++ b/tests/test_multipolygon.py
@@ -110,3 +110,16 @@ def test_getitem_deprecated():
           [((0.25, 0.25), (0.25, 0.5), (0.5, 0.5), (0.5, 0.25))])])
     with pytest.warns(ShapelyDeprecationWarning, match="__getitem__"):
         part = geom[0]
+
+
+@shapely20_deprecated
+@pytest.mark.filterwarnings("error:An exception was ignored")  # NumPy 1.21
+def test_numpy_object_array():
+    np = pytest.importorskip("numpy")
+
+    geom = MultiPolygon(
+        [(((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0)),
+          [((0.25, 0.25), (0.25, 0.5), (0.5, 0.5), (0.5, 0.25))])])
+    ar = np.empty(1, object)
+    ar[:] = [geom]
+    assert ar[0] == geom
diff --git a/tests/test_point.py b/tests/test_point.py
index 264118c..4a80564 100644
--- a/tests/test_point.py
+++ b/tests/test_point.py
@@ -262,6 +262,7 @@ def test_point_array_interface_numpy_deprecated():
         np.array(p)
 
 
+@pytest.mark.filterwarnings("error:An exception was ignored")  # NumPy 1.21
 def test_numpy_empty_point_coords():
     np = pytest.importorskip("numpy")
 
@@ -270,3 +271,14 @@ def test_numpy_empty_point_coords():
     # Access the coords
     a = np.asarray(pe.coords)
     assert a.shape == (0,)
+
+
+@shapely20_deprecated
+@pytest.mark.filterwarnings("error:An exception was ignored")  # NumPy 1.21
+def test_numpy_object_array():
+    np = pytest.importorskip("numpy")
+
+    geom = Point(3.0, 4.0)
+    ar = np.empty(1, object)
+    ar[:] = [geom]
+    assert ar[0] == geom
diff --git a/tests/test_polygon.py b/tests/test_polygon.py
index 863b623..62b2f96 100644
--- a/tests/test_polygon.py
+++ b/tests/test_polygon.py
@@ -109,6 +109,7 @@ def test_numpy_linearring_coords():
     assert_array_equal(ra, expected)
 
 
+@pytest.mark.filterwarnings("error:An exception was ignored")  # NumPy 1.21
 def test_numpy_empty_linearring_coords():
     np = pytest.importorskip("numpy")
 
@@ -116,6 +117,16 @@ def test_numpy_empty_linearring_coords():
     assert np.asarray(ring.coords).shape == (0,)
 
 
+@pytest.mark.filterwarnings("error:An exception was ignored")  # NumPy 1.21
+def test_numpy_object_array():
+    np = pytest.importorskip("numpy")
+
+    geom = Polygon([(0.0, 0.0), (0.0, 1.0), (1.0, 1.0)])
+    ar = np.empty(1, object)
+    ar[:] = [geom]
+    assert ar[0] == geom
+
+
 def test_polygon_from_coordinate_sequence():
     coords = [(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (0.0, 0.0)]
 
diff --git a/tests/test_strtree.py b/tests/test_strtree.py
index 55b61c4..e1402cc 100644
--- a/tests/test_strtree.py
+++ b/tests/test_strtree.py
@@ -7,14 +7,19 @@ import sys
 import pytest
 
 from shapely.errors import ShapelyDeprecationWarning
-from shapely.geometry import Point, Polygon
+from shapely.geometry import Point, Polygon, box
 from shapely.geos import geos_version
 from shapely import strtree
 from shapely.strtree import STRtree
+from shapely import wkt
 
 from .conftest import requires_geos_342
 
 
+point = Point(2, 3)
+empty = wkt.loads("GEOMETRYCOLLECTION EMPTY")
+
+
 @requires_geos_342
 @pytest.mark.parametrize("geoms", [[Point(i, i) for i in range(10)]])
 @pytest.mark.parametrize(
@@ -37,11 +42,43 @@ def test_query(geoms, query_geom, num_results):
 def test_query_enumeration_idx(geoms, query_geom, expected):
     """Store enumeration idx"""
     with pytest.warns(ShapelyDeprecationWarning):
-        tree = STRtree((g, i) for i, g in enumerate(geoms))
+        tree = STRtree(geoms, range(len(geoms)))
     results = tree.query_items(query_geom)
     assert sorted(results) == sorted(expected)
 
 
+@requires_geos_342
+@pytest.mark.parametrize("geoms", [[Point(i, i) for i in range(5)]])
+@pytest.mark.parametrize("items", [None, list(range(1, 6)), list("abcde")])
+@pytest.mark.parametrize(
+    "query_geom,expected",
+    [(Point(2, 2).buffer(0.99), [2]), (Point(2, 2).buffer(1.0), [1, 2, 3])],
+)
+def test_query_items(geoms, items, query_geom, expected):
+    """Store enumeration idx"""
+    with pytest.warns(ShapelyDeprecationWarning):
+        tree = STRtree(geoms, items)
+    results = tree.query_items(query_geom)
+    expected = [items[idx] for idx in expected] if items is not None else expected
+    assert sorted(results) == sorted(expected)
+
+
+@pytest.mark.parametrize(
+    "tree_geometry, geometry,expected",
+    [
+        ([point], box(0, 0, 10, 10), [0]),
+        # None/empty is ignored in the tree, but the index of the valid geometry
+        # should be retained.
+        ([None, point], box(0, 0, 10, 10), [1]),
+        ([None, empty, point], box(0, 0, 10, 10), [2]),
+    ],
+)
+def test_query_items_with_empty(tree_geometry, geometry, expected):
+    with pytest.warns(ShapelyDeprecationWarning):
+        tree = STRtree(tree_geometry)
+    assert tree.query_items(geometry) == expected
+
+
 @requires_geos_342
 def test_insert_empty_geometry():
     """
@@ -112,7 +149,7 @@ def test_pickle_persistence():
     Don't crash trying to use unpickled GEOS handle.
     """
     with pytest.warns(ShapelyDeprecationWarning):
-        tree = STRtree([(Point(i, i).buffer(0.1), i) for i in range(3)])
+        tree = STRtree([Point(i, i).buffer(0.1) for i in range(3)], range(3))
 
     pickled_strtree = pickle.dumps(tree)
     unpickle_script_file_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), "unpickle-strtree.py")
@@ -165,14 +202,14 @@ def test_nearest_item(geoms, items, query_geom):
     assert tree.nearest_item(query_geom) == items[0]
 
 
-@pytest.mark.parametrize(["geoms", "items"], [(None, None), ([], None)])
+@pytest.mark.parametrize(["geoms", "items"], [([], None), ([], [])])
 def test_nearest_empty(geoms, items):
     with pytest.warns(ShapelyDeprecationWarning):
         tree = STRtree(geoms, items)
     assert tree.nearest_item(None) is None
 
 
-@pytest.mark.parametrize(["geoms", "items"], [(None, None), ([], None)])
+@pytest.mark.parametrize(["geoms", "items"], [([], None), ([], [])])
 def test_nearest_items(geoms, items):
     with pytest.warns(ShapelyDeprecationWarning):
         tree = STRtree(geoms, items)