Reject non-monotonic 1D coords in proximity/allocation/direction

docs/source/reference/proximity.rst

@@ -18,6 +18,13 @@ Proximity
    into a single chunk.  Set a finite ``max_distance`` to keep memory
    bounded.
 
+.. caution::
+
+   ``proximity()``, ``allocation()``, and ``direction()`` require monotonic
+   1D ``x`` and ``y`` coordinates (strictly increasing or strictly
+   decreasing). A non-monotonic axis raises a ``ValueError``; sort the raster
+   along that axis first.
+
 Allocation
 ==========
 .. autosummary::

xrspatial/proximity.py

@@ -254,6 +254,36 @@ def _distance(x1, x2, y1, y2, metric):
     return np.float32(d)
 
 
+def _check_monotonic_coords(x_coords, y_coords, x, y):
+    """Reject non-monotonic 1D coordinates.
+
+    Every backend in this module assumes the 1D axis coordinates are
+    monotonic: ``max_possible_distance`` is taken from the endpoints, the
+    dask halo and the NumPy line-sweep treat array adjacency as spatial
+    adjacency, and the tiled KDTree convergence check lower-bounds the
+    out-of-region distance with chunk-boundary coordinate gaps. None of
+    those hold when a coordinate axis is not monotonic, so a non-monotonic
+    axis silently yields wrong proximity/allocation/direction. Reject it up
+    front with a clear message instead.
+
+    A single-element axis has no order to violate and is allowed.
+    """
+    for coords, name in ((x_coords, x), (y_coords, y)):
+        if len(coords) < 2:
+            continue
+        diffs = np.diff(coords)
+        ascending = np.all(diffs > 0)
+        descending = np.all(diffs < 0)
+        if not (ascending or descending):
+            raise ValueError(
+                "proximity/allocation/direction require strictly monotonic "
+                "(strictly increasing or strictly decreasing, no duplicate or "
+                "NaN values) 1D coordinates, but the {0!r} axis does not "
+                "qualify. Sort the raster along {0!r} before calling.".format(
+                    name)
+            )
+
+
 def _halo_depth(x_coords, y_coords, max_distance, distance_metric):
     """Overlap depth in pixels for the bounded dask map_overlap call.
 
@@ -1226,6 +1256,10 @@ def _process(
     if da is not None and isinstance(y_coords, da.Array):
         y_coords = y_coords.compute()
 
+    # The endpoint-based max distance, the dask halo, the NumPy line-sweep,
+    # and the tiled KDTree convergence check all assume monotonic 1D coords.
+    _check_monotonic_coords(x_coords, y_coords, x, y)
+
     # Compute max_possible_distance using coordinate endpoints directly
     max_possible_distance = _distance(
         x_coords[0], x_coords[-1], y_coords[0], y_coords[-1], distance_metric
@@ -1588,6 +1622,9 @@ def proximity(
         2D array image with `raster.shape` = (height, width).
         If a Dataset is passed, the function is applied to each
         data variable independently, returning a Dataset.
+        The 1D ``x`` and ``y`` coordinates must be monotonic (strictly
+        increasing or strictly decreasing); a non-monotonic axis raises
+        a ValueError.
 
     x : str, default='x'
         Name of x-coordinates.
@@ -1744,6 +1781,9 @@ def allocation(
         2D array of target data.
         If a Dataset is passed, the function is applied to each
         data variable independently, returning a Dataset.
+        The 1D ``x`` and ``y`` coordinates must be monotonic (strictly
+        increasing or strictly decreasing); a non-monotonic axis raises
+        a ValueError.
 
     x : str, default='x'
         Name of x-coordinates.
@@ -1902,6 +1942,9 @@ def direction(
         2D array image with `raster.shape` = (height, width).
         If a Dataset is passed, the function is applied to each
         data variable independently, returning a Dataset.
+        The 1D ``x`` and ``y`` coordinates must be monotonic (strictly
+        increasing or strictly decreasing); a non-monotonic axis raises
+        a ValueError.
 
     x : str, default='x'
         Name of x-coordinates.

xrspatial/tests/test_proximity.py

@@ -1709,6 +1709,103 @@ def test_great_circle_numpy_off_by_more_than_a_metre_is_fixed(
     assert np.nanmax(np.abs(result - expected)) < 1.0
 
 
+# --- non-monotonic 1D coordinate rejection (issue #2851) ---
+
+
+def _nonmonotonic_raster(backend, axis):
+    """Build a raster whose `axis` ('lon' or 'lat') is non-monotonic.
+
+    All other inputs are valid so the only reason _process can fail is the
+    coordinate check.
+    """
+    data = np.asarray([[0., 0., 1., 0.],
+                       [0., 0., 0., 0.],
+                       [2., 0., 0., 0.],
+                       [0., 0., 0., 3.]])
+    lon = np.array([-20., -10., 0., 10.])
+    lat = np.array([20., 10., 0., -10.])
+    if axis == 'lon':
+        lon = np.array([-20., 0., -10., 10.])  # not sorted
+    else:
+        lat = np.array([20., 0., 10., -10.])  # not sorted
+
+    raster = xr.DataArray(data, dims=['lat', 'lon'])
+    raster['lon'] = lon
+    raster['lat'] = lat
+    if has_cuda_and_cupy() and 'cupy' in backend:
+        import cupy
+        raster.data = cupy.asarray(data)
+    if 'dask' in backend and da is not None:
+        raster.data = da.from_array(raster.data, chunks=(2, 2))
+    return raster
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy', 'cupy', 'dask+cupy'])
+@pytest.mark.parametrize("func", [proximity, allocation, direction])
+@pytest.mark.parametrize("axis", ['lon', 'lat'])
+def test_nonmonotonic_coords_raise(backend, func, axis):
+    if 'cupy' in backend and not has_cuda_and_cupy():
+        pytest.skip("cupy not available")
+    if 'dask' in backend and da is None:
+        pytest.skip("dask not available")
+
+    raster = _nonmonotonic_raster(backend, axis)
+    with pytest.raises(ValueError, match="monotonic"):
+        func(raster, x='lon', y='lat')
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy', 'cupy', 'dask+cupy'])
+def test_descending_coords_allowed(backend, result_default_proximity):
+    """A strictly decreasing axis is monotonic and must be accepted.
+
+    The default test raster already uses a descending lat axis; assert the
+    standard backends still produce the reference proximity.
+    """
+    if 'cupy' in backend and not has_cuda_and_cupy():
+        pytest.skip("cupy not available")
+    if 'dask' in backend and da is None:
+        pytest.skip("dask not available")
+
+    data = np.asarray([[0., 0., 0., 0., 0., 2.],
+                       [0., 0., 1., 0., 0., 0.],
+                       [0., np.inf, 3., 0., 0., 0.],
+                       [4., 0., 0., 0., np.nan, 0.]])
+    raster = xr.DataArray(data, dims=['lat', 'lon'])
+    raster['lon'] = np.linspace(-20, 20, 6)   # ascending
+    raster['lat'] = np.linspace(20, -20, 4)   # descending
+    if has_cuda_and_cupy() and 'cupy' in backend:
+        import cupy
+        raster.data = cupy.asarray(data)
+    if 'dask' in backend and da is not None:
+        raster.data = da.from_array(raster.data, chunks=(4, 3))
+
+    result = proximity(raster, x='lon', y='lat')
+    general_output_checks(raster, result, result_default_proximity)
+
+
+@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy', 'cupy', 'dask+cupy'])
+def test_single_element_axis_allowed(backend):
+    """A length-1 axis has no order to violate and must not be rejected."""
+    if 'cupy' in backend and not has_cuda_and_cupy():
+        pytest.skip("cupy not available")
+    if 'dask' in backend and da is None:
+        pytest.skip("dask not available")
+
+    data = np.asarray([[0., 1., 0., 2.]])
+    raster = xr.DataArray(data, dims=['lat', 'lon'])
+    raster['lon'] = np.linspace(0, 30, 4)
+    raster['lat'] = np.array([0.])
+    if has_cuda_and_cupy() and 'cupy' in backend:
+        import cupy
+        raster.data = cupy.asarray(data)
+    if 'dask' in backend and da is not None:
+        raster.data = da.from_array(raster.data, chunks=(1, 2))
+
+    result = proximity(raster, x='lon', y='lat')
+    # no exception; finite proximity at the target columns
+    assert result.shape == (1, 4)
+
+
 # ---------------------------------------------------------------------------
 # Regression: unbounded dask fallback must not mutate the caller's input
 # (issue #2847). _process_dask used to do raster.data = raster.data.rechunk(),