diff --git a/doc/source/api.md b/doc/source/api.md index 95de2807..cf4df17d 100644 --- a/doc/source/api.md +++ b/doc/source/api.md @@ -24,7 +24,7 @@ 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:: @@ -32,6 +32,17 @@ Below, we only repeat the core attributes and methods of GeoUtils, see DEM DEM.info + DEM.save +``` + +### Plotting a DEM + +```{eval-rst} +.. autosummary:: + :toctree: gen_modules/ + + DEM + DEM.plot ``` ### Create from an array diff --git a/doc/source/dem_class.md b/doc/source/dem_class.md index 49856b15..0e159015 100644 --- a/doc/source/dem_class.md +++ b/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 \ No newline at end of file diff --git a/doc/source/index.md b/doc/source/index.md index cd6fbc5a..c74caf4f 100644 --- a/doc/source/index.md +++ b/doc/source/index.md @@ -104,7 +104,7 @@ terrain coregistration biascorr comparison -spatialstats +uncertainty ``` ```{toctree} diff --git a/doc/source/spatialstats.md b/doc/source/uncertainty.md similarity index 96% rename from doc/source/spatialstats.md rename to doc/source/uncertainty.md index 497d4c39..d38158bb 100644 --- a/doc/source/spatialstats.md +++ b/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,