Merge 3019a1b into 0772979

docs/release/2.x.rst

@@ -6,9 +6,16 @@ Version 2.x
 Version 2.1.0 (unreleased)
 --------------------------
 
+API changes:
+
+- Equality of geometries (``geom1 == geom2``) now considers NaN coordinate values
+  in the same location to be equal (#1775). It is recommended however to ensure
+  geometries don't have NaN values in the first place, for which you can now use
+  the ``handle_nan`` parameter in construction functions.
+
 Improvements:
 
-- Require GEOS >= 3.7, NumPy >= 1.16, and Python >= 3.8 (#1802)
+- Require GEOS >= 3.8, NumPy >= 1.16, and Python >= 3.8 (#1802, #1885)
 - 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).

setup.py

@@ -167,6 +167,7 @@ def finalize_options(self):
                 "src/geos.c",
                 "src/lib.c",
                 "src/pygeom.c",
+                "src/pygeos.c",
                 "src/strtree.c",
                 "src/ufuncs.c",
                 "src/vector.c",

shapely/geometry/base.py

@@ -197,23 +197,6 @@ def __sub__(self, other):
     def __xor__(self, other):
         return self.symmetric_difference(other)
 
-    def __eq__(self, other):
-        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) == type(self) and np.array_equal(
-            self.coords, other.coords  # , equal_nan=False
-        )
-
-    def __ne__(self, other):
-        if not isinstance(other, BaseGeometry):
-            return NotImplemented
-        return not self.__eq__(other)
-
-    def __hash__(self):
-        return super().__hash__()
-
     # Coordinate access
     # -----------------
 

shapely/geometry/polygon.py

@@ -257,35 +257,6 @@ def coords(self):
             "Component rings have coordinate sequences, but the polygon does not"
         )
 
-    def __eq__(self, other):
-        if not isinstance(other, BaseGeometry):
-            return NotImplemented
-        if not isinstance(other, Polygon):
-            return False
-        check_empty = (self.is_empty, other.is_empty)
-        if all(check_empty):
-            return True
-        elif any(check_empty):
-            return False
-        my_coords = [self.exterior.coords] + [
-            interior.coords for interior in self.interiors
-        ]
-        other_coords = [other.exterior.coords] + [
-            interior.coords for interior in other.interiors
-        ]
-        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)
-                for left, right in zip(my_coords, other_coords)
-            ]
-        )
-
-    def __hash__(self):
-        return super().__hash__()
-
     @property
     def __geo_interface__(self):
         if self.exterior == LinearRing():

shapely/predicates.py

@@ -34,6 +34,7 @@
     "touches",
     "within",
     "equals_exact",
+    "equals_identical",
     "relate",
     "relate_pattern",
 ]
@@ -1015,6 +1016,46 @@ def equals_exact(a, b, tolerance=0.0, normalize=False, **kwargs):
     return lib.equals_exact(a, b, tolerance, **kwargs)
 
 
+@multithreading_enabled
+def equals_identical(a, b, **kwargs):
+    """Returns True if the geometries are pointwise equivalent by checking
+    that the structure, ordering, and values of all vertices are identical
+    in all dimensions.
+
+    Similarly to :func:`equals_exact`, this function uses exact coordinate
+    equality and requires coordinates to be in the same order for all
+    components (vertices, rings, or parts) of a geometry. However, in
+    contrast :func:`equals_exact`, this function does not allow to specify
+    a tolerance, but does require all dimensions to be the same
+    (:func:`equals_exact` ignores the z-dimension), and NaN values are
+    considered to be equal to other NaN values.
+
+    This function is the vectorized equivalent of scalar equality of
+    geometry objects (``a == b``, i.e. ``__eq__``).
+
+    Parameters
+    ----------
+    a, b : Geometry or array_like
+    **kwargs
+        See :ref:`NumPy ufunc docs <ufuncs.kwargs>` for other keyword arguments.
+
+    See Also
+    --------
+    equals_exact : Check if geometries are structurally equal given a specified
+        tolerance.
+    equals : Check if geometries are spatially (topologically) equal.
+
+    Examples
+    --------
+    >>> from shapely import Point
+    >>> point1 = Point(50, 50, 1)
+    >>> point2 = Point(50, 50, 1)
+    >>> equals_identical(point1, point2)
+    True
+    """
+    return lib.equals_identical(a, b, **kwargs)
+
+
 def relate(a, b, **kwargs):
     """
     Returns a string representation of the DE-9IM intersection matrix.

shapely/tests/geometry/test_collection.py

@@ -1,6 +1,7 @@
 import numpy as np
 import pytest
 
+import shapely
 from shapely import GeometryCollection, LineString, Point, wkt
 from shapely.geometry import shape
 
@@ -36,7 +37,11 @@ def test_empty_subgeoms():
     assert geom.geom_type == "GeometryCollection"
     assert geom.is_empty
     assert len(geom.geoms) == 2
-    assert list(geom.geoms) == [Point(), LineString()]
+    parts = list(geom.geoms)
+    if shapely.geos_version < (3, 9, 0):
+        # the accessed empty 2D point has a 3D coordseq on GEOS 3.8
+        parts[0] = shapely.force_2d(parts[0])
+    assert parts == [Point(), LineString()]
 
 
 def test_child_with_deleted_parent():

shapely/tests/geometry/test_equality.py

@@ -14,6 +14,14 @@
 def test_equality(geom):
     assert geom == geom
     transformed = shapely.transform(geom, lambda x: x, include_z=True)
+    if (
+        shapely.geos_version < (3, 9, 0)
+        and isinstance(geom, Point)
+        and geom.is_empty
+        and not geom.has_z
+    ):
+        # the transformed empty 2D point has become 3D on GEOS 3.8
+        transformed = shapely.force_2d(geom)
     assert geom == transformed
     assert not (geom != transformed)
 
@@ -75,11 +83,8 @@ def test_equality_false(left, right):
 
 @pytest.mark.parametrize("left, right", cases1)
 def test_equality_with_nan(left, right):
-    # TODO currently those evaluate as not equal, but we are considering to change this
-    # assert left == right
-    assert not (left == right)
-    # assert not (left != right)
-    assert left != right
+    assert left == right
+    assert not (left != right)
 
 
 with ignore_invalid():
@@ -97,13 +102,8 @@ def test_equality_with_nan(left, right):
 
 @pytest.mark.parametrize("left, right", cases2)
 def test_equality_with_nan_z(left, right):
-    # TODO: those are currently considered equal because z dimension is ignored
-    if shapely.geos_version < (3, 12, 0):
-        assert left == right
-        assert not (left != right)
-    else:
-        # on GEOS main z dimension is not ignored -> NaNs cause inequality
-        assert left != right
+    assert left == right
+    assert not (left != right)
 
 
 with ignore_invalid():
@@ -127,9 +127,7 @@ def test_equality_with_nan_z_false():
 
     if shapely.geos_version < (3, 10, 0):
         # GEOS <= 3.9 fill the NaN with 0, so the z dimension is different
-        # assert left != right
-        # however, has_z still returns False, so z dimension is ignored in .coords
-        assert left == right
+        assert left != right
     elif shapely.geos_version < (3, 12, 0):
         # GEOS 3.10-3.11 ignore NaN for Z also when explicitly created with 3D
         # and so the geometries are considered as 2D (and thus z dimension is ignored)

shapely/tests/test_geometry.py

@@ -247,11 +247,11 @@ def test_set_unique(geom):
 
 
 def test_set_nan():
-    # As NaN != NaN, you can have multiple "NaN" points in a set
-    # set([float("nan"), float("nan")]) also returns a set with 2 elements
+    # Although NaN != NaN, you cannot have multiple "NaN" points in a set
+    # This is because "NaN" coordinates in a geometry are considered as equal.
     with ignore_invalid():
         a = set(shapely.linestrings([[[np.nan, np.nan], [np.nan, np.nan]]] * 10))
-    assert len(a) == 10  # different objects: NaN != NaN
+    assert len(a) == 1  # same objects: NaN == NaN (as geometry coordinates)
 
 
 def test_set_nan_same_objects():

shapely/tests/test_predicates.py

@@ -58,6 +58,7 @@
     ),
     shapely.equals,
     shapely.equals_exact,
+    shapely.equals_identical,
 )
 
 BINARY_PREPARED_PREDICATES = BINARY_PREDICATES[:-2]
@@ -210,6 +211,27 @@ def test_equals_exact_normalize():
     assert shapely.equals_exact(l1, l2, normalize=True)
 
 
+def test_equals_identical():
+    # more elaborate tests are done at the Geometry.__eq__ level
+    # requires same order of coordinates
+    l1 = LineString([(0, 0), (1, 1)])
+    l2 = LineString([(1, 1), (0, 0)])
+    assert not shapely.equals_identical(l1, l2)
+
+    # checks z-dimension (in contrast to equals_exact)
+    l1 = LineString([(0, 0, 0), (1, 1, 0)])
+    l2 = LineString([(0, 0, 1), (1, 1, 1)])
+    assert not shapely.equals_identical(l1, l2)
+    assert shapely.equals_exact(l1, l2)
+
+    # NaNs in same place are equal (in contrast to equals_exact)
+    with ignore_invalid():
+        l1 = LineString([(0, np.nan), (1, 1)])
+        l2 = LineString([(0, np.nan), (1, 1)])
+    assert shapely.equals_identical(l1, l2)
+    assert not shapely.equals_exact(l1, l2)
+
+
 @pytest.mark.skipif(shapely.geos_version < (3, 10, 0), reason="GEOS < 3.10")
 def test_dwithin():
     p1 = shapely.points(50, 4)

src/coords.c

@@ -25,7 +25,7 @@ position `cursor` in the array `out`. Increases the cursor correspondingly.
 Returns 0 on error, 1 on success */
 static char get_coordinates_simple(GEOSContextHandle_t ctx, GEOSGeometry* geom, int type,
                                    PyArrayObject* out, npy_intp* cursor, int include_z) {
-  unsigned int n, dims;
+  unsigned int n;
   double *buf;
   const GEOSCoordSequence* seq;
   char is_empty;
@@ -50,9 +50,8 @@ static char get_coordinates_simple(GEOSContextHandle_t ctx, GEOSGeometry* geom,
     return 0;
   }
 
-  dims = (include_z) ? 3 : 2;
   buf = PyArray_GETPTR2(out, *cursor, 0);
-  if (!coordseq_to_buffer(ctx, seq, buf, n, dims)) {
+  if (!coordseq_to_buffer(ctx, seq, buf, n, include_z, 0)) {
     return 0;
   }
   *cursor += n;

src/geos.c

@@ -1223,16 +1223,17 @@ enum ShapelyErrorCode coordseq_from_buffer(GEOSContextHandle_t ctx, const double
  * Returns 0 on error, 1 on success.
  */
 int coordseq_to_buffer(GEOSContextHandle_t ctx, const GEOSCoordSequence* coord_seq,
-                       double* buf, unsigned int size, unsigned int dims) {
+                       double* buf, unsigned int size, int hasZ, int hasM) {
+
 #if GEOS_SINCE_3_10_0
 
-  int hasZ = dims == 3;
-  return GEOSCoordSeq_copyToBuffer_r(ctx, coord_seq, buf, hasZ, 0);
+  return GEOSCoordSeq_copyToBuffer_r(ctx, coord_seq, buf, hasZ, hasM);
 
 #else
 
   char *cp1, *cp2;
   unsigned int i, j;
+  unsigned int dims = 2 + hasZ + hasM;
 
   cp1 = (char*)buf;
   for (i = 0; i < size; i++, cp1 += 8 * dims) {

src/geos.h

@@ -203,6 +203,6 @@ extern enum ShapelyErrorCode coordseq_from_buffer(GEOSContextHandle_t ctx,
                                                   npy_intp cs2,
                                 GEOSCoordSequence** coord_seq);
 extern int coordseq_to_buffer(GEOSContextHandle_t ctx, const GEOSCoordSequence* coord_seq,
-                              double* buf, unsigned int size, unsigned int dims);
+                              double* buf, unsigned int size, int hasZ, int hasM);
 
 #endif  // _GEOS_H

src/pygeom.c

@@ -7,6 +7,7 @@
 #include <structmember.h>
 
 #include "geos.h"
+#include "pygeos.h"
 
 /* This initializes a global geometry type registry */
 PyObject* geom_registry[1] = {NULL};
@@ -292,11 +293,11 @@ static PyObject* GeometryObject_richcompare(GeometryObject* self, PyObject* othe
         break;
       case Py_EQ:
         result =
-            GEOSEqualsExact_r(ctx, self->ptr, other_geom->ptr, 0) ? Py_True : Py_False;
+            PyGEOSEqualsIdentical(ctx, self->ptr, other_geom->ptr) ? Py_True : Py_False;
         break;
       case Py_NE:
         result =
-            GEOSEqualsExact_r(ctx, self->ptr, other_geom->ptr, 0) ? Py_False : Py_True;
+            PyGEOSEqualsIdentical(ctx, self->ptr, other_geom->ptr) ? Py_False : Py_True;
         break;
       case Py_GT:
         result = Py_NotImplemented;