slope: cellsize_x, cellsize_y from res attribute of input xarray

xrspatial/slope.py

@@ -6,7 +6,7 @@
 
 
 @ngjit
-def _horn_slope(data, cellsize):
+def _horn_slope(data, cellsize_x, cellsize_y):
     out = np.zeros_like(data)
     rows, cols = data.shape
     for y in range(1, rows-1):
@@ -19,8 +19,8 @@ def _horn_slope(data, cellsize):
             g = data[y-1, x-1]
             h = data[y-1, x]
             i = data[y-1, x+1]
-            dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * cellsize)
-            dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * cellsize)
+            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)
             p = (dz_dx * dz_dx + dz_dy * dz_dy) ** .5
             out[y, x] = np.arctan(p) * 57.29578
     return out
@@ -29,34 +29,45 @@ def _horn_slope(data, cellsize):
 # TODO: add optional name parameter `name='slope'`
 def slope(agg):
     """Returns slope of input aggregate in degrees.
-
     Parameters
     ----------
     agg : DataArray
-
     Returns
     -------
     data: DataArray
-
     Notes:
     ------
     Algorithm References:
      - http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-slope-works.htm
-     - Burrough, P. A., and McDonell, R. A., 1998. Principles of Geographical Information Systems (Oxford University Press, New York), pp 406
+     - Burrough, P. A., and McDonell, R. A., 1998.
+      Principles of Geographical Information Systems
+      (Oxford University Press, New York), pp 406
     """
 
     if not isinstance(agg, DataArray):
         raise TypeError("agg must be instance of DataArray")
 
     if not agg.attrs.get('res'):
         #TODO: maybe monkey-patch a "res" attribute valueing unity is reasonable
-        raise ValueError('input xarray must have numeric `res` attr.')
+        raise ValueError('input xarray must have `res` attr.')
 
-    slope_agg = _horn_slope(agg.data, agg.attrs['res'])
+    # get cellsize out from 'res' attribute
+    cellsize = agg.attrs.get('res')
+    if isinstance(cellsize, tuple) and len(cellsize) == 2 \
+            and isinstance(cellsize[0], (int, float)) \
+            and isinstance(cellsize[1], (int, float)):
+        cellsize_x, cellsize_y = cellsize
+    elif isinstance(cellsize, (int, float)):
+        cellsize_x = cellsize
+        cellsize_y = cellsize
+    else:
+        raise ValueError('`res` attr of input xarray must be a numeric'
+                         ' or a tuple of numeric values.')
+
+    slope_agg = _horn_slope(agg.data, cellsize_x, cellsize_y)
 
     return DataArray(slope_agg,
                      name='slope',
                      coords=agg.coords,
                      dims=agg.dims,
                      attrs=agg.attrs)
-

xrspatial/tests/test_slope.py

@@ -1,3 +1,4 @@
+import pytest
 import xarray as xr
 import numpy as np
 
@@ -33,6 +34,31 @@ def _do_gaussian_array():
 data_gaussian = _do_gaussian_array()
 
 
+def test_invalid_res_attr():
+    """
+    Assert 'res' attribute of input xarray
+    """
+    da = xr.DataArray(data_gaussian, attrs={})
+    with pytest.raises(ValueError) as e_info:
+        da_slope = slope(da)
+        assert e_info
+
+    da = xr.DataArray(data_gaussian, attrs={'res': 'any_string'})
+    with pytest.raises(ValueError) as e_info:
+        da_slope = slope(da)
+        assert e_info
+
+    da = xr.DataArray(data_gaussian, attrs={'res': ('str_tuple', 'str_tuple')})
+    with pytest.raises(ValueError) as e_info:
+        da_slope = slope(da)
+        assert e_info
+
+    da = xr.DataArray(data_gaussian, attrs={'res': (1, 2, 3)})
+    with pytest.raises(ValueError) as e_info:
+        da_slope = slope(da)
+        assert e_info
+
+
 def test_slope_transfer_function():
     """
     Assert slope transfer function
@@ -50,3 +76,18 @@ def test_slope_transfer_function():
     # And there the slope must be zero.
     _imax = np.where(da == da.max())
     assert da_slope[_imax] == 0
+
+    # same result when cellsize_x = cellsize_y = 1
+    da = xr.DataArray(data_gaussian, attrs={'res': (1.0, 1.0)})
+    da_slope = slope(da)
+    assert da_slope.dims == da.dims
+    assert da_slope.coords == da.coords
+    assert da_slope.attrs == da.attrs
+    assert da.shape == da_slope.shape
+
+    assert da_slope.sum() > 0
+
+    # In the middle of the array, there is the maximum of the gaussian;
+    # And there the slope must be zero.
+    _imax = np.where(da == da.max())
+    assert da_slope[_imax] == 0