Explode polygons into rings: get_rings (#342)

CHANGELOG.rst

@@ -12,6 +12,8 @@ Version 0.10 (unreleased)
   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).
+* Added the ``get_rings`` function, similar as ``get_parts`` but specifically
+  to extract the rings of Polygon geometries (#342).
 * Updated ``box`` ufunc to use internal C function for creating polygon
   (about 2x faster) and added ``ccw`` parameter to create polygon in
   counterclockwise (default) or clockwise direction (#308).

pygeos/_geometry.pyx

@@ -28,7 +28,9 @@ from pygeos._geos cimport (
     GEOSGeom_destroy_r,
     GEOSGeometry,
     GEOSGeomTypeId_r,
+    GEOSGetExteriorRing_r,
     GEOSGetGeometryN_r,
+    GEOSGetInteriorRingN_r,
     get_geos_handle,
 )
 from pygeos._pygeos_api cimport (
@@ -145,17 +147,33 @@ def simple_geometries_1d(object coordinates, object indices, int geometry_type):
     return result
 
 
+
+cdef const GEOSGeometry* GetRingN(GEOSContextHandle_t handle, GEOSGeometry* polygon, int n):
+    if n == 0:
+        return GEOSGetExteriorRing_r(handle, polygon)
+    else:
+        return GEOSGetInteriorRingN_r(handle, polygon, n - 1)
+
+
+
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def get_parts(object[:] array):
+def get_parts(object[:] array, bint extract_rings=0):
     cdef Py_ssize_t geom_idx = 0
     cdef Py_ssize_t part_idx = 0
     cdef Py_ssize_t idx = 0
+    cdef Py_ssize_t count
     cdef GEOSGeometry *geom = NULL
     cdef const GEOSGeometry *part = NULL
 
-    counts = pygeos.get_num_geometries(array)
-    cdef Py_ssize_t count = counts.sum()
+    if extract_rings:
+        counts = pygeos.get_num_interior_rings(array)
+        is_polygon = (pygeos.get_type_id(array) == 3) & (~pygeos.is_empty(array))
+        counts += is_polygon
+        count = counts.sum()
+    else:
+        counts = pygeos.get_num_geometries(array)
+        count = counts.sum()
 
     if count == 0:
         # return immediately if there are no geometries to return
@@ -186,7 +204,11 @@ 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 extract_rings:
+                    part = GetRingN(geos_handle, geom, part_idx)
+                else:
+                    part = GEOSGetGeometryN_r(geos_handle, geom, part_idx)
                 if part == NULL:
                     return  # GEOSException is raised by get_geos_handle
 

pygeos/_geos.pxd

@@ -22,6 +22,8 @@ cdef extern from "geos_c.h":
 
     # Geometry functions
     const GEOSGeometry* GEOSGetGeometryN_r(GEOSContextHandle_t handle, const GEOSGeometry* g, int n) nogil
+    const GEOSGeometry* GEOSGetExteriorRing_r(GEOSContextHandle_t handle, const GEOSGeometry* g) nogil
+    const GEOSGeometry* GEOSGetInteriorRingN_r(GEOSContextHandle_t handle, const GEOSGeometry* g, int n) nogil
     int GEOSGeomTypeId_r(GEOSContextHandle_t handle, GEOSGeometry* g) nogil
 
     # Geometry creation / destruction

pygeos/geometry.py

@@ -25,6 +25,7 @@
     "get_interior_ring",
     "get_geometry",
     "get_parts",
+    "get_rings",
     "get_precision",
     "set_precision",
 ]
@@ -494,7 +495,7 @@ def get_geometry(geometry, index, **kwargs):
 
     See also
     --------
-    get_num_geometries
+    get_num_geometries, get_parts
 
     Examples
     --------
@@ -525,13 +526,17 @@ def get_parts(geometry, return_index=False):
     ----------
     geometry : Geometry or array_like
     return_index : bool, optional (default: False)
-        If True, will return a tuple of ndarrys of (parts, indexes), where indexes
+        If True, will return a tuple of ndarrays of (parts, indexes), where indexes
         are the indexes of the original geometries in the source array.
 
     Returns
     -------
     ndarray of parts or tuple of (parts, indexes)
 
+    See also
+    --------
+    get_geometry, get_rings
+
     Examples
     --------
     >>> get_parts(Geometry("MULTIPOINT (0 1, 2 3)")).tolist()
@@ -554,6 +559,61 @@ def get_parts(geometry, return_index=False):
     return _geometry.get_parts(geometry)[0]
 
 
+def get_rings(geometry, return_index=False):
+    """Gets rings of Polygon geometry object.
+
+    For each Polygon, the first returned ring is always the exterior ring
+    and potential subsequent rings are interior rings.
+
+    If the geometry is not a Polygon, nothing is returned (empty array for
+    scalar geometry input or no element in output array for array input).
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    return_index : bool, optional (default: False)
+        If True, will return a tuple of ndarrays of (rings, indexes), where
+        indexes are the indexes of the original geometries in the source array.
+
+    Returns
+    -------
+    ndarray of rings or tuple of (rings, indexes)
+
+    See also
+    --------
+    get_exterior_ring, get_interior_ring, get_parts
+
+    Examples
+    --------
+    >>> polygon_with_hole = Geometry("POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0), \
+(2 2, 2 4, 4 4, 4 2, 2 2))")
+    >>> get_rings(polygon_with_hole).tolist()
+    [<pygeos.Geometry LINEARRING (0 0, 0 10, 10 10, 10 0, 0 0)>,
+     <pygeos.Geometry LINEARRING (2 2, 2 4, 4 4, 4 2, 2 2)>]
+
+    With ``return_index=True``:
+
+    >>> polygon = Geometry("POLYGON ((0 0, 2 0, 2 2, 0 2, 0 0))")
+    >>> rings, index = get_rings([polygon, polygon_with_hole], return_index=True)
+    >>> rings.tolist()
+    [<pygeos.Geometry LINEARRING (0 0, 2 0, 2 2, 0 2, 0 0)>,
+     <pygeos.Geometry LINEARRING (0 0, 0 10, 10 10, 10 0, 0 0)>,
+     <pygeos.Geometry LINEARRING (2 2, 2 4, 4 4, 4 2, 2 2)>]
+    >>> index.tolist()
+    [0, 1, 1]
+    """
+    geometry = np.asarray(geometry, dtype=np.object_)
+    geometry = np.atleast_1d(geometry)
+
+    if geometry.ndim != 1:
+        raise ValueError("Array should be one dimensional")
+
+    if return_index:
+        return _geometry.get_parts(geometry, extract_rings=True)
+
+    return _geometry.get_parts(geometry, extract_rings=True)[0]
+
+
 @multithreading_enabled
 def get_num_geometries(geometry, **kwargs):
     """Returns number of geometries in a collection.

pygeos/test/test_geometry.py

@@ -358,6 +358,68 @@ def test_get_parts_invalid_geometry(geom):
         pygeos.get_parts(geom)
 
 
+@pytest.mark.parametrize(
+    "geom",
+    [
+        point,
+        multi_point,
+        line_string,
+        multi_line_string,
+        polygon,
+        multi_polygon,
+        geometry_collection,
+        empty_point,
+        empty_line_string,
+        empty_polygon,
+        empty_geometry_collection,
+        None,
+    ],
+)
+def test_get_rings(geom):
+    if (pygeos.get_type_id(geom) != pygeos.GeometryType.POLYGON) or pygeos.is_empty(
+        geom
+    ):
+        rings = pygeos.get_rings(geom)
+        assert len(rings) == 0
+    else:
+        rings = pygeos.get_rings(geom)
+        assert len(rings) == 1
+        assert rings[0] == pygeos.get_exterior_ring(geom)
+
+
+def test_get_rings_holes():
+    rings = pygeos.get_rings(polygon_with_hole)
+    assert len(rings) == 2
+    assert rings[0] == pygeos.get_exterior_ring(polygon_with_hole)
+    assert rings[1] == pygeos.get_interior_ring(polygon_with_hole, 0)
+
+
+def test_get_rings_return_index():
+    geom = np.array([polygon, None, empty_polygon, polygon_with_hole])
+    expected_parts = []
+    expected_index = []
+    for i, g in enumerate(geom):
+        if g is None or pygeos.is_empty(g):
+            continue
+        expected_parts.append(pygeos.get_exterior_ring(g))
+        expected_index.append(i)
+        for j in range(0, pygeos.get_num_interior_rings(g)):
+            expected_parts.append(pygeos.get_interior_ring(g, j))
+            expected_index.append(i)
+
+    parts, index = pygeos.get_rings(geom, return_index=True)
+    assert len(parts) == len(expected_parts)
+    assert np.all(pygeos.equals_exact(parts, expected_parts))
+    assert np.array_equal(index, expected_index)
+
+
+@pytest.mark.parametrize("geom", [[[None]], [[polygon]]])
+def test_get_rings_invalid_dimensions(geom):
+    """Only 1D inputs are supported"""
+    with pytest.raises(ValueError, match="Array should be one dimensional"):
+        pygeos.get_parts(geom)
+
+
 @pytest.mark.skipif(pygeos.geos_version < (3, 6, 0), reason="GEOS < 3.6")
 def test_get_precision():
     geometries = all_types + (point_z, empty_point, empty_line_string, empty_polygon)

setup.cfg

@@ -11,3 +11,6 @@ ignore = E203, E266, E501, W503
 [tool:isort]
 profile = black
 force_alphabetical_sort_within_sections = true
+
+[tool:pytest]
+doctest_optionflags = NORMALIZE_WHITESPACE