From 94860bc11919378f0ecc03367cfafbee5b712d31 Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Mon, 1 Jun 2026 10:32:14 -0700
Subject: [PATCH] Warn at runtime on degree coords with meter elevation in
 slope planar path (#2764)

The planar slope path uses coordinate spacing directly as the cell size.
With degree (lat/lon) coordinates and meter elevation values, the result
is wrong by orders of magnitude and nothing told the caller. Wire the
existing warn_if_unit_mismatch helper into the planar branch so a
UserWarning fires when the mismatch is detected.

Add tests asserting the warning is emitted for degree coordinates and not
emitted for projected/meter coordinates. Document the runtime warning in
the slope docstring.
---
 xrspatial/slope.py            | 12 +++++++++-
 xrspatial/tests/test_slope.py | 41 +++++++++++++++++++++++++++++++++++
 2 files changed, 52 insertions(+), 1 deletion(-)

diff --git a/xrspatial/slope.py b/xrspatial/slope.py
index 56b41df3..0f3f34fc 100644
--- a/xrspatial/slope.py
+++ b/xrspatial/slope.py
@@ -27,7 +27,8 @@ class cupy(object):
                                 _cpu_geodesic_slope, _run_gpu_geodesic_slope)
 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)
+                             _validate_raster, cuda_args, get_dataarray_resolution, ngjit,
+                             warn_if_unit_mismatch)
 
 
 def _geodesic_cuda_dims(shape):
@@ -339,6 +340,14 @@ def slope(agg: xr.DataArray,
         2D array of slope 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
     ----------
         - arcgis: http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-slope-works.htm # noqa
@@ -376,6 +385,7 @@ def slope(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,
diff --git a/xrspatial/tests/test_slope.py b/xrspatial/tests/test_slope.py
index ed90af00..5e68fc99 100644
--- a/xrspatial/tests/test_slope.py
+++ b/xrspatial/tests/test_slope.py
@@ -1,3 +1,5 @@
+import warnings
+
 import numpy as np
 import pytest
 import xarray as xr
@@ -274,3 +276,42 @@ def test_degenerate_shape_geodesic(shape):
     general_output_checks(raster, result)
     assert result.shape == shape
     assert np.all(np.isnan(result.data))
+
+
+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"):
+        slope(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)
+        slope(raster)  # must not raise: no unit-mismatch warning expected