diff --git a/xrspatial/slope.py b/xrspatial/slope.py index 23906c4a9..1b229dbf9 100644 --- a/xrspatial/slope.py +++ b/xrspatial/slope.py @@ -360,7 +360,9 @@ def slope(agg: xr.DataArray, size. If the coordinates are in degrees (lat/lon) but the elevation values are in meters, the result is wrong by orders of magnitude. When this mismatch is detected, a ``UserWarning`` is emitted suggesting you reproject - to a projected CRS or use ``method='geodesic'``. + to a projected CRS or use ``method='geodesic'``. The ``'planar'`` method also + raises a ``ValueError`` if the cell size on either axis is zero, negative, or + non-finite, since those values cannot produce a valid slope. References ---------- @@ -401,6 +403,17 @@ def slope(agg: xr.DataArray, if method == 'planar': warn_if_unit_mismatch(agg) cellsize_x, cellsize_y = get_dataarray_resolution(agg) + # Reject negatives too (curvature() only checks == 0): a negative cell + # size would silently flip the slope sign rather than error out. + if (not np.isfinite(cellsize_x) or not np.isfinite(cellsize_y) + or cellsize_x <= 0 or cellsize_y <= 0): + raise ValueError( + "slope(method='planar') requires a positive, finite cell size " + f"on both axes; got cellsize_x={cellsize_x!r}, " + f"cellsize_y={cellsize_y!r}. " + "Set agg.attrs['res'] to a (x, y) tuple of positive floats, " + "or attach numeric x/y coordinates to the DataArray." + ) mapper = ArrayTypeFunctionMapping( numpy_func=_run_numpy, cupy_func=_run_cupy, diff --git a/xrspatial/tests/test_slope.py b/xrspatial/tests/test_slope.py index 7768fb8f6..b705f3f13 100644 --- a/xrspatial/tests/test_slope.py +++ b/xrspatial/tests/test_slope.py @@ -473,6 +473,94 @@ def test_center_nan_propagates_dask_cupy(): verify_attrs=False) +# --------------------------------------------------------------------------- +# Planar cell-size validation (issue #2897) +# --------------------------------------------------------------------------- + +INVALID_RES = [ + (0, 0), # zero on both axes + (0, 1), # zero on x + (1, 0), # zero on y + (-1, 1), # negative on x + (1, -1), # negative on y + (np.inf, 1), # inf on x + (1, np.inf), # inf on y + (np.nan, 1), # nan on x + (1, np.nan), # nan on y +] + + +@pytest.mark.parametrize("res", INVALID_RES) +def test_planar_invalid_cellsize_numpy(res): + """slope(method='planar') must reject zero, negative, or non-finite cell + sizes with a clear ValueError instead of a raw ZeroDivisionError, silent + NaN interior, or plausible-but-wrong values. + """ + data = np.zeros((5, 5), dtype=np.float32) + agg = create_test_raster(data, backend='numpy', attrs={'res': res}) + with pytest.raises(ValueError, match="positive, finite cell size"): + slope(agg) + + +@dask_array_available +@pytest.mark.parametrize("res", INVALID_RES) +def test_planar_invalid_cellsize_dask_numpy(res): + data = np.zeros((5, 5), dtype=np.float32) + agg = create_test_raster(data, backend='dask+numpy', + attrs={'res': res}, chunks=(3, 3)) + with pytest.raises(ValueError, match="positive, finite cell size"): + slope(agg) + + +@cuda_and_cupy_available +@pytest.mark.parametrize("res", INVALID_RES) +def test_planar_invalid_cellsize_cupy(res): + data = np.zeros((5, 5), dtype=np.float32) + agg = create_test_raster(data, backend='cupy', attrs={'res': res}) + with pytest.raises(ValueError, match="positive, finite cell size"): + slope(agg) + + +@dask_array_available +@cuda_and_cupy_available +@pytest.mark.parametrize("res", INVALID_RES) +def test_planar_invalid_cellsize_dask_cupy(res): + data = np.zeros((5, 5), dtype=np.float32) + agg = create_test_raster(data, backend='dask+cupy', + attrs={'res': res}, chunks=(3, 3)) + with pytest.raises(ValueError, match="positive, finite cell size"): + slope(agg) + + +def test_planar_valid_cellsize_still_works(): + """A normal positive cell size must not trip the new guard.""" + data = np.zeros((5, 5), dtype=np.float32) + agg = create_test_raster(data, backend='numpy', attrs={'res': (10, 10)}) + result = slope(agg) + # slope() sets attrs['units'] = 'degrees' (#2894), so the output attrs + # legitimately differ from the input; the units contract is covered + # separately by the units tests below. + general_output_checks(agg, result, verify_attrs=False) + + +def test_geodesic_unaffected_by_invalid_res_attr(): + """The geodesic path takes its scale from lat/lon coords, so a bogus 'res' + attribute must not trigger the planar cell-size guard. + """ + lat = np.linspace(40.0, 40.4, 5) + lon = np.linspace(-105.0, -104.6, 5) + data = np.zeros((5, 5), dtype=np.float64) + agg = xr.DataArray( + data, + coords={'y': lat, 'x': lon}, + dims=['y', 'x'], + attrs={'res': (0, 0)}, + ) + # Must not raise the planar guard; geodesic uses lat/lon coords. + result = slope(agg, method='geodesic') + general_output_checks(agg, result, verify_attrs=False) + + # slope is reported in degrees, so the output must advertise units='degrees' # rather than leaking the input elevation units (e.g. 'm'). Regression for #2894. def test_units_overridden_to_degrees_numpy():