Apply cell-size correction to planar aspect (#2760)

xrspatial/aspect.py

@@ -19,7 +19,8 @@
                                 _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, has_dask_array, ngjit)
+                             _validate_raster, cuda_args, get_dataarray_resolution, has_dask_array,
+                             ngjit)
 
 
 def _geodesic_cuda_dims(shape):
@@ -47,7 +48,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 +65,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 +84,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 +107,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,44 +130,57 @@ 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
     if (i-di >= 0 and
         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),
                            meta=np.array((), dtype=np.float32))
     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),
@@ -349,7 +366,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.
@@ -407,13 +425,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:

xrspatial/tests/test_aspect.py

@@ -252,6 +252,63 @@ def test_degenerate_shape_dask_cupy(func, 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)
+
+
 # ---- Independent analytic oracle for planar aspect ----
 #
 # The other planar tests only prove the four backends agree with each
@@ -340,14 +397,8 @@ def test_planar_aspect_oracle_square(backend, gx, gy, cellsize):
     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,
-)
+# Rectangular cells (xres != yres). aspect() now reads cell size from the
+# coordinate spacing (#2760), so the East/North gradient ratio is correct.
 @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)])
@@ -375,10 +426,6 @@ def test_planar_aspect_oracle_square_cupy(backend, gx, gy, cellsize):
 
 
 @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)])