From af1a713efcdfe57270f282b098aeb6310fa17f47 Mon Sep 17 00:00:00 2001 From: thuydotm Date: Fri, 12 Jun 2020 17:18:01 -0400 Subject: [PATCH 1/2] _horn_slope: accept cellsize_x and cellsize_y arguments, validate 'res' attr of input xarray --- xrspatial/slope.py | 43 ++++++++++++++++------------------- xrspatial/tests/test_slope.py | 41 +++++++++++++++++++++++++++++++++ 2 files changed, 60 insertions(+), 24 deletions(-) diff --git a/xrspatial/slope.py b/xrspatial/slope.py index 91d616c4..42946b24 100644 --- a/xrspatial/slope.py +++ b/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 @@ -28,35 +28,30 @@ 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 - """ - 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.') - - slope_agg = _horn_slope(agg.data, agg.attrs['res']) + raise ValueError('input xarray must have `res` attr.') + + # 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) - diff --git a/xrspatial/tests/test_slope.py b/xrspatial/tests/test_slope.py index bb9ed352..d4103914 100644 --- a/xrspatial/tests/test_slope.py +++ b/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 From a1bf568d63bfca293dfb9796ced8de8d287225f2 Mon Sep 17 00:00:00 2001 From: thuydotm Date: Sat, 13 Jun 2020 19:39:39 -0400 Subject: [PATCH 2/2] slope: docstring --- xrspatial/slope.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/xrspatial/slope.py b/xrspatial/slope.py index 42946b24..0cf6ef44 100644 --- a/xrspatial/slope.py +++ b/xrspatial/slope.py @@ -28,6 +28,22 @@ def _horn_slope(data, cellsize_x, cellsize_y): # 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 + """ + if not isinstance(agg, DataArray): raise TypeError("agg must be instance of DataArray")