[Done] Fix segfault when getting coordinates from empty points (#415)

.github/workflows/test-pip.yml

@@ -14,7 +14,7 @@ jobs:
       matrix:
         os: [ubuntu-latest, macos-latest, windows-2019]
         architecture: [x64]
-        geos: [3.6.4, 3.7.3, 3.8.1, 3.9.1, 3.10.0, main]
+        geos: [3.6.4, 3.7.3, 3.8.0, 3.9.1, 3.10.0, main]
         include:
           # 2017
           - python: 3.6
@@ -26,7 +26,7 @@ jobs:
             numpy: 1.15.4
           # 2019
           - python: 3.8
-            geos: 3.8.1
+            geos: 3.8.0
             numpy: 1.17.5
           # 2020
           - python: 3.9

CHANGELOG.rst

@@ -39,6 +39,7 @@ Version 0.11 (unreleased)
 * 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).
+* Fixed a segfault when getting coordinates from empty points in GEOS 3.8.0.
 
 **Acknowledgments**
 

pygeos/tests/test_constructive.py

@@ -276,7 +276,7 @@ def test_offset_curve_join_style_invalid():
             pygeos.Geometry("MULTILINESTRING ((0 0, 1 2), (3 3, 4 4))"),
             pygeos.Geometry("MULTILINESTRING ((1 2, 0 0), (4 4, 3 3))"),
             marks=pytest.mark.skipif(
-                pygeos.geos_version < (3, 8, 0), reason="GEOS < 3.8"
+                pygeos.geos_version < (3, 8, 1), reason="GEOS < 3.8.1"
             ),
         ),
         (

pygeos/tests/test_coordinates.py

@@ -125,6 +125,7 @@ def test_get_coords_index_multidim(order):
         ([line_string_z], [0, 1, 1], [0, 0, 1], [4, 4, 4]),
         ([polygon_z], [0, 2, 2, 0, 0], [0, 0, 2, 2, 0], [4, 4, 4, 4, 4]),
         ([geometry_collection_z], [2, 0, 1, 1], [3, 0, 0, 1], [4, 4, 4, 4]),
+        ([point, empty_point], [2], [3], [np.nan]),
     ],
 )  # fmt: on
 def test_get_coords_3d(geoms, x, y, z, include_z):

src/coords.c

@@ -23,11 +23,26 @@ static void* set_coordinates(GEOSContextHandle_t, GEOSGeometry*, PyArrayObject*,
 /* Get coordinates from a point, linestring or linearring and puts them at
 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,
+static char get_coordinates_simple(GEOSContextHandle_t ctx, GEOSGeometry* geom, int type,
                                    PyArrayObject* out, npy_intp* cursor, int include_z) {
   unsigned int n, i;
   double *x, *y, *z;
-  const GEOSCoordSequence* seq = GEOSGeom_getCoordSeq_r(ctx, geom);
+  const GEOSCoordSequence* seq;
+  char is_empty;
+
+  /* For points, directly check if they are empty. This is because empty points
+   * internally have a coordinate sequence of length 1, but we did not count them in
+   * the allocation of the coordinate array */
+  if (type == GEOS_POINT) {
+    is_empty = GEOSisEmpty_r(ctx, geom);
+    if (is_empty == 2) {
+      return 0;
+    } else if (is_empty == 1) {
+      return 1;
+    }
+  }
+
+  seq = GEOSGeom_getCoordSeq_r(ctx, geom);
   if (seq == NULL) {
     return 0;
   }
@@ -67,7 +82,7 @@ static char get_coordinates_polygon(GEOSContextHandle_t ctx, GEOSGeometry* geom,
   if (ring == NULL) {
     return 0;
   }
-  if (!get_coordinates_simple(ctx, ring, out, cursor, include_z)) {
+  if (!get_coordinates_simple(ctx, ring, GEOS_LINEARRING, out, cursor, include_z)) {
     return 0;
   }
 
@@ -80,7 +95,7 @@ static char get_coordinates_polygon(GEOSContextHandle_t ctx, GEOSGeometry* geom,
     if (ring == NULL) {
       return 0;
     }
-    if (!get_coordinates_simple(ctx, ring, out, cursor, include_z)) {
+    if (!get_coordinates_simple(ctx, ring, GEOS_LINEARRING, out, cursor, include_z)) {
       return 0;
     }
   }
@@ -119,7 +134,7 @@ static char get_coordinates(GEOSContextHandle_t ctx, GEOSGeometry* geom,
                             PyArrayObject* out, npy_intp* cursor, int include_z) {
   int type = GEOSGeomTypeId_r(ctx, geom);
   if ((type == 0) || (type == 1) || (type == 2)) {
-    return get_coordinates_simple(ctx, geom, out, cursor, include_z);
+    return get_coordinates_simple(ctx, geom, type, out, cursor, include_z);
   } else if (type == 3) {
     return get_coordinates_polygon(ctx, geom, out, cursor, include_z);
   } else if ((type >= 4) && (type <= 7)) {
@@ -421,7 +436,7 @@ PyObject* GetCoords(PyArrayObject* arr, int include_z, int return_index) {
   }
 
   /* Handle zero-sized arrays specially */
-  if (PyArray_SIZE(arr) == 0) {
+  if (size == 0) {
     if (return_index) {
       PyObject* result_tpl = PyTuple_New(2);
       PyTuple_SET_ITEM(result_tpl, 0, (PyObject*)result);