Merge 0b12930 into 42392bb

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/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

@@ -10,6 +10,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)
 
@@ -68,11 +76,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():
@@ -90,13 +95,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():
@@ -120,9 +120,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

@@ -239,11 +239,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():

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

@@ -1093,16 +1093,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

@@ -193,6 +193,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};
@@ -272,11 +273,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;