diff --git a/xrspatial/tests/test_utils.py b/xrspatial/tests/test_utils.py index 6cf0ac65..5ebdcfde 100644 --- a/xrspatial/tests/test_utils.py +++ b/xrspatial/tests/test_utils.py @@ -178,3 +178,96 @@ def test_validate_arrays_mixed_eager_numpy_and_cupy_rejected(): b = xr.DataArray(cupy.zeros((6, 6))) with pytest.raises(ValueError, match="numpy.*cupy"): validate_arrays(a, b) + + +# --------------------------------------------------------------------------- +# calc_res irregular-spacing warning (issue #2766) +# --------------------------------------------------------------------------- + + +def _grid(x, y): + x = np.asarray(x, dtype=float) + y = np.asarray(y, dtype=float) + return xr.DataArray( + np.zeros((y.size, x.size)), + dims=("y", "x"), + coords={"y": y, "x": x}, + ) + + +def test_calc_res_regular_grid_no_warning(): + raster = _grid([0, 1, 2, 3, 4], [0, 1, 2]) + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + xres, yres = utils.calc_res(raster) + assert (xres, yres) == (1.0, 1.0) + assert len(w) == 0 + + +def test_calc_res_irregular_x_warns(): + raster = _grid([0, 1, 2, 4, 8], [0, 1, 2]) + with pytest.warns(UserWarning, match="'x' coordinate is not evenly spaced"): + xres, yres = utils.calc_res(raster) + # averaged span: (8 - 0) / (5 - 1) == 2.0 + assert xres == 2.0 + assert yres == 1.0 + + +def test_calc_res_irregular_y_warns(): + raster = _grid([0, 1, 2], [0, 1, 3, 7]) + with pytest.warns(UserWarning, match="'y' coordinate is not evenly spaced"): + utils.calc_res(raster) + + +def test_calc_res_descending_axis_no_warning(): + # north-up rasters have a descending y axis: regular but negative steps + raster = _grid([0, 1, 2, 3], [3, 2, 1, 0]) + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + xres, yres = utils.calc_res(raster) + assert (xres, yres) == (1.0, 1.0) + assert len(w) == 0 + + +def test_calc_res_no_coords_no_warning(): + raster = xr.DataArray(np.zeros((3, 5))) + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + xres, yres = utils.calc_res(raster) + assert (xres, yres) == (1.0, 1.0) + assert len(w) == 0 + + +def test_calc_res_float_regular_grid_no_warning(): + # floating-point regular spacing must not trip the relative tolerance + x = np.linspace(-74.93, -74.9275, 10) + y = np.linspace(5.0, 5.0025, 10) + raster = _grid(x, y) + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + utils.calc_res(raster) + assert len(w) == 0 + + +def test_get_dataarray_resolution_irregular_with_res_attr_no_warning(): + # attrs['res'] is honored before calc_res, so no averaging warning fires + raster = _grid([0, 1, 2, 4, 8], [0, 1, 2]) + raster.attrs["res"] = (1.0, 1.0) + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + cx, cy = utils.get_dataarray_resolution(raster) + assert (cx, cy) == (1.0, 1.0) + assert len(w) == 0 + + +def test_slope_irregular_coords_warns(): + # the user-facing symptom: planar slope averages cell size silently + from xrspatial import slope + + x = [0, 1, 2, 4, 8] + data = np.tile(np.asarray(x, dtype=float), (4, 1)) + raster = xr.DataArray( + data, dims=("y", "x"), coords={"y": [0, 1, 2, 3], "x": x} + ) + with pytest.warns(UserWarning, match="'x' coordinate is not evenly spaced"): + slope(raster) diff --git a/xrspatial/utils.py b/xrspatial/utils.py index 1c9c7cd8..3da402cf 100644 --- a/xrspatial/utils.py +++ b/xrspatial/utils.py @@ -446,13 +446,65 @@ def calc_res(raster, xdim=None, ydim=None): Tuple of (x-resolution, y-resolution). """ + if ydim is None: + ydim = raster.dims[-2] + if xdim is None: + xdim = raster.dims[-1] + h, w = raster.shape[-2:] xrange, yrange = get_xy_range(raster, xdim, ydim) xres = (xrange[-1] - xrange[0]) / (w - 1) yres = (yrange[-1] - yrange[0]) / (h - 1) + + _warn_if_irregular_spacing(raster, xdim, xres, "x") + _warn_if_irregular_spacing(raster, ydim, yres, "y") + return xres, yres +def _warn_if_irregular_spacing(raster, dim, res, axis_label): + """Warn when a 1-D coordinate on `dim` is not evenly spaced. + + `calc_res` reduces the coordinate to a single average cell size + (full span divided by ``n - 1``). On an irregular grid that average + misrepresents every cell, and the caller gets no signal. Emit a + ``UserWarning`` so the averaging is visible and point at + ``attrs['res']`` for an explicit override (which + ``get_dataarray_resolution`` honors before it reaches `calc_res`). + """ + coord = raster.coords.get(dim, None) + # A 2-point axis has a single step that always equals the averaged + # step, so it cannot be "irregular"; only check axes with >= 3 points. + if coord is None or coord.ndim != 1 or coord.size < 3: + return + + values = np.asarray(coord.values) + if not np.issubdtype(values.dtype, np.number): + return + + diffs = np.diff(values) + if not np.all(np.isfinite(diffs)): + return + + # Compare each step magnitude against the averaged step. `res` comes + # from min/max span so it is non-negative regardless of axis + # direction; a descending (north-up) axis has negative diffs but is + # still regular, so compare absolute values. The relative tolerance + # keeps floating-point jitter from tripping the warning. + if np.allclose(np.abs(diffs), abs(res), rtol=1e-5, atol=0): + return + + warnings.warn( + f"xrspatial: '{dim}' coordinate is not evenly spaced; " + f"using an averaged {axis_label}-resolution of {res}. " + "Per-cell spacing varies, so distance-based results may be " + "inaccurate. Set attrs['res'] to an explicit resolution to " + "silence this warning.", + UserWarning, + stacklevel=3, + ) + + def get_dataarray_resolution( agg: xr.DataArray, xdim: str = None,