Fixes on warnings, wrong edits when accounting for comments

docs/source/intro_accuracy_precision.rst

@@ -22,21 +22,34 @@ calculations consistent, reproducible, and easy.
 Accuracy and precision
 ----------------------
 
-Both `accuracy and precision <https://en.wikipedia.org/wiki/Accuracy_and_precision>`_ are important factors to account
-for when analyzing DEMs:
+`Accuracy and precision <https://en.wikipedia.org/wiki/Accuracy_and_precision>`_ describe random and systematic errors,
+respectively.
 
-- the **accuracy** (systematic error) of a DEM describes how close a DEM is to the true location of measured elevations on the Earth's surface,
-- the **precision** (random error) of a DEM describes the typical spread of its error in measurement, independently of a possible bias from the true positioning.
-
-*Note: sometimes "accuracy" is also used to describe both types of errors, and "trueness" systematic errors, following
-`ISO 5725-1 <https://www.iso.org/obp/ui/#iso:std:iso:5725:-1:ed-1:v1:en>`_. Here, we used accuracy for systematic errors
-as, to our knowledge, it is more common in remote sensing applications.*
+*Note: sometimes "accuracy" is also used to describe both types of errors, and "trueness" systematic errors, as defined
+in* `ISO 5725-1 <https://www.iso.org/obp/ui/#iso:std:iso:5725:-1:ed-1:v1:en>`_ *. Here, we used accuracy for systematic
+errors as, to our knowledge, it is a more commonly used terminology in remote sensing applications.*
 
 .. figure:: imgs/precision_accuracy.png
     :width: 80%
 
     Source: `antarcticglaciers.org <http://www.antarcticglaciers.org/glacial-geology/dating-glacial-sediments2/precision-and-accuracy-glacial-geology/>`_, accessed 29.06.21.
 
+
+For DEMs, we thus have:
+
+- **DEM accuracy** (systematic error) describes how close a DEM is to the true location of measured elevations on the Earth's surface,
+- **DEM precision** (random error) of a DEM describes the typical spread of its error in measurement, independently of a possible bias from the true positioning.
+
+The spatial structure of DEMs complexifies the notion of accuracy and precision, however. Spatially structured
+systematic errors are often related to the gridded nature of DEMs, creating **affine biases** while other, **specific
+biases** exist at the pixel scale. For random errors, a variability in error magnitude or **heteroscedasticity** exists
+across the DEM, while spatially structured patterns of errors are linked to **spatial correlations**.
+
+.. figure:: https://github.com/rhugonnet/dem_error_study/blob/main/figures/fig_2.png?raw=true
+    :width: 100%
+
+    Source: `Hugonnet et al. (2022) <https://doi.org/10.1109/jstars.2022.3188922>`_.
+
 Absolute or relative accuracy
 -----------------------------
 
@@ -50,10 +63,6 @@ TODO: Add another little schematic!
 Optimizing DEM absolute accuracy
 --------------------------------
 
-.. figure:: https://github.com/rhugonnet/dem_error_study/blob/main/figures/fig_2.png?raw=true
-    :width: 100%
-
-    Source: `Hugonnet et al. (2022) <https://doi.org/10.1109/jstars.2022.3188922>`_.
 
 
 Shifts due to poor absolute accuracy are common in elevation datasets, and can be easily corrected by performing a DEM

examples/plot_heterosc_estimation_modelling.py

@@ -177,9 +177,10 @@
 # We also identify that, steep slopes (> 40°) only correspond to high curvature, while the opposite is not true, hence
 # the importance of mapping the variability in two dimensions.
 #
-# Now we need to account for the non-stationarities identified. For this, the simplest approach is a numerical
-# approximation i.e. a piecewise linear interpolation/extrapolation based on the binning results.
-# To ensure that only robust statistic values are used in the interpolation, we set a ``min_count`` value at 30 samples.
+# Now we need to account for the heteroscedasticity identified. For this, the simplest approach is a numerical
+# approximation i.e. a piecewise linear interpolation/extrapolation based on the binning results available through
+# the function :func:`xdem.spatialstats.interp_nd_binning`. To ensure that only robust statistic values are used
+# in the interpolation, we set a ``min_count`` value at 30 samples.
 
 unscaled_dh_err_fun = xdem.spatialstats.interp_nd_binning(df, list_var_names=['slope', 'maxc'],
                                                            statistic='nmad', min_count=30)
@@ -197,7 +198,7 @@
 
 # %%
 # Thus, we rescale the function to exactly match the spread on stable terrain using the
-# ``xdem.spatialstats.two_step_standardization`` function, and get our final error function.
+# :func:`xdem.spatialstats.two_step_standardization` function, and get our final error function.
 
 zscores, dh_err_fun = xdem.spatialstats.two_step_standardization(dh_arr, list_var=[slope_arr, maxc_arr],
                                                            unscaled_error_fun=unscaled_dh_err_fun)
@@ -209,23 +210,6 @@
 # %%
 # This function can be used to estimate the spatial distribution of the elevation error on the extent of our DEMs:
 maxc = np.maximum(np.abs(profc), np.abs(planc))
-errors = dh_err_fun((slope, maxc))
-
-plt.figure(figsize=(8, 5))
-plt_extent = [
-    ref_dem.bounds.left,
-    ref_dem.bounds.right,
-    ref_dem.bounds.bottom,
-    ref_dem.bounds.top,
-]
-plt.imshow(errors.squeeze(), cmap="Reds", vmin=2, vmax=8, extent=plt_extent)
-cbar = plt.colorbar()
-cbar.set_label('Elevation error ($1\sigma$, m)')
-plt.show()
-
-# %%
-# These 3 steps can be done in one go using the ``xdem.spatialstats.estimate_model_heteroscedasticity`` function, which
-# wraps those three funtions, expecting ``np.ndarray`` inputs subset to stable terrain.
-
-df, dh_err_fun = xdem.spatialstats.estimate_model_heteroscedasticity(dvalues=dh_arr, list_var=[slope_arr, maxc_arr],
-                                                                     list_var_names=['slope', 'maxc'])
+errors = dh.copy(new_array= dh_err_fun((slope, maxc)))
+
+errors.show(cmap='Reds', vmin=2, vmax=8, cb_title='Elevation error ($1\sigma$, m)')

examples/plot_variogram_estimation_modelling.py

@@ -33,12 +33,12 @@
 
 # %%
 # We exclude values on glacier terrain in order to isolate stable terrain, our proxy for elevation errors.
-dh.data[mask_glacier] = np.nan
+dh.set_mask(mask_glacier)
 
 # %%
 # We estimate the average per-pixel elevation error on stable terrain, using both the standard deviation
 # and normalized median absolute deviation. For this example, we do not account for elevation heteroscedasticity.
-print('STD: {:.2f} meters.'.format(np.nanstd(dh.data)))
+print('STD: {:.2f} meters.'.format(np.ma.std(dh.data)))
 print('NMAD: {:.2f} meters.'.format(xdem.spatialstats.nmad(dh.data)))
 
 # %%

tests/test_spatialstats.py

@@ -517,7 +517,7 @@ def test_number_effective_samples(self):
     def test_spatial_error_propagation(self):
         """Test that the spatial error propagation wrapper function runs properly"""
 
-        ref, diff, vector_glacier = load_ref_and_diff()
+        ref, diff, _ , vector_glacier = load_ref_and_diff()
 
         # Get the error map and variogram model with standardization
         slope, maxc = xdem.terrain.get_terrain_attribute(ref, attribute=['slope', 'maximum_curvature'])

xdem/spatialstats.py

@@ -38,7 +38,7 @@ def nmad(data: np.ndarray, nfact: float = 1.4826) -> float:
     Calculate the normalized median absolute deviation (NMAD) of an array.
     Default scaling factor is 1.4826 to scale the median absolute deviation (MAD) to the dispersion of a normal
     distribution (see https://en.wikipedia.org/wiki/Median_absolute_deviation#Relation_to_standard_deviation, and
-    e.g. http://dx.doi.org/10.1016/j.isprsjprs.2009.02.003)
+    e.g. Höhle and Höhle (2009), http://dx.doi.org/10.1016/j.isprsjprs.2009.02.003)
 
     :param data: Input data
     :param nfact: Normalization factor for the data