Add terrain attributes (TRI, TPI, Roughness, Rugosity), add methods o…

…f slope, aspect, hillshade and test robustness against GDAL and richDEM (#227)

* fix deramp residual and optimizer

* add terrain references in function doc

* add tpi, tri, roughness and tests

* add maximum curvature

* add rugosity draft

* incremental commit

* fix documentation

* fix hillshade, refine tests comparing attributes to GDALDEM

* Add slope methods from Horn (1981) and slope_method and tri_method arguments

* fix spatial_tools package

* Fix scipy.ndimage.morphology deprecation for binary dilation

* Update doctest quadric coefficients

* Fix quadratic test on finite values now that default edge method had changed

* Add rugosity wrapper and references

* Fix doctests

* Add test skipping for coreg failures, issue opened

* Remove test on finite values for terrain attributes because default fill_method and edge_method do not fill gaps anymore

* Update value of patches method test with new NaNs of dataset diff

* Fix minor sphinx errors

* Fix tests after update of the diff example data

* Catch DeprecationWarning to exclude them from fatal test errors

* Fix quantile calculation for data with NaNs

* Fix the zeroing of profile and planform curvature to avoid NaNs on flat surfaces, and re-add a test of finite values

* Add warning for rugosity computation and window sizes

* Add or refine documentation and gallery example for all terrain attributes

* Move richDEM wrapper functions, move terrain ruggedness to maintain consistent order of roughness indexes, and fix minor docstrings issues

* Capitalize TRI and TPI attributes

* Correct error in profile curvature calculation, add docstring examples for all attributes

* Add missing return to line in docstring test

* Fix bug in rugosity calculation

* Add test for rugosity calculation

* Fix typo and clarify

* Adapt colorbar

* Fix and add rugosity tests

* Move rugosity to quadratic_coefficient due to fixed window size and need of resolution, and fix doctests

* Fix doctests

* Add richDEM methods, tests based on richDEM against curvatures and slope/aspect, tests on errors raised, changed default edge and fill methods to None and adjusted doctests

* Add richdem dev dependency for testing terrain attributes

* Fix doctests signs and spacing

* Fix doctest spacing issue

* Try to fix RuntimeWarning on NaN addition

* Add slope method example, add multiplot to terrain attribute example, modify thumbnails and link those examples to the documentation chapter

* Add comment to describe regex escape in test to catch raised errors

* Add fractal roughness, fix edge detection of get_windowed_indexes

* Fix doctest

* Add fractal roughness to documentation

* Fix addresses of scrapped gallery images

* Refine doc descriptions of individiual terrain attributes

* Clarify attribute units in documentation chapter

* Filter warning in aspect calculation and reenable test failure from warnings

dev-environment.yml

@@ -13,6 +13,7 @@ dependencies:
   - scipy
   - tqdm
   - scikit-image
+  - richdem
   - proj-data
   - scikit-gstat>=0.6.8
   - pytransform3d

docs/source/api.rst

@@ -14,7 +14,6 @@ Full information about xdem's functionality is provided on this page.
         coreg
         dem
         filters
-        spatial_tools
         spatialstats
         terrain
         volume

docs/source/robuststats.rst

@@ -58,7 +58,8 @@ The half difference between 84\ :sup:`th` and 16\ :sup:`th` percentiles, or the
 can also be used as a robust dispersion measure equivalent to the standard deviation.
 
 .. code-block:: python
-        nmad = xdem.spatialstats.nmad
+
+    nmad = xdem.spatialstats.nmad
 
 When working with weighted data, the difference between the 84\ :sup:`th` and 16\ :sup:`th` `weighted percentile <https://en.wikipedia.org
 /wiki/Percentile#Weighted_percentile>`_, or the absolute 68\ :sup:`th` weighted percentile can be used as a robust measure of dispersion.

docs/source/terrain.rst

@@ -13,19 +13,36 @@ Slope
 -----
 :func:`xdem.terrain.slope`
 
-The slope map of a DEM describes the tilt (or gradient) of each pixel in relation to its neighbours.
+The slope of a DEM describes the tilt, or gradient, of each pixel in relation to its neighbours.
 It is most often described in degrees, where a flat surface is 0° and a vertical cliff is 90°.
 No tilt direction is stored in the slope map; a 45° tilt westward is identical to a 45° tilt eastward.
 
+The slope can be computed either by the method of `Horn (1981) <http://dx.doi.org/10.1109/PROC.1981.11918>`_ (default)
+based on a refined gradient formulation on a 3x3 pixel window, or by the method of `Zevenbergen and Thorne (1987)
+<http://dx.doi.org/10.1002/esp.3290120107>`_ based on a plane fit on a 3x3 pixel window.
+
+The differences between methods are illustrated in the :ref:`sphx_glr_auto_examples_plot_terrain_attributes.py`
+example.
+
+.. image:: auto_examples/images/sphx_glr_plot_terrain_attributes_001.png
+  :width: 600
+
 .. minigallery:: xdem.terrain.slope
         :add-heading:
 
 Aspect
 ------
 :func:`xdem.terrain.aspect`
 
-The aspect map of a DEM describes which direction the slope is tilting of each pixel in relation to its neighbours.
-It is most often described in degrees, where a slope tilting straight north would have an aspect of 0°, an eastern aspect is 90°, south is 180° and west is 270°.
+The aspect describes the orientation of strongest slope.
+It is often reported in degrees, where a slope tilting straight north corresponds to an aspect of 0°, and an eastern
+aspect is 90°, south is 180° and west is 270°. By default, a flat slope is given an arbitrary aspect of 180°.
+
+As the aspect is directly based on the slope, it varies between the method of `Horn (1981) <http://dx.doi.org/10.
+1109/PROC.1981.11918>`_ (default) and that of `Zevenbergen and Thorne (1987) <http://dx.doi.org/10.1002/esp.3290120107>`_.
+
+.. image:: auto_examples/images/sphx_glr_plot_terrain_attributes_002.png
+  :width: 600
 
 .. minigallery:: xdem.terrain.aspect
         :add-heading:
@@ -41,8 +58,16 @@ Hillshades are therefore often preferable for visualizing DEMs.
 With a westerly azimuth (a simulated sun coming from the west), all eastern slopes are slightly darker.
 This mode of shading the slopes often generates a map that is much more easily interpreted than the slope map.
 
+
+As the hillshade is directly based on the slope and aspect, it varies between the method of `Horn (1981) <http://dx.doi
+.org/10.1109/PROC.1981.11918>`_ (default) and that of `Zevenbergen and Thorne (1987) <http://dx.doi.org/10.1002/esp.
+3290120107>`_.
+
 Note, however, that the hillshade is not a shadow map; no occlusion is taken into account so it does not represent "true" shading.
-It therefore has little analytic purpose, but it still a great tool for visualization.
+It therefore has little analytic purpose, but it still constitutes a great visualization tool.
+
+.. image:: auto_examples/images/sphx_glr_plot_terrain_attributes_003.png
+  :width: 600
 
 .. minigallery:: xdem.terrain.hillshade
         :add-heading:
@@ -51,13 +76,15 @@ Curvature
 ---------
 :func:`xdem.terrain.curvature`
 
-The curvature map is the second derivative of elevation.
-It highlights the convexity or convavity of the terrain, and has many both analytic and visual purposes.
+The curvature map is the second derivative of elevation, which highlights the convexity or concavity of the terrain.
 If a surface is convex (like a mountain peak), it will have positive curvature.
 If a surface is concave (like a trough or a valley bottom), it will have negative curvature.
+The curvature values in units of m\ :sup:`-1` are quite small, so they are by convention multiplied by 100.
+
+The curvature is based on the method of `Zevenbergen and Thorne (1987) <http://dx.doi.org/10.1002/esp.3290120107>`_.
 
-Usually, the curvature values are quite small, so they are by convention multiplied by 100.
-For analytic purposes, it may therefore be worth considering dividing the output by 100.
+.. image:: auto_examples/images/sphx_glr_plot_terrain_attributes_004.png
+  :width: 600
 
 .. minigallery:: xdem.terrain.curvature
         :add-heading:
@@ -66,7 +93,12 @@ Planform curvature
 ------------------
 :func:`xdem.terrain.planform_curvature`
 
-TODO: Add text.
+The planform curvature is the curvature perpendicular to the direction of slope, reported in 100 m\ :sup:`-1`.
+
+It is based on the method of `Zevenbergen and Thorne (1987) <http://dx.doi.org/10.1002/esp.3290120107>`_.
+
+.. image:: auto_examples/images/sphx_glr_plot_terrain_attributes_005.png
+  :width: 600
 
 .. minigallery:: xdem.terrain.planform_curvature
         :add-heading:
@@ -75,16 +107,101 @@ Profile curvature
 -----------------
 :func:`xdem.terrain.profile_curvature`
 
-TODO: Add text.
+The profile curvature is the curvature parallel to the direction of slope, reported in 100 m\ :sup:`-1`..
+
+It is based on the method of `Zevenbergen and Thorne (1987) <http://dx.doi.org/10.1002/esp.3290120107>`_.
+
+.. image:: auto_examples/images/sphx_glr_plot_terrain_attributes_006.png
+  :width: 600
 
 .. minigallery:: xdem.terrain.profile_curvature
         :add-heading:
 
+Topographic Position Index
+--------------------------
+:func:`xdem.terrain.topographic_position_index`
+
+The Topographic Position Index (TPI) is a metric of slope position, based on the method of `Weiss (2001) <http://www
+.jennessent.com/downloads/TPI-poster-TNC_18x22.pdf>`_ that corresponds to the difference of the elevation of a central
+pixel with the average of that of neighbouring pixels. Its unit is that of the DEM (typically meters) and it can be
+computed for any window size (default 3x3 pixels).
+
+.. image:: auto_examples/images/sphx_glr_plot_terrain_attributes_007.png
+  :width: 600
+
+.. minigallery:: xdem.terrain.topographic_position_index
+        :add-heading:
+
+Terrain Ruggedness Index
+------------------------
+:func:`xdem.terrain.terrain_ruggedness_index`
+
+The Terrain Ruggedness Index (TRI) is a metric of terrain ruggedness, based on cumulated differences in elevation between
+a central pixel and its surroundings. Its unit is that of the DEM (typically meters) and it can be computed for any
+window size (default 3x3 pixels).
+
+For topography (default), the method of `Riley et al. (1999) <http://download.osgeo.org/qgis/doc/reference-docs/Terrain_
+Ruggedness_Index.pdf>`_ is generally used, where the TRI is computed by the squareroot of squared differences with
+neighbouring pixels.
+
+For bathymetry, the method of `Wilson et al. (2007) <http://dx.doi.org/10.1080/01490410701295962>`_ is generally used,
+where the TRI is defined by the mean absolute difference with neighbouring pixels
+
+.. image:: auto_examples/images/sphx_glr_plot_terrain_attributes_008.png
+  :width: 600
+
+.. minigallery:: xdem.terrain.terrain_ruggedness_index
+        :add-heading:
+
+Roughness
+---------
+:func:`xdem.terrain.roughness`
+
+The roughness is a metric of terrain ruggedness, based on the maximum difference in elevation in the surroundings.
+The roughness is based on the method of `Dartnell (2000) <http://dx.doi.org/10.14358/PERS.70.9.
+1081>`_. Its unit is that of the DEM (typically meters) and it can be computed for any window size (default 3x3 pixels).
+
+.. image:: auto_examples/images/sphx_glr_plot_terrain_attributes_009.png
+  :width: 600
+
+.. minigallery:: xdem.terrain.roughness
+        :add-heading:
+
+Rugosity
+--------
+:func:`xdem.terrain.rugosity`
+
+The rugosity is a metric of terrain ruggedness, based on the ratio between planimetric and real surface area. The
+rugosity is based on the method of `Jenness (2004) <https://doi.org/10.2193/0091-7648(2004)032[0829:CLSAFD]2.0.CO;2>`_.
+It is unitless, and is only supported for a 3x3 window size.
+
+.. image:: auto_examples/images/sphx_glr_plot_terrain_attributes_010.png
+  :width: 600
+
+.. minigallery:: xdem.terrain.rugosity
+        :add-heading:
+
+Fractal roughness
+-----------------
+:func:`xdem.terrain.fractal_roughness`
+
+The fractal roughness is a metric of terrain ruggedness based on the local fractal dimension estimated by the volume
+box-counting method of `Taud and Parrot (2005) <https://doi.org/10.4000/geomorphologie.622>`_.
+The fractal roughness is computed by estimating the fractal dimension in 3D space, for a local window centered on the
+DEM pixels. Its unit is that of a dimension, and is always between 1 (dimension of a line in 3D space) and 3
+(dimension of a cube in 3D space). It can only be computed on window sizes larger than 5x5 pixels, and defaults to 13x13.
+
+.. image:: auto_examples/images/sphx_glr_plot_terrain_attributes_011.png
+  :width: 600
+
+.. minigallery:: xdem.terrain.fractal_roughness
+        :add-heading:
+
 Generating multiple attributes at once
 --------------------------------------
 
 Often, one may seek more terrain attributes than one, e.g. both the slope and the aspect.
-Since both are dependent on the gradient of the DEM, calculating them separately is unneccesarily repetitive.
+Since both are dependent on the gradient of the DEM, calculating them separately is computationally repetitive.
 Multiple terrain attributes can be calculated from the same gradient using the :func:`xdem.terrain.get_terrain_attribute` function.
 
 .. minigallery:: xdem.terrain.get_terrain_attribute

examples/plot_nonstationary_error.py

@@ -155,13 +155,13 @@
 #
 # If we use custom quantiles for both binning variables, and adjust the plot scale:
 
-custom_bin_slope = np.unique(np.concatenate([np.quantile(slope_arr,np.linspace(0,0.95,20)),
-                                             np.quantile(slope_arr,np.linspace(0.96,0.99,5)),
-                                             np.quantile(slope_arr,np.linspace(0.991,1,10))]))
+custom_bin_slope = np.unique(np.concatenate([np.nanquantile(slope_arr,np.linspace(0,0.95,20)),
+                                             np.nanquantile(slope_arr,np.linspace(0.96,0.99,5)),
+                                             np.nanquantile(slope_arr,np.linspace(0.991,1,10))]))
 
-custom_bin_curvature = np.unique(np.concatenate([np.quantile(maxc_arr,np.linspace(0,0.95,20)),
-                                             np.quantile(maxc_arr,np.linspace(0.96,0.99,5)),
-                                             np.quantile(maxc_arr,np.linspace(0.991,1,10))]))
+custom_bin_curvature = np.unique(np.concatenate([np.nanquantile(maxc_arr,np.linspace(0,0.95,20)),
+                                             np.nanquantile(maxc_arr,np.linspace(0.96,0.99,5)),
+                                             np.nanquantile(maxc_arr,np.linspace(0.991,1,10))]))
 
 df = xdem.spatialstats.nd_binning(values=dh_arr, list_var=[slope_arr, maxc_arr], list_var_names=['slope', 'maxc'],
                                   statistics=['count', np.nanmedian, xdem.spatialstats.nmad],

examples/plot_slope_methods.py

@@ -0,0 +1,111 @@
+"""
+Slope and aspect methods
+========================
+
+Terrain slope and aspect can be estimated using different methods.
+Here is an example of how to generate the two with each method, and understand their differences.
+
+For more information, see the :ref:`terrain_attributes` chapter and the
+:ref:`sphx_glr_auto_examples_plot_terrain_attributes.py` example.
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+
+import xdem
+
+# %%
+# **Example data**
+
+dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
+
+
+def plot_attribute(attribute, cmap, label=None, vlim=None):
+    plt.figure(figsize=(8, 5))
+
+    if vlim is not None:
+        if isinstance(vlim, (int, float)):
+            vlims = {"vmin": -vlim, "vmax": vlim}
+        elif len(vlim) == 2:
+            vlims = {"vmin": vlim[0], "vmax": vlim[1]}
+    else:
+        vlims = {}
+
+    plt.imshow(
+        attribute.squeeze(),
+        cmap=cmap,
+        extent=[dem.bounds.left, dem.bounds.right, dem.bounds.bottom, dem.bounds.top],
+        **vlims
+    )
+    if label is not None:
+        cbar = plt.colorbar()
+        cbar.set_label(label)
+
+    plt.xticks([])
+    plt.yticks([])
+    plt.tight_layout()
+
+    plt.show()
+
+
+# %%
+# Slope with method of `Horn (1981) <http://dx.doi.org/10.1109/PROC.1981.11918>`_ (GDAL default), based on a refined
+# approximation of the gradient  (page 18, bottom left, and pages 20-21).
+
+slope_horn = xdem.terrain.slope(dem.data, resolution=dem.res)
+
+plot_attribute(slope_horn, "Reds", "Slope of Horn (1981) (°)")
+
+# %%
+# Slope with method of `Zevenbergen and Thorne (1987) <http://dx.doi.org/10.1002/esp.3290120107>`_, Equation 13.
+
+slope_zevenberg = xdem.terrain.slope(dem.data, resolution=dem.res, method='ZevenbergThorne')
+
+plot_attribute(slope_zevenberg, "Reds", "Slope of Zevenberg and Thorne (1987) (°)")
+
+# %%
+# We compute the difference between the slopes computed with each method.
+
+diff_slope = slope_horn - slope_zevenberg
+
+plot_attribute(diff_slope, "RdYlBu", "Slope of Horn (1981) minus\n slope of Zevenberg and Thorne (1987) (°)", vlim=3)
+
+# %%
+# The differences are negative, implying that the method of Horn always provides flatter slopes.
+# Additionally, they seem to occur in places of high curvatures. We verify this by plotting the maximum curvature.
+
+maxc = xdem.terrain.maximum_curvature(dem.data, resolution=dem.res)
+
+plot_attribute(maxc, "RdYlBu", "Maximum curvature (100 m $^{-1}$)", vlim=2)
+
+# %%
+# We quantify the relationship by computing the median of slope differences in bins of curvatures, and plot the
+# result. We define custom bins for curvature, due to its skewed distribution.
+
+df_bin = xdem.spatialstats.nd_binning(values=diff_slope[:],list_var=[maxc[:]], list_var_names=["maxc"],
+                                      list_var_bins=30, statistics=[np.nanmedian, "count"])
+
+xdem.spatialstats.plot_1d_binning(df_bin, var_name="maxc", statistic_name="nanmedian",
+                                  label_var="Maximum absolute curvature (100 m$^{-1}$)",
+                                  label_statistic="Slope of Horn (1981) minus\n "
+                                                  "slope of Zevenberg and Thorne (1987) (°)")
+
+
+# %%
+# We perform the same exercise to analyze the differences in terrain aspect. We compute the difference modulo 360°,
+# to account for the circularity of aspect.
+
+aspect_horn = xdem.terrain.aspect(dem.data)
+aspect_zevenberg = xdem.terrain.aspect(dem.data, method='ZevenbergThorne')
+
+diff_aspect = aspect_horn - aspect_zevenberg
+diff_aspect_mod = np.minimum(np.mod(diff_aspect, 360), 360 - np.mod(diff_aspect, 360))
+
+plot_attribute(diff_aspect_mod, "Spectral", "Aspect of Horn (1981) minus\n aspect of Zevenberg and Thorne (1987) (°)",
+               vlim=[0, 90])
+
+# %%
+# Same as for slope, differences in aspect seem to coincide with high curvature areas. We observe also observe large
+# differences for areas with nearly flat slopes, owing to the high sensitivity of orientation estimation
+# for flat terrain.
+
+# Note: the default aspect for a 0° slope is 180°, as in GDAL.

examples/plot_standardization.py

@@ -49,13 +49,13 @@
 dh_arr[np.abs(dh_arr) > 4 * xdem.spatialstats.nmad(dh_arr)] = np.nan
 
 # Define bins for 2D binning
-custom_bin_slope = np.unique(np.concatenate([np.quantile(slope_arr,np.linspace(0,0.95,20)),
-                                             np.quantile(slope_arr,np.linspace(0.96,0.99,5)),
-                                             np.quantile(slope_arr,np.linspace(0.991,1,10))]))
+custom_bin_slope = np.unique(np.concatenate([np.nanquantile(slope_arr,np.linspace(0,0.95,20)),
+                                             np.nanquantile(slope_arr,np.linspace(0.96,0.99,5)),
+                                             np.nanquantile(slope_arr,np.linspace(0.991,1,10))]))
 
-custom_bin_curvature = np.unique(np.concatenate([np.quantile(maxc_arr,np.linspace(0,0.95,20)),
-                                             np.quantile(maxc_arr,np.linspace(0.96,0.99,5)),
-                                             np.quantile(maxc_arr,np.linspace(0.991,1,10))]))
+custom_bin_curvature = np.unique(np.concatenate([np.nanquantile(maxc_arr,np.linspace(0,0.95,20)),
+                                             np.nanquantile(maxc_arr,np.linspace(0.96,0.99,5)),
+                                             np.nanquantile(maxc_arr,np.linspace(0.991,1,10))]))
 
 # Perform 2D binning to estimate the measurement error with slope and maximum curvature
 df = xdem.spatialstats.nd_binning(values=dh_arr, list_var=[slope_arr, maxc_arr], list_var_names=['slope', 'maxc'],