From bd8f1cdad5f8e4713cae0d6527e80b285537503d Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 1 Jun 2026 10:33:44 -0700 Subject: [PATCH 1/2] aspect: add independent rectangular-cell oracle for planar aspect (#2768) The planar tests only check that the four backends agree. They share the same _cpu/_gpu kernels, so a shared defect like the missing-cellsize bug (#2760) passes anyway. The QGIS fixture is a real reference but uses square cells, so it never exercises rectangular cells or checks that cell size is read from the coordinates. Add an analytic oracle: build a plane z = gx*X + gy*Y on a coordinate grid with no 'res' attr (so resolution comes from the x/y spacing) and compare aspect() against the closed-form Horn compass value instead of against another backend. - Square-cell oracle passes unconditionally on numpy, dask, cupy, dask+cupy. - Rectangular-cell oracle is xfail(strict=False) pending #2760; it flips to XPASS automatically once that fix lands, flagging the marker for removal. --- xrspatial/tests/test_aspect.py | 135 +++++++++++++++++++++++++++++++++ 1 file changed, 135 insertions(+) diff --git a/xrspatial/tests/test_aspect.py b/xrspatial/tests/test_aspect.py index 32d0f66a..5493f6df 100644 --- a/xrspatial/tests/test_aspect.py +++ b/xrspatial/tests/test_aspect.py @@ -1,5 +1,6 @@ import numpy as np import pytest +import xarray as xr try: import dask.array as da @@ -249,3 +250,137 @@ def test_degenerate_shape_dask_cupy(func, shape): result = func(agg) assert result.shape == shape assert np.all(np.isnan(_to_numpy(result))) + + +# ---- Independent analytic oracle for planar aspect ---- +# +# The other planar tests only prove the four backends agree with each +# other. They all share the same _cpu/_gpu kernels, so a defect common to +# every backend (like the missing-cellsize bug in #2760) passes anyway. +# The QGIS fixture is a real external reference but uses square cells, so +# it never checks rectangular cells or that cell size is read from the +# coordinates. +# +# These tests build a surface with a known constant gradient and compare +# aspect() against the closed-form Horn value computed from the real cell +# sizes, not against another backend. The raster carries no 'res' attr, so +# resolution must come from the x/y coordinate spacing. + +_RADIAN = 180 / np.pi + + +def _aspect_oracle(gx, gy): + """Closed-form compass aspect for a plane z = gx*X + gy*Y (X east, Y + south-by-row). The cell-size-correct Horn gradients of such a plane are + exactly gx and gy, independent of cell size, so the oracle only needs + the gradients.""" + if gx == 0 and gy == 0: + return -1.0 + a = np.arctan2(gy, -gx) * _RADIAN + if a < 0: + return 90.0 - a + elif a > 90.0: + return 360.0 - a + 90.0 + else: + return 90.0 - a + + +def _gradient_plane(gx, gy, cellsize_x, cellsize_y, n=7): + """A DataArray of z = gx*X + gy*Y on a grid with the given cell sizes. + + Coordinates encode the resolution; no 'res' attr is set so cell size + must be recovered from the coordinate spacing. + """ + xs = np.arange(n) * cellsize_x + ys = np.arange(n) * cellsize_y + z = gx * xs[None, :] + gy * ys[:, None] + return xr.DataArray( + z.astype(np.float64), dims=['y', 'x'], coords={'y': ys, 'x': xs}, + name='gradient_plane') + + +def _to_backend(agg, backend, chunks=(3, 3)): + data = agg.data + if 'cupy' in backend: + import cupy + data = cupy.asarray(data) + if 'dask' in backend: + data = da.from_array(data, chunks=chunks) + return xr.DataArray(data, dims=agg.dims, coords=agg.coords, + attrs=agg.attrs, name=agg.name) + + +_ORACLE_BACKENDS = ['numpy'] +if has_dask_array(): + _ORACLE_BACKENDS.append('dask+numpy') + + +def _aspect_interior(agg): + """Compute aspect and return the realized numpy interior (drop the + NaN edge ring).""" + result = aspect(agg, method='planar') + return _to_numpy(result)[1:-1, 1:-1] + + +# Square cells: aspect is correct today, so this oracle must pass +# unconditionally on every backend. +@pytest.mark.parametrize("backend", _ORACLE_BACKENDS) +@pytest.mark.parametrize("gx,gy", [(1.0, 1.0), (2.0, -3.0), (-1.0, 0.0)]) +@pytest.mark.parametrize("cellsize", [1.0, 0.5, 10.0]) +def test_planar_aspect_oracle_square(backend, gx, gy, cellsize): + agg = _to_backend(_gradient_plane(gx, gy, cellsize, cellsize), backend) + interior = _aspect_interior(agg) + expected = _aspect_oracle(gx, gy) + np.testing.assert_allclose(interior, expected, rtol=1e-4, atol=1e-3) + + +# Rectangular cells (xres != yres). aspect() ignores cell size today +# (#2760), so the East/North gradient ratio is wrong and these fail until +# that fix lands. xfail(strict=False) flips to XPASS automatically once +# #2760 merges, signalling the marker can be removed. +@pytest.mark.xfail( + reason="planar aspect ignores cellsize_x/cellsize_y; pending #2760", + strict=False, +) +@pytest.mark.parametrize("backend", _ORACLE_BACKENDS) +@pytest.mark.parametrize("gx,gy", [(1.0, 1.0), (2.0, -3.0)]) +@pytest.mark.parametrize("cellsize_x,cellsize_y", [(10.0, 1.0), (1.0, 4.0)]) +def test_planar_aspect_oracle_rectangular( + backend, gx, gy, cellsize_x, cellsize_y): + agg = _to_backend(_gradient_plane(gx, gy, cellsize_x, cellsize_y), backend) + interior = _aspect_interior(agg) + expected = _aspect_oracle(gx, gy) + np.testing.assert_allclose(interior, expected, rtol=1e-4, atol=1e-3) + + +# Cupy and dask+cupy oracle (separate functions so the cuda/cupy skip +# marker applies). Same assertion target as the numpy/dask versions. +@cuda_and_cupy_available +@pytest.mark.parametrize("backend", ['cupy', 'dask+cupy']) +@pytest.mark.parametrize("gx,gy", [(1.0, 1.0), (2.0, -3.0), (-1.0, 0.0)]) +@pytest.mark.parametrize("cellsize", [1.0, 0.5, 10.0]) +def test_planar_aspect_oracle_square_cupy(backend, gx, gy, cellsize): + if backend == 'dask+cupy' and not has_dask_array(): + pytest.skip("Requires dask.Array") + agg = _to_backend(_gradient_plane(gx, gy, cellsize, cellsize), backend) + interior = _aspect_interior(agg) + expected = _aspect_oracle(gx, gy) + np.testing.assert_allclose(interior, expected, rtol=1e-4, atol=1e-3) + + +@cuda_and_cupy_available +@pytest.mark.xfail( + reason="planar aspect ignores cellsize_x/cellsize_y; pending #2760", + strict=False, +) +@pytest.mark.parametrize("backend", ['cupy', 'dask+cupy']) +@pytest.mark.parametrize("gx,gy", [(1.0, 1.0), (2.0, -3.0)]) +@pytest.mark.parametrize("cellsize_x,cellsize_y", [(10.0, 1.0), (1.0, 4.0)]) +def test_planar_aspect_oracle_rectangular_cupy( + backend, gx, gy, cellsize_x, cellsize_y): + if backend == 'dask+cupy' and not has_dask_array(): + pytest.skip("Requires dask.Array") + agg = _to_backend(_gradient_plane(gx, gy, cellsize_x, cellsize_y), backend) + interior = _aspect_interior(agg) + expected = _aspect_oracle(gx, gy) + np.testing.assert_allclose(interior, expected, rtol=1e-4, atol=1e-3) From 12be9683daf1ec4e2aa0241b24fb368707d80f0a Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 1 Jun 2026 10:35:49 -0700 Subject: [PATCH 2/2] Address review: exercise flat-surface oracle case, note duplication (#2768) - Add (0.0, 0.0) to the square oracle params so the flat-surface -1 branch in _aspect_oracle is actually tested (it was dead code before). - Note in the docstring that the arctan2/compass conversion is duplicated from aspect._cpu on purpose, with the square tests guarding against drift. --- xrspatial/tests/test_aspect.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/xrspatial/tests/test_aspect.py b/xrspatial/tests/test_aspect.py index 5493f6df..eadf6c99 100644 --- a/xrspatial/tests/test_aspect.py +++ b/xrspatial/tests/test_aspect.py @@ -273,7 +273,13 @@ def _aspect_oracle(gx, gy): """Closed-form compass aspect for a plane z = gx*X + gy*Y (X east, Y south-by-row). The cell-size-correct Horn gradients of such a plane are exactly gx and gy, independent of cell size, so the oracle only needs - the gradients.""" + the gradients. + + The arctan2/compass conversion below is duplicated from aspect._cpu on + purpose: an oracle that called aspect() would not be independent. The + square-cell tests re-check this against the real function every run, so + drift is caught. + """ if gx == 0 and gy == 0: return -1.0 a = np.arctan2(gy, -gx) * _RADIAN @@ -325,7 +331,7 @@ def _aspect_interior(agg): # Square cells: aspect is correct today, so this oracle must pass # unconditionally on every backend. @pytest.mark.parametrize("backend", _ORACLE_BACKENDS) -@pytest.mark.parametrize("gx,gy", [(1.0, 1.0), (2.0, -3.0), (-1.0, 0.0)]) +@pytest.mark.parametrize("gx,gy", [(1.0, 1.0), (2.0, -3.0), (-1.0, 0.0), (0.0, 0.0)]) @pytest.mark.parametrize("cellsize", [1.0, 0.5, 10.0]) def test_planar_aspect_oracle_square(backend, gx, gy, cellsize): agg = _to_backend(_gradient_plane(gx, gy, cellsize, cellsize), backend)