Set slope output units to 'degrees' instead of leaking input units (#2894)

slope() returned degrees but copied the input DataArray's attrs verbatim,
so a DEM with units='m' produced a slope array still labelled 'm'. That
stale label made _infer_vertical_unit_type misclassify slope output as
elevation. Override attrs['units'] to 'degrees' in the shared public
wrapper (covers numpy, cupy, dask+numpy, dask+cupy and Dataset inputs),
preserving all other attrs.

Add a verify_attrs passthrough to the cross-backend equality helpers and
use it in the slope tests whose attrs now legitimately differ. Add tests
asserting the degrees override on every backend and the geodesic path.

Author: brendancol

xrspatial/slope.py

@@ -351,8 +351,8 @@ def slope(agg: xr.DataArray,
         If `agg` is a DataArray, returns a DataArray of the same type.
         If `agg` is a Dataset, returns a Dataset with slope computed
         for each data variable.
-        2D array of slope values.
-        All other input attributes are preserved.
+        2D array of slope values in degrees. The output ``attrs['units']``
+        is set to ``'degrees'``; all other input attributes are preserved.
 
     Notes
     -----
@@ -429,11 +429,17 @@ def slope(agg: xr.DataArray,
         )
         out = mapper(agg)(agg.data, lat_2d, lon_2d, WGS84_A2, WGS84_B2, z_factor, boundary)
 
+    # slope is reported in degrees, so override the input elevation 'units'
+    # (e.g. 'm') with the angle unit the result actually carries. Other input
+    # attrs are preserved.
+    attrs = dict(agg.attrs)
+    attrs['units'] = 'degrees'
+
     result = xr.DataArray(out,
                           name=name,
                           coords=agg.coords,
                           dims=agg.dims,
-                          attrs=agg.attrs)
+                          attrs=attrs)
     # On dask backends, xr.DataArray keeps the dask array's internal graph
     # token as .name when name=None, so reset it post-construction to match
     # the numpy/cupy backends. (Same fix as zonal #2611, focal #2733.)

xrspatial/tests/general_checks.py

@@ -174,34 +174,36 @@ def assert_boundary_mode_correctness(numpy_agg, dask_agg, func, depth=1, rtol=1e
             )
 
 
-def assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, func, nan_edges=True):
+def assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, func, nan_edges=True,
+                                   verify_attrs=True):
     numpy_result = func(numpy_agg)
     if nan_edges:
         assert_nan_edges_effect(numpy_result)
 
     dask_result = func(dask_agg)
-    general_output_checks(dask_agg, dask_result)
+    general_output_checks(dask_agg, dask_result, verify_attrs=verify_attrs)
     np.testing.assert_allclose(numpy_result.data, dask_result.data.compute(), equal_nan=True)
 
 
-def assert_numpy_equals_cupy(numpy_agg, cupy_agg, func, nan_edges=True, atol=0, rtol=1e-7):
+def assert_numpy_equals_cupy(numpy_agg, cupy_agg, func, nan_edges=True, atol=0, rtol=1e-7,
+                             verify_attrs=True):
     numpy_result = func(numpy_agg)
     if nan_edges:
         assert_nan_edges_effect(numpy_result)
 
     cupy_result = func(cupy_agg)
-    general_output_checks(cupy_agg, cupy_result)
+    general_output_checks(cupy_agg, cupy_result, verify_attrs=verify_attrs)
     np.testing.assert_allclose(
         numpy_result.data, cupy_result.data.get(), equal_nan=True, atol=atol, rtol=rtol)
 
 
 def assert_numpy_equals_dask_cupy(numpy_agg, dask_cupy_agg, func,
-                                  nan_edges=True, atol=0, rtol=1e-7):
+                                  nan_edges=True, atol=0, rtol=1e-7, verify_attrs=True):
     numpy_result = func(numpy_agg)
     if nan_edges:
         assert_nan_edges_effect(numpy_result)
 
     dask_cupy_result = func(dask_cupy_agg)
-    general_output_checks(dask_cupy_agg, dask_cupy_result)
+    general_output_checks(dask_cupy_agg, dask_cupy_result, verify_attrs=verify_attrs)
     np.testing.assert_allclose(numpy_result.data, dask_cupy_result.data.compute().get(),
                                equal_nan=True, atol=atol, rtol=rtol)

xrspatial/tests/test_slope.py

@@ -42,7 +42,8 @@ def test_numpy_equals_qgis(elevation_raster, qgis_slope):
     # slope by xrspatial
     numpy_agg = input_data(elevation_raster, backend='numpy')
     xrspatial_slope_numpy = slope(numpy_agg, name='slope_numpy')
-    general_output_checks(numpy_agg, xrspatial_slope_numpy)
+    # slope overrides attrs['units'] to 'degrees', so skip exact-attrs check
+    general_output_checks(numpy_agg, xrspatial_slope_numpy, verify_attrs=False)
     assert xrspatial_slope_numpy.name == 'slope_numpy'
     print('numpy_agg', numpy_agg)
     print('xrspatial_slope_numpy', xrspatial_slope_numpy)
@@ -61,15 +62,15 @@ def test_numpy_equals_dask_qgis_data(elevation_raster):
     # compare using the data run through QGIS
     numpy_agg = input_data(elevation_raster, 'numpy')
     dask_agg = input_data(elevation_raster, 'dask+numpy')
-    assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, slope)
+    assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, slope, verify_attrs=False)
 
 
 @cuda_and_cupy_available
 def test_numpy_equals_cupy_qgis_data(elevation_raster):
     # compare using the data run through QGIS
     numpy_agg = input_data(elevation_raster, 'numpy')
     cupy_agg = input_data(elevation_raster, 'cupy')
-    assert_numpy_equals_cupy(numpy_agg, cupy_agg, slope)
+    assert_numpy_equals_cupy(numpy_agg, cupy_agg, slope, verify_attrs=False)
 
 
 @dask_array_available
@@ -80,7 +81,8 @@ def test_numpy_equals_cupy_qgis_data(elevation_raster):
 def test_numpy_equals_dask_cupy_random_data(random_data):
     numpy_agg = create_test_raster(random_data, backend='numpy')
     dask_cupy_agg = create_test_raster(random_data, backend='dask+cupy')
-    assert_numpy_equals_dask_cupy(numpy_agg, dask_cupy_agg, slope, atol=1e-6, rtol=1e-6)
+    assert_numpy_equals_dask_cupy(numpy_agg, dask_cupy_agg, slope, atol=1e-6, rtol=1e-6,
+                                  verify_attrs=False)
 
 
 @dask_array_available
@@ -330,7 +332,7 @@ def test_degenerate_shape_numpy(shape):
     agg = create_test_raster(_degenerate_data(shape), backend='numpy',
                              attrs={'res': (1, 1)})
     result = slope(agg)
-    general_output_checks(agg, result)
+    general_output_checks(agg, result, verify_attrs=False)
     assert result.shape == shape
     assert np.all(np.isnan(result.data))
 
@@ -342,7 +344,8 @@ def test_degenerate_shape_dask_numpy(shape):
     numpy_agg = create_test_raster(data, backend='numpy', attrs={'res': (1, 1)})
     dask_agg = create_test_raster(data, backend='dask+numpy',
                                   attrs={'res': (1, 1)}, chunks=shape)
-    assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, slope, nan_edges=False)
+    assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, slope, nan_edges=False,
+                                   verify_attrs=False)
 
 
 @cuda_and_cupy_available
@@ -351,7 +354,8 @@ def test_degenerate_shape_cupy(shape):
     data = _degenerate_data(shape)
     numpy_agg = create_test_raster(data, backend='numpy', attrs={'res': (1, 1)})
     cupy_agg = create_test_raster(data, backend='cupy', attrs={'res': (1, 1)})
-    assert_numpy_equals_cupy(numpy_agg, cupy_agg, slope, nan_edges=False)
+    assert_numpy_equals_cupy(numpy_agg, cupy_agg, slope, nan_edges=False,
+                             verify_attrs=False)
 
 
 @dask_array_available
@@ -362,7 +366,8 @@ def test_degenerate_shape_dask_cupy(shape):
     numpy_agg = create_test_raster(data, backend='numpy', attrs={'res': (1, 1)})
     dask_cupy_agg = create_test_raster(data, backend='dask+cupy',
                                        attrs={'res': (1, 1)}, chunks=shape)
-    assert_numpy_equals_dask_cupy(numpy_agg, dask_cupy_agg, slope, nan_edges=False)
+    assert_numpy_equals_dask_cupy(numpy_agg, dask_cupy_agg, slope, nan_edges=False,
+                                  verify_attrs=False)
 
 
 @pytest.mark.parametrize("shape", DEGENERATE_SHAPES)
@@ -376,7 +381,7 @@ def test_degenerate_shape_geodesic(shape):
         coords={'lat': lat, 'lon': lon},
     )
     result = slope(raster, method='geodesic')
-    general_output_checks(raster, result)
+    general_output_checks(raster, result, verify_attrs=False)
     assert result.shape == shape
     assert np.all(np.isnan(result.data))
 
@@ -432,7 +437,7 @@ def _center_nan_data():
 def test_center_nan_propagates_numpy():
     agg = create_test_raster(_center_nan_data(), backend='numpy', attrs={'res': (1, 1)})
     result = slope(agg)
-    general_output_checks(agg, result)
+    general_output_checks(agg, result, verify_attrs=False)
     assert np.isnan(result.data[2, 2])
 
 
@@ -442,7 +447,8 @@ def test_center_nan_propagates_dask_numpy():
     numpy_agg = create_test_raster(data, backend='numpy', attrs={'res': (1, 1)})
     dask_agg = create_test_raster(data, backend='dask+numpy',
                                   attrs={'res': (1, 1)}, chunks=(3, 3))
-    assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, slope, nan_edges=False)
+    assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, slope, nan_edges=False,
+                                   verify_attrs=False)
     assert np.isnan(slope(dask_agg).data.compute()[2, 2])
 
 
@@ -451,7 +457,8 @@ def test_center_nan_propagates_cupy():
     data = _center_nan_data()
     numpy_agg = create_test_raster(data, backend='numpy', attrs={'res': (1, 1)})
     cupy_agg = create_test_raster(data, backend='cupy', attrs={'res': (1, 1)})
-    assert_numpy_equals_cupy(numpy_agg, cupy_agg, slope, nan_edges=False)
+    assert_numpy_equals_cupy(numpy_agg, cupy_agg, slope, nan_edges=False,
+                             verify_attrs=False)
     assert np.isnan(slope(cupy_agg).data.get()[2, 2])
 
 
@@ -462,4 +469,67 @@ def test_center_nan_propagates_dask_cupy():
     numpy_agg = create_test_raster(data, backend='numpy', attrs={'res': (1, 1)})
     dask_cupy_agg = create_test_raster(data, backend='dask+cupy',
                                        attrs={'res': (1, 1)}, chunks=(3, 3))
-    assert_numpy_equals_dask_cupy(numpy_agg, dask_cupy_agg, slope, nan_edges=False)
+    assert_numpy_equals_dask_cupy(numpy_agg, dask_cupy_agg, slope, nan_edges=False,
+                                  verify_attrs=False)
+
+
+# slope is reported in degrees, so the output must advertise units='degrees'
+# rather than leaking the input elevation units (e.g. 'm'). Regression for #2894.
+def test_units_overridden_to_degrees_numpy():
+    agg = create_test_raster(_degenerate_data((5, 5)), backend='numpy',
+                             attrs={'res': (1, 1), 'units': 'm', 'crs': 'EPSG:5070'})
+    result = slope(agg)
+    assert result.attrs['units'] == 'degrees'
+    # other attrs are preserved untouched
+    assert result.attrs['res'] == (1, 1)
+    assert result.attrs['crs'] == 'EPSG:5070'
+    # input is not mutated
+    assert agg.attrs['units'] == 'm'
+
+
+@dask_array_available
+def test_units_overridden_to_degrees_dask_numpy():
+    agg = create_test_raster(_degenerate_data((5, 5)), backend='dask+numpy',
+                             attrs={'res': (1, 1), 'units': 'm'}, chunks=(3, 3))
+    result = slope(agg)
+    assert result.attrs['units'] == 'degrees'
+    assert agg.attrs['units'] == 'm'
+
+
+@cuda_and_cupy_available
+def test_units_overridden_to_degrees_cupy():
+    agg = create_test_raster(_degenerate_data((5, 5)), backend='cupy',
+                             attrs={'res': (1, 1), 'units': 'm'})
+    result = slope(agg)
+    assert result.attrs['units'] == 'degrees'
+    assert agg.attrs['units'] == 'm'
+
+
+@dask_array_available
+@cuda_and_cupy_available
+def test_units_overridden_to_degrees_dask_cupy():
+    agg = create_test_raster(_degenerate_data((5, 5)), backend='dask+cupy',
+                             attrs={'res': (1, 1), 'units': 'm'}, chunks=(3, 3))
+    result = slope(agg)
+    assert result.attrs['units'] == 'degrees'
+    assert agg.attrs['units'] == 'm'
+
+
+def test_units_set_to_degrees_when_input_has_no_units():
+    agg = create_test_raster(_degenerate_data((5, 5)), backend='numpy',
+                             attrs={'res': (1, 1)})
+    result = slope(agg)
+    assert result.attrs['units'] == 'degrees'
+
+
+def test_units_overridden_to_degrees_geodesic():
+    H, W = 5, 5
+    raster = xr.DataArray(
+        _degenerate_data((H, W)),
+        dims=['lat', 'lon'],
+        coords={'lat': np.linspace(40.0, 41.0, H), 'lon': np.linspace(10.0, 11.0, W)},
+        attrs={'units': 'm'},
+    )
+    result = slope(raster, method='geodesic')
+    assert result.attrs['units'] == 'degrees'
+    assert raster.attrs['units'] == 'm'