Repository: shapely/shapely (python)

Issue:
REGR: equality (__eq__) now ignores the z dimension
With shapely 2.0, we switched the implementation of `Geometry.__eq__` to use the "equals_exact" predicate from GEOS, instead of manually comparing the coordinates as we did before (`geom1.coords == geom2.coords`), in the assumption this is equivalent. But it seems that `GEOSEqualsExact` only checks the X and Y coordinates (https://github.com/libgeos/geos/issues/809), thus introducing a regression in comparing 3D geometries:

With shapely 1.8
```python
>>> Point(1, 2, 3) == Point(1, 2, 4)
False
```

With shapely 2.0:

```python
>>> Point(1, 2, 3) == Point(1, 2, 4)
True
```

Pull request description:
API: use equals_identical for __eq__
Closes #1732, closes https://github.com/shapely/shapely/issues/1775

This might be a bit overkill for _just_ implementing `__eq__`, but this should also allow to actually expose "equals_identical" as a ufunc for all GEOS versions (and not only for GEOS>=3.12). And it is also faster than the easier python alternative (but not sure how significant this will be in real-world cases):

```
geom = shapely.LineString(np.random.randn(1000, 2))

In [19]: %timeit geom == geom
10.8 µs ± 198 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

In [20]: %timeit np.allclose(shapely.get_coordinates(geom), shapely.get_coordinates(geom), rtol=0, equal_nan=True)
76.5 µs ± 599 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
```

Patch:
diff --git a/docs/release/2.x.rst b/docs/release/2.x.rst
index 37a9369..e53f7b4 100644
--- a/docs/release/2.x.rst
+++ b/docs/release/2.x.rst
@@ -6,6 +6,13 @@ 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.9, NumPy >= 1.20, and Python >= 3.9 (#1802, #1885, #2124)
diff --git a/setup.py b/setup.py
index 0574923..6995f84 100644
--- a/setup.py
+++ b/setup.py
@@ -165,6 +165,7 @@ else:
                 "src/geos.c",
                 "src/lib.c",
                 "src/pygeom.c",
+                "src/pygeos.c",
                 "src/strtree.c",
                 "src/ufuncs.c",
                 "src/vector.c",
diff --git a/shapely/geometry/base.py b/shapely/geometry/base.py
index aa7ca0e..d131c14 100644
--- a/shapely/geometry/base.py
+++ b/shapely/geometry/base.py
@@ -210,24 +210,6 @@ class BaseGeometry(shapely.Geometry):
         """Return the symmetric difference of the geometries."""
         return self.symmetric_difference(other)
 
-    def __eq__(self, other):
-        """Return True if geometries are equal, else False."""
-        if not isinstance(other, BaseGeometry):
-            return NotImplemented
-        return type(other) is type(self) and np.array_equal(
-            self.coords, other.coords, equal_nan=False
-        )
-
-    def __ne__(self, other):
-        """Return True if geometries are not equal, else False."""
-        if not isinstance(other, BaseGeometry):
-            return NotImplemented
-        return not self.__eq__(other)
-
-    def __hash__(self):
-        """Return the hash value of the geometry."""
-        return super().__hash__()
-
     # Coordinate access
     # -----------------
 
diff --git a/shapely/geometry/polygon.py b/shapely/geometry/polygon.py
index 071225a..add56e1 100644
--- a/shapely/geometry/polygon.py
+++ b/shapely/geometry/polygon.py
@@ -264,36 +264,6 @@ class Polygon(BaseGeometry):
             "Component rings have coordinate sequences, but the polygon does not"
         )
 
-    def __eq__(self, other):
-        """Test for equality with another geometry."""
-        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
-        return np.all(
-            [
-                np.array_equal(left, right, equal_nan=False)
-                for left, right in zip(my_coords, other_coords)
-            ]
-        )
-
-    def __hash__(self):
-        """Hash the polygon."""
-        return super().__hash__()
-
     @property
     def __geo_interface__(self):
         """Return a GeoJSON-like mapping of the Polygon geometry."""
diff --git a/shapely/predicates.py b/shapely/predicates.py
index 213c6e5..316fe15 100644
--- a/shapely/predicates.py
+++ b/shapely/predicates.py
@@ -36,6 +36,7 @@ __all__ = [
     "touches",
     "within",
     "equals_exact",
+    "equals_identical",
     "relate",
     "relate_pattern",
 ]
@@ -1074,6 +1075,51 @@ 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):
+    """Return True if the geometries are identical.
+
+    This function verifies whether 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 and M dimensions), 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__``).
+
+    .. versionadded:: 2.1
+
+    Parameters
+    ----------
+    a, b : Geometry or array_like
+        Geometry or geometries to check.
+    **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
+    >>> equals_identical(Point(1, 2, 3), Point(1, 2, 3))
+    True
+    >>> equals_identical(Point(1, 2, 3), Point(1, 2, 0))
+    False
+    """
+    return lib.equals_identical(a, b, **kwargs)
+
+
 def relate(a, b, **kwargs):
     """Return a string representation of the DE-9IM intersection matrix.
 
diff --git a/src/coords.c b/src/coords.c
index 6a17598..ddc0e82 100644
--- a/src/coords.c
+++ b/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;
diff --git a/src/geos.c b/src/geos.c
index 97a96f8..7a338a7 100644
--- a/src/geos.c
+++ b/src/geos.c
@@ -1156,16 +1156,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 has_z, int has_m) {
+
 #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, has_z, has_m);
 
 #else
 
   char *cp1, *cp2;
   unsigned int i, j;
+  unsigned int dims = 2 + has_z + has_m;
 
   cp1 = (char*)buf;
   for (i = 0; i < size; i++, cp1 += 8 * dims) {
diff --git a/src/geos.h b/src/geos.h
index bc2f8e1..9615a6a 100644
--- a/src/geos.h
+++ b/src/geos.h
@@ -201,6 +201,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 has_z, int has_m);
 
 #endif  // _GEOS_H
diff --git a/src/pygeom.c b/src/pygeom.c
index e6294fa..812c8ac 100644
--- a/src/pygeom.c
+++ b/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};
@@ -271,11 +272,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;
diff --git a/src/pygeos.c b/src/pygeos.c
new file mode 100644
index 0000000..f6adefe
--- /dev/null
+++ b/src/pygeos.c
@@ -0,0 +1,253 @@
+#define PY_SSIZE_T_CLEAN
+#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
+
+#include <Python.h>
+
+#define NO_IMPORT_ARRAY
+#define PY_ARRAY_UNIQUE_SYMBOL shapely_ARRAY_API
+
+#include <numpy/ndarraytypes.h>
+#include <numpy/npy_math.h>
+
+#include "geos.h"
+
+
+/* Check simple geometries (Point, LineString, LinearRing; backed by a single CoordSequence).
+Returns 1 on true, 0 on false, 2 on exception */
+char equals_identical_simple(GEOSContextHandle_t ctx, const GEOSGeometry* geom1,
+                             const GEOSGeometry* geom2) {
+  char hasZ1, hasZ2;
+  const GEOSCoordSequence *coord_seq1 = NULL;
+  const GEOSCoordSequence *coord_seq2 = NULL;
+  unsigned int n1, n2, dims1, dims2;
+  double* buf1 = NULL;
+  double* buf2 = NULL;
+  npy_intp i;
+  char ret;
+
+  hasZ1 = GEOSHasZ_r(ctx, geom1);
+  if (hasZ1 == 2) {
+    return 2;
+  }
+  hasZ2 = GEOSHasZ_r(ctx, geom2);
+  if (hasZ2 == 2) {
+    return 2;
+  }
+  if (hasZ1 != hasZ2) {
+    return 0;
+  }
+
+  coord_seq1 = GEOSGeom_getCoordSeq_r(ctx, geom1);
+  if (coord_seq1 == NULL) {
+    return 2;
+  }
+  coord_seq2 = GEOSGeom_getCoordSeq_r(ctx, geom2);
+  if (coord_seq2 == NULL) {
+    return 2;
+  }
+
+  if (GEOSCoordSeq_getSize_r(ctx, coord_seq1, &n1) == 0) {
+    return 2;
+  }
+  if (GEOSCoordSeq_getSize_r(ctx, coord_seq2, &n2) == 0) {
+    return 2;
+  }
+  if (n1 != n2) {
+    return 0;
+  }
+
+  if (GEOSCoordSeq_getDimensions_r(ctx, coord_seq1, &dims1) == 0) {
+    return 2;
+  }
+  if (GEOSCoordSeq_getDimensions_r(ctx, coord_seq2, &dims2) == 0) {
+    return 2;
+  }
+  if (dims1 != dims2) {
+    return 0;
+  }
+
+  int hasM = 0;
+  if ((dims1 == 3) & (hasZ1 == 0)) {
+    hasM = 1;
+  }
+  if (dims1 == 4) {
+    hasM = 1;
+  }
+
+  buf1 = malloc(sizeof(double) * n1 * dims1);
+  if (coordseq_to_buffer(ctx, coord_seq1, buf1, n1, hasZ1, hasM) == 0) {
+    free(buf1);
+    return 2;
+  }
+  buf2 = malloc(sizeof(double) * n1 * dims1);
+  if (coordseq_to_buffer(ctx, coord_seq2, buf2, n1, hasZ1, hasM) == 0) {
+    free(buf1);
+    free(buf2);
+    return 2;
+  }
+
+  ret = 1;
+  for (i = 0; i < n1 * dims1; i++) {
+    const double a = buf1[i];
+    const double b = buf2[i];
+    // NaN coordinates in same position are considered equal
+    if (a != b && !(npy_isnan(a) && npy_isnan(b))) {
+      ret = 0;
+      break;
+    }
+  }
+
+  free(buf1);
+  free(buf2);
+  return ret;
+}
+
+/* Check the exterior/interior rings of Polygons.
+Returns 1 on true, 0 on false, 2 on exception */
+char equals_identical_polygon(GEOSContextHandle_t ctx, const GEOSGeometry* geom1,
+                              const GEOSGeometry* geom2) {
+  int i, n1, n2;
+  const GEOSGeometry *ring1 = NULL;
+  const GEOSGeometry *ring2 = NULL;
+  char ret;
+
+  n1 = GEOSGetNumInteriorRings_r(ctx, geom1);
+  if (n1 == -1) {
+    return 2;
+  }
+  n2 = GEOSGetNumInteriorRings_r(ctx, geom2);
+  if (n2 == -1) {
+    return 2;
+  }
+  if (n1 != n2) {
+    return 0;
+  }
+
+  /* validate the exterior ring */
+  ring1 = GEOSGetExteriorRing_r(ctx, geom1);
+  if (ring1 == NULL) {
+    return 2;
+  }
+  ring2 = GEOSGetExteriorRing_r(ctx, geom2);
+  if (ring2 == NULL) {
+    return 2;
+  }
+  ret = equals_identical_simple(ctx, ring1, ring2);
+  if (ret != 1) {
+    return ret;
+  }
+  // we don't own the shells, so don't need to destroy them
+
+  /* validate the interior rings */
+  for (i = 0; i < n1; i++) {
+    ring1 = GEOSGetInteriorRingN_r(ctx, geom1, i);
+    if (ring1 == NULL) {
+      return 2;
+    }
+    ring2 = GEOSGetInteriorRingN_r(ctx, geom2, i);
+    if (ring2 == NULL) {
+      return 2;
+    }
+    ret = equals_identical_simple(ctx, ring1, ring2);
+    if (ret != 1) {
+      return ret;
+    }
+  }
+  return 1;
+}
+
+/* Note: forward declaration of the final function that is defined further below,
+so that it can be used in equals_identical_collection just below (due to the
+recursive nature of the implementation). */
+char equals_identical(GEOSContextHandle_t ctx, const GEOSGeometry* geom1,
+                      const GEOSGeometry* geom2);
+
+/* Check the sub-geometries of Multi-part geometries.
+Returns 1 on true, 0 on false, 2 on exception */
+char equals_identical_collection(GEOSContextHandle_t ctx, const GEOSGeometry* geom1,
+                                 const GEOSGeometry* geom2) {
+  int i, n1, n2;
+  const GEOSGeometry *sub_geom1 = NULL;
+  const GEOSGeometry *sub_geom2 = NULL;
+  char ret;
+
+  n1 = GEOSGetNumGeometries_r(ctx, geom1);
+  if (n1 == -1) {
+    return 2;
+  }
+  n2 = GEOSGetNumGeometries_r(ctx, geom2);
+  if (n2 == -1) {
+    return 2;
+  }
+  if (n1 != n2) {
+    return 0;
+  }
+
+  /* validate the sub geometries */
+  for (i = 0; i < n1; i++) {
+    sub_geom1 = GEOSGetGeometryN_r(ctx, geom1, i);
+    if (sub_geom1 == NULL) {
+      return 2;
+    }
+    sub_geom2 = GEOSGetGeometryN_r(ctx, geom2, i);
+    if (sub_geom2 == NULL) {
+      return 2;
+    }
+    ret = equals_identical(ctx, sub_geom1, sub_geom2);
+    if (ret != 1) {
+      return ret;
+    }
+  }
+  return 1;
+}
+
+/* Check geometries for being identical.
+Returns 1 on true, 0 on false, 2 on exception */
+char equals_identical(GEOSContextHandle_t ctx, const GEOSGeometry* geom1,
+                      const GEOSGeometry* geom2) {
+  int type1, type2;
+
+  type1 = GEOSGeomTypeId_r(ctx, geom1);
+  if (type1 == -1) {
+    return 2;
+  }
+  type2 = GEOSGeomTypeId_r(ctx, geom2);
+  if (type2 == -1) {
+    return 2;
+  }
+  if (type1 != type2) {
+    return 0;
+  }
+
+  switch (type1) {
+    case GEOS_POINT:
+    case GEOS_LINESTRING:
+    case GEOS_LINEARRING:
+      return equals_identical_simple(ctx, geom1, geom2);
+    case GEOS_POLYGON:
+      return equals_identical_polygon(ctx, geom1, geom2);
+    case GEOS_MULTIPOINT:
+    case GEOS_MULTILINESTRING:
+    case GEOS_MULTIPOLYGON:
+    case GEOS_GEOMETRYCOLLECTION:
+      return equals_identical_collection(ctx, geom1, geom2);
+    default:
+      return 2;
+  }
+}
+
+/*
+ * Determine pointwise equivalence of two geometries by checking
+ * that the structure, ordering, and values of all vertices are
+ * identical in all dimensions. NaN values are considered to be
+ * equal to other NaN values in the same position within coordinate sequences.
+ * \returns 1 on true, 0 on false, 2 on exception
+ */
+char PyGEOSEqualsIdentical(GEOSContextHandle_t ctx, const GEOSGeometry* geom1,
+                           const GEOSGeometry* geom2) {
+#if GEOS_SINCE_3_12_0
+  return GEOSEqualsIdentical_r(ctx, geom1, geom2);
+#else
+  return equals_identical(ctx, geom1, geom2);
+#endif
+}
diff --git a/src/pygeos.h b/src/pygeos.h
new file mode 100644
index 0000000..b0dde29
--- /dev/null
+++ b/src/pygeos.h
@@ -0,0 +1,8 @@
+#ifndef _PYGEOS_H
+#define _PYGEOS_H
+
+#include "geos.h"
+
+char PyGEOSEqualsIdentical(GEOSContextHandle_t ctx, const GEOSGeometry* geom1, const GEOSGeometry* geom2);
+
+#endif  // _PYGEOS_H
diff --git a/src/ufuncs.c b/src/ufuncs.c
index 6e31da5..0c96bb5 100644
--- a/src/ufuncs.c
+++ b/src/ufuncs.c
@@ -14,6 +14,7 @@
 
 #include "fast_loop_macros.h"
 #include "geos.h"
+#include "pygeos.h"
 #include "pygeom.h"
 
 /* This initializes a global value for interrupt checking */
@@ -252,6 +253,7 @@ static PyUFuncGenericFunction O_b_funcs[1] = {&O_b_func};
 
 /* Define the geom, geom -> bool functions (YY_b) */
 static void* equals_data[1] = {GEOSEquals_r};
+static void* equals_identical_data[1] = {PyGEOSEqualsIdentical};
 typedef char FuncGEOS_YY_b(void* context, void* a, void* b);
 static char YY_b_dtypes[3] = {NPY_OBJECT, NPY_OBJECT, NPY_BOOL};
 static void YY_b_func(char** args, const npy_intp* dimensions, const npy_intp* steps, void* data) {
@@ -3685,6 +3687,7 @@ int init_ufuncs(PyObject* m, PyObject* d) {
   DEFINE_YY_b_p(contains_properly);
   DEFINE_YY_b_p(overlaps);
   DEFINE_YY_b(equals);
+  DEFINE_YY_b(equals_identical);
   DEFINE_YY_b_p(covers);
   DEFINE_YY_b_p(covered_by);
   DEFINE_Ydd_b_p(contains_xy);


Test patch:
diff --git a/shapely/tests/geometry/test_collection.py b/shapely/tests/geometry/test_collection.py
index 606a667..95813fc 100644
--- a/shapely/tests/geometry/test_collection.py
+++ b/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():
diff --git a/shapely/tests/geometry/test_equality.py b/shapely/tests/geometry/test_equality.py
index da3bf27..26704a1 100644
--- a/shapely/tests/geometry/test_equality.py
+++ b/shapely/tests/geometry/test_equality.py
@@ -10,10 +10,19 @@ all_non_empty_types = np.array(all_types + all_types_z)[
 ]
 
 
+# TODO add all_types_m and all_types_zm once tranform supports M coordinates
 @pytest.mark.parametrize("geom", all_types + all_types_z)
 def test_equality(geom):
     assert geom == geom  # noqa: PLR0124
     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)
 
@@ -36,6 +45,21 @@ def test_equality(geom):
             MultiLineString([[(1, 1), (2, 2)], [(2, 2), (3, 3)]]),
             MultiLineString([[(2, 2), (3, 3)], [(1, 1), (2, 2)]]),
         ),
+        # M coordinates (don't work yet with automated cases)
+        pytest.param(
+            shapely.from_wkt("POINT M (0 0 0)"),
+            shapely.from_wkt("POINT M (0 0 1)"),
+            marks=pytest.mark.skipif(
+                shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12"
+            ),
+        ),
+        pytest.param(
+            shapely.from_wkt("POINT ZM (0 0 0 0)"),
+            shapely.from_wkt("POINT ZM (0 0 0 1)"),
+            marks=pytest.mark.skipif(
+                shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12"
+            ),
+        ),
     ],
 )
 def test_equality_false(left, right):
@@ -75,11 +99,8 @@ with ignore_invalid():
 
 @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 +118,8 @@ with ignore_invalid():
 
 @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():
@@ -112,6 +128,20 @@ with ignore_invalid():
         (LineString([(0, 1), (2, np.nan)]), LineString([(0, 1), (2, 3)])),
         (LineString([(0, 1, np.nan), (2, 3, 4)]), LineString([(0, 1, 2), (2, 3, 4)])),
         (LineString([(0, 1, 2), (2, 3, np.nan)]), LineString([(0, 1, 2), (2, 3, 4)])),
+        pytest.param(
+            shapely.from_wkt("POINT M (0 0 0)"),
+            shapely.from_wkt("POINT M (0 0 NaN)"),
+            marks=pytest.mark.skipif(
+                shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12"
+            ),
+        ),
+        pytest.param(
+            shapely.from_wkt("POINT ZM (0 0 0 0)"),
+            shapely.from_wkt("POINT ZM (0 0 0 NaN)"),
+            marks=pytest.mark.skipif(
+                shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12"
+            ),
+        ),
     ]
 
 
@@ -127,9 +157,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)
diff --git a/shapely/tests/test_geometry.py b/shapely/tests/test_geometry.py
index db0eea3..f4de7f0 100644
--- a/shapely/tests/test_geometry.py
+++ b/shapely/tests/test_geometry.py
@@ -246,11 +246,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():
diff --git a/shapely/tests/test_predicates.py b/shapely/tests/test_predicates.py
index 8ca08cb..99b5ccd 100644
--- a/shapely/tests/test_predicates.py
+++ b/shapely/tests/test_predicates.py
@@ -58,6 +58,7 @@ BINARY_PREDICATES = (
     ),
     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)