ENH: add minimum_bounding_circle and minimum_bounding_radius (pygeos/…

…pygeos#315)

CHANGELOG.rst

@@ -38,6 +38,7 @@ Version 0.10 (unreleased)
 * Addition of a ``polygonize`` function (#275)
 * Addition of a ``polygonize_full`` function (#298)
 * Addition of a ``segmentize`` function for GEOS >= 3.10 (#299)
+* Addition of ``minimum_bounding_circle`` and ``minimum_bounding_radius`` functions for GEOS >= 3.8 (#315)
 
 **Bug fixes**
 
@@ -55,6 +56,7 @@ People with a "+" by their names contributed a patch for the first time.
 * Brendan Ward
 * Casper van der Wel
 * Joris Van den Bossche
+* Martin Fleischmann
 * 0phoff +
 
 

pygeos/constructive.py

@@ -28,6 +28,7 @@
     "simplify",
     "snap",
     "voronoi_polygons",
+    "minimum_bounding_circle",
 ]
 
 
@@ -821,3 +822,35 @@ def voronoi_polygons(
     <pygeos.Geometry GEOMETRYCOLLECTION EMPTY>
     """
     return lib.voronoi_polygons(geometry, tolerance, extend_to, only_edges, **kwargs)
+
+
+@requires_geos("3.8.0")
+@multithreading_enabled
+def minimum_bounding_circle(geometry, **kwargs):
+    """Computes the minimum bounding circle that encloses an input geometry.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> minimum_bounding_circle(Geometry("POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0))"))
+    <pygeos.Geometry POLYGON ((12.1 5, 11.9 3.62, 11.5 2.29, 10.9 1.07, 10 4.44e...>
+    >>> minimum_bounding_circle(Geometry("LINESTRING (1 1, 10 10)"))
+    <pygeos.Geometry POLYGON ((11.9 5.5, 11.7 4.26, 11.4 3.06, 10.8 1.96, 10 1, ...>
+    >>> minimum_bounding_circle(Geometry("MULTIPOINT (2 2, 4 2)"))
+    <pygeos.Geometry POLYGON ((4 2, 3.98 1.8, 3.92 1.62, 3.83 1.44, 3.71 1.29, 3...>
+    >>> minimum_bounding_circle(Geometry("POINT (0 1)"))
+    <pygeos.Geometry POINT (0 1)>
+    >>> minimum_bounding_circle(Geometry("GEOMETRYCOLLECTION EMPTY"))
+    <pygeos.Geometry POLYGON EMPTY>
+
+    See also
+    --------
+    minimum_bounding_radius
+    """
+    return lib.minimum_bounding_circle(geometry, **kwargs)

pygeos/measurement.py

@@ -13,6 +13,7 @@
     "hausdorff_distance",
     "frechet_distance",
     "minimum_clearance",
+    "minimum_bounding_radius",
 ]
 
 
@@ -277,3 +278,36 @@ def minimum_clearance(geometry, **kwargs):
     nan
     """
     return lib.minimum_clearance(geometry, **kwargs)
+
+
+@requires_geos("3.8.0")
+@multithreading_enabled
+def minimum_bounding_radius(geometry, **kwargs):
+    """Computes the radius of the minimum bounding circle that encloses an input geometry.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+
+    Examples
+    --------
+    >>> minimum_bounding_radius(Geometry("POLYGON ((0 5, 5 10, 10 5, 5 0, 0 5))"))
+    5.0
+    >>> minimum_bounding_radius(Geometry("LINESTRING (1 1, 1 10)"))
+    4.5
+    >>> minimum_bounding_radius(Geometry("MULTIPOINT (2 2, 4 2)"))
+    1.0
+    >>> minimum_bounding_radius(Geometry("POINT (0 1)"))
+    0.0
+    >>> minimum_bounding_radius(Geometry("GEOMETRYCOLLECTION EMPTY"))
+    0.0
+
+    See also
+    --------
+    minimum_bounding_circle
+    """
+    return lib.minimum_bounding_radius(geometry, **kwargs)

pygeos/test/test_constructive.py

@@ -687,3 +687,45 @@ def test_segmentize_none():
 def test_segmentize(geometry, tolerance, expected):
     actual = pygeos.segmentize(geometry, tolerance)
     assert pygeos.equals(actual, geometry).all()
+
+
+@pytest.mark.skipif(pygeos.geos_version < (3, 8, 0), reason="GEOS < 3.8")
+@pytest.mark.parametrize("geometry", all_types)
+def test_minimum_bounding_circle_all_types(geometry):
+    actual = pygeos.minimum_bounding_circle([geometry, geometry])
+    assert actual.shape == (2,)
+    assert actual[0] is None or isinstance(actual[0], Geometry)
+
+    actual = pygeos.minimum_bounding_circle(None)
+    assert actual is None
+
+
+@pytest.mark.skipif(pygeos.geos_version < (3, 8, 0), reason="GEOS < 3.8")
+@pytest.mark.parametrize(
+    "geometry, expected",
+    [
+        (
+            pygeos.Geometry("POLYGON ((0 5, 5 10, 10 5, 5 0, 0 5))"),
+            pygeos.buffer(pygeos.Geometry("POINT (5 5)"), 5),
+        ),
+        (
+            pygeos.Geometry("LINESTRING (1 0, 1 10)"),
+            pygeos.buffer(pygeos.Geometry("POINT (1 5)"), 5),
+        ),
+        (
+            pygeos.Geometry("MULTIPOINT (2 2, 4 2)"),
+            pygeos.buffer(pygeos.Geometry("POINT (3 2)"), 1),
+        ),
+        (
+            pygeos.Geometry("POINT (2 2)"),
+            pygeos.Geometry("POINT (2 2)"),
+        ),
+        (
+            pygeos.Geometry("GEOMETRYCOLLECTION EMPTY"),
+            pygeos.Geometry("POLYGON EMPTY"),
+        ),
+    ],
+)
+def test_minimum_bounding_circle(geometry, expected):
+    actual = pygeos.minimum_bounding_circle(geometry)
+    assert pygeos.equals(actual, expected).all()

pygeos/test/test_measurement.py

@@ -303,3 +303,34 @@ def test_minimum_clearance_nonexistent():
 def test_minimum_clearance_missing():
     actual = pygeos.minimum_clearance(None)
     assert np.isnan(actual)
+
+
+@pytest.mark.skipif(pygeos.geos_version < (3, 8, 0), reason="GEOS < 3.8")
+@pytest.mark.parametrize(
+    "geometry, expected",
+    [
+        (
+            pygeos.Geometry("POLYGON ((0 5, 5 10, 10 5, 5 0, 0 5))"),
+            5,
+        ),
+        (
+            pygeos.Geometry("LINESTRING (1 0, 1 10)"),
+            5,
+        ),
+        (
+            pygeos.Geometry("MULTIPOINT (2 2, 4 2)"),
+            1,
+        ),
+        (
+            pygeos.Geometry("POINT (2 2)"),
+            0,
+        ),
+        (
+            pygeos.Geometry("GEOMETRYCOLLECTION EMPTY"),
+            0,
+        ),
+    ],
+)
+def test_minimum_bounding_radius(geometry, expected):
+    actual = pygeos.minimum_bounding_radius(geometry)
+    assert actual == pytest.approx(expected, abs=1e-12)

src/ufuncs.c

@@ -390,6 +390,17 @@ static void* normalize_data[1] = {GEOSNormalize_r_with_clone};
 static void* build_area_data[1] = {GEOSBuildArea_r};
 static void* make_valid_data[1] = {GEOSMakeValid_r};
 static void* coverage_union_data[1] = {GEOSCoverageUnion_r};
+static void* GEOSMinimumBoundingCircleWithReturn(void* context, void* geom) {
+  GEOSGeometry* center = NULL;
+  double radius;
+  GEOSGeometry* ret = GEOSMinimumBoundingCircle_r(context, geom, &radius, &center);
+  if (ret == NULL) {
+    return NULL;
+  }
+  GEOSGeom_destroy_r(context, center);
+  return ret;
+}
+static void* minimum_bounding_circle_data[1] = {GEOSMinimumBoundingCircleWithReturn};
 #endif
 #if GEOS_SINCE_3_7_0
 static void* reverse_data[1] = {GEOSReverse_r};
@@ -921,6 +932,19 @@ static int MinimumClearance(void* context, void* a, double* b) {
 }
 static void* minimum_clearance_data[1] = {MinimumClearance};
 #endif
+#if GEOS_SINCE_3_8_0
+static int GEOSMinimumBoundingRadius(void* context, GEOSGeometry* geom, double* radius) {
+  GEOSGeometry* center = NULL;
+  GEOSGeometry* ret = GEOSMinimumBoundingCircle_r(context, geom, radius, &center);
+  if (ret == NULL) {
+    return 0;  // exception code
+  }
+  GEOSGeom_destroy_r(context, center);
+  GEOSGeom_destroy_r(context, ret);
+  return 1;  // success code
+}
+static void* minimum_bounding_radius_data[1] = {GEOSMinimumBoundingRadius};
+#endif
 typedef int FuncGEOS_Y_d(void* context, void* a, double* b);
 static char Y_d_dtypes[2] = {NPY_OBJECT, NPY_DOUBLE};
 static void Y_d_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
@@ -3080,6 +3104,8 @@ int init_ufuncs(PyObject* m, PyObject* d) {
   DEFINE_Y_Y(make_valid);
   DEFINE_Y_Y(build_area);
   DEFINE_Y_Y(coverage_union);
+  DEFINE_Y_Y(minimum_bounding_circle);
+  DEFINE_Y_d(minimum_bounding_radius);
 #endif
 
 #if GEOS_SINCE_3_9_0