Warn on degree/meter unit mismatch in planar aspect()

xrspatial/aspect.py

@@ -21,7 +21,7 @@
 from xrspatial.utils import (Z_UNITS, ArrayTypeFunctionMapping, _boundary_to_dask,
                              _extract_latlon_coords, _pad_array, _validate_boundary,
                              _validate_raster, cuda_args, get_dataarray_resolution,
-                             ngjit)
+                             ngjit, warn_if_unit_mismatch)
 
 
 def _geodesic_cuda_dims(shape):
@@ -400,6 +400,14 @@ def aspect(agg: xr.DataArray,
         2D aggregate array of calculated aspect values.
         All other input attributes are preserved.
 
+    Notes
+    -----
+    The ``'planar'`` method uses the coordinate spacing directly as the cell
+    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'``.
+
     References
     ----------
 
@@ -434,6 +442,7 @@ def aspect(agg: xr.DataArray,
     _validate_boundary(boundary)
 
     if method == 'planar':
+        warn_if_unit_mismatch(agg)
         cellsize_x, cellsize_y = get_dataarray_resolution(agg)
         mapper = ArrayTypeFunctionMapping(
             numpy_func=_run_numpy,

xrspatial/tests/test_aspect.py

@@ -1,3 +1,5 @@
+import warnings
+
 import numpy as np
 import pytest
 import xarray as xr
@@ -477,6 +479,45 @@ def test_near_360_aspect_cupy_matches_numpy(backend):
         gpu_interior, np_interior, rtol=1e-6, atol=1e-3)
 
 
+def _degree_coord_meter_elevation_raster():
+    # degree-like coords (lon/lat) with meter-scale elevation values
+    data = np.linspace(0, 999, 10 * 10, dtype=float).reshape(10, 10)
+    y = np.linspace(5.0, 5.0025, 10)
+    x = np.linspace(-74.93, -74.9275, 10)
+    return xr.DataArray(
+        data,
+        dims=('y', 'x'),
+        coords={'y': y, 'x': x},
+        attrs={'units': 'm'},
+    )
+
+
+def _projected_meter_raster():
+    # projected coords in meters (UTM-ish) with meter-scale elevation values
+    data = np.linspace(0, 999, 10 * 10, dtype=float).reshape(10, 10)
+    y = np.arange(10) * 30.0
+    x = 500_000.0 + np.arange(10) * 30.0
+    return xr.DataArray(
+        data,
+        dims=('y', 'x'),
+        coords={'y': y, 'x': x},
+        attrs={'units': 'm'},
+    )
+
+
+def test_planar_warns_on_degree_coords_meter_elevation():
+    raster = _degree_coord_meter_elevation_raster()
+    with pytest.warns(UserWarning, match="appears to have coordinates in degrees"):
+        aspect(raster)  # method='planar' is the default
+
+
+def test_planar_no_warn_on_projected_meter_coords():
+    raster = _projected_meter_raster()
+    with warnings.catch_warnings():
+        warnings.simplefilter("error", UserWarning)
+        aspect(raster)  # must not raise: no unit-mismatch warning expected
+
+
 # The output .name must agree across backends. xr.DataArray keeps a dask
 # array's internal graph token (e.g. '_trim-<hash>' planar, 'getitem-<hash>'
 # geodesic) as .name when name=None, so the dask backends used to disagree