Fix 3D empty WKT serialization (#403)

CHANGELOG.rst

@@ -28,6 +28,7 @@ Version 0.11 (unreleased)
 
 * 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 single part 3D empty geometries for GEOS >= 3.9.0 (#402).
 * Fixed the WKT serialization of multipoints with empty points for GEOS >= 3.9.0 (#392).
 
 **Acknowledgments**

pygeos/geometry.py

@@ -776,6 +776,9 @@ def force_3d(geometry, z=0.0, **kwargs):
     2D geometries will get the provided Z coordinate; Z coordinates of 3D geometries
     are unchanged (unless they are nan).
 
+    Note that for empty geometries, 3D is only supported since GEOS 3.9 and then
+    still only for simple geometries (non-collections).
+
     Parameters
     ----------
     geometry : Geometry or array_like

pygeos/io.py

@@ -99,6 +99,13 @@ def to_wkt(
     The Well-known Text format is defined in the `OGC Simple Features
     Specification for SQL <https://www.opengeospatial.org/standards/sfs>`__.
 
+    The following limitations apply to WKT serialization:
+
+    - for GEOS <= 3.8 a multipoint with an empty sub-geometry will raise an exception
+    - for GEOS <= 3.8 empty geometries are always serialized to 2D
+    - for GEOS >= 3.9 only simple empty geometries can be 3D, collections are still
+      always 2D
+
     Parameters
     ----------
     geometry : Geometry or array_like

pygeos/tests/test_io.py

@@ -276,6 +276,23 @@ def test_to_wkt_point_empty():
     assert pygeos.to_wkt(empty_point) == "POINT EMPTY"
 
 
+@pytest.mark.skipif(
+    pygeos.geos_version < (3, 9, 0),
+    reason="Empty geometries have no dimensionality on GEOS < 3.9",
+)
+@pytest.mark.parametrize(
+    "wkt",
+    [
+        "POINT Z EMPTY",
+        "LINESTRING Z EMPTY",
+        "LINEARRING Z EMPTY",
+        "POLYGON Z EMPTY",
+    ],
+)
+def test_to_wkt_empty_z(wkt):
+    assert pygeos.to_wkt(pygeos.Geometry(wkt)) == wkt
+
+
 def test_to_wkt_geometrycollection_with_point_empty():
     collection = pygeos.geometrycollections([empty_point, point])
     # do not check the full value as some GEOS versions give
@@ -325,6 +342,14 @@ def test_repr_multipoint_with_point_empty():
     assert repr(geom) == "<pygeos.Geometry Exception in WKT writer>"
 
 
+@pytest.mark.skipif(
+    pygeos.geos_version < (3, 9, 0),
+    reason="Empty geometries have no dimensionality on GEOS < 3.9",
+)
+def test_repr_point_z_empty():
+    assert repr(empty_point_z) == "<pygeos.Geometry POINT Z EMPTY>"
+
+
 def test_to_wkb():
     point = pygeos.points(1, 1)
     actual = pygeos.to_wkb(point, byte_order=1)

src/geos.c

@@ -291,6 +291,74 @@ char check_to_wkt_compatible(GEOSContextHandle_t ctx, GEOSGeometry* geom) {
 
 #endif  // !GEOS_SINCE_3_9_0
 
+#if GEOS_SINCE_3_9_0
+
+/* Checks whether the geometry is a 3D empty geometry and, if so, create the WKT string
+ *
+ * GEOS 3.9.* is able to distiguish 2D and 3D simple geometries (non-collections). But the
+ * but the WKT serialization never writes a 3D empty geometry. This function fixes that.
+ * It only makes sense to use this for GEOS versions >= 3.9.
+ *
+ * Pending GEOS ticket: https://trac.osgeo.org/geos/ticket/1129
+ *
+ * If the geometry is a 3D empty, then the (char**) wkt argument is filled with the
+ * correct WKT string. Else, wkt becomes NULL and the GEOS WKT writer should be used.
+ *
+ * The return value is one of:
+ * - PGERR_SUCCESS
+ * - PGERR_GEOS_EXCEPTION
+ */
+char wkt_empty_3d_geometry(GEOSContextHandle_t ctx, GEOSGeometry* geom, char** wkt) {
+  char is_empty;
+  int geom_type;
+
+  is_empty = GEOSisEmpty_r(ctx, geom);
+  if (is_empty == 2) {
+    return PGERR_GEOS_EXCEPTION;
+  } else if (is_empty == 0) {
+    *wkt = NULL;
+    return PGERR_SUCCESS;
+  }
+  if (GEOSGeom_getCoordinateDimension_r(ctx, geom) == 2) {
+    *wkt = NULL;
+    return PGERR_SUCCESS;
+  }
+  geom_type = GEOSGeomTypeId_r(ctx, geom);
+  switch (geom_type) {
+    case GEOS_POINT:
+      *wkt = "POINT Z EMPTY";
+      return PGERR_SUCCESS;
+    case GEOS_LINESTRING:
+      *wkt = "LINESTRING Z EMPTY";
+      break;
+    case GEOS_LINEARRING:
+      *wkt = "LINEARRING Z EMPTY";
+      break;
+    case GEOS_POLYGON:
+      *wkt = "POLYGON Z EMPTY";
+      break;
+    // Note: Empty collections cannot be 3D in GEOS.
+    // We do include the options in case of future support.
+    case GEOS_MULTIPOINT:
+      *wkt = "MULTIPOINT Z EMPTY";
+      break;
+    case GEOS_MULTILINESTRING:
+      *wkt = "MULTILINESTRING Z EMPTY";
+      break;
+    case GEOS_MULTIPOLYGON:
+      *wkt = "MULTIPOLYGON Z EMPTY";
+      break;
+    case GEOS_GEOMETRYCOLLECTION:
+      *wkt = "GEOMETRYCOLLECTION Z EMPTY";
+      break;
+    default:
+      return PGERR_GEOS_EXCEPTION;
+  }
+  return PGERR_SUCCESS;
+}
+
+#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

@@ -152,6 +152,10 @@ extern GEOSGeometry* point_empty_to_nan_all_geoms(GEOSContextHandle_t ctx,
                                                   GEOSGeometry* geom);
 extern char check_to_wkt_compatible(GEOSContextHandle_t ctx, GEOSGeometry* geom);
 #endif  // !GEOS_SINCE_3_9_0
+#if GEOS_SINCE_3_9_0
+extern char wkt_empty_3d_geometry(GEOSContextHandle_t ctx, GEOSGeometry* geom,
+                                  char** wkt);
+#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

@@ -65,20 +65,31 @@ static PyMemberDef GeometryObject_members[] = {
 static PyObject* GeometryObject_ToWKT(GeometryObject* obj) {
   char* wkt;
   PyObject* result;
-  if (obj->ptr == NULL) {
+  GEOSGeometry* geom = obj->ptr;
+
+  if (geom == NULL) {
     Py_INCREF(Py_None);
     return Py_None;
   }
 
   GEOS_INIT;
 
+#if GEOS_SINCE_3_9_0
+  errstate = wkt_empty_3d_geometry(ctx, geom, &wkt);
+  if (errstate != PGERR_SUCCESS) {
+    goto finish;
+  }
+  if (wkt != NULL) {
+    result = PyUnicode_FromString(wkt);
+    goto finish;
+  }
+#else
   // 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);
+  errstate = check_to_wkt_compatible(ctx, geom);
   if (errstate != PGERR_SUCCESS) {
     goto finish;
   }
-  #endif
+#endif
 
   GEOSWKTWriter* writer = GEOSWKTWriter_create_r(ctx);
   if (writer == NULL) {
@@ -101,7 +112,7 @@ static PyObject* GeometryObject_ToWKT(GeometryObject* obj) {
     goto finish;
   }
 
-  wkt = GEOSWKTWriter_write_r(ctx, writer, obj->ptr);
+  wkt = GEOSWKTWriter_write_r(ctx, writer, geom);
   result = PyUnicode_FromString(wkt);
   GEOSFree_r(ctx, wkt);
   GEOSWKTWriter_destroy_r(ctx, writer);

src/ufuncs.c

@@ -2940,8 +2940,18 @@ static void to_wkt_func(char** args, npy_intp* dimensions, npy_intp* steps, void
       Py_INCREF(Py_None);
       *out = Py_None;
     } else {
-// Before GEOS 3.9.0, there was as segfault on e.g. MULTIPOINT (1 1, EMPTY)
-#if !GEOS_SINCE_3_9_0
+#if GEOS_SINCE_3_9_0
+      errstate = wkt_empty_3d_geometry(ctx, in1, &wkt);
+      if (errstate != PGERR_SUCCESS) {
+        goto finish;
+      }
+      if (wkt != NULL) {
+        *out = PyUnicode_FromString(wkt);
+        goto finish;
+      }
+
+#else
+      // Before GEOS 3.9.0, there was as segfault on e.g. MULTIPOINT (1 1, EMPTY)
       errstate = check_to_wkt_compatible(ctx, in1);
       if (errstate != PGERR_SUCCESS) {
         goto finish;