PERF: Use GEOSCoordSeq_copyFromBuffer for creating simple geometries (#…

…436)

CHANGELOG.rst

@@ -22,6 +22,8 @@ Version 0.12 (unreleased)
   to find geometries within a search distance for GEOS >= 3.10 (#425).
 * Added GeoJSON input/output capabilities (``pygeos.from_geojson``,
   ``pygeos.to_geojson``) for GEOS >= 3.10 (#413).
+* Performance improvement in constructing LineStrings or LinearRings from
+  numpy arrays for GEOS >= 3.10 (#436)
 
 **API Changes**
 

pygeos/_geometry.pyx

@@ -12,11 +12,6 @@ import pygeos
 
 from pygeos._geos cimport (
     GEOSContextHandle_t,
-    GEOSCoordSeq_create_r,
-    GEOSCoordSeq_destroy_r,
-    GEOSCoordSeq_setX_r,
-    GEOSCoordSeq_setY_r,
-    GEOSCoordSeq_setZ_r,
     GEOSCoordSequence,
     GEOSGeom_clone_r,
     GEOSGeom_createCollection_r,
@@ -35,6 +30,7 @@ from pygeos._geos cimport (
 )
 from pygeos._pygeos_api cimport (
     import_pygeos_c_api,
+    PyGEOS_CoordSeq_FromBuffer,
     PyGEOS_CreateGeometry,
     PyGEOS_GetGEOSGeometry,
 )
@@ -43,18 +39,6 @@ from pygeos._pygeos_api cimport (
 import_pygeos_c_api()
 
 
-cdef char _set_xyz(GEOSContextHandle_t geos_handle, GEOSCoordSequence *seq, unsigned int coord_idx,
-                   unsigned int dims, double[:, :] coord_view, Py_ssize_t idx):
-    if GEOSCoordSeq_setX_r(geos_handle, seq, coord_idx, coord_view[idx, 0]) == 0:
-        return 0
-    if GEOSCoordSeq_setY_r(geos_handle, seq, coord_idx, coord_view[idx, 1]) == 0:
-        return 0
-    if dims == 3:
-        if GEOSCoordSeq_setZ_r(geos_handle, seq, coord_idx, coord_view[idx, 2]) == 0:
-            return 0
-    return 1
-
-
 def _check_out_array(object out, Py_ssize_t size):
     if out is None:
         return np.empty(shape=(size, ), dtype=object)
@@ -79,12 +63,11 @@ def simple_geometries_1d(object coordinates, object indices, int geometry_type,
     cdef Py_ssize_t geom_idx = 0
     cdef unsigned int geom_size = 0
     cdef unsigned int ring_closure = 0
-    cdef Py_ssize_t coll_geom_idx = 0
     cdef GEOSGeometry *geom = NULL
     cdef GEOSCoordSequence *seq = NULL
 
     # Cast input arrays and define memoryviews for later usage
-    coordinates = np.asarray(coordinates, dtype=np.float64)
+    coordinates = np.asarray(coordinates, dtype=np.float64, order="C")
     if coordinates.ndim != 2:
         raise TypeError("coordinates must be a two-dimensional array.")
 
@@ -149,22 +132,16 @@ def simple_geometries_1d(object coordinates, object indices, int geometry_type,
                     # the error equals PGERR_LINEARRING_NCOORDS (in pygeos/src/geos.h)
                     raise ValueError("A linearring requires at least 4 coordinates.")
 
-            seq = GEOSCoordSeq_create_r(geos_handle, geom_size + ring_closure, dims)
-            for coord_idx in range(geom_size):
-                if _set_xyz(geos_handle, seq, coord_idx, dims, coord_view, idx) == 0:
-                    GEOSCoordSeq_destroy_r(geos_handle, seq)
-                    return  # GEOSException is raised by get_geos_handle
-                idx += 1
+            seq = PyGEOS_CoordSeq_FromBuffer(geos_handle, &coord_view[idx, 0], geom_size, dims, ring_closure)
+            if seq == NULL:
+                return  # GEOSException is raised by get_geos_handle
+            idx += geom_size
 
             if geometry_type == 0:
                 geom = GEOSGeom_createPoint_r(geos_handle, seq)
             elif geometry_type == 1:
                 geom = GEOSGeom_createLineString_r(geos_handle, seq)
             elif geometry_type == 2:
-                if ring_closure == 1:
-                    if _set_xyz(geos_handle, seq, geom_size, dims, coord_view, idx - geom_size) == 0:
-                        GEOSCoordSeq_destroy_r(geos_handle, seq)
-                        return  # GEOSException is raised by get_geos_handle
                 geom = GEOSGeom_createLinearRing_r(geos_handle, seq)
 
             if geom == NULL:

pygeos/_pygeos_api.pxd

@@ -16,7 +16,7 @@ Segfaults will occur if the C API is not imported properly.
 cimport numpy as np
 from cpython.ref cimport PyObject
 
-from pygeos._geos cimport GEOSContextHandle_t, GEOSGeometry
+from pygeos._geos cimport GEOSContextHandle_t, GEOSCoordSequence, GEOSGeometry
 
 
 cdef extern from "c_api.h":
@@ -30,3 +30,6 @@ cdef extern from "c_api.h":
     # they are declared that way in the header file).
     object PyGEOS_CreateGeometry(GEOSGeometry *ptr, GEOSContextHandle_t ctx)
     char PyGEOS_GetGEOSGeometry(PyObject *obj, GEOSGeometry **out) nogil
+    GEOSCoordSequence* PyGEOS_CoordSeq_FromBuffer(GEOSContextHandle_t ctx, const double* buf,
+                                                 unsigned int size, unsigned int dims,
+                                                 char ring_closure) nogil

pygeos/tests/test_creation.py

@@ -163,6 +163,25 @@ def test_linearrings_all_nan():
         pygeos.linearrings(coords)
 
 
+@pytest.mark.parametrize("dim", [2, 3])
+@pytest.mark.parametrize("order", ["C", "F"])
+def test_linearrings_buffer(dim, order):
+    coords1 = np.random.randn(10, 4, dim)
+    coords1 = np.asarray(coords1, order=order)
+    result1 = pygeos.linearrings(coords1)
+
+    # with manual closure -> can directly copy from buffer if C order
+    coords2 = np.hstack((coords1, coords1[:, [0], :]))
+    coords2 = np.asarray(coords2, order=order)
+    result2 = pygeos.linearrings(coords2)
+    assert_geometries_equal(result1, result2)
+
+    # create scalar -> can also directly copy from buffer if F order
+    coords3 = np.asarray(coords2[0], order=order)
+    result3 = pygeos.linearrings(coords3)
+    assert_geometries_equal(result3, result1[0])
+
+
 def test_polygon_from_linearring():
     actual = pygeos.polygons(pygeos.linearrings(box_tpl(0, 0, 1, 1)))
     assert_geometries_equal(

pygeos/tests/test_creation_indices.py

@@ -222,6 +222,22 @@ def test_linearrings_out(indices, expected):
     assert actual is out
 
 
+@pytest.mark.parametrize("dim", [2, 3])
+@pytest.mark.parametrize("order", ["C", "F"])
+def test_linearrings_buffer(dim, order):
+    coords = np.random.randn(10, 4, dim)
+    coords1 = np.asarray(coords.reshape(10 * 4, dim), order=order)
+    indices1 = np.repeat(range(10), 4)
+    result1 = pygeos.linearrings(coords1, indices=indices1)
+
+    # with manual closure -> can directly copy from buffer if C order
+    coords2 = np.hstack((coords, coords[:, [0], :]))
+    coords2 = np.asarray(coords2.reshape(10 * 5, dim), order=order)
+    indices2 = np.repeat(range(10), 5)
+    result2 = pygeos.linearrings(coords2, indices=indices2)
+    assert_geometries_equal(result1, result2)
+
+
 hole_1 = pygeos.linearrings([(0.2, 0.2), (0.2, 0.4), (0.4, 0.4)])
 hole_2 = pygeos.linearrings([(0.6, 0.6), (0.6, 0.8), (0.8, 0.8)])
 poly = pygeos.polygons(linear_ring)

src/c_api.c

@@ -24,3 +24,8 @@ extern PyObject* PyGEOS_CreateGeometry(GEOSGeometry *ptr, GEOSContextHandle_t ct
 extern char PyGEOS_GetGEOSGeometry(PyObject *obj, GEOSGeometry **out) {
     return get_geom((GeometryObject*)obj, out);
 }
+
+extern GEOSCoordSequence* PyGEOS_CoordSeq_FromBuffer(GEOSContextHandle_t ctx, const double* buf,
+                                                    unsigned int size, unsigned int dims, char ring_closure) {
+    return coordseq_from_buffer(ctx, buf, size, dims, ring_closure, dims * 8, 8);
+}

src/c_api.h

@@ -32,8 +32,14 @@
 #define PyGEOS_GetGEOSGeometry_RETURN char
 #define PyGEOS_GetGEOSGeometry_PROTO (PyObject * obj, GEOSGeometry * *out)
 
+/* GEOSCoordSequence* PyGEOS_CoordSeq_FromBuffer(GEOSContextHandle_t ctx, const double* buf,
+                                                unsigned int size, unsigned int dims, char ring_closure)*/
+#define PyGEOS_CoordSeq_FromBuffer_NUM 2
+#define PyGEOS_CoordSeq_FromBuffer_RETURN GEOSCoordSequence*
+#define PyGEOS_CoordSeq_FromBuffer_PROTO (GEOSContextHandle_t ctx, const double* buf, unsigned int size, unsigned int dims, char ring_closure)
+
 /* Total number of C API pointers */
-#define PyGEOS_API_num_pointers 2
+#define PyGEOS_API_num_pointers 3
 
 #ifdef PyGEOS_API_Module
 /* This section is used when compiling pygeos.lib C extension.
@@ -42,6 +48,7 @@
 
 extern PyGEOS_CreateGeometry_RETURN PyGEOS_CreateGeometry PyGEOS_CreateGeometry_PROTO;
 extern PyGEOS_GetGEOSGeometry_RETURN PyGEOS_GetGEOSGeometry PyGEOS_GetGEOSGeometry_PROTO;
+extern PyGEOS_CoordSeq_FromBuffer_RETURN PyGEOS_CoordSeq_FromBuffer PyGEOS_CoordSeq_FromBuffer_PROTO;
 
 #else
 /* This section is used in modules that use the PyGEOS C API
@@ -57,6 +64,9 @@ static void **PyGEOS_API;
 #define PyGEOS_GetGEOSGeometry \
     (*(PyGEOS_GetGEOSGeometry_RETURN(*) PyGEOS_GetGEOSGeometry_PROTO)PyGEOS_API[PyGEOS_GetGEOSGeometry_NUM])
 
+#define PyGEOS_CoordSeq_FromBuffer \
+    (*(PyGEOS_CoordSeq_FromBuffer_RETURN(*) PyGEOS_CoordSeq_FromBuffer_PROTO)PyGEOS_API[PyGEOS_CoordSeq_FromBuffer_NUM])
+
 /* Dynamically load C API from PyCapsule.
  * This MUST be called prior to using C API functions in other modules; otherwise
  * segfaults will occur when the PyGEOS C API functions are called.

src/geos.c

@@ -4,6 +4,7 @@
 #include "geos.h"
 
 #include <Python.h>
+#include <numpy/ndarraytypes.h>
 #include <numpy/npy_math.h>
 #include <structmember.h>
 
@@ -880,3 +881,60 @@ GEOSGeometry* PyGEOSForce2D(GEOSContextHandle_t ctx, GEOSGeometry* geom) {
 GEOSGeometry* PyGEOSForce3D(GEOSContextHandle_t ctx, GEOSGeometry* geom, double z) {
   return force_dims(ctx, geom, 3, z);
 }
+
+GEOSCoordSequence* coordseq_from_buffer(GEOSContextHandle_t ctx, const double* buf,
+                                        unsigned int size, unsigned int dims, char ring_closure,
+                                        npy_intp cs1, npy_intp cs2) {
+  GEOSCoordSequence* coord_seq;
+  char *cp1, *cp2;
+  unsigned int i, j;
+  double first_coord;
+
+#if GEOS_SINCE_3_10_0
+
+  if (!ring_closure) {
+    if ((cs1 == dims * 8) && (cs2 == 8)) {
+      /* C-contiguous memory */
+      int hasZ = dims == 3;
+      coord_seq = GEOSCoordSeq_copyFromBuffer_r(ctx, buf, size, hasZ, 0);
+      return coord_seq;
+    }
+    else if ((cs1 == 8) && (cs2 == size * 8)) {
+      /* F-contiguous memory (note: this for the subset, so we don't necessarily
+      end up here if the full array is F-contiguous) */
+      const double* x = buf;
+      const double* y =  (double*)((char*)buf + cs2);
+      const double* z = (dims == 3) ? (double*)((char*)buf + 2 * cs2) : NULL;
+      coord_seq = GEOSCoordSeq_copyFromArrays_r(ctx, x, y, z, NULL, size);
+      return coord_seq;
+    }
+  }
+
+#endif
+
+  coord_seq = GEOSCoordSeq_create_r(ctx, size + ring_closure, dims);
+  if (coord_seq == NULL) {
+    return NULL;
+  }
+  cp1 = (char*)buf;
+  for (i = 0; i < size; i++, cp1 += cs1) {
+    cp2 = cp1;
+    for (j = 0; j < dims; j++, cp2 += cs2) {
+      if (!GEOSCoordSeq_setOrdinate_r(ctx, coord_seq, i, j, *(double*)cp2)) {
+        GEOSCoordSeq_destroy_r(ctx, coord_seq);
+        return NULL;
+      }
+    }
+  }
+  /* add the closing coordinate if necessary */
+  if (ring_closure) {
+    for (j = 0; j < dims; j++){
+      first_coord = *(double*)((char*)buf + j * cs2);
+      if (!GEOSCoordSeq_setOrdinate_r(ctx, coord_seq, size, j, first_coord)) {
+        GEOSCoordSeq_destroy_r(ctx, coord_seq);
+        return NULL;
+      }
+    }
+  }
+  return coord_seq;
+}

src/geos.h

@@ -2,6 +2,7 @@
 #define _GEOS_H
 
 #include <Python.h>
+#include <numpy/ndarraytypes.h>
 
 /* To avoid accidental use of non reentrant GEOS API. */
 #ifndef GEOS_USE_ONLY_R_API
@@ -176,4 +177,8 @@ GEOSGeometry* create_point(GEOSContextHandle_t ctx, double x, double y);
 GEOSGeometry* PyGEOSForce2D(GEOSContextHandle_t ctx, GEOSGeometry* geom);
 GEOSGeometry* PyGEOSForce3D(GEOSContextHandle_t ctx, GEOSGeometry* geom, double z);
 
+GEOSCoordSequence* coordseq_from_buffer(GEOSContextHandle_t ctx, const double* buf,
+                                        unsigned int size, unsigned int dims, char ring_closure,
+                                        npy_intp cs1, npy_intp cs2);
+
 #endif  // _GEOS_H

src/lib.c

@@ -89,6 +89,7 @@ PyMODINIT_FUNC PyInit_lib(void) {
   /* Initialize the C API pointer array */
   PyGEOS_API[PyGEOS_CreateGeometry_NUM] = (void*)PyGEOS_CreateGeometry;
   PyGEOS_API[PyGEOS_GetGEOSGeometry_NUM] = (void*)PyGEOS_GetGEOSGeometry;
+  PyGEOS_API[PyGEOS_CoordSeq_FromBuffer_NUM] = (void*)PyGEOS_CoordSeq_FromBuffer;
 
   /* Create a Capsule containing the API pointer array's address */
   c_api_object = PyCapsule_New((void*)PyGEOS_API, "pygeos.lib._C_API", NULL);

src/ufuncs.c

@@ -2198,22 +2198,12 @@ static void linestrings_func(char** args, npy_intp* dimensions, npy_intp* steps,
   GEOS_INIT_THREADS;
 
   DOUBLE_COREDIM_LOOP_OUTER {
-    coord_seq = GEOSCoordSeq_create_r(ctx, n_c1, n_c2);
+    coord_seq = coordseq_from_buffer(ctx, (double*)ip1, n_c1, n_c2, 0, cs1, cs2);
     if (coord_seq == NULL) {
       errstate = PGERR_GEOS_EXCEPTION;
       destroy_geom_arr(ctx, geom_arr, i - 1);
       goto finish;
     }
-    DOUBLE_COREDIM_LOOP_INNER_1 {
-      DOUBLE_COREDIM_LOOP_INNER_2 {
-        if (!GEOSCoordSeq_setOrdinate_r(ctx, coord_seq, i_c1, i_c2, *(double*)cp2)) {
-          errstate = PGERR_GEOS_EXCEPTION;
-          GEOSCoordSeq_destroy_r(ctx, coord_seq);
-          destroy_geom_arr(ctx, geom_arr, i - 1);
-          goto finish;
-        }
-      }
-    }
     geom_arr[i] = GEOSGeom_createLineString_r(ctx, coord_seq);
     // Note: coordinate sequence is owned by linestring; if linestring fails to construct,
     // it will automatically clean up the coordinate sequence
@@ -2271,34 +2261,12 @@ static void linearrings_func(char** args, npy_intp* dimensions, npy_intp* steps,
       goto finish;
     }
     /* fill the coordinate sequence */
-    coord_seq = GEOSCoordSeq_create_r(ctx, n_c1 + ring_closure, n_c2);
+    coord_seq = coordseq_from_buffer(ctx, (double*)ip1, n_c1, n_c2, ring_closure, cs1, cs2);
     if (coord_seq == NULL) {
       errstate = PGERR_GEOS_EXCEPTION;
       destroy_geom_arr(ctx, geom_arr, i - 1);
       goto finish;
     }
-    DOUBLE_COREDIM_LOOP_INNER_1 {
-      DOUBLE_COREDIM_LOOP_INNER_2 {
-        if (!GEOSCoordSeq_setOrdinate_r(ctx, coord_seq, i_c1, i_c2, *(double*)cp2)) {
-          errstate = PGERR_GEOS_EXCEPTION;
-          GEOSCoordSeq_destroy_r(ctx, coord_seq);
-          destroy_geom_arr(ctx, geom_arr, i - 1);
-          goto finish;
-        }
-      }
-    }
-    /* add the closing coordinate if necessary */
-    if (ring_closure) {
-      DOUBLE_COREDIM_LOOP_INNER_2 {
-        first_coord = *(double*)(ip1 + i_c2 * cs2);
-        if (!GEOSCoordSeq_setOrdinate_r(ctx, coord_seq, n_c1, i_c2, first_coord)) {
-          errstate = PGERR_GEOS_EXCEPTION;
-          GEOSCoordSeq_destroy_r(ctx, coord_seq);
-          destroy_geom_arr(ctx, geom_arr, i - 1);
-          goto finish;
-        }
-      }
-    }
     geom_arr[i] = GEOSGeom_createLinearRing_r(ctx, coord_seq);
     // Note: coordinate sequence is owned by linearring; if linearring fails to construct,
     // it will automatically clean up the coordinate sequence