Stubbed out function for Analytics module (#621)

* stubbed out summarize_terrain function in analytics module #620

* added docs

* renamed  parameter to

* update tests

* add example

* create_test_raster(): added name param

* test analytics: use random_data pytest fixture

* flake8

* add a new line to avoid sphinx warnings

Co-authored-by: thuydotm <[email protected]>

Author: brendancol

xrspatial/analytics.py

@@ -0,0 +1,89 @@
+from xrspatial import slope
+from xrspatial import curvature
+from xrspatial import aspect
+
+import xarray as xr
+
+
+def summarize_terrain(terrain: xr.DataArray):
+    """
+    Calculates slope, aspect, and curvature of an elevation terrain and return a dataset
+    of the computed data.
+
+    Parameters
+    ----------
+    terrain: xarray.DataArray
+        2D NumPy, CuPy, or Dask with NumPy-backed xarray DataArray of elevation values.
+
+    Returns
+    -------
+    summarized_terrain: xarray.Dataset
+        Dataset with slope, aspect, curvature variables with a naming convention of
+        `terrain.name-variable_name`
+
+    Examples
+    --------
+    .. sourcecode:: python
+
+        >>> import numpy as np
+        >>> import xarray as xr
+        >>> from xrspatial.analytics import summarize_terrain
+        >>> data = np.array([
+            [0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 1, 0, 0, -1, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0]], dtype=np.float64)
+        >>> raster = xr.DataArray(data, name='myraster', attrs={'res': (1, 1)})
+        >>> summarized_terrain = summarize_terrain(raster)
+        >>> summarized_terrain
+        <xarray.Dataset>
+        Dimensions:             (dim_0: 5, dim_1: 8)
+        Dimensions without coordinates: dim_0, dim_1
+        Data variables:
+            myraster            (dim_0, dim_1) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
+            myraster-slope      (dim_0, dim_1) float32 nan nan nan nan ... nan nan nan
+            myraster-curvature  (dim_0, dim_1) float64 nan nan nan nan ... nan nan nan
+            myraster-aspect     (dim_0, dim_1) float32 nan nan nan nan ... nan nan nan
+        >>> summarized_terrain['myraster-slope']
+        <xarray.DataArray 'myraster-slope' (dim_0: 5, dim_1: 8)>
+        array([[      nan,       nan,       nan,       nan,       nan,       nan,       nan,   nan],
+               [      nan, 10.024988, 14.036243, 10.024988, 10.024988, 14.036243, 10.024988,   nan],
+               [      nan, 14.036243,  0.      , 14.036243, 14.036243,  0.      , 14.036243,   nan],
+               [      nan, 10.024988, 14.036243, 10.024988, 10.024988, 14.036243, 10.024988,   nan],
+               [      nan,       nan,       nan,       nan,       nan,       nan,       nan,   nan]], dtype=float32)  # noqa
+        Dimensions without coordinates: dim_0, dim_1
+        Attributes:
+            res:      (1, 1)
+
+        >>> summarized_terrain['myraster-curvature']
+        <xarray.DataArray 'myraster-curvature' (dim_0: 5, dim_1: 8)>
+        array([[  nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan],
+               [  nan,   -0., -100.,   -0.,   -0.,  100.,   -0.,   nan],
+               [  nan, -100.,  400., -100.,  100., -400.,  100.,   nan],
+               [  nan,   -0., -100.,   -0.,   -0.,  100.,   -0.,   nan],
+               [  nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan]])
+        Dimensions without coordinates: dim_0, dim_1
+        Attributes:
+            res:      (1, 1)
+
+        >>> summarized_terrain['myraster-aspect']
+        <xarray.DataArray 'myraster-aspect' (dim_0: 5, dim_1: 8)>
+        array([[ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan],
+               [ nan, 315.,   0.,  45., 135., 180., 225.,  nan],
+               [ nan, 270.,  -1.,  90.,  90.,  -1., 270.,  nan],
+               [ nan, 225., 180., 135.,  45.,   0., 315.,  nan],
+               [ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan]], dtype=float32)
+        Dimensions without coordinates: dim_0, dim_1
+        Attributes:
+            res:      (1, 1)
+    """
+
+    if terrain.name is None:
+        raise NameError('Requires xr.DataArray.name property to be set')
+
+    ds = terrain.to_dataset()
+    ds[f'{terrain.name}-slope'] = slope(terrain)
+    ds[f'{terrain.name}-curvature'] = curvature(terrain)
+    ds[f'{terrain.name}-aspect'] = aspect(terrain)
+    return ds

xrspatial/tests/general_checks.py

@@ -6,8 +6,10 @@
 from xrspatial.utils import ArrayTypeFunctionMapping
 
 
-def create_test_raster(data, backend='numpy', dims=['y', 'x'], attrs=None, chunks=(3, 3)):
-    raster = xr.DataArray(data, dims=dims, attrs=attrs)
+def create_test_raster(
+    data, backend='numpy', name='myraster', dims=['y', 'x'], attrs=None, chunks=(3, 3)
+):
+    raster = xr.DataArray(data, name=name, dims=dims, attrs=attrs)
     # set coords for test raster
     for i, dim in enumerate(dims):
         raster[dim] = np.linspace(0, data.shape[i] - 1, data.shape[i])

xrspatial/tests/test_analytics.py

@@ -0,0 +1,28 @@
+import pytest
+import numpy as np
+
+from xrspatial import aspect
+from xrspatial import curvature
+from xrspatial import slope
+from xrspatial.analytics import summarize_terrain
+
+from xrspatial.tests.general_checks import create_test_raster
+
+
+@pytest.mark.parametrize("size", [(2, 4), (100, 150)])
+@pytest.mark.parametrize(
+    "dtype", [np.int32, np.int64, np.uint32, np.uint64, np.float32, np.float64]
+)
+def test_summarize_terrain(random_data):
+    test_terrain = create_test_raster(random_data, name='myterrain')
+    summarized_ds = summarize_terrain(test_terrain)
+    variables = [v for v in summarized_ds]
+    should_have = ['myterrain',
+                   'myterrain-slope',
+                   'myterrain-curvature',
+                   'myterrain-aspect']
+    assert variables == should_have
+
+    np.testing.assert_allclose(summarized_ds['myterrain-slope'], slope(test_terrain))
+    np.testing.assert_allclose(summarized_ds['myterrain-curvature'], curvature(test_terrain))
+    np.testing.assert_allclose(summarized_ds['myterrain-aspect'], aspect(test_terrain))