Incremental commit on doc

doc/source/api.md

@@ -24,14 +24,25 @@ Below, we only repeat the core attributes and methods of GeoUtils, see
 [the Raster API in GeoUtils](https://geoutils.readthedocs.io/en/latest/api.html#raster) for the full list.
 ```
 
-### Opening a DEM
+### Opening or saving a DEM
 
 ```{eval-rst}
 .. autosummary::
     :toctree: gen_modules/
 
     DEM
     DEM.info
+    DEM.save
+```
+
+### Plotting a DEM
+
+```{eval-rst}
+.. autosummary::
+    :toctree: gen_modules/
+
+    DEM
+    DEM.plot
 ```
 
 ### Create from an array

doc/source/dem_class.md

@@ -24,18 +24,17 @@ A {class}`~xdem.DEM` is a {class}`~geoutils.Raster` with an additional georefere
 It inherits the **four main attributes** of {class}`~geoutils.Raster` which are {attr}`~xdem.DEM.data`,
 {attr}`~xdem.DEM.transform`, {attr}`~xdem.DEM.crs` and {attr}`~xdem.DEM.nodata`.
 
-Many other useful raster attributes, such as {attr}`~xdem.DEM.bounds` and {attr}`~xdem.DEM.res` and raster methods
-such {attr}`~xdem.DEM.reproject` or {attr}`~xdem.DEM.crop` are available through the {class}`~geoutils.Raster` object.
+Many other useful raster attributes and methods are available through the {class}`~geoutils.Raster` object, such as 
+{attr}`~geoutils.Raster.bounds`, {attr}`~geoutils.Raster.res`, {func}`~xdem.DEM.reproject` and {func}`~xdem.DEM.crop` .
 
-```{important}
-Below, we only cover a few core aspects linked to the {class}`~geoutils.Raster` object. For more details, see [GeoUtils' Raster documentation page](https://geoutils.readthedocs.io/en/stable/raster_class.html).
-The complete list of {class}`~geoutils.Raster` attributes and methods in [GeoUtils' API](https://geoutils.readthedocs.io/en/stable/api.html#raster).
+```{tip}
+The complete list of {class}`~geoutils.Raster` attributes and methods can be found in [GeoUtils' API](https://geoutils.readthedocs.io/en/stable/api.html#raster) and more info on rasters on [GeoUtils' Raster documentation page](https://geoutils.readthedocs.io/en/stable/raster_class.html).
 ```
 
 ## Open and save
 
-A {class}`~xdem.DEM` is opened, as for a {class}`~geoutils.Raster`, by instantiating with either a {class}`str`, a {class}`pathlib.Path`, a {class}`rasterio.io.DatasetReader` or a
-{class}`rasterio.io.MemoryFile`.
+A {class}`~xdem.DEM` is opened by instantiating with either a {class}`str`, a {class}`pathlib.Path`, a {class}`rasterio.io.DatasetReader` or a
+{class}`rasterio.io.MemoryFile`, as for a {class}`~geoutils.Raster`.
 
 
 ```{code-cell} ipython3
@@ -54,8 +53,10 @@ Detailed information on the {class}`~xdem.DEM` is printed using {func}`~geoutils
 dem.info(stats=True)
 ```
 
-```{note}
-Calling {class}`~xdem.DEM.info()` with `stats=True` automatically loads the array in-memory, like any other operation calling {attr}`~xdem.DEM.data`.
+```{important}
+The {class}`~xdem.DEM` data array remains implicitly unloaded until {attr}`~xdem.DEM.data` is called. For instance, here, calling {class}`~xdem.DEM.info()` with `stats=True` automatically loads the array in-memory.
+
+The georeferencing metadata ({attr}`~xdem.DEM.transform`, {attr}`~xdem.DEM.crs`, {attr}`~xdem.DEM.nodata`), however, is always loaded. This allows to pass it effortlessly to other objects requiring it for geospatial operations (reproject-match, rasterizing a vector, etc).
 ```
 
 A {class}`~xdem.DEM` is saved to file by calling {func}`~xdem.DEM.save` with a {class}`str` or a {class}`pathlib.Path`.
@@ -69,3 +70,97 @@ dem.save("mydem.tif")
 import os
 os.remove("mydem.tif")
 ```
+
+## Plotting
+
+Plotting a DEM is done using {func}`~xdem.DEM.plot`, and can be done alongside a vector file.
+
+```{code-cell} ipython3
+# Open a vector file of glacier outlines near the DEM
+import geoutils as gu
+fn_glacier_outlines = xdem.examples.get_path("longyearbyen_glacier_outlines")
+vect_gla = gu.Vector(fn_glacier_outlines)
+
+# Plot the DEM and the vector file
+dem.plot(cmap="terrain", cbar_title="Elevation (m)")
+vect_gla.plot(dem)  # We pass the DEM as reference for the plot CRS/extent
+```
+
+## Vertical referencing
+
+The vertical reference of a {class}`~xdem.DEM` is stored in {attr}`~xdem.DEM.vcrs`, and derived either from its 
+{attr}`~xdem.DEM.crs` (if 3D) or assessed from the DEM product name during instantiation.
+
+```{code-cell} ipython3
+# Check vertical CRS of dataset
+dem.vcrs
+```
+
+In this case, the DEM has no defined vertical CRS, which is quite common. To set the vertical CRS manually, 
+use {class}`~xdem.DEM.set_vcrs`. Then, to transform into another vertical CRS, use {class}`~xdem.DEM.to_vcrs`.
+
+```{code-cell} ipython3
+# Define the vertical CRS as the 3D ellipsoid of the 2D CRS
+dem.set_vcrs("Ellipsoid")
+# Transform to the EGM96 geoid
+dem.to_vcrs("EGM96")
+```
+
+```{note}
+For more details on vertical referencing, see the {ref}`vertical-ref` page.
+```
+
+## Terrain attributes
+
+A wide range of terrain attributes can be derived from a {class}`~xdem.DEM`, with several methods and options available, 
+by calling the function corresponding to the attribute name such as {func}`~xdem.DEM.slope`.
+
+```{code-cell} ipython3
+# Derive slope using the Zevenberg and Thorne (1987) method
+slope = dem.slope(method="ZevenbergThorne")
+slope.plot(cmap="Reds", cbar_title="Slope (degrees)")
+```
+
+```{note}
+For the full list of terrain attributes, see the {ref}`terrain-attributes` page.
+```
+
+## Coregistration
+
+3D coregistration is performed with {func}`~xdem.DEM.coregister_3d`, which aligns the 
+{class}`~xdem.DEM` to another DEM using a pipeline defined with a {class}`~xdem.coreg.Coreg` 
+object (defaults to horizontal and vertical shifts).
+
+```{code-cell} ipython3
+# Another DEM to-be-aligned to the first one
+filename_tba = xdem.examples.get_path("longyearbyen_tba_dem")
+dem_tba = xdem.DEM(filename_tba)
+
+# Coregister (horizontal and vertical shifts)
+dem_tba_coreg = dem_tba.coregister_3d(dem)
+
+# Plot the elevation change of the DEM due to coregistration
+dh_tba = dem_tba - dem_tba_coreg.reproject(dem_tba)
+dh_tba.plot(cmap="Spectral", cbar_title="Elevation change due to coreg (m)")
+```
+
+```{note}
+For more details on building coregistration pipelines and accessing metadata, see the {ref}`coregistration` page.
+```
+
+## Uncertainty analysis
+
+Estimation of DEM-related uncertainty can be performed with {func}`~xdem.DEM.estimate_uncertainty`, which estimates both 
+**an error map** considering variable per-pixel errors and **the spatial correlation of errors**. The workflow relies 
+on an independent elevation data object that is **assumed to have either much finer or similar precision**, and uses 
+stable terrain as a proxy.
+
+```{code-cell} ipython3
+# Estimate elevation uncertainty assuming both DEMs have similar precision
+sig_dem, rho_sig = dem.estimate_uncertainty(dem_tba_coreg, precision_of_other="same")
+
+# The error map variability is estimated from slope and curvature by default
+sig_dem.plot(cmap="Purples", cbar_title=r"Error in elevation (1$\sigma$, m)")
+
+# The spatial correlation function represents how much errors are correlated at a certain distance
+rho_sig(1000)  # Correlation at 1 km

doc/source/index.md

@@ -104,7 +104,7 @@ terrain
 coregistration
 biascorr
 comparison
-spatialstats
+uncertainty
 ```
 
 ```{toctree}

doc/source/spatialstats.md → doc/source/uncertainty.md

@@ -10,21 +10,19 @@ kernelspec:
   language: python
   name: xdem
 ---
-(spatialstats)=
+(uncertainty)=
 
-# Spatial statistics
+# Uncertainty analysis
 
-Spatial statistics, also referred to as [geostatistics](https://en.wikipedia.org/wiki/Geostatistics), are essential
-for the analysis of observations distributed in space.
-To analyze DEMs, xDEM integrates spatial statistics tools specific to DEMs described in recent literature,
+To analyze DEMs, xDEM integrates spatial uncertainty analysis tools from the recent DEM literature,
 in particular in [Hugonnet et al. (2022)](https://doi.org/10.1109/jstars.2022.3188922) and
 [Rolstad et al. (2009)](https://doi.org/10.3189/002214309789470950). The implementation of these methods relies
-partially on the package [scikit-gstat](https://mmaelicke.github.io/scikit-gstat/index.html).
+partially on the package [scikit-gstat](https://mmaelicke.github.io/scikit-gstat/index.html) for spatial statistics.
 
-The spatial statistics tools can be used to assess the precision of DEMs (see the definition of precision in {ref}`intro`).
-In particular, these tools help to:
+The uncertainty analysis tools can be used to assess the precision of DEMs (see the definition of precision in {ref}`intro`).
+In particular, they help to:
 
-> - account for elevation heteroscedasticity (e.g., varying precision with terrain slope),
+> - account for elevation heteroscedasticity (e.g., varying precision such as with terrain slope or stereo-correlation),
 > - quantify the spatial correlation of errors in DEMs (e.g., native spatial resolution, instrument noise),
 > - estimate robust errors for observations analyzed in space (e.g., average or sum of elevation, or of elevation changes),
 > - propagate errors between spatial ensembles at different scales (e.g., sum of glacier volume changes).
@@ -35,7 +33,10 @@ In particular, these tools help to:
 
 ### Assumptions for statistical inference in spatial statistics
 
-Spatial statistics are valid if the variable of interest verifies [the assumption of second-order stationarity](https://www.aspexit.com/en/fundamental-assumptions-of-the-variogram-second-order-stationarity-intrinsic-stationarity-what-is-this-all-about/).
+
+Spatial statistics, also referred to as [geostatistics](https://en.wikipedia.org/wiki/Geostatistics), are essential
+for the analysis of observations distributed in space. Spatial statistics are valid if the variable of interest 
+verifies [the assumption of second-order stationarity](https://www.aspexit.com/en/fundamental-assumptions-of-the-variogram-second-order-stationarity-intrinsic-stationarity-what-is-this-all-about/).
 That is, if the three following assumptions are verified:
 
 > 1. The mean of the variable of interest is stationary in space, i.e. constant over sufficiently large areas,