[Done] dwithin for GEOS 3.10.0 (#417)

CHANGELOG.rst

@@ -7,6 +7,7 @@ Version 0.12 (unreleased)
 
 **Major enhancements**
 
+* Added ``pygeos.dwithin`` for GEOS >= 3.10 (#417).
 * Added GeoJSON input/output capabilities (``pygeos.from_geojson``, 
   ``pygeos.to_geojson``) (#413).
 

pygeos/predicates.py

@@ -23,6 +23,7 @@
     "covered_by",
     "covers",
     "disjoint",
+    "dwithin",
     "equals",
     "intersects",
     "overlaps",
@@ -1001,3 +1002,41 @@ def relate_pattern(a, b, pattern, **kwargs):
     True
     """
     return lib.relate_pattern(a, b, pattern, **kwargs)
+
+
+@multithreading_enabled
+@requires_geos("3.10.0")
+def dwithin(a, b, distance, **kwargs):
+    """
+    Returns True if the geometries are within a given distance.
+
+    Using this function is more efficient than computing the distance and
+    comparing the result.
+
+    Parameters
+    ----------
+    a, b : Geometry or array_like
+    distance : float
+        Negative distances always return False.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    distance : compute the actual distance between A and B
+    prepare : improve performance by preparing ``a`` (the first argument)
+
+    Examples
+    --------
+    >>> point = Geometry("POINT (0.5 0.5)")
+    >>> dwithin(point, Geometry("POINT (2 0.5)"), 2)
+    True
+    >>> dwithin(point, Geometry("POINT (2 0.5)"), [2, 1.5, 1]).tolist()
+    [True, True, False]
+    >>> dwithin(point, Geometry("POINT (0.5 0.5)"), 0)
+    True
+    >>> dwithin(point, None, 100)
+    False
+    """
+    return lib.dwithin(a, b, distance, **kwargs)

pygeos/tests/test_predicates.py

@@ -1,3 +1,5 @@
+from functools import partial
+
 import numpy as np
 import pytest
 
@@ -39,15 +41,19 @@
     pygeos.contains,
     pygeos.contains_properly,
     pygeos.overlaps,
-    pygeos.equals,
     pygeos.covers,
     pygeos.covered_by,
+    pytest.param(
+        partial(pygeos.dwithin, distance=1.0),
+        marks=pytest.mark.skipif(
+            pygeos.geos_version < (3, 10, 0), reason="GEOS < 3.10"
+        ),
+    ),
+    pygeos.equals,
     pygeos.equals_exact,
 )
 
-BINARY_PREPARED_PREDICATES = tuple(
-    set(BINARY_PREDICATES) - {pygeos.equals, pygeos.equals_exact}
-)
+BINARY_PREPARED_PREDICATES = BINARY_PREDICATES[:-2]
 
 
 @pytest.mark.parametrize("geometry", all_types)
@@ -111,6 +117,26 @@ def test_equals_exact_tolerance():
     assert pygeos.equals_exact(p1, p1).item() is True
     assert pygeos.equals_exact(p1, p2).item() is False
 
+    # an array of tolerances
+    actual = pygeos.equals_exact(p1, p2, tolerance=[0.05, 0.2, np.nan])
+    np.testing.assert_allclose(actual, [False, True, False])
+
+
+@pytest.mark.skipif(pygeos.geos_version < (3, 10, 0), reason="GEOS < 3.10")
+def test_dwithin():
+    p1 = pygeos.points(50, 4)
+    p2 = pygeos.points(50.1, 4.1)
+    actual = pygeos.dwithin([p1, p2, None], p1, distance=0.05)
+    np.testing.assert_equal(actual, [True, False, False])
+    assert actual.dtype == np.bool_
+    actual = pygeos.dwithin([p1, p2, None], p1, distance=0.2)
+    np.testing.assert_allclose(actual, [True, True, False])
+    assert actual.dtype == np.bool_
+
+    # an array of distances
+    actual = pygeos.dwithin(p1, p2, distance=[0.05, 0.2, np.nan])
+    np.testing.assert_allclose(actual, [False, True, False])
+
 
 @pytest.mark.parametrize(
     "geometry,expected",

src/ufuncs.c

@@ -1560,6 +1560,7 @@ static char equals_exact_dtypes[4] = {NPY_OBJECT, NPY_OBJECT, NPY_DOUBLE, NPY_BO
 static void equals_exact_func(char** args, npy_intp* dimensions, npy_intp* steps,
                               void* data) {
   GEOSGeometry *in1 = NULL, *in2 = NULL;
+  double in3;
   npy_bool ret;
 
   GEOS_INIT_THREADS;
@@ -1574,7 +1575,7 @@ static void equals_exact_func(char** args, npy_intp* dimensions, npy_intp* steps
       errstate = PGERR_NOT_A_GEOMETRY;
       goto finish;
     }
-    double in3 = *(double*)ip3;
+    in3 = *(double*)ip3;
     if ((in1 == NULL) || (in2 == NULL) || npy_isnan(in3)) {
       /* return 0 (False) for missing values */
       ret = 0;
@@ -1593,6 +1594,56 @@ static void equals_exact_func(char** args, npy_intp* dimensions, npy_intp* steps
 }
 static PyUFuncGenericFunction equals_exact_funcs[1] = {&equals_exact_func};
 
+#if GEOS_SINCE_3_10_0
+
+static char dwithin_dtypes[4] = {NPY_OBJECT, NPY_OBJECT, NPY_DOUBLE, NPY_BOOL};
+static void dwithin_func(char** args, npy_intp* dimensions, npy_intp* steps, void* data) {
+  GEOSGeometry *in1 = NULL, *in2 = NULL;
+  GEOSPreparedGeometry* in1_prepared = NULL;
+  double in3;
+  npy_bool ret;
+
+  GEOS_INIT_THREADS;
+
+  TERNARY_LOOP {
+    /* get the geometries: return on error */
+    if (!get_geom_with_prepared(*(GeometryObject**)ip1, &in1, &in1_prepared)) {
+      errstate = PGERR_NOT_A_GEOMETRY;
+      goto finish;
+    }
+    if (!get_geom(*(GeometryObject**)ip2, &in2)) {
+      errstate = PGERR_NOT_A_GEOMETRY;
+      goto finish;
+    }
+    in3 = *(double*)ip3;
+    if ((in1 == NULL) || (in2 == NULL) || npy_isnan(in3)) {
+      /* in case of a missing value: return 0 (False) */
+      ret = 0;
+    } else {
+      if (in1_prepared == NULL) {
+        /* call the GEOS function */
+        ret = GEOSDistanceWithin_r(ctx, in1, in2, in3);
+      } else {
+        /* call the prepared GEOS function */
+        ret = GEOSPreparedDistanceWithin_r(ctx, in1_prepared, in2, in3);
+      }
+      /* return for illegal values */
+      if (ret == 2) {
+        errstate = PGERR_GEOS_EXCEPTION;
+        goto finish;
+      }
+    }
+    *(npy_bool*)op1 = ret;
+  }
+
+finish:
+
+  GEOS_FINISH_THREADS;
+}
+static PyUFuncGenericFunction dwithin_funcs[1] = {&dwithin_func};
+
+#endif  // GEOS_SINCE_3_10_0
+
 static char delaunay_triangles_dtypes[4] = {NPY_OBJECT, NPY_DOUBLE, NPY_BOOL, NPY_OBJECT};
 static void delaunay_triangles_func(char** args, npy_intp* dimensions, npy_intp* steps,
                                     void* data) {
@@ -3385,6 +3436,7 @@ int init_ufuncs(PyObject* m, PyObject* d) {
 
 #if GEOS_SINCE_3_10_0
   DEFINE_Yd_Y(segmentize);
+  DEFINE_CUSTOM(dwithin, 3);
   DEFINE_CUSTOM(from_geojson, 2);
   DEFINE_CUSTOM(to_geojson, 2);
 #endif