Spatial and robust stats, intro documentation + Update of `spatialsta…

…ts.py` with new `scikit-gstat` features (#159)

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit on spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* accouting for pr comments

* add accuracy/precision image

* incremental commit for spatialstats documentation

* try fix for image display

* add doc image locally

* try static dir in conf

* add source link

* larger size image

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* incremental commit for spatialstats documentation

* testing image width as percentage

* fix testing image width as percentage

* add plot 1d and 2d functions for binning

* remove remaining download example lines

* first draft example nonstationarity

* playing with sphinx gallery

* nonstationarity gallery example

* fix warnings and improve prints

* fix curv unit

* fix min_count

* fix slope variable name

* shorten title

* reformulation

* text fixes

* fix text

* fix text

* use subsample_raster, simplify with skgstat count property

* incremental commit for text

* updating variogram sampling with skgstat custom binning

* incremental commit for adapting to new skgstat

* pr comments + refactor name

* refactor spstats into spatialstats

* incremental commit on documentation

* incremental commit on documentation

* incremental commit on documentation

* incremental commit on documentation

* incremental commit on documentation

* incremental commit on documentation

* fix spatial variogram function with new metricspace

* deal with nodata and simplify quantile binning

* reorganize variogram sampling with skgstat update, pass kwargs down, pass random states for robust testing

* more robust variogram tests, fix small issues

* other merge commit left behind

* reduce test length

* improve plot_vgm

* finish skeleton of gallery example

* remove old plotting

* fix test

* improve plot_vgm and patches_method

* finalize plot_vgm gallery example

* refine plot_vgm example

* fix patches_method test

* update spatial stats doc page with new functions

* improve modularity of patches_method, random_state, add more tests

* fixes with amaurys comments

* fix histogram bug with log scale for plot_vgm

* polish text

* refine plot_vgm xmin

* refactor nruns into n_variograms and nproc into n_jobs (like scikit)

* fix gallery example display

* Fix spatialstats code with pytest.

* add random state, polish fit_sum_variogram

* pass down child random states to get independent variogram runs with fixed parent random state

* change nuth and kaab reference to documentation to avoid sphinx gallery referencing

* streamline gallery examples

* refactor nmax argument of patches method into n_patches

* incremental commit on documentation

* incremental commit on documentation

* incremental commit on documentation

* incremental commit on documentation

* incremental commit on documentation

* incremental commit on documentation

* incremental commit on documentation

* refactor variogram function names

* add subsections

* incremental commit on documentation

* add plot to illustrate non stationarity

* fix doc tests

* draft standardization gallery example

* finalize standardization gallery example

* add standardizing doc plot

* incremental commit on documentation

* fix test_docs

* incremental commit on documentation

* raise warning and return nan dataframe when no valid patch is found

* nan is clearly than nodata

* accelerate patch sampling for large nan matrices

* update minimal skgstat version

* try previous skgstat version for CI

* eriks comments on gallery examples

Co-authored-by: Erik Mannerfelt <[email protected]>

dev-environment.yml

@@ -14,7 +14,7 @@ dependencies:
   - tqdm
   - scikit-image
   - proj-data
-  - scikit-gstat
+  - scikit-gstat>=0.6.7
   - pytransform3d
   - geoutils
 

docs/source/api.rst

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

docs/source/biascorr.rst

@@ -0,0 +1,16 @@
+.. _biascorr:
+
+Bias corrections
+================
+
+Bias corrections correspond to transformations that cannot be described as a 3-dimensional affine function (see :ref:`coregistration`).
+
+Directional biases
+------------------
+
+TODO
+
+Terrain biases
+--------------
+
+TODO

docs/source/code/spatialstats.py

@@ -0,0 +1,41 @@
+"""Code example for spatial statistics"""
+import xdem
+import geoutils as gu
+import numpy as np
+
+
+# Load data
+dh = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_ddem"))
+ref_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
+glacier_mask = gu.geovector.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
+mask = glacier_mask.create_mask(dh)
+slope = xdem.terrain.get_terrain_attribute(ref_dem.data, resolution=ref_dem.res[0], attribute=['slope'])
+
+# Keep only stable terrain data
+dh.data[mask] = np.nan
+
+# Estimate the measurement error by bin of slope, using the NMAD as robust estimator
+df_ns = xdem.spatialstats.nd_binning(dh.data.ravel(), list_var=[slope.ravel()], list_var_names=['slope'],
+                                     statistics=['count', xdem.spatialstats.nmad])
+
+# Derive a numerical function of the measurement error
+err_dh = xdem.spatialstats.interp_nd_binning(df_ns, list_var_names=['slope'])
+
+# Standardize the data
+z_dh = dh.data.ravel() / err_dh(slope.ravel())
+
+# Sample empirical variogram
+df_vgm = xdem.spatialstats.sample_empirical_variogram(values=dh.data, gsd=dh.res[0], subsample=50,
+                                                      random_state=42, runs=10)
+# Fit sum of double-range spherical model
+fun, coefs = xdem.spatialstats.fit_sum_model_variogram(list_model=['Sph', 'Sph'], empirical_variogram=df_vgm)
+
+# Calculate the area-averaged uncertainty with these models
+list_vgm = []
+for i in range(2):
+    list_vgm.append((coefs[2 * i], "Sph", coefs[2 * i + 1]))
+area = 1000
+neff = xdem.spatialstats.neff_circ(area, list_vgm)
+
+
+

docs/source/code/spatialstats_nonstationarity_slope.py

@@ -0,0 +1,22 @@
+"""Code example for spatial statistics"""
+import xdem
+import geoutils as gu
+import numpy as np
+
+# Load data
+dh = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_ddem"))
+ref_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
+glacier_mask = gu.geovector.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
+mask = glacier_mask.create_mask(dh)
+
+# Get slope for non-stationarity
+slope = xdem.terrain.get_terrain_attribute(dem=ref_dem.data, resolution=dh.res, attribute=['slope'])
+
+# Keep only stable terrain data
+dh.data[mask] = np.nan
+
+# Estimate the measurement error by bin of slope, using the NMAD as robust estimator
+df_ns = xdem.spatialstats.nd_binning(dh.data.ravel(), list_var=[slope.ravel()], list_var_names=['slope'],
+                                     statistics=['count', xdem.spatialstats.nmad], list_var_bins=30)
+
+xdem.spatialstats.plot_1d_binning(df_ns, 'slope', 'nmad', 'Slope (degrees)', 'Elevation measurement error (m)')

docs/source/code/spatialstats_standardizing.py

@@ -0,0 +1,44 @@
+"""Documentation plot illustrating standardization of a distribution"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Example x vector
+mu = 15
+sig = 5
+np.random.seed(42)
+y = np.random.normal(mu, sig, size=300)
+
+fig, ax1 = plt.subplots(figsize=(8,3))
+
+# Original histogram
+ax1.hist(y, color='tab:blue', edgecolor='white', linewidth=0.5, alpha=0.7)
+ax1.vlines(mu, ymin=0, ymax=90, color='tab:blue', linestyle='dashed', lw=2)
+ax1.vlines([mu-2*sig, mu+2*sig], ymin=0, ymax=90, colors=['tab:blue','tab:blue'], linestyles='dotted', lw=2)
+ax1.annotate('Original\ndata $x$\n$\\mu_{x} = 15$\n$\\sigma_{x} = 5$', xy=(mu+0.5, 85), xytext=(mu+5, 110), arrowprops=dict(color='tab:blue', width=0.5, headwidth=8),
+             color='tab:blue', fontweight='bold', ha='left')
+ax1.spines['right'].set_visible(False)
+ax1.spines['top'].set_visible(False)
+ax1.spines['left'].set_visible(False)
+ax1.set_yticks([])
+ax1.set_ylim((0,130))
+
+# Standardized histogram
+ax1.hist((y-mu)/sig, color='tab:olive', edgecolor='white', linewidth=0.5, alpha=0.7)
+ax1.vlines(0, ymin=0, ymax=90, color='tab:olive', linestyle='dashed', lw=2)
+ax1.vlines([-2, 2], ymin=0, ymax=90, colors=['tab:olive','tab:olive'], linestyles='dotted', lw=2)
+ax1.annotate('Standardized\ndata $z$\n$\\mu_{z} = 0$\n$\\sigma_{z} = 1$', xy=(-0.3, 85), xytext=(-5, 110), arrowprops=dict(color='tab:olive', width=0.5, headwidth=8),
+             color='tab:olive', fontweight='bold', ha='left')
+ax1.spines['right'].set_visible(False)
+ax1.spines['top'].set_visible(False)
+ax1.spines['left'].set_visible(False)
+ax1.set_yticks([])
+ax1.set_ylim((0,130))
+
+ax1.annotate('', xy=(0, 65), xytext=(mu, 65), arrowprops=dict(arrowstyle="-|>", connectionstyle="arc3,rad=0.2", fc="w"),
+             color='black')
+ax1.text(mu/2, 90, 'Standardization:\n$z = \\frac{x - \\mu}{\\sigma}$', color='black', ha='center', fontsize=14,
+         fontweight='bold')
+ax1.plot([], [], color='tab:gray', linestyle='dashed', label='Mean')
+ax1.plot([], [], color='tab:gray', linestyle='dotted', label='Standard\ndeviation (2$\\sigma$)')
+ax1.legend(loc='center right')
+

docs/source/code/spatialstats_stationarity_assumption.py

@@ -0,0 +1,71 @@
+"""Documentation plot illustrating stationarity of mean and variance"""
+import matplotlib.pyplot as plt
+import numpy as np
+import xdem
+
+# Example x vector
+x = np.linspace(0,1,200)
+
+sig = 0.2
+np.random.seed(42)
+y_rand1 = np.random.normal(0, sig, size=len(x))
+y_rand2 = np.random.normal(0, sig, size=len(x))
+y_rand3 = np.random.normal(0, sig, size=len(x))
+
+
+y_mean = np.array([0.5*xval - 0.25 if xval >0.5 else 0.5*(1-xval) - 0.25 for xval in x])
+
+fac_y_std = 0.5 + 2*x
+
+
+fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8,4))
+
+# Stationary mean and variance
+ax1.plot(x, y_rand1, color='tab:blue', linewidth=0.5)
+ax1.hlines(0, xmin=0, xmax=1, color='black', label='Mean', linestyle='dashed')
+ax1.hlines([-2*sig, 2*sig], xmin=0, xmax=1, colors=['tab:gray','tab:gray'], label='Standard deviation',linestyles='dashed')
+ax1.set_xlim((0,1))
+ax1.set_title('Stationary mean\nStationary variance')
+# ax1.legend()
+ax1.spines['right'].set_visible(False)
+ax1.spines['top'].set_visible(False)
+ax1.set_ylim((-1,1))
+ax1.set_xticks([])
+ax1.set_yticks([])
+ax1.plot(1, 0, ">k", transform=ax1.transAxes, clip_on=False)
+ax1.plot(0, 1, "^k", transform=ax1.transAxes, clip_on=False)
+
+# Non-stationary mean and stationary variance
+ax2.plot(x, y_rand2 + y_mean, color='tab:olive', linewidth=0.5)
+ax2.plot(x, y_mean, color='black', label='Mean', linestyle='dashed')
+ax2.plot(x, y_mean + 2*sig, color='tab:gray', label='Dispersion (2$\\sigma$)', linestyle='dashed')
+ax2.plot(x, y_mean - 2*sig, color='tab:gray', linestyle='dashed')
+ax2.set_xlim((0,1))
+ax2.set_title('Non-stationary mean\nStationary variance')
+ax2.legend(loc='lower center')
+ax2.spines['right'].set_visible(False)
+ax2.spines['top'].set_visible(False)
+ax2.set_xticks([])
+ax2.set_yticks([])
+ax2.set_ylim((-1,1))
+ax2.plot(1, 0, ">k", transform=ax2.transAxes, clip_on=False)
+ax2.plot(0, 1, "^k", transform=ax2.transAxes, clip_on=False)
+
+# Stationary mean and non-stationary variance
+ax3.plot(x, y_rand3 * fac_y_std, color='tab:orange', linewidth=0.5)
+ax3.hlines(0, xmin=0, xmax=1, color='black', label='Mean', linestyle='dashed')
+ax3.plot(x, 2*sig*fac_y_std, color='tab:gray', linestyle='dashed')
+ax3.plot(x, -2*sig*fac_y_std, color='tab:gray', linestyle='dashed')
+ax3.set_xlim((0,1))
+ax3.set_title('Stationary mean\nNon-stationary variance')
+# ax1.legend()
+ax3.spines['right'].set_visible(False)
+ax3.spines['top'].set_visible(False)
+ax3.set_xticks([])
+ax3.set_yticks([])
+ax3.set_ylim((-1,1))
+ax3.plot(1, 0, ">k", transform=ax3.transAxes, clip_on=False)
+ax3.plot(0, 1, "^k", transform=ax3.transAxes, clip_on=False)
+
+plt.tight_layout()
+plt.show()

docs/source/code/spatialstats_variogram_covariance.py

@@ -0,0 +1,26 @@
+"""Documentation plot illustrating the link between variogram and covariance"""
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Example variogram function
+def variogram_exp(h):
+    from xdem.spatialstats import vgm
+    val = vgm(h, 15, model='Exp', psill=10)
+    return val
+
+fig, ax = plt.subplots()
+x = np.linspace(0,100,100)
+ax.plot(x, variogram_exp(x), color='tab:blue', linewidth=2)
+ax.plot(x, 10 - variogram_exp(x), color='black', linewidth=2)
+ax.hlines(10, xmin=0, xmax=100, linestyles='dashed', colors='tab:red')
+ax.text(75, variogram_exp(75)-1, 'Semi-variogram $\\gamma(l)$', ha='center', va='top', color='tab:blue')
+ax.text(75, 10 - variogram_exp(75) + 1, 'Covariance $C(l) = \\sigma^{2} - \\gamma(l)$', ha='center', va='bottom', color='black')
+ax.text(75, 11, 'Variance $\\sigma^{2}$', ha='center', va='bottom', color='tab:red')
+ax.set_xlim((0, 100))
+ax.set_ylim((0, 12))
+ax.set_xlabel('Spatial lag $l$')
+ax.set_ylabel('Variance of elevation differences (m²)')
+ax.spines['right'].set_visible(False)
+ax.spines['top'].set_visible(False)
+plt.tight_layout()
+plt.show()

docs/source/comparison.rst

@@ -1,4 +1,4 @@
-DEM subtraction and volume change
+Differencing and volume change
 =================================
 
 .. contents:: Contents 

docs/source/conf.py

@@ -70,7 +70,7 @@
     },
      # directory where function/class granular galleries are stored
     "backreferences_dir"  : "gen_modules/backreferences",
-    "doc_module": ("xdem", "geoutils")  # I honestly don't know what this is.
+    "doc_module": ("xdem", "geoutils")  # which function/class levels are used to create galleries
 }
 
 # Add any paths that contain templates here, relative to this directory.
@@ -87,7 +87,7 @@
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
-# html_static_path = ['_static']  # Commented out as we have no custom static data
+html_static_path = ['imgs']  # Commented out as we have no custom static data
 
 
 exclude_patterns = [

docs/source/coregistration.rst

@@ -1,7 +1,16 @@
 .. _coregistration:
 
-DEM Coregistration
-==================
+Coregistration
+===============
+
+Coregistration between DEMs correspond to aligning the digital elevation models in three dimension.
+
+Transformations that can be described by a 3-dimensional `affine <https://en.wikipedia.org/wiki/Affine_transformation>`_ function are included in coregistration methods.
+Those transformations include for instance:
+
+- vertical and horizontal translations,
+- rotations, reflections,
+- scalings.
 
 .. contents:: Contents 
    :local:
@@ -208,11 +217,10 @@ For larger rotations, ICP is the only reliable approach (but does not outperform
 
         coreg.ICP() + coreg.NuthKaab()
 
-
 For large biases, rotations and high amounts of noise:
 
 .. code-block:: python
 
-        coreg.BiasCorr() + coreg.ICP() + coreg.NuthKaab()
+        coreg.VerticalShift() + coreg.ICP() + coreg.NuthKaab()
         
 

docs/source/filters.rst

@@ -0,0 +1,6 @@
+.. _filters:
+
+Filtering
+=========
+
+In construction

docs/source/index.rst

@@ -5,7 +5,7 @@
 
 Welcome to xdem's documentation!
 ================================
-xdem aims to make Digital Elevation Model (DEM) comparisons easy.
+``xdem`` aims to make Digital Elevation Model (DEM) analysis easy.
 Coregistration, subtraction (and volume measurements), and error statistics should be available to anyone with the correct input data.
 
 
@@ -27,9 +27,13 @@ Simple usage
    :maxdepth: 2
    :caption: Contents:
 
+   intro
    coregistration
+   biascorr
+   filters
    comparison
-   spatial_stats
+   spatialstats
+   robuststats
    terrain
    auto_examples/index.rst
    api.rst