From 48ebb332de5b843509558949421c94a1de78d3f4 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Tue, 2 Jun 2026 14:52:00 -0700 Subject: [PATCH 1/2] Reject non-monotonic 1D coords in proximity/allocation/direction (#2851) The endpoint-based max distance, the dask map_overlap halo, the NumPy line-sweep, and the tiled KDTree convergence check all assume monotonic 1D coordinates. Non-monotonic coords were accepted silently and returned wrong results. Validate both axes once in the shared _process entrypoint and raise a clear ValueError naming the offending axis. A single-element axis is still allowed. Document the requirement in the docstrings and reference page, and add cross-backend tests for the raised error, the descending-axis accept path, and the single-element accept path. --- docs/source/reference/proximity.rst | 7 +++ xrspatial/proximity.py | 42 +++++++++++++ xrspatial/tests/test_proximity.py | 92 +++++++++++++++++++++++++++++ 3 files changed, 141 insertions(+) diff --git a/docs/source/reference/proximity.rst b/docs/source/reference/proximity.rst index e9b4530c2..f1caa1bdc 100644 --- a/docs/source/reference/proximity.rst +++ b/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:: diff --git a/xrspatial/proximity.py b/xrspatial/proximity.py index 35e331d80..01c1c8a73 100644 --- a/xrspatial/proximity.py +++ b/xrspatial/proximity.py @@ -254,6 +254,35 @@ 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 monotonic " + "(strictly increasing or strictly decreasing) 1D " + "coordinates, but the {0!r} axis is not monotonic. 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. @@ -1194,6 +1223,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 @@ -1537,6 +1570,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. @@ -1686,6 +1722,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. @@ -1837,6 +1876,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. diff --git a/xrspatial/tests/test_proximity.py b/xrspatial/tests/test_proximity.py index 726e82f89..40290e8ab 100644 --- a/xrspatial/tests/test_proximity.py +++ b/xrspatial/tests/test_proximity.py @@ -1621,3 +1621,95 @@ def test_great_circle_numpy_off_by_more_than_a_metre_is_fixed( raster, x='lon', y='lat', distance_metric='GREAT_CIRCLE').data 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']) +def test_single_element_axis_allowed(backend): + """A length-1 axis has no order to violate and must not be rejected.""" + 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 '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) From 5c0cd237e3b4c05f7e9d4d3f149e4c2a6253f847 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Tue, 2 Jun 2026 14:53:52 -0700 Subject: [PATCH 2/2] Address review: clarify error message and add cupy single-element case (#2851) - Spell out that duplicate and NaN coordinate values also fail the strict monotonic check, so the message is accurate for those cases. - Parametrize test_single_element_axis_allowed over cupy and dask+cupy to match the other new tests. --- xrspatial/proximity.py | 9 +++++---- xrspatial/tests/test_proximity.py | 7 ++++++- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/xrspatial/proximity.py b/xrspatial/proximity.py index 01c1c8a73..afdc997c6 100644 --- a/xrspatial/proximity.py +++ b/xrspatial/proximity.py @@ -276,10 +276,11 @@ def _check_monotonic_coords(x_coords, y_coords, x, y): descending = np.all(diffs < 0) if not (ascending or descending): raise ValueError( - "proximity/allocation/direction require monotonic " - "(strictly increasing or strictly decreasing) 1D " - "coordinates, but the {0!r} axis is not monotonic. Sort the " - "raster along {0!r} before calling.".format(name) + "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) ) diff --git a/xrspatial/tests/test_proximity.py b/xrspatial/tests/test_proximity.py index 40290e8ab..09224e55a 100644 --- a/xrspatial/tests/test_proximity.py +++ b/xrspatial/tests/test_proximity.py @@ -1697,9 +1697,11 @@ def test_descending_coords_allowed(backend, result_default_proximity): general_output_checks(raster, result, result_default_proximity) -@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy']) +@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") @@ -1707,6 +1709,9 @@ def test_single_element_axis_allowed(backend): 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))