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 fc88bf549..99841eecc 100644 --- a/xrspatial/proximity.py +++ b/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. diff --git a/xrspatial/tests/test_proximity.py b/xrspatial/tests/test_proximity.py index f628a9ec3..b8a20ad74 100644 --- a/xrspatial/tests/test_proximity.py +++ b/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(),