Create simple geometries from coordinates and indices (#322)

CHANGELOG.rst

@@ -8,8 +8,11 @@ Version 0.10 (unreleased)
 
 * Addition of ``nearest`` and ``nearest_all`` functions to ``STRtree`` for
   GEOS >= 3.6 to find the nearest neighbors (#272).
+* Enable bulk construction of geometries with different number of coordinates
+  by optionally taking index arrays in the constructors ``multipoints``,
+  ``points``, ``linestrings``, and ``linearrings`` (#322).
 * Enable bulk construction of collections with different number of geometries
-  by optionally taking an index arrays in the constructors ``multipoints``,
+  by optionally taking index arrays in the constructors ``multipoints``,
   ``multilinestrings``, ``multipolygons``, and ``geometrycollections`` (#290).
 * Released GIL for ``points``, ``linestrings``, ``linearrings``, and 
   ``polygons`` (without holes) (#310).

pygeos/_geometry.pyx

@@ -10,12 +10,22 @@ import pygeos
 
 from pygeos._geos cimport (
     GEOSGeometry,
+    GEOSCoordSequence,
+    GEOSContextHandle_t,
     GEOSGeom_clone_r,
     GEOSGetGeometryN_r,
     get_geos_handle,
     GEOSGeom_destroy_r,
     GEOSGeom_createCollection_r,
     GEOSGeomTypeId_r,
+    GEOSCoordSeq_create_r,
+    GEOSCoordSeq_destroy_r,
+    GEOSCoordSeq_setX_r,
+    GEOSCoordSeq_setY_r,
+    GEOSCoordSeq_setZ_r,
+    GEOSGeom_createPoint_r,
+    GEOSGeom_createLineString_r,
+    GEOSGeom_createLinearRing_r
 )
 from pygeos._pygeos_api cimport (
     import_pygeos_c_api,
@@ -27,6 +37,110 @@ 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
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def simple_geometries_1d(object coordinates, object indices, int geometry_type):
+    cdef Py_ssize_t idx = 0
+    cdef unsigned int coord_idx = 0
+    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)
+    if coordinates.ndim != 2:
+        raise TypeError("coordinates is not a two-dimensional array.")
+
+    indices = np.asarray(indices, dtype=np.intp)
+    if indices.ndim != 1:
+        raise TypeError("indices is not a one-dimensional array.")
+
+    if coordinates.shape[0] != indices.shape[0]:
+        raise ValueError("geometries and indices do not have equal size.")
+
+    cdef unsigned int dims = coordinates.shape[1]
+    if dims not in {2, 3}:
+        raise ValueError("coordinates should N by 2 or N by 3.")
+
+    if geometry_type not in {0, 1, 2}:
+        raise ValueError(f"Invalid geometry_type: {geometry_type}.")
+
+    if coordinates.shape[0] == 0:
+        # return immediately if there are no geometries to return
+        return np.empty(shape=(0, ), dtype=np.object_)
+
+    if np.any(indices[1:] < indices[:indices.shape[0] - 1]):
+        raise ValueError("The indices must be sorted.")  
+
+    cdef double[:, :] coord_view = coordinates
+    cdef np.intp_t[:] index_view = indices
+
+    # get the geometry count per collection
+    cdef unsigned int[:] coord_counts = np.bincount(indices).astype(np.uint32)
+
+    # The final target array
+    cdef Py_ssize_t n_geoms = coord_counts.shape[0]
+    result = np.empty(shape=(n_geoms, ), dtype=object)
+    cdef object[:] result_view = result
+
+    with get_geos_handle() as geos_handle:
+        for geom_idx in range(n_geoms):
+            geom_size = coord_counts[geom_idx]
+
+            # insert None if there are no coordinates
+            if geom_size == 0:
+                result_view[geom_idx] = PyGEOS_CreateGeometry(NULL, geos_handle)
+                continue
+
+            # check if we need to close a linearring
+            if geometry_type == 2:
+                ring_closure = 0
+                for coord_idx in range(dims):
+                    if coord_view[idx, coord_idx] != coord_view[idx + geom_size - 1, coord_idx]:
+                        ring_closure = 1
+                        break
+
+            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
+
+            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:
+                return  # GEOSException is raised by get_geos_handle
+
+            result_view[geom_idx] = PyGEOS_CreateGeometry(geom, geos_handle)
+
+    return result
+
+
 @cython.boundscheck(False)
 @cython.wraparound(False)
 def get_parts(object[:] array):
@@ -42,11 +156,11 @@ def get_parts(object[:] array):
     if count == 0:
         # return immediately if there are no geometries to return
         return (
-            np.empty(shape=(0, ), dtype=np.object_),
+            np.empty(shape=(0, ), dtype=object),
             np.empty(shape=(0, ), dtype=np.intp)
         )
 
-    parts = np.empty(shape=(count, ), dtype=np.object_)
+    parts = np.empty(shape=(count, ), dtype=object)
     index = np.empty(shape=(count, ), dtype=np.intp)
 
     cdef int[:] counts_view = counts
@@ -69,9 +183,13 @@ def get_parts(object[:] array):
             for part_idx in range(counts_view[geom_idx]):
                 index_view[idx] = geom_idx
                 part = GEOSGetGeometryN_r(geos_handle, geom, part_idx)
+                if part == NULL:
+                    return  # GEOSException is raised by get_geos_handle
 
                 # clone the geometry to keep it separate from the inputs
                 part = GEOSGeom_clone_r(geos_handle, part)
+                if part == NULL:
+                    return  # GEOSException is raised by get_geos_handle
                 # cast part back to <GEOSGeometry> to discard const qualifier
                 # pending issue #227
                 parts_view[idx] = PyGEOS_CreateGeometry(<GEOSGeometry *>part, geos_handle)
@@ -81,6 +199,17 @@ def get_parts(object[:] array):
     return parts, index
 
 
+cdef _deallocate_arr(void* handle, np.intp_t[:] arr, Py_ssize_t last_geom_i):
+    """Deallocate a temporary geometry array to prevent memory leaks"""
+    cdef Py_ssize_t i = 0
+    cdef GEOSGeometry *g
+
+    for i in range(last_geom_i):
+        g = <GEOSGeometry *>arr[i]
+        if g != NULL:
+            GEOSGeom_destroy_r(handle, <GEOSGeometry *>arr[i])
+
+
 @cython.boundscheck(False)
 @cython.wraparound(False)
 def collections_1d(object geometries, object indices, int geometry_type = 7):
@@ -107,7 +236,7 @@ def collections_1d(object geometries, object indices, int geometry_type = 7):
         raise ValueError(f"Invalid geometry_type: {geometry_type}.")
 
     # Cast input arrays and define memoryviews for later usage
-    geometries = np.asarray(geometries, dtype=np.object)
+    geometries = np.asarray(geometries, dtype=object)
     if geometries.ndim != 1:
         raise TypeError("geometries is not a one-dimensional array.")
 
@@ -120,7 +249,7 @@ def collections_1d(object geometries, object indices, int geometry_type = 7):
 
     if geometries.shape[0] == 0:
         # return immediately if there are no geometries to return
-        return np.empty(shape=(0, ), dtype=np.object_)
+        return np.empty(shape=(0, ), dtype=object)
 
     if np.any(indices[1:] < indices[:indices.shape[0] - 1]):
         raise ValueError("The indices should be sorted.")  
@@ -138,7 +267,7 @@ def collections_1d(object geometries, object indices, int geometry_type = 7):
 
     # The final target array
     cdef Py_ssize_t n_colls = collection_size.shape[0]
-    result = np.empty(shape=(n_colls, ), dtype=np.object_)
+    result = np.empty(shape=(n_colls, ), dtype=object)
     cdef object[:] result_view = result
 
     with get_geos_handle() as geos_handle:
@@ -148,9 +277,7 @@ def collections_1d(object geometries, object indices, int geometry_type = 7):
             # fill the temporary array with geometries belonging to this collection
             for coll_geom_idx in range(collection_size[coll_idx]):
                 if PyGEOS_GetGEOSGeometry(<PyObject *>geometries_view[geom_idx_1 + coll_geom_idx], &geom) == 0:
-                    # deallocate previous temp geometries (preventing memory leaks)
-                    for coll_geom_idx in range(coll_size):
-                        GEOSGeom_destroy_r(geos_handle, <GEOSGeometry *>temp_geoms_view[coll_geom_idx])
+                    _deallocate_arr(geos_handle, temp_geoms_view, coll_size)
                     raise TypeError(
                         "One of the arguments is of incorrect type. Please provide only Geometry objects."
                     )
@@ -162,16 +289,21 @@ def collections_1d(object geometries, object indices, int geometry_type = 7):
                 # Check geometry subtype for non-geometrycollections
                 if geometry_type != 7:
                     curr_type = GEOSGeomTypeId_r(geos_handle, geom)
+                    if curr_type == -1:
+                        _deallocate_arr(geos_handle, temp_geoms_view, coll_size)
+                        return  # GEOSException is raised by get_geos_handle
                     if curr_type != expected_type and curr_type != expected_type_alt:
-                        # deallocate previous temp geometries (preventing memory leaks)
-                        for coll_geom_idx in range(coll_size):
-                            GEOSGeom_destroy_r(geos_handle, <GEOSGeometry *>temp_geoms_view[coll_geom_idx])
+                        _deallocate_arr(geos_handle, temp_geoms_view, coll_size)
                         raise TypeError(
                             f"One of the arguments has unexpected geometry type {curr_type}."
                         )
 
-                # assign to the temporary geometry array                   
-                temp_geoms_view[coll_size] = <np.intp_t>GEOSGeom_clone_r(geos_handle, geom)
+                # assign to the temporary geometry array  
+                geom = GEOSGeom_clone_r(geos_handle, geom)
+                if geom == NULL:
+                    _deallocate_arr(geos_handle, temp_geoms_view, coll_size)
+                    return  # GEOSException is raised by get_geos_handle           
+                temp_geoms_view[coll_size] = <np.intp_t>geom
                 coll_size += 1
 
             # create the collection
@@ -181,6 +313,8 @@ def collections_1d(object geometries, object indices, int geometry_type = 7):
                 <GEOSGeometry**> &temp_geoms_view[0],
                 coll_size
             )
+            if coll == NULL:
+                return  # GEOSException is raised by get_geos_handle
 
             result_view[coll_idx] = PyGEOS_CreateGeometry(coll, geos_handle)
 

pygeos/_geos.pxd

@@ -8,30 +8,40 @@ Example:
 """
 
 cdef extern from "geos_c.h":
+    # Types
     ctypedef void *GEOSContextHandle_t
+    ctypedef struct GEOSGeometry
+    ctypedef struct GEOSCoordSequence
+    ctypedef void (*GEOSMessageHandler_r)(const char *message, void *userdata)
 
+    # GEOS Context & Messaging
     GEOSContextHandle_t GEOS_init_r() nogil
     void GEOS_finish_r(GEOSContextHandle_t handle) nogil
+    void GEOSContext_setErrorMessageHandler_r(GEOSContextHandle_t handle, GEOSMessageHandler_r ef, void* userData) nogil
+    void GEOSContext_setNoticeMessageHandler_r(GEOSContextHandle_t handle, GEOSMessageHandler_r nf, void* userData) nogil
 
-    ctypedef struct GEOSGeometry
+    # Geometry functions
+    const GEOSGeometry* GEOSGetGeometryN_r(GEOSContextHandle_t handle, const GEOSGeometry* g, int n) nogil
+    int GEOSGeomTypeId_r(GEOSContextHandle_t handle, GEOSGeometry* g) nogil
 
-    const GEOSGeometry* GEOSGetGeometryN_r(GEOSContextHandle_t handle,
-                                           const GEOSGeometry* g,
-                                           int n) nogil except NULL
-    int GEOSGeomTypeId_r(GEOSContextHandle_t handle, GEOSGeometry* g) nogil except -1
-    void GEOSGeom_destroy_r(GEOSContextHandle_t handle,
-                                    GEOSGeometry* g) nogil
-    GEOSGeometry* GEOSGeom_clone_r(GEOSContextHandle_t handle,
-                                   const GEOSGeometry* g) nogil except NULL
+    # Geometry creation / destruction
+    GEOSGeometry* GEOSGeom_clone_r(GEOSContextHandle_t handle, const GEOSGeometry* g) nogil
+    GEOSGeometry* GEOSGeom_createPoint_r(GEOSContextHandle_t handle, GEOSCoordSequence* s) nogil
+    GEOSGeometry* GEOSGeom_createLineString_r(GEOSContextHandle_t handle, GEOSCoordSequence* s) nogil
+    GEOSGeometry* GEOSGeom_createLinearRing_r(GEOSContextHandle_t handle, GEOSCoordSequence* s) nogil
+    GEOSGeometry* GEOSGeom_createCollection_r(GEOSContextHandle_t handle, int type, GEOSGeometry** geoms, unsigned int ngeoms) nogil
+    void GEOSGeom_destroy_r(GEOSContextHandle_t handle, GEOSGeometry* g) nogil
 
-    GEOSGeometry* GEOSGeom_createCollection_r(
-        GEOSContextHandle_t handle,
-        int type,
-        GEOSGeometry** geoms,
-        unsigned int ngeoms
-    ) nogil except NULL
+    # Coordinate sequences
+    GEOSCoordSequence* GEOSCoordSeq_create_r(GEOSContextHandle_t handle, unsigned int size, unsigned int dims) nogil
+    void GEOSCoordSeq_destroy_r(GEOSContextHandle_t handle, GEOSCoordSequence* s) nogil
+    int GEOSCoordSeq_setX_r(GEOSContextHandle_t handle, GEOSCoordSequence* s, unsigned int idx, double val) nogil
+    int GEOSCoordSeq_setY_r(GEOSContextHandle_t handle, GEOSCoordSequence* s, unsigned int idx, double val) nogil
+    int GEOSCoordSeq_setZ_r(GEOSContextHandle_t handle, GEOSCoordSequence* s, unsigned int idx, double val) nogil
 
 
 cdef class get_geos_handle:
     cdef GEOSContextHandle_t handle
+    cdef char* last_error
+    cdef char* last_warning
     cdef GEOSContextHandle_t __enter__(self)

pygeos/_geos.pyx

@@ -1,4 +1,14 @@
 # distutils: define_macros=GEOS_USE_ONLY_R_API
+#from pygeos import GEOSException
+from libc.stdio cimport snprintf
+from libc.stdlib cimport malloc, free
+import warnings
+from pygeos import GEOSException
+
+
+cdef void geos_message_handler(const char* message, void* userdata):
+    snprintf(<char *>userdata, 1024, "%s", message)
+
 
 cdef class get_geos_handle:
     '''This class provides a context manager that wraps the GEOS context handle.
@@ -10,7 +20,25 @@ cdef class get_geos_handle:
     '''
     cdef GEOSContextHandle_t __enter__(self):
         self.handle = GEOS_init_r()
+        self.last_error = <char *> malloc((1025) * sizeof(char))
+        self.last_error[0] = 0
+        self.last_warning = <char *> malloc((1025) * sizeof(char))
+        self.last_warning[0] = 0
+        GEOSContext_setErrorMessageHandler_r(
+            self.handle, &geos_message_handler, self.last_error
+        )
+        GEOSContext_setNoticeMessageHandler_r(
+            self.handle, &geos_message_handler, self.last_warning
+        )
         return self.handle
 
     def __exit__(self, type, value, traceback):
-        GEOS_finish_r(self.handle)
+        try:
+            if self.last_error[0] != 0:
+                raise GEOSException(self.last_error)
+            if self.last_warning[0] != 0:
+                warnings.warn(self.last_warning)
+        finally:
+            GEOS_finish_r(self.handle)
+            free(self.last_error)
+            free(self.last_warning)