Fix WKB/WKT of empty (multi)points on GEOS 3.9 (#392)

CHANGELOG.rst

@@ -24,6 +24,8 @@ Version 0.11 (unreleased)
 **Bug fixes**
 
 * Return True instead of False for LINEARRING geometries in ``is_closed`` (#379).
+* Fixed the WKB serialization of 3D empty points for GEOS >= 3.9.0 (#392).
+* Fixed the WKT serialization of multipoints with empty points for GEOS >= 3.9.0 (#392).
 
 **Acknowledgments**
 

docs/geometry.rst

@@ -35,7 +35,7 @@ Geometries can be serialized using pickle:
   >>> pickle.loads(point_1)
   <pygeos.Geometry POINT (5.2 52.1)>
 
-.. warning:: Pickling will convert linearrings to linestrings and it may convert empty points.
+.. warning:: Pickling will convert linearrings to linestrings.
              See :func:`pygeos.io.to_wkb` for a complete list of limitations.
 
 Hashing
@@ -53,7 +53,7 @@ Therefore, geometries are equal if and only if their WKB representations are equ
   {<pygeos.Geometry POINT (5.2 52.1)>, <pygeos.Geometry POINT (1 1)>}
 
 .. warning:: Due to limitations of WKB, linearrings will equal linestrings if they contain the exact same points.
-             Also, different kind of empty points will be considered equal. See :func:`pygeos.io.to_wkb`.
+             See :func:`pygeos.io.to_wkb`.
 
 Comparing two geometries directly is also supported.
 This is the same as using :func:`pygeos.predicates.equals_exact` with a ``tolerance`` value of zero.

pygeos/io.py

@@ -174,8 +174,9 @@ def to_wkb(
 
     - linearrings will be converted to linestrings
     - a point with only NaN coordinates is converted to an empty point
-    - empty points are transformed to 3D in GEOS < 3.8
-    - empty points are transformed to 2D in GEOS 3.8
+    - for GEOS <= 3.7, empty points are always serialized to 3D if
+      output_dimension=3, and to 2D if output_dimension=2
+    - for GEOS == 3.8, empty points are always serialized to 2D
 
     Parameters
     ----------

pygeos/tests/common.py

@@ -34,6 +34,7 @@
     "POLYGON((0 0, 0 10, 10 10, 10 0, 0 0), (2 2, 2 4, 4 4, 4 2, 2 2))"
 )
 empty_point = pygeos.Geometry("POINT EMPTY")
+empty_point_z = pygeos.Geometry("POINT Z EMPTY")
 empty_line_string = pygeos.Geometry("LINESTRING EMPTY")
 empty_polygon = pygeos.Geometry("POLYGON EMPTY")
 empty = pygeos.Geometry("GEOMETRYCOLLECTION EMPTY")

pygeos/tests/test_io.py

@@ -7,7 +7,7 @@
 
 import pygeos
 
-from .common import all_types, empty_point, point, point_z
+from .common import all_types, empty_point, empty_point_z, point, point_z
 
 # fmt: off
 POINT11_WKB = b"\x01\x01\x00\x00\x00" + struct.pack("<2d", 1.0, 1.0)
@@ -283,8 +283,21 @@ def test_to_wkt_geometrycollection_with_point_empty():
     assert pygeos.to_wkt(collection).endswith("(POINT EMPTY, POINT (2 3))")
 
 
+@pytest.mark.skipif(
+    pygeos.geos_version < (3, 9, 0),
+    reason="MULTIPOINT (EMPTY, 2 3) only works for GEOS >= 3.9",
+)
+def test_to_wkt_multipoint_with_point_empty():
+    geom = pygeos.multipoints([empty_point, point])
+    assert pygeos.to_wkt(geom) == "MULTIPOINT (EMPTY, 2 3)"
+
+
+@pytest.mark.skipif(
+    pygeos.geos_version >= (3, 9, 0),
+    reason="MULTIPOINT (EMPTY, 2 3) gives ValueError on GEOS < 3.9",
+)
 def test_to_wkt_multipoint_with_point_empty_errors():
-    # Test if segfault is prevented
+    # test if segfault is prevented
     geom = pygeos.multipoints([empty_point, point])
     with pytest.raises(ValueError):
         pygeos.to_wkt(geom)
@@ -302,6 +315,10 @@ def test_repr_max_length():
     assert representation.endswith("...>")
 
 
+@pytest.mark.skipif(
+    pygeos.geos_version >= (3, 9, 0),
+    reason="MULTIPOINT (EMPTY, 2 3) gives Exception on GEOS < 3.9",
+)
 def test_repr_multipoint_with_point_empty():
     # Test if segfault is prevented
     geom = pygeos.multipoints([point, empty_point])
@@ -374,88 +391,116 @@ def test_to_wkb_srid():
     assert np.frombuffer(result[5:9], "<u4").item() == 4326
 
 
-@pytest.mark.skipif(
-    pygeos.geos_version >= (3, 8, 0), reason="Pre GEOS 3.8.0 has 3D empty points"
-)
 @pytest.mark.parametrize(
-    "geom,dims,expected",
+    "geom,expected",
     [
-        (empty_point, 2, POINT_NAN_WKB),
-        (empty_point, 3, POINTZ_NAN_WKB),
-        (pygeos.multipoints([empty_point]), 2, MULTIPOINT_NAN_WKB),
-        (pygeos.multipoints([empty_point]), 3, MULTIPOINTZ_NAN_WKB),
-        (pygeos.geometrycollections([empty_point]), 2, GEOMETRYCOLLECTION_NAN_WKB),
-        (pygeos.geometrycollections([empty_point]), 3, GEOMETRYCOLLECTIONZ_NAN_WKB),
+        (empty_point, POINT_NAN_WKB),
+        (empty_point_z, POINT_NAN_WKB),
+        (pygeos.multipoints([empty_point]), MULTIPOINT_NAN_WKB),
+        (pygeos.multipoints([empty_point_z]), MULTIPOINT_NAN_WKB),
+        (pygeos.geometrycollections([empty_point]), GEOMETRYCOLLECTION_NAN_WKB),
+        (pygeos.geometrycollections([empty_point_z]), GEOMETRYCOLLECTION_NAN_WKB),
         (
             pygeos.geometrycollections([pygeos.multipoints([empty_point])]),
-            2,
             NESTED_COLLECTION_NAN_WKB,
         ),
         (
-            pygeos.geometrycollections([pygeos.multipoints([empty_point])]),
-            3,
-            NESTED_COLLECTIONZ_NAN_WKB,
+            pygeos.geometrycollections([pygeos.multipoints([empty_point_z])]),
+            NESTED_COLLECTION_NAN_WKB,
         ),
     ],
 )
-def test_to_wkb_point_empty_pre_geos38(geom, dims, expected):
-    actual = pygeos.to_wkb(geom, output_dimension=dims, byte_order=1)
-    # Use numpy.isnan; there are many byte representations for NaN
-    assert actual[: -dims * 8] == expected[: -dims * 8]
-    assert np.isnan(struct.unpack("<{}d".format(dims), actual[-dims * 8 :])).all()
-
-
-@pytest.mark.skipif(
-    pygeos.geos_version < (3, 8, 0), reason="Post GEOS 3.8.0 has 2D empty points"
+def test_to_wkb_point_empty_2d(geom, expected):
+    actual = pygeos.to_wkb(geom, output_dimension=2, byte_order=1)
+    # Split 'actual' into header and coordinates
+    coordinate_length = 16
+    header_length = len(expected) - coordinate_length
+    # Check the total length (this checks the correct dimensionality)
+    assert len(actual) == header_length + coordinate_length
+    # Check the header
+    assert actual[:header_length] == expected[:header_length]
+    # Check the coordinates (using numpy.isnan; there are many byte representations for NaN)
+    assert np.isnan(struct.unpack("<2d", actual[header_length:])).all()
+
+
+@pytest.mark.xfail(
+    pygeos.geos_version[:2] == (3, 8), reason="GEOS==3.8 never outputs 3D empty points"
 )
 @pytest.mark.parametrize(
-    "geom,dims,expected",
+    "geom,expected",
     [
-        (empty_point, 2, POINT_NAN_WKB),
-        (empty_point, 3, POINT_NAN_WKB),
-        (pygeos.multipoints([empty_point]), 2, MULTIPOINT_NAN_WKB),
-        (pygeos.multipoints([empty_point]), 3, MULTIPOINT_NAN_WKB),
-        (pygeos.geometrycollections([empty_point]), 2, GEOMETRYCOLLECTION_NAN_WKB),
-        (pygeos.geometrycollections([empty_point]), 3, GEOMETRYCOLLECTION_NAN_WKB),
+        (empty_point_z, POINTZ_NAN_WKB),
+        (pygeos.multipoints([empty_point_z]), MULTIPOINTZ_NAN_WKB),
+        (pygeos.geometrycollections([empty_point_z]), GEOMETRYCOLLECTIONZ_NAN_WKB),
         (
-            pygeos.geometrycollections([pygeos.multipoints([empty_point])]),
-            2,
-            NESTED_COLLECTION_NAN_WKB,
+            pygeos.geometrycollections([pygeos.multipoints([empty_point_z])]),
+            NESTED_COLLECTIONZ_NAN_WKB,
         ),
+    ],
+)
+def test_to_wkb_point_empty_3d(geom, expected):
+    actual = pygeos.to_wkb(geom, output_dimension=3, byte_order=1)
+    # Split 'actual' into header and coordinates
+    coordinate_length = 24
+    header_length = len(expected) - coordinate_length
+    # Check the total length (this checks the correct dimensionality)
+    assert len(actual) == header_length + coordinate_length
+    # Check the header
+    assert actual[:header_length] == expected[:header_length]
+    # Check the coordinates (using numpy.isnan; there are many byte representations for NaN)
+    assert np.isnan(struct.unpack("<3d", actual[header_length:])).all()
+
+
+@pytest.mark.xfail(
+    pygeos.geos_version < (3, 8, 0),
+    reason="GEOS<3.8 always outputs 3D empty points if output_dimension=3",
+)
+@pytest.mark.parametrize(
+    "geom,expected",
+    [
+        (empty_point, POINT_NAN_WKB),
+        (pygeos.multipoints([empty_point]), MULTIPOINT_NAN_WKB),
+        (pygeos.geometrycollections([empty_point]), GEOMETRYCOLLECTION_NAN_WKB),
         (
             pygeos.geometrycollections([pygeos.multipoints([empty_point])]),
-            3,
             NESTED_COLLECTION_NAN_WKB,
         ),
     ],
 )
-def test_to_wkb_point_empty_post_geos38(geom, dims, expected):
-    # Post GEOS 3.8: empty point is 2D
-    actual = pygeos.to_wkb(geom, output_dimension=dims, byte_order=1)
-    # Use numpy.isnan; there are many byte representations for NaN
-    assert actual[: -2 * 8] == expected[: -2 * 8]
-    assert np.isnan(struct.unpack("<2d", actual[-2 * 8 :])).all()
+def test_to_wkb_point_empty_2d_output_dim_3(geom, expected):
+    actual = pygeos.to_wkb(geom, output_dimension=3, byte_order=1)
+    # Split 'actual' into header and coordinates
+    coordinate_length = 16
+    header_length = len(expected) - coordinate_length
+    # Check the total length (this checks the correct dimensionality)
+    assert len(actual) == header_length + coordinate_length
+    # Check the header
+    assert actual[:header_length] == expected[:header_length]
+    # Check the coordinates (using numpy.isnan; there are many byte representations for NaN)
+    assert np.isnan(struct.unpack("<2d", actual[header_length:])).all()
 
 
 @pytest.mark.parametrize(
-    "wkb,expected_type",
+    "wkb,expected_type,expected_dim",
     [
-        (POINT_NAN_WKB, 0),
-        (POINTZ_NAN_WKB, 0),
-        (MULTIPOINT_NAN_WKB, 4),
-        (MULTIPOINTZ_NAN_WKB, 4),
-        (GEOMETRYCOLLECTION_NAN_WKB, 7),
-        (GEOMETRYCOLLECTIONZ_NAN_WKB, 7),
-        (NESTED_COLLECTION_NAN_WKB, 7),
-        (NESTED_COLLECTIONZ_NAN_WKB, 7),
+        (POINT_NAN_WKB, 0, 2),
+        (POINTZ_NAN_WKB, 0, 3),
+        (MULTIPOINT_NAN_WKB, 4, 2),
+        (MULTIPOINTZ_NAN_WKB, 4, 3),
+        (GEOMETRYCOLLECTION_NAN_WKB, 7, 2),
+        (GEOMETRYCOLLECTIONZ_NAN_WKB, 7, 3),
+        (NESTED_COLLECTION_NAN_WKB, 7, 2),
+        (NESTED_COLLECTIONZ_NAN_WKB, 7, 3),
     ],
 )
-def test_from_wkb_point_empty(wkb, expected_type):
+def test_from_wkb_point_empty(wkb, expected_type, expected_dim):
     geom = pygeos.from_wkb(wkb)
     # POINT (nan nan) transforms to an empty point
-    # Note that the dimensionality (2D/3D) is GEOS-version dependent
     assert pygeos.is_empty(geom)
     assert pygeos.get_type_id(geom) == expected_type
+    # The dimensionality (2D/3D) is only read correctly for GEOS >= 3.9.0
+    if pygeos.geos_version >= (3, 9, 0):
+        assert pygeos.get_coordinate_dimension(geom) == expected_dim
 
 
 def test_to_wkb_point_empty_srid():

src/geos.c

@@ -25,6 +25,12 @@ void destroy_geom_arr(void* context, GEOSGeometry** array, int length) {
   }
 }
 
+/* These functions are used to workaround two Pre-GEOS 3.9.0 issues:
+ * - POINT EMPTY was not handled correctly (we do it ourselves)
+ * - MULTIPOINT (EMPTY) resulted in segfault (we check for it and raise)
+ */
+#if !GEOS_SINCE_3_9_0
+
 /* Returns 1 if a multipoint has an empty point, 0 otherwise, 2 on error.
  */
 char multipoint_has_point_empty(GEOSContextHandle_t ctx, GEOSGeometry* geom) {
@@ -50,10 +56,6 @@ char multipoint_has_point_empty(GEOSContextHandle_t ctx, GEOSGeometry* geom) {
   return 0;
 }
 
-// POINT EMPTY is converted to POINT (nan nan)
-// by GEOS >= 3.10.0. Before that, we do it ourselves here.
-#if !GEOS_SINCE_3_10_0
-
 /* Returns 1 if geometry is an empty point, 0 otherwise, 2 on error.
  */
 char is_point_empty(GEOSContextHandle_t ctx, GEOSGeometry* geom) {
@@ -255,8 +257,6 @@ GEOSGeometry* point_empty_to_nan_all_geoms(GEOSContextHandle_t ctx, GEOSGeometry
   return result;
 }
 
-#endif  // !GEOS_SINCE_3_10_0
-
 /* Checks whether the geometry is a multipoint with an empty point in it
  *
  * According to https://github.com/libgeos/geos/issues/305, this check is not
@@ -289,6 +289,8 @@ char check_to_wkt_compatible(GEOSContextHandle_t ctx, GEOSGeometry* geom) {
   }
 }
 
+#endif  // !GEOS_SINCE_3_9_0
+
 /* GEOSInterpolate_r and GEOSInterpolateNormalized_r segfault on empty
  * geometries and also on collections with the first geometry empty.
  *

src/geos.h

@@ -141,12 +141,12 @@ extern PyObject* geos_exception[1];
 extern void geos_error_handler(const char* message, void* userdata);
 extern void geos_notice_handler(const char* message, void* userdata);
 extern void destroy_geom_arr(void* context, GEOSGeometry** array, int length);
-#if !GEOS_SINCE_3_10_0
+#if !GEOS_SINCE_3_9_0
 extern char has_point_empty(GEOSContextHandle_t ctx, GEOSGeometry* geom);
 extern GEOSGeometry* point_empty_to_nan_all_geoms(GEOSContextHandle_t ctx,
                                                   GEOSGeometry* geom);
-#endif  // !GEOS_SINCE_3_10_0
 extern char check_to_wkt_compatible(GEOSContextHandle_t ctx, GEOSGeometry* geom);
+#endif  // !GEOS_SINCE_3_9_0
 extern char geos_interpolate_checker(GEOSContextHandle_t ctx, GEOSGeometry* geom);
 
 extern int init_geos(PyObject* m);

src/pygeom.c

@@ -72,10 +72,13 @@ static PyObject* GeometryObject_ToWKT(GeometryObject* obj) {
 
   GEOS_INIT;
 
+  // Before GEOS 3.9.0, there was as segfault on e.g. MULTIPOINT (1 1, EMPTY)
+  #if !GEOS_SINCE_3_9_0
   errstate = check_to_wkt_compatible(ctx, obj->ptr);
   if (errstate != PGERR_SUCCESS) {
     goto finish;
   }
+  #endif
 
   GEOSWKTWriter* writer = GEOSWKTWriter_create_r(ctx);
   if (writer == NULL) {
@@ -126,8 +129,8 @@ static PyObject* GeometryObject_ToWKB(GeometryObject* obj) {
 
   GEOS_INIT;
 
-#if !GEOS_SINCE_3_10_0
-  // WKB Does not allow empty points in GEOS < 3.10.
+#if !GEOS_SINCE_3_9_0
+  // WKB Does not allow empty points in GEOS < 3.9.
   // We check for that and patch the POINT EMPTY if necessary
   has_empty = has_point_empty(ctx, obj->ptr);
   if (has_empty == 2) {
@@ -141,7 +144,7 @@ static PyObject* GeometryObject_ToWKB(GeometryObject* obj) {
   }
 #else
   geom = obj->ptr;
-#endif  // !GEOS_SINCE_3_10_0
+#endif  // !GEOS_SINCE_3_9_0
 
   /* Create the WKB writer */
   writer = GEOSWKBWriter_create_r(ctx);