Skip to content

Commit

Permalink
Merge polygons_1d into collections_1d + finish implementation for rls…
Browse files Browse the repository at this point in the history 
… 0.10 (#346)
  • Loading branch information
@caspervdw
caspervdw committed May 6, 2021
1 parent 9a959d1 commit 8a50f86
Show file tree
Hide file tree
Showing 7 changed files with 312 additions and 250 deletions.
4 changes: 3 additions & 1 deletion CHANGELOG.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ 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 all creation functions (#230, #322, #326).
by optionally taking index arrays in all creation functions (#230, #322, #326, #346).
* 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
Expand All @@ -31,8 +31,10 @@ Version 0.10 (unreleased)
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).
* ``polygons`` will now return an empty polygon for `None` inputs (#346).
* Removed compatibility with Python 3.5 (#341).


**Added GEOS functions**

* Addition of a ``contains_properly`` function (#267)
Expand Down
199 changes: 43 additions & 156 deletions pygeos/_geometry.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ from pygeos._geos cimport (
GEOSCoordSequence,
GEOSGeom_clone_r,
GEOSGeom_createCollection_r,
GEOSGeom_createEmptyPolygon_r,
GEOSGeom_createLinearRing_r,
GEOSGeom_createLineString_r,
GEOSGeom_createPoint_r,
Expand Down Expand Up @@ -69,7 +70,7 @@ def simple_geometries_1d(object coordinates, object indices, int geometry_type):
if coordinates.ndim != 2:
raise TypeError("coordinates is not a two-dimensional array.")

indices = np.asarray(indices, dtype=np.intp)
indices = np.asarray(indices, dtype=np.intp) # intp is what bincount takes
if indices.ndim != 1:
raise TypeError("indices is not a one-dimensional array.")

Expand All @@ -91,7 +92,6 @@ def simple_geometries_1d(object coordinates, object indices, int geometry_type):
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 (this raises on negative indices)
cdef unsigned int[:] coord_counts = np.bincount(indices).astype(np.uint32)
Expand All @@ -105,10 +105,11 @@ def simple_geometries_1d(object coordinates, object indices, int geometry_type):
for geom_idx in range(n_geoms):
geom_size = coord_counts[geom_idx]

# insert None if there are no coordinates
# for now, raise if there are indices missing (decision on this in GH345)
if geom_size == 0:
result_view[geom_idx] = PyGEOS_CreateGeometry(NULL, geos_handle)
continue
raise ValueError(
f"Index {geom_idx} is missing from the input indices."
)

# check if we need to close a linearring
if geometry_type == 2:
Expand Down Expand Up @@ -217,6 +218,15 @@ cdef _deallocate_arr(void* handle, np.intp_t[:] arr, Py_ssize_t last_geom_i):
@cython.boundscheck(False)
@cython.wraparound(False)
def collections_1d(object geometries, object indices, int geometry_type = 7):
"""Converts geometries + indices to collections
Allowed geometry type conversions are:
- linearrings to polygons
- points to multipoints
- linestrings/linearrings to multilinestrings
- polygons to multipolygons
- any to geometrycollections
"""
cdef Py_ssize_t geom_idx_1 = 0
cdef Py_ssize_t coll_idx = 0
cdef unsigned int coll_size = 0
Expand All @@ -227,7 +237,9 @@ def collections_1d(object geometries, object indices, int geometry_type = 7):
cdef int expected_type_alt = -1
cdef int curr_type = -1

if geometry_type == 4: # MULTIPOINT
if geometry_type == 3: # POLYGON
expected_type = 2
elif geometry_type == 4: # MULTIPOINT
expected_type = 0
elif geometry_type == 5: # MULTILINESTRING
expected_type = 1
Expand All @@ -244,7 +256,7 @@ def collections_1d(object geometries, object indices, int geometry_type = 7):
if geometries.ndim != 1:
raise TypeError("geometries is not a one-dimensional array.")

indices = np.asarray(indices, dtype=np.int32)
indices = np.asarray(indices, dtype=np.intp) # intp is what bincount takes
if indices.ndim != 1:
raise TypeError("indices is not a one-dimensional array.")

Expand All @@ -258,9 +270,6 @@ def collections_1d(object geometries, object indices, int geometry_type = 7):
if np.any(indices[1:] < indices[:indices.shape[0] - 1]):
raise ValueError("The indices should be sorted.")

cdef object[:] geometries_view = geometries
cdef int[:] indices_view = indices

# get the geometry count per collection (this raises on negative indices)
cdef int[:] collection_size = np.bincount(indices).astype(np.int32)

Expand All @@ -276,11 +285,15 @@ def collections_1d(object geometries, object indices, int geometry_type = 7):

with get_geos_handle() as geos_handle:
for coll_idx in range(n_colls):
if collection_size[coll_idx] == 0:
raise ValueError(
f"Index {coll_idx} is missing from the input indices."
)
coll_size = 0

# 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:
if PyGEOS_GetGEOSGeometry(<PyObject *>geometries[geom_idx_1 + coll_geom_idx], &geom) == 0:
_deallocate_arr(geos_handle, temp_geoms_view, coll_size)
raise TypeError(
"One of the arguments is of incorrect type. Please provide only Geometry objects."
Expand Down Expand Up @@ -311,156 +324,30 @@ def collections_1d(object geometries, object indices, int geometry_type = 7):
coll_size += 1

# create the collection
coll = GEOSGeom_createCollection_r(
geos_handle,
geometry_type,
<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)

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."
if geometry_type != 3: # Collection
coll = GEOSGeom_createCollection_r(
geos_handle,
geometry_type,
<GEOSGeometry**> &temp_geoms_view[0],
coll_size
)

# 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}."
elif coll_size != 0: # Polygon, non-empty
coll = GEOSGeom_createPolygon_r(
geos_handle,
<GEOSGeometry*> temp_geoms_view[0],
NULL if coll_size <= 1 else <GEOSGeometry**> &temp_geoms_view[1],
coll_size - 1
)
else: # Polygon, empty
coll = GEOSGeom_createEmptyPolygon_r(
geos_handle
)

# 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.
if coll == NULL:
return # GEOSException is raised by get_geos_handle

result_view[poly_idx] = PyGEOS_CreateGeometry(poly, geos_handle)
result_view[coll_idx] = PyGEOS_CreateGeometry(coll, geos_handle)

hole_idx_1 += hole_count[poly_idx]
geom_idx_1 += collection_size[coll_idx]

return result
1 change: 1 addition & 0 deletions pygeos/_geos.pxd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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_createEmptyPolygon_r(GEOSContextHandle_t handle) 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
Expand Down

0 comments on commit 8a50f86

Please sign in to comment.