From 0fb342b45bdc2db9c60280cd0afd9b1c6d2452f4 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 1 Jun 2026 10:33:21 -0700 Subject: [PATCH 1/2] Apply cell-size correction to planar aspect kernels (#2760) Planar aspect divided both Horn gradients by 8 without accounting for cellsize_x / cellsize_y, so it returned wrong directions for non-square cells. slope already did this correction; aspect did not. Thread cellsize_x / cellsize_y through the planar path the same way slope does: extract resolution in the public aspect(), pass it into the numpy/cupy/dask wrappers, and divide dz_dx by 8*cellsize_x and dz_dy by 8*cellsize_y in both the CPU and GPU kernels. Square cells are unchanged because the factor cancels in the atan2 ratio. northness() and eastness() are fixed for free since they derive from aspect(). Add a regression test: a constant-gradient raster with xres=10, yres=1 now yields aspect 315 (was 275.7106), with matching northness/eastness. --- xrspatial/aspect.py | 57 ++++++++++++++++++++++------------ xrspatial/tests/test_aspect.py | 57 ++++++++++++++++++++++++++++++++++ 2 files changed, 94 insertions(+), 20 deletions(-) diff --git a/xrspatial/aspect.py b/xrspatial/aspect.py index 3025c0e7..f618f189 100644 --- a/xrspatial/aspect.py +++ b/xrspatial/aspect.py @@ -19,7 +19,7 @@ _cpu_geodesic_aspect, _run_gpu_geodesic_aspect) from xrspatial.utils import (Z_UNITS, ArrayTypeFunctionMapping, _boundary_to_dask, _extract_latlon_coords, _pad_array, _validate_boundary, - _validate_raster, cuda_args, ngjit) + _validate_raster, cuda_args, get_dataarray_resolution, ngjit) def _geodesic_cuda_dims(shape): @@ -47,7 +47,7 @@ class cupy(object): # ===================================================================== @ngjit -def _cpu(data: np.ndarray): +def _cpu(data: np.ndarray, cellsize_x, cellsize_y): data = data.astype(np.float32) out = np.empty(data.shape, dtype=np.float32) out[:] = np.nan @@ -64,8 +64,8 @@ def _cpu(data: np.ndarray): h = data[y+1, x] i = data[y+1, x+1] - dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) / 8 - dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / 8 + dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * cellsize_x) + dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * cellsize_y) if dz_dx == 0 and dz_dy == 0: # flat surface, slope = 0, thus invalid aspect @@ -83,16 +83,19 @@ def _cpu(data: np.ndarray): return out -def _run_numpy(data: np.ndarray, boundary: str = 'nan') -> np.ndarray: +def _run_numpy(data: np.ndarray, + cellsize_x, + cellsize_y, + boundary: str = 'nan') -> np.ndarray: if boundary == 'nan': - return _cpu(data) + return _cpu(data, cellsize_x, cellsize_y) padded = _pad_array(data, 1, boundary) - result = _cpu(padded) + result = _cpu(padded, cellsize_x, cellsize_y) return result[1:-1, 1:-1] @cuda.jit(device=True) -def _gpu(arr): +def _gpu(arr, cellsize_x, cellsize_y): a = arr[0, 0] b = arr[0, 1] @@ -103,8 +106,8 @@ def _gpu(arr): h = arr[2, 1] i = arr[2, 2] - dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) / 8 - dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / 8 + dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * cellsize_x[0]) + dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * cellsize_y[0]) if dz_dx == 0 and dz_dy == 0: # flat surface, slope = 0, thus invalid aspect @@ -126,7 +129,7 @@ def _gpu(arr): @cuda.jit -def _run_gpu(arr, out): +def _run_gpu(arr, cellsize_x_arr, cellsize_y_arr, out): i, j = cuda.grid(2) di = 1 dj = 1 @@ -134,26 +137,36 @@ def _run_gpu(arr, out): i+di < out.shape[0] and j-dj >= 0 and j+dj < out.shape[1]): - out[i, j] = _gpu(arr[i-di:i+di+1, j-dj:j+dj+1]) + out[i, j] = _gpu(arr[i-di:i+di+1, j-dj:j+dj+1], + cellsize_x_arr, + cellsize_y_arr) -def _run_cupy(data: cupy.ndarray, boundary: str = 'nan') -> cupy.ndarray: +def _run_cupy(data: cupy.ndarray, + cellsize_x, + cellsize_y, + boundary: str = 'nan') -> cupy.ndarray: if boundary != 'nan': padded = _pad_array(data, 1, boundary) - result = _run_cupy(padded) + result = _run_cupy(padded, cellsize_x, cellsize_y) return result[1:-1, 1:-1] + cellsize_x_arr = cupy.array([float(cellsize_x)], dtype='f4') + cellsize_y_arr = cupy.array([float(cellsize_y)], dtype='f4') data = data.astype(cupy.float32) griddim, blockdim = cuda_args(data.shape) out = cupy.empty(data.shape, dtype='f4') out[:] = cupy.nan - _run_gpu[griddim, blockdim](data, out) + _run_gpu[griddim, blockdim](data, cellsize_x_arr, cellsize_y_arr, out) return out -def _run_dask_numpy(data: da.Array, boundary: str = 'nan') -> da.Array: +def _run_dask_numpy(data: da.Array, + cellsize_x, + cellsize_y, + boundary: str = 'nan') -> da.Array: data = data.astype(np.float32) - _func = partial(_cpu) + _func = partial(_cpu, cellsize_x=cellsize_x, cellsize_y=cellsize_y) out = data.map_overlap(_func, depth=(1, 1), boundary=_boundary_to_dask(boundary), @@ -161,9 +174,12 @@ def _run_dask_numpy(data: da.Array, boundary: str = 'nan') -> da.Array: return out -def _run_dask_cupy(data: da.Array, boundary: str = 'nan') -> da.Array: +def _run_dask_cupy(data: da.Array, + cellsize_x, + cellsize_y, + boundary: str = 'nan') -> da.Array: data = data.astype(cupy.float32) - _func = partial(_run_cupy) + _func = partial(_run_cupy, cellsize_x=cellsize_x, cellsize_y=cellsize_y) out = data.map_overlap(_func, depth=(1, 1), boundary=_boundary_to_dask(boundary, is_cupy=True), @@ -407,13 +423,14 @@ def aspect(agg: xr.DataArray, _validate_boundary(boundary) if method == 'planar': + cellsize_x, cellsize_y = get_dataarray_resolution(agg) mapper = ArrayTypeFunctionMapping( numpy_func=_run_numpy, dask_func=_run_dask_numpy, cupy_func=_run_cupy, dask_cupy_func=_run_dask_cupy, ) - out = mapper(agg)(agg.data, boundary=boundary) + out = mapper(agg)(agg.data, cellsize_x, cellsize_y, boundary=boundary) else: # geodesic if z_unit not in Z_UNITS: diff --git a/xrspatial/tests/test_aspect.py b/xrspatial/tests/test_aspect.py index 32d0f66a..1ba90f9a 100644 --- a/xrspatial/tests/test_aspect.py +++ b/xrspatial/tests/test_aspect.py @@ -249,3 +249,60 @@ def test_degenerate_shape_dask_cupy(func, shape): result = func(agg) assert result.shape == shape assert np.all(np.isnan(_to_numpy(result))) + + +# ---- Rectangular-cell correctness (issue #2760) ---- +# +# Planar aspect must honour cellsize_x / cellsize_y. For a raster whose +# elevation increases by one unit per map-unit in both x and y, the +# East and North gradients are equal in map units, so the downslope +# points southwest and aspect is 315 degrees regardless of cell aspect +# ratio. With the divide-by-8-only kernels (no cell-size correction) +# this case returned 275.7106. See issue #2760. + +def _rectangular_cell_raster(backend='numpy', xres=10.0, yres=1.0): + ny, nx = 6, 6 + yc = np.arange(ny) * yres + xc = np.arange(nx) * xres + X, Y = np.meshgrid(xc, yc) + data = (X + Y).astype(np.float64) + return create_test_raster( + data, backend=backend, + attrs={'res': (xres, yres), 'crs': 'EPSG: 5070'}) + + +@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy']) +def test_rectangular_cell_aspect(backend): + if backend.startswith('dask') and not dask_array_available: + pytest.skip("dask not available") + agg = _rectangular_cell_raster(backend=backend) + result = aspect(agg) + interior = _to_numpy(result)[1:-1, 1:-1] + np.testing.assert_allclose(interior, 315.0, rtol=1e-4) + + +@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy']) +def test_rectangular_cell_northness_eastness(backend): + if backend.startswith('dask') and not dask_array_available: + pytest.skip("dask not available") + agg = _rectangular_cell_raster(backend=backend) + north = _to_numpy(northness(agg))[1:-1, 1:-1] + east = _to_numpy(eastness(agg))[1:-1, 1:-1] + # aspect == 315: cos(315) = +sqrt(2)/2, sin(315) = -sqrt(2)/2 + np.testing.assert_allclose(north, np.cos(np.deg2rad(315.0)), atol=1e-5) + np.testing.assert_allclose(east, np.sin(np.deg2rad(315.0)), atol=1e-5) + + +@cuda_and_cupy_available +def test_rectangular_cell_aspect_cupy(): + agg = _rectangular_cell_raster(backend='cupy') + interior = _to_numpy(aspect(agg))[1:-1, 1:-1] + np.testing.assert_allclose(interior, 315.0, rtol=1e-4) + + +@dask_array_available +@cuda_and_cupy_available +def test_rectangular_cell_aspect_dask_cupy(): + agg = _rectangular_cell_raster(backend='dask+cupy') + interior = _to_numpy(aspect(agg))[1:-1, 1:-1] + np.testing.assert_allclose(interior, 315.0, rtol=1e-4) From 8ae02af9ee79dfafe44bac4e789d4a13efa9fa7d Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 1 Jun 2026 10:35:14 -0700 Subject: [PATCH 2/2] Update aspect docstring for cell-size handling (#2760) --- xrspatial/aspect.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/xrspatial/aspect.py b/xrspatial/aspect.py index f618f189..aff6801c 100644 --- a/xrspatial/aspect.py +++ b/xrspatial/aspect.py @@ -365,7 +365,8 @@ def aspect(agg: xr.DataArray, name : str, default='aspect' Name of ouput DataArray. method : str, default='planar' - ``'planar'`` uses the classic Horn algorithm with uniform cell size. + ``'planar'`` uses the classic Horn algorithm, scaling the gradients + by the x and y cell sizes so non-square cells are handled correctly. ``'geodesic'`` converts cells to Earth-Centered Earth-Fixed (ECEF) coordinates and fits a 3D plane, yielding accurate results for geographic (lat/lon) coordinate systems.