ENH: Add shortest_line (nearest_points) ufunc (#334)

CHANGELOG.rst

@@ -45,6 +45,7 @@ Version 0.10 (unreleased)
 * Addition of a ``segmentize`` function for GEOS >= 3.10 (#299)
 * Addition of ``oriented_envelope`` and ``minimum_rotated_rectangle`` functions (#314)
 * Addition of ``minimum_bounding_circle`` and ``minimum_bounding_radius`` functions for GEOS >= 3.8 (#315)
+* Addition of a ``shortest_line`` ("nearest points") function (#334)
 
 **Bug fixes**
 

pygeos/linear.py

@@ -4,7 +4,13 @@
 from . import lib
 from .decorators import multithreading_enabled
 
-__all__ = ["line_interpolate_point", "line_locate_point", "line_merge", "shared_paths"]
+__all__ = [
+    "line_interpolate_point",
+    "line_locate_point",
+    "line_merge",
+    "shared_paths",
+    "shortest_line",
+]
 
 
 @multithreading_enabled
@@ -140,3 +146,36 @@ def shared_paths(a, b, **kwargs):
     <pygeos.Geometry GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING ...>
     """
     return lib.shared_paths(a, b, **kwargs)
+
+
+@multithreading_enabled
+def shortest_line(a, b, **kwargs):
+    """
+    Returns the shortest line between two geometries.
+
+    The resulting line consists of two points, representing the nearest
+    points between the geometry pair. The line always starts in the first
+    geometry `a` and ends in he second geometry `b`. The endpoints of the
+    line will not necessarily be existing vertices of the input geometries
+    `a` and `b`, but can also be a point along a line segment.
+
+    Parameters
+    ----------
+    a : Geometry or array_like
+    b : 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>`_.
+
+    See also
+    --------
+    prepare : improve performance by preparing ``a`` (the first argument) (for GEOS>=3.9)
+
+    Examples
+    --------
+    >>> geom1 = Geometry("LINESTRING (0 0, 1 0, 1 1, 0 1, 0 0)")
+    >>> geom2 = Geometry("LINESTRING (0 3, 3 0, 5 3)")
+    >>> shortest_line(geom1, geom2)
+    <pygeos.Geometry LINESTRING (1 1, 1.5 1.5)>
+    """
+    return lib.shortest_line(a, b, **kwargs)

pygeos/test/test_linear.py

@@ -154,3 +154,38 @@ def test_shared_paths_non_linestring():
     g2 = pygeos.points(0, 1)
     with pytest.raises(pygeos.GEOSException):
         pygeos.shared_paths(g1, g2)
+
+
+def _prepare_input(geometry, prepare):
+    """Prepare without modifying inplace"""
+    if prepare:
+        geometry = pygeos.apply(geometry, lambda x: x)  # makes a copy
+        pygeos.prepare(geometry)
+        return geometry
+    else:
+        return geometry
+
+
+@pytest.mark.parametrize("prepare", [True, False])
+def test_shortest_line(prepare):
+    g1 = pygeos.linestrings([(0, 0), (1, 0), (1, 1)])
+    g2 = pygeos.linestrings([(0, 3), (3, 0)])
+    actual = pygeos.shortest_line(_prepare_input(g1, prepare), g2)
+    expected = pygeos.linestrings([(1, 1), (1.5, 1.5)])
+    assert pygeos.equals(actual, expected)
+
+
+@pytest.mark.parametrize("prepare", [True, False])
+def test_shortest_line_none(prepare):
+    assert pygeos.shortest_line(_prepare_input(line_string, prepare), None) is None
+    assert pygeos.shortest_line(None, line_string) is None
+    assert pygeos.shortest_line(None, None) is None
+
+
+@pytest.mark.parametrize("prepare", [True, False])
+def test_shortest_line_empty(prepare):
+    g1 = _prepare_input(line_string, prepare)
+    assert pygeos.shortest_line(g1, empty_line_string) is None
+    g1_empty = _prepare_input(empty_line_string, prepare)
+    assert pygeos.shortest_line(g1_empty, line_string) is None
+    assert pygeos.shortest_line(g1_empty, empty_line_string) is None

src/ufuncs.c

@@ -1945,6 +1945,78 @@ static void polygonize_full_func(char** args, npy_intp* dimensions, npy_intp* st
 }
 static PyUFuncGenericFunction polygonize_full_funcs[1] = {&polygonize_full_func};
 
+static char shortest_line_dtypes[3] = {NPY_OBJECT, NPY_OBJECT, NPY_OBJECT};
+static void shortest_line_func(char** args, npy_intp* dimensions, npy_intp* steps,
+                                void* data) {
+  GEOSGeometry* in1 = NULL;
+  GEOSGeometry* in2 = NULL;
+  GEOSPreparedGeometry* in1_prepared = NULL;
+  GEOSGeometry** geom_arr;
+  GEOSCoordSequence* coord_seq = NULL;
+
+  CHECK_NO_INPLACE_OUTPUT(2);
+
+  // allocate a temporary array to store output GEOSGeometry objects
+  geom_arr = malloc(sizeof(void*) * dimensions[0]);
+  CHECK_ALLOC(geom_arr);
+
+  GEOS_INIT_THREADS;
+
+  BINARY_LOOP {
+    /* get the geometries: return on error */
+    if (!get_geom_with_prepared(*(GeometryObject**)ip1, &in1, &in1_prepared)) {
+      errstate = PGERR_NOT_A_GEOMETRY;
+      destroy_geom_arr(ctx, geom_arr, i - 1);
+      break;
+    }
+    if (!get_geom(*(GeometryObject**)ip2, &in2)) {
+      errstate = PGERR_NOT_A_GEOMETRY;
+      destroy_geom_arr(ctx, geom_arr, i - 1);
+      break;
+    }
+
+    if ((in1 == NULL) || (in2 == NULL) || GEOSisEmpty_r(ctx, in1) ||
+        GEOSisEmpty_r(ctx, in2)) {
+      // in case of a missing value or empty geometry: return NULL (None)
+      // GEOSNearestPoints_r returns NULL for empty geometries
+      // but this is not distinguishable from an actual error, so we handle this ourselves
+      geom_arr[i] = NULL;
+      continue;
+    }
+#if GEOS_SINCE_3_9_0
+    if (in1_prepared != NULL) {
+      coord_seq = GEOSPreparedNearestPoints_r(ctx, in1_prepared, in2);
+    } else {
+      coord_seq = GEOSNearestPoints_r(ctx, in1, in2);
+    }
+#else
+    coord_seq = GEOSNearestPoints_r(ctx, in1, in2);
+#endif
+    if (coord_seq == NULL) {
+      errstate = PGERR_GEOS_EXCEPTION;
+      destroy_geom_arr(ctx, geom_arr, i - 1);
+      break;
+    }
+    geom_arr[i] = GEOSGeom_createLineString_r(ctx, coord_seq);
+    // Note: coordinate sequence is owned by linestring; if linestring fails to
+    // construct, it will automatically clean up the coordinate sequence
+    if (geom_arr[i] == NULL) {
+      errstate = PGERR_GEOS_EXCEPTION;
+      destroy_geom_arr(ctx, geom_arr, i - 1);
+      break;
+    }
+  }
+
+  GEOS_FINISH_THREADS;
+
+  // fill the numpy array with PyObjects while holding the GIL
+  if (errstate == PGERR_SUCCESS) {
+    geom_arr_to_npy(geom_arr, args[2], steps[2], dimensions[0]);
+  }
+  free(geom_arr);
+}
+static PyUFuncGenericFunction shortest_line_funcs[1] = {&shortest_line_func};
+
 #if GEOS_SINCE_3_6_0
 static char set_precision_dtypes[4] = {NPY_OBJECT, NPY_DOUBLE, NPY_BOOL, NPY_OBJECT};
 static void set_precision_func(char** args, npy_intp* dimensions, npy_intp* steps,
@@ -3080,6 +3152,7 @@ int init_ufuncs(PyObject* m, PyObject* d) {
   DEFINE_CUSTOM(relate_pattern, 3);
   DEFINE_GENERALIZED(polygonize, 1, "(d)->()");
   DEFINE_GENERALIZED_NOUT4(polygonize_full, 1, "(d)->(),(),(),()");
+  DEFINE_CUSTOM(shortest_line, 2);
 
   DEFINE_GENERALIZED(points, 1, "(d)->()");
   DEFINE_GENERALIZED(linestrings, 1, "(i, d)->()");