[ENH] Add bounds function (#69)

pygeos/measurement.py

@@ -1,10 +1,12 @@
+import numpy as np
+
 from . import lib
 from . import Geometry  # NOQA
 
-__all__ = ["area", "distance", "length", "hausdorff_distance"]
+__all__ = ["area", "distance", "bounds", "length", "hausdorff_distance"]
 
 
-def area(geometry):
+def area(geometry, **kwargs):
     """Computes the area of a (multi)polygon.
 
     Parameters
@@ -22,10 +24,10 @@ def area(geometry):
     >>> area(None)
     nan
     """
-    return lib.area(geometry)
+    return lib.area(geometry, **kwargs)
 
 
-def distance(a, b):
+def distance(a, b, **kwargs):
     """Computes the Cartesian distance between two geometries.
 
     Parameters
@@ -46,10 +48,37 @@ def distance(a, b):
     >>> distance(None, point)
     nan
     """
-    return lib.distance(a, b)
+    return lib.distance(a, b, **kwargs)
+
+
+def bounds(geometry, **kwargs):
+    """Computes the bounds (extent) of a geometry.
+
+    For each geometry these 4 numbers are returned: min x, min y, max x, max y.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+
+    Examples
+    --------
+    >>> bounds(Geometry("POINT (2 3)")).tolist()
+    [2.0, 3.0, 2.0, 3.0]
+    >>> bounds(Geometry("LINESTRING (0 0, 0 2, 3 2)")).tolist()
+    [0.0, 0.0, 3.0, 2.0]
+    >>> bounds(Geometry("POLYGON EMPTY")).tolist()
+    [nan, nan, nan, nan]
+    >>> bounds(None).tolist()
+    [nan, nan, nan, nan]
+    """
+    # We need to provide the `out` argument here for compatibility with
+    # numpy < 1.16. See https://github.com/numpy/numpy/issues/14949
+    geometry_arr = np.asarray(geometry, dtype=np.object)
+    out = np.empty(geometry_arr.shape + (4,), dtype="float64")
+    return lib.bounds(geometry_arr, out=out, **kwargs)
 
 
-def length(geometry):
+def length(geometry, **kwargs):
     """Computes the length of a (multi)linestring or polygon perimeter.
 
     Parameters
@@ -69,10 +98,10 @@ def length(geometry):
     >>> length(None)
     nan
     """
-    return lib.length(geometry)
+    return lib.length(geometry, **kwargs)
 
 
-def hausdorff_distance(a, b, densify=None):
+def hausdorff_distance(a, b, densify=None, **kwargs):
     """Compute the discrete Haussdorf distance between two geometries.
 
     The Haussdorf distance is a measure of similarity: it is the greatest
@@ -101,6 +130,6 @@ def hausdorff_distance(a, b, densify=None):
     nan
     """
     if densify is None:
-        return lib.hausdorff_distance(a, b)
+        return lib.hausdorff_distance(a, b, **kwargs)
     else:
-        return lib.haussdorf_distance_densify(a, b, densify)
+        return lib.haussdorf_distance_densify(a, b, densify, **kwargs)

pygeos/test/test_measurement.py

@@ -46,6 +46,37 @@ def test_distance_missing():
     assert np.isnan(actual)
 
 
+@pytest.mark.parametrize(
+    "geom,expected",
+    [
+        (point, [2.0, 3.0, 2.0, 3.0]),
+        (pygeos.linestrings([[0, 0], [0, 1]]), [0.0, 0.0, 0.0, 1.0]),
+        (pygeos.linestrings([[0, 0], [1, 0]]), [0.0, 0.0, 1.0, 0.0]),
+        (multi_point, [0.0, 0.0, 1.0, 2.0]),
+        (multi_polygon, [0.0, 0.0, 2.2, 2.2]),
+        (geometry_collection, [49.0, -1.0, 52.0, 2.0]),
+    ],
+)
+def test_bounds(geom, expected):
+    actual = pygeos.bounds(geom)
+    assert actual.tolist() == expected
+
+
+def test_bounds_array():
+    actual = pygeos.bounds([[point, multi_point], [polygon, None]])
+    assert actual.shape == (2, 2, 4)
+
+
+def test_bounds_missing():
+    actual = pygeos.bounds(None)
+    assert np.isnan(actual).all()
+
+
+def test_bounds_empty():
+    actual = pygeos.bounds(empty)
+    assert np.isnan(actual).all()
+
+
 def test_length():
     actual = pygeos.length(
         [

src/ufuncs.c

@@ -934,6 +934,71 @@ static void create_collection_func(char **args, npy_intp *dimensions,
 static PyUFuncGenericFunction create_collection_funcs[1] = {&create_collection_func};
 
 
+static char bounds_dtypes[2] = {NPY_OBJECT, NPY_DOUBLE};
+static void bounds_func(char **args, npy_intp *dimensions,
+                        npy_intp *steps, void *data)
+{
+    void *context = geos_context[0];
+    GEOSGeometry *envelope, *in1;
+    const GEOSGeometry *ring;
+    const GEOSCoordSequence *coord_seq;
+    int size;
+    char *ip1 = args[0], *op1 = args[1];
+    double *x1, *y1, *x2, *y2;
+
+    npy_intp is1 = steps[0], os1 = steps[1], cs1 = steps[2];
+    npy_intp n = dimensions[0], i;
+
+    for(i = 0; i < n; i++, ip1 += is1, op1 += os1) {
+        if (!get_geom(*(GeometryObject **)ip1, &in1)) { return; }
+
+        /* get the 4 (pointers to) the bbox values from the "core stride 1" (cs1) */
+        x1 = (double *)(op1);
+        y1 = (double *)(op1 + cs1);
+        x2 = (double *)(op1 + 2 * cs1);
+        y2 = (double *)(op1 + 3 * cs1);
+
+        if (in1 == NULL) {  /* no geometry => bbox becomes (nan, nan, nan, nan) */
+            *x1 = *y1 = *x2 = *y2 = NPY_NAN;
+        } else {
+            /* construct the envelope */
+            envelope = GEOSEnvelope_r(context, in1);
+            if (envelope == NULL) { return; }
+            size = GEOSGetNumCoordinates_r(context, envelope);
+
+            /* get the bbox depending on the number of coordinates in the envelope */
+            if (size == 0) {  /* Envelope is empty */
+                *x1 = *y1 = *x2 = *y2 = NPY_NAN;
+            } else if (size == 1) {  /* Envelope is a point */
+                if (!GEOSGeomGetX_r(context, envelope, x1)) { goto fail; }
+                if (!GEOSGeomGetY_r(context, envelope, y1)) { goto fail; }
+                *x2 = *x1;
+                *y2 = *y1;
+            } else if (size == 5) {  /* Envelope is a box */
+                ring = GEOSGetExteriorRing_r(context, envelope);
+                if (ring == NULL) { goto fail; }
+                coord_seq = GEOSGeom_getCoordSeq_r(context, ring);
+                if (coord_seq == NULL) { goto fail; }
+                if (!GEOSCoordSeq_getX_r(context, coord_seq, 0, x1)) { goto fail; }
+                if (!GEOSCoordSeq_getY_r(context, coord_seq, 0, y1)) { goto fail; }
+                if (!GEOSCoordSeq_getX_r(context, coord_seq, 2, x2)) { goto fail; }
+                if (!GEOSCoordSeq_getY_r(context, coord_seq, 2, y2)) { goto fail; }
+            } else {
+                PyErr_Format(PyExc_ValueError, "Could not determine bounds from an envelope with %d coordinate pairs", size);
+                goto fail;
+            }
+            GEOSGeom_destroy_r(context, envelope);
+        }
+    }
+
+    return;
+    fail:
+        GEOSGeom_destroy_r(context, envelope);
+        return;
+}
+static PyUFuncGenericFunction bounds_funcs[1] = {&bounds_func};
+
+
 
 /* Define the object -> geom functions (O_Y) */
 
@@ -1323,6 +1388,7 @@ int init_ufuncs(PyObject *m, PyObject *d)
     DEFINE_GENERALIZED(points, 1, "(d)->()");
     DEFINE_GENERALIZED(linestrings, 1, "(i, d)->()");
     DEFINE_GENERALIZED(linearrings, 1, "(i, d)->()");
+    DEFINE_GENERALIZED(bounds, 1, "()->(n)");
     DEFINE_Y_Y (polygons_without_holes);
     DEFINE_GENERALIZED(polygons_with_holes, 2, "(),(i)->()");
     DEFINE_GENERALIZED(create_collection, 2, "(i),()->()");