added horizontal vertical unit mismatch heuristic to help address issue #840

xrspatial/aspect.py

@@ -14,7 +14,10 @@
 import xarray as xr
 from numba import cuda
 
-from xrspatial.utils import ArrayTypeFunctionMapping, cuda_args, ngjit
+from xrspatial.utils import ArrayTypeFunctionMapping
+from xrspatial.utils import cuda_args
+from xrspatial.utils import ngjit
+from xrspatial.utils import warn_if_unit_mismatch
 
 # 3rd-party
 try:
@@ -262,6 +265,9 @@ def aspect(agg: xr.DataArray,
               dtype=float32)
         Dimensions without coordinates: y, x
     """
+
+    warn_if_unit_mismatch(agg)
+
     mapper = ArrayTypeFunctionMapping(
         numpy_func=_run_numpy,
         dask_func=_run_dask_numpy,

xrspatial/curvature.py

@@ -21,7 +21,11 @@ class cupy(object):
 from numba import cuda
 
 # local modules
-from xrspatial.utils import (ArrayTypeFunctionMapping, cuda_args, get_dataarray_resolution, ngjit)
+from xrspatial.utils import ArrayTypeFunctionMapping
+from xrspatial.utils import cuda_args
+from xrspatial.utils import get_dataarray_resolution
+from xrspatial.utils import ngjit
+from xrspatial.utils import warn_if_unit_mismatch
 
 
 @ngjit
@@ -220,6 +224,7 @@ def curvature(agg: xr.DataArray,
         Attributes:
             res:      (10, 10)
     """
+    warn_if_unit_mismatch(agg)
 
     cellsize_x, cellsize_y = get_dataarray_resolution(agg)
     cellsize = (cellsize_x + cellsize_y) / 2

xrspatial/slope.py

@@ -22,7 +22,11 @@ class cupy(object):
 from numba import cuda
 
 # local modules
-from xrspatial.utils import (ArrayTypeFunctionMapping, cuda_args, get_dataarray_resolution, ngjit)
+from xrspatial.utils import ArrayTypeFunctionMapping
+from xrspatial.utils import cuda_args
+from xrspatial.utils import get_dataarray_resolution
+from xrspatial.utils import ngjit
+from xrspatial.utils import warn_if_unit_mismatch
 
 
 @ngjit
@@ -178,6 +182,9 @@ def slope(agg: xr.DataArray,
         Dimensions without coordinates: dim_0, dim_1
     """
 
+    # warn if we strongly suspect degrees + meters mismatch
+    warn_if_unit_mismatch(agg)
+
     cellsize_x, cellsize_y = get_dataarray_resolution(agg)
     mapper = ArrayTypeFunctionMapping(
         numpy_func=_run_numpy,

xrspatial/tests/general_checks.py

@@ -32,7 +32,7 @@ def create_test_raster(
         backend='numpy',
         name='myraster',
         dims=['y', 'x'],
-        attrs={'res': (0.5, 0.5), 'crs': 'EPSG: 4326'},
+        attrs={'res': (0.5, 0.5), 'crs': 'EPSG: 5070'},
         chunks=(3, 3)
 ):
     raster = xr.DataArray(data, name=name, dims=dims, attrs=attrs)
@@ -42,12 +42,14 @@ def create_test_raster(
     if attrs is not None:
         if 'res' in attrs:
             res = attrs['res']
+
     # set coords for test raster, 2D coords only
     raster[dims[0]] = np.linspace((data.shape[0] - 1) * res[0], 0, data.shape[0])
     raster[dims[1]] = np.linspace(0, (data.shape[1] - 1) * res[1], data.shape[1])
 
-    raster[dims[0]] = np.linspace((data.shape[0] - 1)/2, 0, data.shape[0])
-    raster[dims[1]] = np.linspace(0, (data.shape[1] - 1)/2, data.shape[1])
+    # assign units to coords
+    raster[dims[0]].attrs["units"] = "m"
+    raster[dims[1]].attrs["units"] = "m"
 
     if has_cuda_and_cupy() and 'cupy' in backend:
         import cupy

xrspatial/tests/test_utils.py

@@ -1,5 +1,11 @@
+import numpy as np
+import xarray as xr
+import pytest
+import warnings
+
+
 from xrspatial.datasets import make_terrain
-from xrspatial.utils import canvas_like
+from xrspatial import utils
 from xrspatial.tests.general_checks import dask_array_available
 
 
@@ -8,5 +14,95 @@ def test_canvas_like():
     # aspect ratio is 1:1
     terrain_shape = (1000, 1000)
     terrain = make_terrain(shape=terrain_shape)
-    terrain_res = canvas_like(terrain, width=50)
+    terrain_res = utils.canvas_like(terrain, width=50)
     assert terrain_res.shape == (50, 50)
+
+
+def test_warn_if_unit_mismatch_degrees_horizontal_elevation_vertical(monkeypatch):
+    """
+    If coordinates look like degrees (lon/lat) and values look like elevation
+    (e.g., meters), warn the user about a likely unit mismatch.
+    """
+    data = np.linspace(0, 999, 10 * 10, dtype=float).reshape(10, 10)
+
+    # Coordinates in degrees (lon/lat-ish)
+    y = np.linspace(5.0, 5.0025, 10)
+    x = np.linspace(-74.93, -74.9275, 10)
+
+    da = xr.DataArray(
+        data,
+        dims=("y", "x"),
+        coords={"y": y, "x": x},
+        attrs={"units": "m"},  # elevation in meters
+    )
+
+    def fake_get_dataarray_resolution(arr):
+        return float(x[1] - x[0]), float(y[1] - y[0])
+
+    monkeypatch.setattr(utils, "get_dataarray_resolution", fake_get_dataarray_resolution)
+
+    # Here we *do* expect a warning
+    with pytest.warns(UserWarning, match="appears to have coordinates in degrees"):
+        utils.warn_if_unit_mismatch(da)
+
+
+def test_warn_if_unit_mismatch_no_warning_for_projected_like_grid(monkeypatch):
+    """
+    If coordinates look like projected linear units (e.g., meters) and values
+    look like elevation, we should NOT warn.
+    """
+    data = np.linspace(0, 999, 10 * 10, dtype=float).reshape(10, 10)
+
+    # Coordinates in meters (projected-looking)
+    y = np.arange(10) * 30.0              # 0, 30, 60, ...
+    x = 500_000.0 + np.arange(10) * 30.0  # UTM-ish eastings
+
+    da = xr.DataArray(
+        data,
+        dims=("y", "x"),
+        coords={"y": y, "x": x},
+        attrs={"units": "m"},  # elevation in meters
+    )
+
+    def fake_get_dataarray_resolution(arr):
+        return float(x[1] - x[0]), float(y[1] - y[0])  # 30 m
+
+    monkeypatch.setattr(utils, "get_dataarray_resolution", fake_get_dataarray_resolution)
+
+    # Capture warnings using the stdlib warnings module
+    with warnings.catch_warnings(record=True) as w:
+        warnings.simplefilter("always")
+        utils.warn_if_unit_mismatch(da)
+
+    assert len(w) == 0, "Expected no warnings for projected-like coordinates"
+
+
+def test_warn_if_unit_mismatch_degrees_but_angle_vertical(monkeypatch):
+    """
+    If coordinates are in degrees but the DataArray itself looks like an angle
+    (e.g., units='degrees'), we should NOT warn; slope/aspect outputs fall into
+    this category.
+    """
+    data = np.linspace(0, 90, 10 * 10, dtype=float).reshape(10, 10)
+
+    # Coordinates in degrees again
+    y = np.linspace(5.0, 5.0025, 10)
+    x = np.linspace(-74.93, -74.9275, 10)
+
+    da = xr.DataArray(
+        data,
+        dims=("y", "x"),
+        coords={"y": y, "x": x},
+        attrs={"units": "degrees"},  # angle, not elevation
+    )
+
+    def fake_get_dataarray_resolution(arr):
+        return float(x[1] - x[0]), float(y[1] - y[0])
+
+    monkeypatch.setattr(utils, "get_dataarray_resolution", fake_get_dataarray_resolution)
+
+    with warnings.catch_warnings(record=True) as w:
+        warnings.simplefilter("always")
+        utils.warn_if_unit_mismatch(da)
+
+    assert len(w) == 0, "Expected no warnings when vertical units are angles"