Comments from Amaury and Erik

doc/source/background.md

@@ -58,7 +58,7 @@ Amaury Dehecq<sup>2</sup> and took place online on November 8, 2020.
 
 The initial core development of xDEM was performed by members of the Glaciology group of the Laboratory of Hydraulics, Hydrology and
 Glaciology (VAW) at ETH Zürich<sup>3</sup>, with contributions by members of the University of Oslo, the University of Washington, and University
-Grenobles Alpes.
+Grenoble Alpes.
 
 We are not software developers but geoscientists, and we try our best to offer tools that can be useful to a larger group,
 documented, reliable and maintained. All development and maintenance is made on a voluntary basis and we welcome

doc/source/biascorr.md

@@ -311,8 +311,8 @@ _ = ax[1].set_yticklabels([])
 
 - **Performs:** Correct biases along a terrain attribute.
 - **Supports weights:** Yes.
-- **Pros:** Useful to correct for instance curvature bias due to different native resolution between elevation data.
-- **Cons:** For curvature biases, only works for elevation data with relatively close natire resolution.
+- **Pros:** Useful to correct for instance curvature-related bias due to different native resolution between elevation data.
+- **Cons:** For curvature-related biases, only works for elevation data with relatively close native resolution.
 
 The default optimizer for terrain biases optimizes a 1D polynomial with an order from 1 to 6,
 and keeps the best performing fit.

doc/source/citation.md

@@ -94,7 +94,7 @@ More details are available on each feature page!
    * - Bilinear
      - N/A
    * - Local and regional hypsometric
-     - [McNabb et al. (2019)](https://tc.copernicus.org/articles/13/895/2019/)
+     - [Arendt et al. (2002)](https://doi.org/10.1126/science.1072497), [McNabb et al. (2019)](https://tc.copernicus.org/articles/13/895/2019/)
 ```
 
 

doc/source/code/intricacies_datatypes.py

@@ -52,7 +52,7 @@
 ax[1, 0].set_yticklabels([])
 ax[1, 0].set_aspect("equal")
 
-plt.title("Types of elevation data")
+plt.suptitle("Types of elevation data")
 
 plt.tight_layout()
 plt.show()

doc/source/coregistration.md

@@ -167,7 +167,8 @@ See coregistration on real data in the **{ref}`examples-basic` and {ref}`example
 
 The [Nuth and Kääb (2011)](https://doi.org/10.5194/tc-5-271-2011) coregistration approach estimates a horizontal
 translation iteratively by solving a cosine equation between the terrain slope, aspect and the elevation differences.
-The iteration stops if it reaches the maximum number of iteration limit or if the tolerance does not improve.
+The iteration stops if it reaches the maximum number of iteration limit, or if the iterative shift amplitude falls 
+below a specified tolerance.
 
 ```{code-cell} ipython3
 :tags: [hide-cell]
@@ -308,7 +309,7 @@ plt.tight_layout()
 - **Cons:** Poor sub-pixel accuracy for horizontal shifts, sensitive to outliers, and runs slowly with large samples.
 
 Iterative Closest Point (ICP) coregistration is an iterative point cloud registration method from [Besl and McKay (1992)](https://doi.org/10.1117/12.57955). It aims at iteratively minimizing the distance between closest neighbours by applying sequential rigid transformations. If DEMs are used as inputs, they are converted to point clouds.
-As for Nuth and Kääb (2011), the iteration stops if it reaches the maximum number of iteration limit or if the tolerance does not improve.
+As for Nuth and Kääb (2011), the iteration stops if it reaches the maximum number of iteration limit or if the iterative transformation amplitude falls below a specified tolerance.
 
 ICP is currently based on [OpenCV's implementation](https://docs.opencv.org/4.x/dc/d9b/classcv_1_1ppf__match__3d_1_1ICP.html) (an optional dependency), which includes outlier removal arguments. This may improve results significantly on outlier-prone data, but care should still be taken, as the risk of landing in [local minima](https://en.wikipedia.org/wiki/Maxima_and_minima) increases.
 

doc/source/dem_class.md

@@ -94,7 +94,7 @@ vect_gla = vect_gla.crop(dem)
 
 # Plot the DEM and the vector file
 dem.plot(cmap="terrain", cbar_title="Elevation (m)")
-vect_gla.plot(dem, ec="k", fc="none")  # We pass the DEM as reference for the plot CRS/extent
+vect_gla.plot(dem, ec="k", fc="none")  # We pass the DEM as reference for the plot CRS
 ```
 
 ## Vertical referencing
@@ -176,3 +176,9 @@ sig_dem.plot(cmap="Purples", cbar_title=r"Random error in elevation (1$\sigma$,
 # The spatial correlation function represents how much errors are correlated at a certain distance
 print("Elevation errors at a distance of 1 km are correlated at {:.2f} %.".format(rho_sig(1000) * 100))
 ```
+
+```{note}
+We use `random_state` to ensure a fixed randomized output. It is **only necessary if you need your results to be exactly reproductible**. 
+
+For more details on quantifying random and structured errors, see the {ref}`uncertainty` page.
+```

doc/source/ecosystem.md

@@ -26,6 +26,11 @@ Complementary Python tools to **analyze elevation data** are for instance:
 If you are working in Julia, the [Geomorphometry](https://github.com/Deltares/Geomorphometry.jl) package provides a
 wide range of terrain analysis for elevation data.
 
+## R
+
+If you are working in R, the [MultiscaleDTM](https://ailich.github.io/MultiscaleDTM/) package provides modular tools 
+for terrain analysis at multiple scales!
+
 ## Other community resources
 
 Whether to retrieve data among their wide range of open datasets, or to dive into their other resources, be sure to check out the

doc/source/elevation_intricacies.md

@@ -28,7 +28,7 @@ See the **{ref}`elevation-objects` features pages** for more details.
 
 
 Additionally, there are a critical differences for elevation point clouds depending on point density:
-- **Sparse elevation point clouds** (e.g., altimetry) are generally be stored as small vector-type datasets (e.g., SHP). Due to their sparsity, for subsequent analysis, they are rarely gridded into a DEM, and instead compared with DEMs at the point cloud coordinates by interpolation of the DEM,
+- **Sparse elevation point clouds** (e.g., altimetry) are generally stored as small vector-type datasets (e.g., SHP). Due to their sparsity, for subsequent analysis, they are rarely gridded into a DEM, and instead compared with DEMs at the point cloud coordinates by interpolation of the DEM,
 - **Dense elevation point clouds** (e.g., lidar) are large datasets generally stored in specific formats (LAS). Due to their high density, they are often gridded into DEMs by triangular interpolation of the point cloud.
 
 ```{note}

doc/source/gapfill.md

@@ -16,13 +16,13 @@ kernelspec:
 xDEM contains routines to gap-fill elevation data or elevation differences depending on the type of terrain.
 
 ```{important}
-Some of the approaches below are application-specific (e.g., glaciers) and might be moved to a separate package
+Most of the approaches below are application-specific (e.g., glaciers) and might be moved to a separate package
 in future releases.
 ```
 
 So far, xDEM has three types of gap-filling methods:
 
-- Linear spatial interpolation,
+- Inverse-distance weighting interpolation,
 - Local hypsometric interpolation (only relevant for elevation differences and glacier applications),
 - Regional hypsometric interpolation (also for glaciers).
 
@@ -79,18 +79,18 @@ We create a difference of DEMs object {class}`xdem.ddem.dDEM` to experiment on:
 ddem = xdem.dDEM(raster=dem_2009 - dem_1990, start_time=dem_1990.datetime, end_time=dem_2009.datetime)
 
 # The example DEMs are void-free, so let's make some random voids.
-# Introduce 50000 nans randomly throughout the dDEM.
+# Introduce a fifth of nans randomly throughout the dDEM.
 mask = np.zeros_like(ddem.data, dtype=bool)
 mask.ravel()[(np.random.choice(ddem.data.size, int(ddem.data.size/5), replace=False))] = True
 ddem.set_mask(mask)
 ```
 
-## Linear spatial interpolation
+## Inverse-distance weighting interpolation
 
-Linear spatial interpolation (also often called bilinear interpolation) of dDEMs is arguably the simplest approach: voids are filled by an average of the surrounding pixels values, weighted by their distance to the void pixel.
+Inverse-distance weighting (IDW) interpolation of elevation differences is arguably the simplest approach: voids are filled by a weighted-mean of the surrounding pixels values, with weight inversely proportional to their distance to the void pixel.
 
 ```{code-cell} ipython3
-ddem_linear = ddem.interpolate(method="linear")
+ddem_idw = ddem.interpolate(method="idw")
 ```
 
 ```{code-cell} ipython3
@@ -99,21 +99,21 @@ ddem_linear = ddem.interpolate(method="linear")
 :  code_prompt_show: "Show plotting code"
 :  code_prompt_hide: "Hide plotting code"
 
-ddem_linear = ddem.copy(new_array=ddem_linear)
+ddem_idw = ddem.copy(new_array=ddem_idw)
 
 # Plot before and after
 f, ax = plt.subplots(1, 2)
-ax[0].set_title("Before linear\ngap-filling")
-ddem.plot(cmap='RdYlBu', vmin=-10, vmax=10, ax=ax[0])
-ax[1].set_title("After linear\ngap-filling")
-ddem_linear.plot(cmap='RdYlBu', vmin=-10, vmax=10, ax=ax[1], cbar_title="Elevation differences (m)")
+ax[0].set_title("Before IDW\ngap-filling")
+ddem.plot(cmap='RdYlBu', vmin=-20, vmax=20, ax=ax[0])
+ax[1].set_title("After IDW\ngap-filling")
+ddem_idw.plot(cmap='RdYlBu', vmin=-20, vmax=20, ax=ax[1], cbar_title="Elevation differences (m)")
 _ = ax[1].set_yticklabels([])
 ```
 
 ## Local hypsometric interpolation
 
 This approach assumes that there is a relationship between the elevation and the elevation change in the dDEM, which is often the case for glaciers.
-Elevation change gradients in late 1900s and 2000s on glaciers often have the signature of large melt in the lower parts, while the upper parts might be less negative, or even positive.
+Elevation change gradients in late 1900s and 2000s on glaciers often have the signature of large thinning in the lower parts, while the upper parts might be less negative, or even positive.
 This relationship is strongly correlated for a specific glacier, and weakly correlated on regional scales.
 With the local (glacier specific) hypsometric approach, elevation change gradients are estimated for each glacier separately.
 This is simply a linear or polynomial model estimated with the dDEM and a reference DEM.
@@ -134,9 +134,9 @@ ddem_localhyps = ddem.copy(new_array=ddem_localhyps)
 # Plot before and after
 f, ax = plt.subplots(1, 2)
 ax[0].set_title("Before local\nhypsometric\ngap-filling")
-ddem.plot(cmap='RdYlBu', vmin=-10, vmax=10, ax=ax[0])
+ddem.plot(cmap='RdYlBu', vmin=-20, vmax=20, ax=ax[0])
 ax[1].set_title("After local\nhypsometric\ngap-filling")
-ddem_localhyps.plot(cmap='RdYlBu', vmin=-10, vmax=10, ax=ax[1], cbar_title="Elevation differences (m)")
+ddem_localhyps.plot(cmap='RdYlBu', vmin=-20, vmax=20, ax=ax[1], cbar_title="Elevation differences (m)")
 _ = ax[1].set_yticklabels([])
 ```
 
@@ -201,10 +201,9 @@ plt.show()
 
 Similarly to local hypsometric interpolation, the elevation change is assumed to be largely elevation-dependent.
 With the regional approach (often also called "global"), elevation change gradients are estimated for all glaciers in an entire region, instead of estimating one by one.
-This is advantageous in respect to areas where voids are frequent, as not even a single dDEM value has to exist on a glacier in order to reconstruct it.
-Of course, the accuracy of such an averaging is much lower than if the local hypsometric approach is used (assuming it is possible).
 
 ```{code-cell} ipython3
+ddem.set_mask(mask)
 ddem_reghyps = ddem.interpolate(method="regional_hypsometric", reference_elevation=dem_2009, mask=glaciers_1990)
 ```
 
@@ -220,7 +219,7 @@ ddem_reghyps = ddem.copy(new_array=ddem_reghyps)
 f, ax = plt.subplots(1, 2)
 ax[0].set_title("Before regional\nhypsometric\ngap-filling")
 ddem.plot(cmap='RdYlBu', vmin=-10, vmax=10, ax=ax[0])
-ax[1].set_title("After regional\hypsometric\ngap-filling")
+ax[1].set_title("After regional\nhypsometric\ngap-filling")
 ddem_reghyps.plot(cmap='RdYlBu', vmin=-10, vmax=10, ax=ax[1], cbar_title="Elevation differences (m)")
 _ = ax[1].set_yticklabels([])
 ```

doc/source/quick_start.md

@@ -116,7 +116,7 @@ os.remove("dh_error.tif")
 ## More examples
 
 To dive into more illustrated code, explore our gallery of examples that is composed of:
-- An {ref}`examples-basic` section on simpler routines (terrain attributes, pre-defined coregistration and uncertainty pipelines),
+- A {ref}`examples-basic` section on simpler routines (terrain attributes, pre-defined coregistration and uncertainty pipelines),
 - An {ref}`examples-advanced` section using advanced pipelines (for in-depth coregistration and uncertainty analysis).
 
 See also the concatenated list of examples below.

doc/source/robust_estimators.md

@@ -5,7 +5,7 @@
 Elevation data often contain outliers that can be traced back to instrument acquisition or processing artefacts, and which hamper further analysis.
 
 In order to mitigate their effect, the analysis of elevation data can integrate [robust statistics](https://en.wikipedia.org/wiki/Robust_statistics) at different levels:
-- **Robust estimators for the central tendency (e.g., mean) and dispersion (e.g., standard deviation)**, to evaluate DEM quality and converge during {ref}`coregistration`,
+- **Robust estimators for the central tendency and statistical dispersion** used during {ref}`coregistration`, {ref}`biascorr` and {ref}`uncertainty`,
 - **Robust estimators for estimating spatial autocorrelation** applied to error propagation in {ref}`uncertainty`,
 - **Robust optimizers for the fitting of parametric models** during {ref}`coregistration` and {ref}`biascorr`.
 

doc/source/terrain.md

@@ -56,8 +56,12 @@ slope = dem.slope()
 slope = xdem.terrain.slope(dem.data, resolution=dem.res)
 ```
 
-If computational performance is key, xDEM can rely on [RichDEM](https://richdem.readthedocs.io/) by specifying
-`use_richdem=True` for speed-up of specific supported attributes (slope, aspect, curvature).
+:::{admonition} Coming soon
+:class: note
+
+We are working on further optimizing the computational performance of certain terrain attributes using convolution.
+:::
+
 
 ## Summary of supported methods
 
@@ -104,6 +108,10 @@ If computational performance is key, xDEM can rely on [RichDEM](https://richdem.
      - [Taud and Parrot (2005)](https://doi.org/10.4000/geomorphologie.622)
 ```
 
+```{note}
+Only grids with **equal pixel size in X and Y** are currently supported. Transform into such a grid with {func}`xdem.DEM.reproject`.
+```
+
 (slope)=
 ## Slope
 

doc/source/uncertainty.md

@@ -109,7 +109,7 @@ The tables below summarize the characteristics of these methods.
 Frequently, in spatial statistics, a single correlation range is considered ("basic" method below).
 However, elevation data often contains errors with correlation ranges spanning different orders of magnitude.
 For this, [Rolstad et al. (2009)](http://dx.doi.org/10.3189/002214309789470950) and
-[Hugonnet et al. (2022)](http://dx.doi.org/10.1109/JSTARS.2022.3188922) considers
+[Hugonnet et al. (2022)](http://dx.doi.org/10.1109/JSTARS.2022.3188922) consider
 potential multiple ranges of spatial correlation (instead of a single one). In addition, [Hugonnet et al. (2022)](http://dx.doi.org/10.1109/JSTARS.2022.3188922)
 considers potential heteroscedasticity or variable errors (instead of homoscedasticity, or constant errors), also common in elevation data.
 
@@ -373,7 +373,7 @@ neff = xdem.spatialstats.number_effective_samples(area=outline_brom, params_vari
 
 **Step 3: Derive final error**
 
-And can now compute our final random error for the mean elevation change in this area of interest:
+And we can now compute our final random error for the mean elevation change in this area of interest:
 
 ```{code-cell} ipython3
 # Compute the standard error

xdem/ddem.py

@@ -2,7 +2,7 @@
 from __future__ import annotations
 
 import warnings
-from typing import Any
+from typing import Any, Literal
 
 import geoutils as gu
 import numpy as np
@@ -163,7 +163,7 @@ def from_array(
 
     def interpolate(
         self,
-        method: str = "linear",
+        method: Literal["idw", "local_hypsometric", "regional_hypsometric"] = "idw",
         reference_elevation: NDArrayf | np.ma.masked_array[Any, np.dtype[np.floating[Any]]] | xdem.DEM = None,
         mask: NDArrayf | xdem.DEM | gu.Vector = None,
     ) -> NDArrayf | None:
@@ -188,8 +188,8 @@ def interpolate(
                 f" different from 'self' ({self.data.shape})"
             )
 
-        if method == "linear":
-            self.filled_data = xdem.volume.linear_interpolation(self.data)
+        if method == "idw":
+            self.filled_data = xdem.volume.idw_interpolation(self.data)
         elif method == "local_hypsometric":
             assert reference_elevation is not None
             assert mask is not None
@@ -231,7 +231,7 @@ def interpolate(
             diff = abs(np.nanmean(interpolated_ddem - self.data))
             assert diff < 0.01, (diff, self.data.mean())
 
-            self.filled_data = xdem.volume.linear_interpolation(interpolated_ddem)
+            self.filled_data = xdem.volume.idw_interpolation(interpolated_ddem)
 
         elif method == "regional_hypsometric":
             assert reference_elevation is not None

xdem/volume.py

@@ -286,7 +286,7 @@ def calculate_hypsometry_area(
     return output
 
 
-def linear_interpolation(
+def idw_interpolation(
     array: NDArrayf | MArrayf,
     max_search_distance: int = 10,
     extrapolate: bool = False,
@@ -532,7 +532,7 @@ def local_hypsometric_interpolation(
     ddem_difference[idealized_ddem == nodata] = np.nan
 
     # Spatially interpolate the difference between these two products.
-    interpolated_ddem_diff = linear_interpolation(np.where(ddem_mask, np.nan, ddem_difference))
+    interpolated_ddem_diff = idw_interpolation(np.where(ddem_mask, np.nan, ddem_difference))
     interpolated_ddem_diff[np.isnan(interpolated_ddem_diff)] = 0
 
     # Correct the idealized dDEM with the difference to the original dDEM.