diff --git a/doc/source/api.md b/doc/source/api.md index 9057234f..cf4df17d 100644 --- a/doc/source/api.md +++ b/doc/source/api.md @@ -11,14 +11,20 @@ documentation. .. currentmodule:: xdem ``` -## DEM - ```{eval-rst} .. minigallery:: xdem.DEM :add-heading: ``` -### Opening a DEM +## DEM + +```{important} +A {class}`~xdem.DEM` inherits all raster methods and attributes from the {class}`~geoutils.Raster` object of GeoUtils. +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 or saving a DEM ```{eval-rst} .. autosummary:: @@ -26,6 +32,17 @@ documentation. DEM DEM.info + DEM.save +``` + +### Plotting a DEM + +```{eval-rst} +.. autosummary:: + :toctree: gen_modules/ + + DEM + DEM.plot ``` ### Create from an array @@ -41,10 +58,20 @@ documentation. ### Unique attributes -```{note} -A {class}`~xdem.DEM` inherits four unique attributes from {class}`~geoutils.Raster`, see [the dedicated section of GeoUtils' API](https://geoutils.readthedocs.io/en/latest/api.html#unique-attributes). +#### Inherited from {class}`~geoutils.Raster` + +```{eval-rst} +.. autosummary:: + :toctree: gen_modules/ + + DEM.data + DEM.crs + DEM.transform + DEM.nodata ``` +#### Specific to {class}`~xdem.DEM` + ```{eval-rst} .. autosummary:: :toctree: gen_modules/ @@ -52,16 +79,19 @@ A {class}`~xdem.DEM` inherits four unique attributes from {class}`~geoutils.Rast DEM.vcrs ``` -### Derived attributes +### Georeferencing + +#### Inherited from {class}`~geoutils.Raster` ```{eval-rst} .. autosummary:: :toctree: gen_modules/ - DEM.ccrs + DEM.reproject + DEM.crop ``` -### Vertical referencing +#### Vertical referencing for {class}`~xdem.DEM` ```{eval-rst} .. autosummary:: @@ -71,6 +101,53 @@ A {class}`~xdem.DEM` inherits four unique attributes from {class}`~geoutils.Rast DEM.to_vcrs ``` +### Raster-vector interface + +```{note} +See the full list of vector methods in [GeoUtils' documentation](https://geoutils.readthedocs.io/en/latest/api.html#vector). +``` + +```{eval-rst} +.. autosummary:: + :toctree: gen_modules/ + + DEM.polygonize + DEM.proximity +``` + +### Coregistration + +```{tip} +To build and pass your coregistration pipeline to {func}`~xdem.DEM.coregister_3d`, see the API of {ref}`api-geo-handle`. +``` + +```{eval-rst} +.. autosummary:: + :toctree: gen_modules/ + + DEM.coregister_3d +``` + +### Terrain attributes + +```{eval-rst} +.. autosummary:: + :toctree: gen_modules/ + + DEM.slope + DEM.aspect + DEM.hillshade + DEM.curvature + DEM.profile_curvature + DEM.planform_curvature + DEM.maximum_curvature + DEM.topographic_position_index + DEM.terrain_ruggedness_index + DEM.roughness + DEM.rugosity + DEM.fractal_roughness +``` + ## dDEM ```{eval-rst} @@ -113,6 +190,25 @@ A {class}`~xdem.DEM` inherits four unique attributes from {class}`~geoutils.Rast xdem.coreg.BlockwiseCoreg ``` +#### Fitting and applying transforms + +```{eval-rst} +.. autosummary:: + :toctree: gen_modules/ + + xdem.coreg.Coreg.fit + xdem.coreg.Coreg.apply +``` + +#### Other functionalities + +```{eval-rst} +.. autosummary:: + :toctree: gen_modules/ + + xdem.coreg.Coreg.residuals +``` + ### Affine coregistration methods diff --git a/doc/source/coregistration.md b/doc/source/coregistration.md index f424b730..5c98d869 100644 --- a/doc/source/coregistration.md +++ b/doc/source/coregistration.md @@ -23,6 +23,14 @@ Those transformations include for instance: - rotations, reflections, - scalings. +## Quick use + +Coregistrations are defined using either a single method or pipeline of {class}`~xdem.coreg.Coreg` methods, that are listed below. + +Performing the coregistration on a pair of DEM is done either by using {func}`xdem.DEM.coregister_3d` from the DEM that will be aligned, or +by specifying the {func}`xdem.coreg.Coreg.fit` and {func}`xdem.coreg.Coreg.apply` steps, which allows array inputs and +to apply the same fitted transformation to several objects (e.g., horizontal shift of both a stereo DEM and its ortho-image). + ## Introduction Coregistration of a DEM is performed when it needs to be compared to a reference, but the DEM does not align with the reference perfectly. diff --git a/doc/source/dem_class.md b/doc/source/dem_class.md new file mode 100644 index 00000000..0e159015 --- /dev/null +++ b/doc/source/dem_class.md @@ -0,0 +1,166 @@ +--- +file_format: mystnb +jupytext: + formats: md:myst + text_representation: + extension: .md + format_name: myst +kernelspec: + display_name: xdem-env + language: python + name: xdem +--- +(dem-class)= + +# The digital elevation model ({class}`~xdem.DEM`) + +Below, a summary of the {class}`~xdem.DEM` object and its methods. + +(dem-obj-def)= + +## Object definition and attributes + +A {class}`~xdem.DEM` is a {class}`~geoutils.Raster` with an additional georeferenced vertical dimension stored in the attribute {attr}`~xdem.DEM.vcrs`. +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 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` . + +```{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 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 +import xdem + +# Instantiate a DEM from a filename on disk +filename_dem = xdem.examples.get_path("longyearbyen_ref_dem") +dem = xdem.DEM(filename_dem) +dem +``` + +Detailed information on the {class}`~xdem.DEM` is printed using {func}`~geoutils.Raster.info`, along with basic statistics using `stats=True`: + +```{code-cell} ipython3 +# Print details of raster +dem.info(stats=True) +``` + +```{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`. + +```{code-cell} ipython3 +# Save raster to disk +dem.save("mydem.tif") +``` +```{code-cell} ipython3 +:tags: [remove-cell] +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/elevation_objects.md b/doc/source/elevation_objects.md new file mode 100644 index 00000000..ff80ae56 --- /dev/null +++ b/doc/source/elevation_objects.md @@ -0,0 +1,12 @@ +(elevation-objects)= +# Elevation data objects + +Elevation data objects of xDEM inherit their characteristics from raster and vector objects of +our sister-package [GeoUtils](https://geoutils.readthedocs.io/en/stable/). + +```{toctree} +:maxdepth: 2 + +dem_class +elevation_point_cloud +``` diff --git a/doc/source/elevation_point_cloud.md b/doc/source/elevation_point_cloud.md new file mode 100644 index 00000000..21b73504 --- /dev/null +++ b/doc/source/elevation_point_cloud.md @@ -0,0 +1,17 @@ +--- +file_format: mystnb +jupytext: + formats: md:myst + text_representation: + extension: .md + format_name: myst +kernelspec: + display_name: xdem-env + language: python + name: xdem +--- +(elevation-point-cloud)= + +# The elevation point cloud ({class}`~xdem.EPC`) + +In construction, planned for 2024. diff --git a/doc/source/index.md b/doc/source/index.md index 0c4c0a83..c74caf4f 100644 --- a/doc/source/index.md +++ b/doc/source/index.md @@ -98,12 +98,13 @@ intro_accuracy_precision :caption: Features :maxdepth: 2 +elevation_objects vertical_ref terrain coregistration biascorr comparison -spatialstats +uncertainty ``` ```{toctree} diff --git a/doc/source/terrain.md b/doc/source/terrain.md index 50391759..3d6d16f7 100644 --- a/doc/source/terrain.md +++ b/doc/source/terrain.md @@ -2,12 +2,21 @@ # Terrain attributes -For analytic and visual purposes, deriving certain attributes of a DEM may be required. -Some are useful for direct analysis, such as a slope map to differentiate features of different angles, while others, like the hillshade, are great tools for visualizing a DEM. +xDEM can derive a wide range of **terrain attributes** from a DEM. + +Attributes are derived in pure Python for modularity (e.g., varying window size) and other uses (e.g., uncertainty), +and tested for consistency against [gdaldem](https://gdal.org/programs/gdaldem.html) and [RichDEM](https://richdem.readthedocs.io/). + +## Quick use + +Terrain attribute methods can either be called directly from a {class}`~xdem.DEM` (e.g., {func}`xdem.DEM.slope`) or +through the {class}`~xdem.terrain` module which allows array input. If computational performance +is key, xDEM can rely on [RichDEM](https://richdem.readthedocs.io/) by specifying `use_richdem=True` for speed-up +of its supported attributes (slope, aspect, curvature). ## Slope -{func}`xdem.terrain.slope` +{func}`xdem.DEM.slope` The slope of a DEM describes the tilt, or gradient, of each pixel in relation to its neighbours. It is most often described in degrees, where a flat surface is 0° and a vertical cliff is 90°. @@ -29,7 +38,7 @@ example. ## Aspect -{func}`xdem.terrain.aspect` +{func}`xdem.DEM.aspect` The aspect describes the orientation of strongest slope. It is often reported in degrees, where a slope tilting straight north corresponds to an aspect of 0°, and an eastern @@ -48,7 +57,7 @@ As the aspect is directly based on the slope, it varies between the method of [H ## Hillshade -{func}`xdem.terrain.hillshade` +{func}`xdem.DEM.hillshade` The hillshade is a slope map, shaded by the aspect of the slope. The slope map is a good tool to visualize terrain, but it does not distinguish between a mountain and a valley. @@ -73,7 +82,7 @@ It therefore has little analytic purpose, but it still constitutes a great visua ## Curvature -{func}`xdem.terrain.curvature` +{func}`xdem.DEM.curvature` The curvature map is the second derivative of elevation, which highlights the convexity or concavity of the terrain. If a surface is convex (like a mountain peak), it will have positive curvature. @@ -93,7 +102,7 @@ The curvature is based on the method of [Zevenbergen and Thorne (1987)](http://d ## Planform curvature -{func}`xdem.terrain.planform_curvature` +{func}`xdem.DEM.planform_curvature` The planform curvature is the curvature perpendicular to the direction of slope, reported in 100 m{sup}`-1`. @@ -110,7 +119,7 @@ It is based on the method of [Zevenbergen and Thorne (1987)](http://dx.doi.org/1 ## Profile curvature -{func}`xdem.terrain.profile_curvature` +{func}`xdem.DEM.profile_curvature` The profile curvature is the curvature parallel to the direction of slope, reported in 100 m{sup}`-1`.. @@ -127,7 +136,7 @@ It is based on the method of [Zevenbergen and Thorne (1987)](http://dx.doi.org/1 ## Topographic Position Index -{func}`xdem.terrain.topographic_position_index` +{func}`xdem.DEM.topographic_position_index` The Topographic Position Index (TPI) is a metric of slope position, based on the method of [Weiss (2001)](http://www.jennessent.com/downloads/TPI-poster-TNC_18x22.pdf) that corresponds to the difference of the elevation of a central pixel with the average of that of neighbouring pixels. Its unit is that of the DEM (typically meters) and it can be @@ -144,7 +153,7 @@ computed for any window size (default 3x3 pixels). ## Terrain Ruggedness Index -{func}`xdem.terrain.terrain_ruggedness_index` +{func}`xdem.DEM.terrain_ruggedness_index` The Terrain Ruggedness Index (TRI) is a metric of terrain ruggedness, based on cumulated differences in elevation between a central pixel and its surroundings. Its unit is that of the DEM (typically meters) and it can be computed for any @@ -167,7 +176,7 @@ where the TRI is defined by the mean absolute difference with neighbouring pixel ## Roughness -{func}`xdem.terrain.roughness` +{func}`xdem.DEM.roughness` The roughness is a metric of terrain ruggedness, based on the maximum difference in elevation in the surroundings. The roughness is based on the method of [Dartnell (2000)](http://dx.doi.org/10.14358/PERS.70.9.1081). Its unit is that of the DEM (typically meters) and it can be computed for any window size (default 3x3 pixels). @@ -183,7 +192,7 @@ The roughness is based on the method of [Dartnell (2000)](http://dx.doi.org/10.1 ## Rugosity -{func}`xdem.terrain.rugosity` +{func}`xdem.DEM.rugosity` The rugosity is a metric of terrain ruggedness, based on the ratio between planimetric and real surface area. The rugosity is based on the method of [Jenness (2004)](). @@ -200,7 +209,7 @@ It is unitless, and is only supported for a 3x3 window size. ## Fractal roughness -{func}`xdem.terrain.fractal_roughness` +{func}`xdem.DEM.fractal_roughness` The fractal roughness is a metric of terrain ruggedness based on the local fractal dimension estimated by the volume box-counting method of [Taud and Parrot (2005)](https://doi.org/10.4000/geomorphologie.622). @@ -221,7 +230,7 @@ DEM pixels. Its unit is that of a dimension, and is always between 1 (dimension Often, one may seek more terrain attributes than one, e.g. both the slope and the aspect. Since both are dependent on the gradient of the DEM, calculating them separately is computationally repetitive. -Multiple terrain attributes can be calculated from the same gradient using the {func}`xdem.terrain.get_terrain_attribute` function. +Multiple terrain attributes can be calculated from the same gradient using the {func}`xdem.DEM.get_terrain_attribute` function. ```{eval-rst} .. minigallery:: xdem.terrain.get_terrain_attribute 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, diff --git a/tests/test_dem.py b/tests/test_dem.py index 0e4224cd..2f894603 100644 --- a/tests/test_dem.py +++ b/tests/test_dem.py @@ -6,8 +6,7 @@ import warnings from typing import Any -import geoutils.raster as gr -import geoutils.raster.satimg as si +import geoutils as gu import numpy as np import pytest import rasterio as rio @@ -36,14 +35,14 @@ def test_init(self) -> None: assert isinstance(dem2, DEM) # From Raster - r = gr.Raster(fn_img) + r = gu.Raster(fn_img) dem3 = DEM(r) assert isinstance(dem3, DEM) # From SatelliteImage with warnings.catch_warnings(): warnings.filterwarnings("ignore", "Parse metadata from file not implemented") - img = si.SatelliteImage(fn_img) + img = gu.SatelliteImage(fn_img) dem4 = DEM(img) assert isinstance(dem4, DEM) @@ -51,7 +50,7 @@ def test_init(self) -> None: # Check all attributes attrs = [at for at in _default_rio_attrs if at not in ["name", "dataset_mask", "driver"]] - all_attrs = attrs + si.satimg_attrs + xdem.dem.dem_attrs + all_attrs = attrs + gu.raster.satimg.satimg_attrs + xdem.dem.dem_attrs for attr in all_attrs: attrs_per_dem = [idem.__getattribute__(attr) for idem in list_dem] assert all(at == attrs_per_dem[0] for at in attrs_per_dem) @@ -193,7 +192,7 @@ def test_copy(self) -> None: # using list directly available in Class attrs = [at for at in _default_rio_attrs if at not in ["name", "dataset_mask", "driver"]] - all_attrs = attrs + si.satimg_attrs + xdem.dem.dem_attrs + all_attrs = attrs + gu.raster.satimg.satimg_attrs + xdem.dem.dem_attrs for attr in all_attrs: assert r.__getattribute__(attr) == r2.__getattribute__(attr) @@ -338,3 +337,44 @@ def test_to_vcrs__grids(self, grid_shifts: dict[str, Any]) -> None: # Check the shift is the one expect within 10% assert z_diff == pytest.approx(grid_shifts["shift"], rel=0.1) + + @pytest.mark.parametrize("terrain_attribute", xdem.terrain.available_attributes) # type: ignore + def test_terrain_attributes_wrappers(self, terrain_attribute: str) -> None: + """Check the terrain attributes corresponds to the ones derived in the terrain module.""" + + fn_dem = xdem.examples.get_path("longyearbyen_ref_dem") + dem = DEM(fn_dem) + + dem_class_attr = getattr(dem, terrain_attribute)() + terrain_module_attr = getattr(xdem.terrain, terrain_attribute)(dem) + + assert dem_class_attr.raster_equal(terrain_module_attr) + + def test_coregister_3d_wrapper(self) -> None: + + fn_ref = xdem.examples.get_path("longyearbyen_ref_dem") + fn_tba = xdem.examples.get_path("longyearbyen_tba_dem") + + dem_ref = DEM(fn_ref) + dem_tba = DEM(fn_tba) + + dem_class_aligned = dem_tba.coregister_3d(dem_ref, random_state=42) + + nk = xdem.coreg.NuthKaab() + nk.fit(dem_ref, dem_tba, random_state=42) + coreg_module_aligned = nk.apply(dem_tba) + + assert dem_class_aligned.raster_equal(coreg_module_aligned) + + def test_estimate_uncertainty(self) -> None: + + fn_ref = xdem.examples.get_path("longyearbyen_ref_dem") + fn_tba = xdem.examples.get_path("longyearbyen_tba_dem") + + dem_ref = DEM(fn_ref) + dem_tba = DEM(fn_tba) + + sig_h, corr_sig = dem_tba.estimate_uncertainty(dem_ref) + + assert isinstance(sig_h, gu.Raster) + assert callable(corr_sig) diff --git a/xdem/dem.py b/xdem/dem.py index 70fb5544..b81cb1ff 100644 --- a/xdem/dem.py +++ b/xdem/dem.py @@ -5,14 +5,22 @@ import warnings from typing import Any, Literal +import numpy as np import rasterio as rio from affine import Affine from geoutils import SatelliteImage -from geoutils.raster import RasterType +from geoutils.raster import Mask, RasterType from pyproj import CRS from pyproj.crs import CompoundCRS, VerticalCRS - -from xdem._typing import MArrayf, NDArrayf +from skgstat import Variogram + +from xdem import coreg, terrain +from xdem._typing import MArrayf, NDArrayb, NDArrayf +from xdem.misc import copy_doc +from xdem.spatialstats import ( + infer_heteroscedasticity_from_stable, + infer_spatial_correlation_from_stable, +) from xdem.vcrs import ( _build_ccrs_from_crs_and_vcrs, _grid_from_user_input, @@ -290,3 +298,182 @@ def to_vcrs( # Update vcrs (which will update ccrs if called) self.set_vcrs(new_vcrs=vcrs) + + @copy_doc(terrain, remove_dem_res_params=True) + def slope( + self, + method: str = "Horn", + degrees: bool = True, + use_richdem: bool = False, + ) -> RasterType: + return terrain.slope(self, method=method, degrees=degrees, use_richdem=use_richdem) + + @copy_doc(terrain, remove_dem_res_params=True) + def aspect( + self, + method: str = "Horn", + degrees: bool = True, + use_richdem: bool = False, + ) -> RasterType: + + return terrain.aspect(self, method=method, degrees=degrees, use_richdem=use_richdem) + + @copy_doc(terrain, remove_dem_res_params=True) + def hillshade( + self, + method: str = "Horn", + azimuth: float = 315.0, + altitude: float = 45.0, + z_factor: float = 1.0, + use_richdem: bool = False, + ) -> RasterType: + + return terrain.hillshade( + self, method=method, azimuth=azimuth, altitude=altitude, z_factor=z_factor, use_richdem=use_richdem + ) + + @copy_doc(terrain, remove_dem_res_params=True) + def curvature( + self, + use_richdem: bool = False, + ) -> RasterType: + + return terrain.curvature(self, use_richdem=use_richdem) + + @copy_doc(terrain, remove_dem_res_params=True) + def planform_curvature( + self, + use_richdem: bool = False, + ) -> RasterType: + + return terrain.planform_curvature(self, use_richdem=use_richdem) + + @copy_doc(terrain, remove_dem_res_params=True) + def profile_curvature( + self, + use_richdem: bool = False, + ) -> RasterType: + + return terrain.profile_curvature(self, use_richdem=use_richdem) + + @copy_doc(terrain, remove_dem_res_params=True) + def maximum_curvature( + self, + use_richdem: bool = False, + ) -> RasterType: + + return terrain.maximum_curvature(self, use_richdem=use_richdem) + + @copy_doc(terrain, remove_dem_res_params=True) + def topographic_position_index(self, window_size: int = 3) -> RasterType: + + return terrain.topographic_position_index(self, window_size=window_size) + + @copy_doc(terrain, remove_dem_res_params=True) + def terrain_ruggedness_index(self, method: str = "Riley", window_size: int = 3) -> RasterType: + + return terrain.terrain_ruggedness_index(self, method=method, window_size=window_size) + + @copy_doc(terrain, remove_dem_res_params=True) + def roughness(self, window_size: int = 3) -> RasterType: + + return terrain.roughness(self, window_size=window_size) + + @copy_doc(terrain, remove_dem_res_params=True) + def rugosity(self) -> RasterType: + + return terrain.rugosity(self) + + @copy_doc(terrain, remove_dem_res_params=True) + def fractal_roughness(self, window_size: int = 13) -> RasterType: + + return terrain.fractal_roughness(self, window_size=window_size) + + @copy_doc(terrain, remove_dem_res_params=True) + def get_terrain_attribute(self, attribute: str | list[str], **kwargs: Any) -> RasterType | list[RasterType]: + return terrain.get_terrain_attribute(self, attribute=attribute, **kwargs) + + def coregister_3d( + self, + reference_dem: DEM, + coreg_method: coreg.Coreg = None, + inlier_mask: Mask | NDArrayb = None, + bias_vars: dict[str, NDArrayf | MArrayf | RasterType] = None, + **kwargs: Any, + ) -> DEM: + """ + Coregister DEM to another DEM in three dimensions. + + Any coregistration method or pipeline can be passed, default is only horizontal and vertical shift of + Nuth and Kääb (2011). + + :param reference_dem: Reference DEM the alignment is made towards. + :param coreg_method: Coregistration method or pipeline. + :param inlier_mask: Optional. 2D boolean array or mask of areas to include in the analysis (inliers=True). + :param bias_vars: Optional, only for some bias correction methods. 2D array or rasters of bias variables used. + :param kwargs: Keyword arguments passed to Coreg.fit(). + + :return: Coregistered DEM. + """ + + if coreg_method is None: + coreg_method = coreg.NuthKaab() + + coreg_method.fit( + reference_dem=reference_dem, dem_to_be_aligned=self, inlier_mask=inlier_mask, bias_vars=bias_vars, **kwargs + ) + return coreg_method.apply(self) + + def estimate_uncertainty( + self, + other_dem: DEM, + stable_terrain: Mask | NDArrayb = None, + precision_of_other: Literal["finer"] | Literal["same"] = "finer", + list_vars: tuple[RasterType | str, ...] = ("slope", "maximum_curvature"), + list_vario_models: str | tuple[str, ...] = ("spherical", "spherical"), + ) -> tuple[RasterType, Variogram]: + """ + Estimate uncertainty of DEM. + + Derives a map of per-pixel errors (based on slope and curvature by default) and a function describing the + spatial correlation of error (between 0 and 1) with spatial lag (distance between observations). + + Uses stable terrain as an error proxy and assumes a higher or similar-precision DEM is used as reference, + see Hugonnet et al. (2022). + + :param other_dem: Other DEM to use for estimation, either of finer or similar precision for reliable estimates. + :param stable_terrain: Mask of stable terrain to use as error proxy. + :param precision_of_other: Whether finer precision (3 times more precise = 95% of estimated error will come from + this DEM) or similar precision (for instance another acquisition of the same DEM). + :param list_vars: Variables to use to predict error variability (= elevation heteroscedasticity). Either rasters + or names of a terrain attributes. Defaults to slope and maximum curvature of the DEM. + :param list_vario_models: Variogram forms to model the spatial correlation of error. A list translates into + a sum of models. + + :return: Uncertainty raster, Variogram of uncertainty correlation. + """ + + # Elevation change + dh = other_dem.reproject(self) - self + + # If the precision of the other DEM is the same, divide the dh values by sqrt(2) + # See Equation 7 and 8 of Hugonnet et al. (2022) + if precision_of_other == "same": + dh /= np.sqrt(2) + + # Derive terrain attributes of DEM if string are passed in the list of variables + list_var_rast = [] + for var in list_vars: + if isinstance(var, str): + list_var_rast.append(getattr(terrain, var)(self)) + else: + list_var_rast.append(var) + + # Estimate per-pixel uncertainty + sig_dh = infer_heteroscedasticity_from_stable(dvalues=dh, list_var=list_var_rast, stable_mask=stable_terrain)[0] + + # Estimate spatial correlation + corr_sig = infer_spatial_correlation_from_stable( + dvalues=dh, list_models=list(list_vario_models), stable_mask=stable_terrain, errors=sig_dh + )[2] + return sig_dh, corr_sig diff --git a/xdem/misc.py b/xdem/misc.py index 7b0eb0b8..b36ad838 100644 --- a/xdem/misc.py +++ b/xdem/misc.py @@ -130,6 +130,54 @@ def new_func(*args: Any, **kwargs: Any) -> Any: return deprecator_func +def copy_doc( + module_to_copy: object, + remove_dem_res_params: bool = False, +) -> Callable: # type: ignore + """ + A decorator to copy docstring from a function to another one while replacing the docstring. + Used for copying xdem.terrain documentation to xdem.DEM. + + :param module_to_copy: Name of module to copy the function from + :param remove_dem_res_params: To remove the parameters dem: and resolution: in the terrain docstring, + as they are useless for the DEM class. + """ + + def decorator(decorated: Callable) -> Callable: # type: ignore + # Get name of decorated object + # If object is a property, get name through fget + try: + decorated_name = decorated.fget.__name__ + # Otherwise, directly with the name attribute + except AttributeError: + decorated_name = decorated.__name__ + + # Get parent doc + other_doc = getattr(module_to_copy, decorated_name).__doc__ + + # Replace argument description of dem and resolution (not used in the DEM class, only in terrain) + if remove_dem_res_params: + + # Find and remove them if they exist + if ":param dem:" in other_doc: + dem_section = "\n :param dem:" + other_doc.split("\n :param dem:")[1].split("\n")[0] + other_doc = other_doc.replace(dem_section, "") + if ":param resolution:" in other_doc: + resolution_section = ( + "\n :param resolution:" + other_doc.split("\n :param resolution:")[1].split("\n")[0] + ) + other_doc = other_doc.replace(resolution_section, "") + + # Remove docstring examples + other_doc = other_doc.split(":examples:")[0] + + decorated.__doc__ = other_doc + + return decorated + + return decorator + + def diff_environment_yml( fn_env: str | dict[str, Any], fn_devenv: str | dict[str, Any], print_dep: str = "both", input_dict: bool = False ) -> None: diff --git a/xdem/terrain.py b/xdem/terrain.py index f0ba2b13..fc7a0f14 100644 --- a/xdem/terrain.py +++ b/xdem/terrain.py @@ -18,6 +18,21 @@ except ImportError: _has_rd = False +available_attributes = [ + "slope", + "aspect", + "hillshade", + "curvature", + "planform_curvature", + "profile_curvature", + "maximum_curvature", + "topographic_position_index", + "terrain_ruggedness_index", + "roughness", + "rugosity", + "fractal_roughness", +] + def _raster_to_rda(rst: RasterType) -> rd.rdarray: """ @@ -1134,7 +1149,7 @@ def slope( :param dem: The DEM to generate a slope map for. :param resolution: The X/Y or (X, Y) resolution of the DEM. :param method: Method to calculate slope: "Horn" or "ZevenbergThorne". - :param degrees: Return a slope map in degrees (False means radians). + :param degrees: Whether to use degrees or radians (False means radians). :param use_richdem: Whether to use RichDEM to compute the attribute. :examples: @@ -1192,7 +1207,7 @@ def aspect( :param dem: The DEM to calculate the aspect from. :param method: Method to calculate aspect: "Horn" or "ZevenbergThorne". - :param degrees: Return an aspect map in degrees (if False, returns radians) + :param degrees: Whether to use degrees or radians (False means radians). :param use_richdem: Whether to use RichDEM to compute the attribute. :examples: