[Done] Create polygons from indices (#326)

CHANGELOG.rst

@@ -9,13 +9,8 @@ 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 index arrays in the constructors ``multipoints``,
-  ``multilinestrings``, ``multipolygons``, and ``geometrycollections`` (#290).
-* Released GIL for ``points``, ``linestrings``, ``linearrings``, and
-  ``polygons`` (without holes) (#310).
+  by optionally taking index arrays in all creation functions (#230, #322, #326).
+* Released the GIL in all geometry creation functions (#310, #326).
 * Added the option to return the geometry index in ``get_coordinates`` (#318).
 * Updated ``box`` ufunc to use internal C function for creating polygon
   (about 2x faster) and added ``ccw`` parameter to create polygon in
@@ -31,9 +26,11 @@ Version 0.10 (unreleased)
   in downstream libraries using the ``pygeos.strtree.BinaryPredicate`` enum.
   This will be removed in a future release.
 * ``points``, ``linestrings``, ``linearrings``, and ``polygons`` now return a ``GEOSException``
-  instead of a ``ValueError`` for invalid input (#310).
+  instead of a ``ValueError`` or ``TypeError`` for invalid input (#310, #326).
 * Addition of ``on_invalid`` parameter to ``from_wkb`` and ``from_wkt`` to
   optionally return invalid WKB geometries as ``None``.
+* Removed the (internal) function ``lib.polygons_without_holes`` and renamed
+  ``lib.polygons_with_holes`` to ``lib.polygons`` (#326).
 
 **Added GEOS functions**
 
@@ -47,6 +44,8 @@ Version 0.10 (unreleased)
 * Fixed portability issue for ARM architecture (#293)
 * Fixed segfault in ``linearrings`` and ``box`` when constructing a geometry with nan
   coordinates (#310).
+* Fixed segfault in ``polygons`` (with holes) when None was provided.
+* Fixed memory leak in ``polygons`` when non-linearring input was provided.
 
 **Acknowledgments**
 

pygeos/_geometry.pyx

@@ -23,6 +23,7 @@ from pygeos._geos cimport (
     GEOSGeom_createLinearRing_r,
     GEOSGeom_createLineString_r,
     GEOSGeom_createPoint_r,
+    GEOSGeom_createPolygon_r,
     GEOSGeom_destroy_r,
     GEOSGeometry,
     GEOSGeomTypeId_r,
@@ -92,7 +93,7 @@ def simple_geometries_1d(object coordinates, object indices, int geometry_type):
     cdef double[:, :] coord_view = coordinates
     cdef np.intp_t[:] index_view = indices
 
-    # get the geometry count per collection
+    # get the geometry count per collection (this raises on negative indices)
     cdef unsigned int[:] coord_counts = np.bincount(indices).astype(np.uint32)
 
     # The final target array
@@ -192,6 +193,7 @@ def get_parts(object[:] array):
                 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)
@@ -259,7 +261,7 @@ def collections_1d(object geometries, object indices, int geometry_type = 7):
     cdef object[:] geometries_view = geometries
     cdef int[:] indices_view = indices
 
-    # get the geometry count per collection
+    # get the geometry count per collection (this raises on negative indices)
     cdef int[:] collection_size = np.bincount(indices).astype(np.int32)
 
     # A temporary array for the geometries that will be given to CreateCollection.
@@ -323,3 +325,142 @@ def collections_1d(object geometries, object indices, int geometry_type = 7):
             geom_idx_1 += collection_size[coll_idx]
 
     return result
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def polygons_1d(object shells, object holes, object indices):
+    cdef Py_ssize_t hole_idx_1 = 0
+    cdef Py_ssize_t poly_idx = 0
+    cdef unsigned int n_holes = 0
+    cdef int geom_type = 0
+    cdef Py_ssize_t poly_hole_idx = 0
+    cdef GEOSGeometry *shell = NULL
+    cdef GEOSGeometry *hole = NULL
+    cdef GEOSGeometry *poly = NULL
+
+    # Cast input arrays and define memoryviews for later usage
+    shells = np.asarray(shells, dtype=object)
+    if shells.ndim != 1:
+        raise TypeError("shells is not a one-dimensional array.")
+
+    # Cast input arrays and define memoryviews for later usage
+    holes = np.asarray(holes, dtype=object)
+    if holes.ndim != 1:
+        raise TypeError("holes is not a one-dimensional array.")
+
+    indices = np.asarray(indices, dtype=np.int32)
+    if indices.ndim != 1:
+        raise TypeError("indices is not a one-dimensional array.")
+
+    if holes.shape[0] != indices.shape[0]:
+        raise ValueError("holes and indices do not have equal size.")
+
+    if shells.shape[0] == 0:
+        # return immediately if there are no geometries to return
+        return np.empty(shape=(0, ), dtype=object)
+
+    if (indices >= shells.shape[0]).any():
+        raise ValueError("Some indices are of bounds of the shells array.")  
+
+    if np.any(indices[1:] < indices[:indices.shape[0] - 1]):
+        raise ValueError("The indices should be sorted.")  
+
+    cdef Py_ssize_t n_poly = shells.shape[0]
+    cdef object[:] shells_view = shells
+    cdef object[:] holes_view = holes
+    cdef int[:] indices_view = indices
+
+    # get the holes count per polygon (this raises on negative indices)
+    cdef int[:] hole_count = np.bincount(indices, minlength=n_poly).astype(np.int32)
+
+    # A temporary array for the holes that will be given to CreatePolygon
+    # Its size equals max(hole_count) to accomodate the largest polygon.
+    temp_holes = np.empty(shape=(np.max(hole_count), ), dtype=np.intp)
+    cdef np.intp_t[:] temp_holes_view = temp_holes
+
+    # The final target array
+    result = np.empty(shape=(n_poly, ), dtype=object)
+    cdef object[:] result_view = result
+
+    with get_geos_handle() as geos_handle:
+        for poly_idx in range(n_poly):
+            n_holes = 0
+
+            # get the shell
+            if PyGEOS_GetGEOSGeometry(<PyObject *>shells_view[poly_idx], &shell) == 0:
+                raise TypeError(
+                    "One of the arguments is of incorrect type. Please provide only Geometry objects."
+                )
+
+            # return None for missing shells (ignore possibly present holes)
+            if shell == NULL:
+                result_view[poly_idx] = PyGEOS_CreateGeometry(NULL, geos_handle)
+                hole_idx_1 += hole_count[poly_idx]
+                continue
+
+            geom_type = GEOSGeomTypeId_r(geos_handle, shell)
+            if geom_type == -1:
+                return  # GEOSException is raised by get_geos_handle
+            elif geom_type != 2:
+                raise TypeError(
+                    f"One of the shells has unexpected geometry type {geom_type}."
+                )
+
+            # fill the temporary array with holes belonging to this polygon
+            for poly_hole_idx in range(hole_count[poly_idx]):
+                if PyGEOS_GetGEOSGeometry(<PyObject *>holes_view[hole_idx_1 + poly_hole_idx], &hole) == 0:
+                    _deallocate_arr(geos_handle, temp_holes_view, n_holes)
+                    raise TypeError(
+                        "One of the arguments is of incorrect type. Please provide only Geometry objects."
+                    )
+
+                # ignore missing holes
+                if hole == NULL:
+                    continue
+
+                # check the type
+                geom_type = GEOSGeomTypeId_r(geos_handle, hole)
+                if geom_type == -1:
+                    _deallocate_arr(geos_handle, temp_holes_view, n_holes)
+                    return  # GEOSException is raised by get_geos_handle
+                elif geom_type != 2:
+                    _deallocate_arr(geos_handle, temp_holes_view, n_holes)
+                    raise TypeError(
+                        f"One of the holes has unexpected geometry type {geom_type}."
+                    )
+
+                # assign to the temporary geometry array  
+                hole = GEOSGeom_clone_r(geos_handle, hole)
+                if hole == NULL:
+                    _deallocate_arr(geos_handle, temp_holes_view, n_holes)
+                    return  # GEOSException is raised by get_geos_handle                 
+                temp_holes_view[n_holes] = <np.intp_t>hole
+                n_holes += 1
+
+            # clone the shell as the polygon will take ownership
+            shell = GEOSGeom_clone_r(geos_handle, shell)
+            if shell == NULL:
+                _deallocate_arr(geos_handle, temp_holes_view, n_holes)
+                return  # GEOSException is raised by get_geos_handle
+
+            # create the polygon
+            poly = GEOSGeom_createPolygon_r(
+                geos_handle,
+                shell, 
+                <GEOSGeometry**> &temp_holes_view[0],
+                n_holes
+            )
+            if poly == NULL:
+                # GEOSGeom_createPolygon_r should take ownership of the input geometries,
+                # but it doesn't in case of an exception. We prefer a memory leak here over
+                # a possible segfault in the future. The pre-emptive type check already covers
+                # all known cases that GEOSGeom_createPolygon_r errors, so this should never
+                # happen anyway. https://trac.osgeo.org/geos/ticket/1111.
+                return  # GEOSException is raised by get_geos_handle
+
+            result_view[poly_idx] = PyGEOS_CreateGeometry(poly, geos_handle)
+
+            hole_idx_1 += hole_count[poly_idx]
+
+    return result

pygeos/_geos.pxd

@@ -29,6 +29,7 @@ cdef extern from "geos_c.h":
     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_createPolygon_r(GEOSContextHandle_t handle, GEOSGeometry* shell, GEOSGeometry** holes, unsigned int nholes) 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
 

pygeos/creation.py

@@ -1,7 +1,7 @@
 import numpy as np
 
 from . import Geometry, GeometryType, lib
-from ._geometry import collections_1d, simple_geometries_1d
+from ._geometry import collections_1d, polygons_1d, simple_geometries_1d
 from .decorators import multithreading_enabled
 
 __all__ = [
@@ -44,10 +44,10 @@ def points(coords, y=None, z=None, indices=None, **kwargs):
     y : array_like
     z : array_like
     indices : array_like or None
-       Indices into the target array where input coordinates belong. If
-       provided, the coords should be 2D with shape (N, 2) or (N, 3) and
-       indices should be 1D with shape (N,). Missing indices will give None
-       values in the output array.
+        Indices into the target array where input coordinates belong. If
+        provided, the coords should be 2D with shape (N, 2) or (N, 3) and
+        indices should be 1D with shape (N,). Missing indices will give None
+        values in the output array.
     **kwargs
         For other keyword-only arguments, see the
         `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
@@ -75,10 +75,10 @@ def linestrings(coords, y=None, z=None, indices=None, **kwargs):
     y : array_like
     z : array_like
     indices : array_like or None
-       Indices into the target array where input coordinates belong. If
-       provided, the coords should be 2D with shape (N, 2) or (N, 3) and
-       indices should be 1D with shape (N,). Missing indices will give None
-       values in the output array.
+        Indices into the target array where input coordinates belong. If
+        provided, the coords should be 2D with shape (N, 2) or (N, 3) and
+        indices should be 1D with shape (N,). Missing indices will give None
+        values in the output array.
     **kwargs
         For other keyword-only arguments, see the
         `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
@@ -108,10 +108,10 @@ def linearrings(coords, y=None, z=None, indices=None, **kwargs):
     y : array_like
     z : array_like
     indices : array_like or None
-       Indices into the target array where input coordinates belong. If
-       provided, the coords should be 2D with shape (N, 2) or (N, 3) and
-       indices should be 1D with shape (N,). Missing indices will give None
-       values in the output array.
+        Indices into the target array where input coordinates belong. If
+        provided, the coords should be 2D with shape (N, 2) or (N, 3) and
+        indices should be 1D with shape (N,). Missing indices will give None
+        values in the output array.
     **kwargs
         For other keyword-only arguments, see the
         `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
@@ -125,30 +125,48 @@ def linearrings(coords, y=None, z=None, indices=None, **kwargs):
 
 
 @multithreading_enabled
-def polygons(shells, holes=None):
+def polygons(shells, holes=None, indices=None, **kwargs):
     """Create an array of polygons.
 
     Parameters
     ----------
     shell : array_like
         An array of linearrings that constitute the out shell of the polygons.
         Coordinates can also be passed, see linearrings.
-    holes : array_like
+    holes : array_like or None
         An array of lists of linearrings that constitute holes for each shell.
+    indices : array_like or None
+        Indices into the shells array where input holes belong. If
+        provided, shells, holes, and indices should all be 1D and the size
+        of holes must equal the size of indices.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+        Ignored if ``indices`` is provided.
     """
     shells = np.asarray(shells)
     if not isinstance(shells, Geometry) and np.issubdtype(shells.dtype, np.number):
         shells = linearrings(shells)
 
-    if holes is None:
-        return lib.polygons_without_holes(shells)
+    if holes is None and indices is not None:
+        raise TypeError("Indices provided without a holes array.")
+    elif holes is None:
+        # no holes provided: initialize an empty holes array matching shells
+        shape = shells.shape + (0,) if isinstance(shells, np.ndarray) else (0,)
+        holes = np.empty(shape, dtype=object)
+    else:
+        holes = np.asarray(holes)
+        # convert holes coordinates into linearrings
+        if np.issubdtype(holes.dtype, np.number):
+            holes = linearrings(holes)
 
-    holes = np.asarray(holes)
-    if not isinstance(holes, Geometry) and np.issubdtype(holes.dtype, np.number):
-        holes = linearrings(holes)
-    return lib.polygons_with_holes(shells, holes)
+    if indices is None:
+        return lib.polygons(shells, holes, **kwargs)
+    else:
+        return polygons_1d(shells, holes, indices)
 
 
+@multithreading_enabled
 def box(xmin, ymin, xmax, ymax, ccw=True, **kwargs):
     """Create box polygons.
 
@@ -178,6 +196,7 @@ def box(xmin, ymin, xmax, ymax, ccw=True, **kwargs):
     return lib.box(xmin, ymin, xmax, ymax, ccw, **kwargs)
 
 
+@multithreading_enabled
 def multipoints(geometries, indices=None, **kwargs):
     """Create multipoints from arrays of points
 
@@ -186,7 +205,7 @@ def multipoints(geometries, indices=None, **kwargs):
     geometries : array_like
         An array of points or coordinates (see points).
     indices : array_like or None
-       Indices into the target array where input geometries belong. If
+        Indices into the target array where input geometries belong. If
         provided, both geometries and indices should be 1D and have matching
         sizes.
     **kwargs
@@ -206,6 +225,7 @@ def multipoints(geometries, indices=None, **kwargs):
         return collections_1d(geometries, indices, typ)
 
 
+@multithreading_enabled
 def multilinestrings(geometries, indices=None, **kwargs):
     """Create multilinestrings from arrays of linestrings
 
@@ -235,6 +255,7 @@ def multilinestrings(geometries, indices=None, **kwargs):
         return collections_1d(geometries, indices, typ)
 
 
+@multithreading_enabled
 def multipolygons(geometries, indices=None, **kwargs):
     """Create multipolygons from arrays of polygons
 
@@ -263,6 +284,7 @@ def multipolygons(geometries, indices=None, **kwargs):
         return collections_1d(geometries, indices, typ)
 
 
+@multithreading_enabled
 def geometrycollections(geometries, indices=None, **kwargs):
     """Create geometrycollections from arrays of geometries