Elevation error map

Digital elevation models have a precision that can vary with terrain and instrument-related variables. Here, we rely on a non-stationary spatial statistics framework to estimate and model this variability in elevation error, using terrain slope and maximum curvature as explanatory variables, with stable terrain as an error proxy for moving terrain.

Reference: Hugonnet et al. (2022), Figs. 4 and S6–S9. Equations 7 or 8 can be used to convert elevation change errors into elevation errors.

import geoutils as gu

import xdem

We load a difference of DEMs at Longyearbyen, already coregistered using Nuth and Kääb (2011) as shown in the Nuth and Kääb coregistration example. We also load the reference DEM to derive terrain attributes and the glacier outlines here corresponding to moving terrain.

dh = xdem.DEM(xdem.examples.get_path("longyearbyen_ddem"))
ref_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
glacier_outlines = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))

We derive the terrain slope and maximum curvature from the reference DEM.

slope, maximum_curvature = xdem.terrain.get_terrain_attribute(ref_dem, attribute=["slope", "maximum_curvature"])

Then, we run the pipeline for inference of elevation heteroscedasticity from stable terrain:

errors, df_binning, error_function = xdem.spatialstats.infer_heteroscedasticity_from_stable(
    dvalues=dh, list_var=[slope, maximum_curvature], list_var_names=["slope", "maxc"], unstable_mask=glacier_outlines
)

The first output corresponds to the error map for the DEM (\(\pm\) 1\(\sigma\) level):

errors.plot(vmin=2, vmax=7, cmap="Reds", cbar_title=r"Elevation error (1$\sigma$, m)")

The second output is the dataframe of 2D binning with slope and maximum curvature:

df_binning

	nd	count	nmad	slope	maxc
0	1	277448.0	2.189584	[0.0, 5.77841004000488)	NaN
1	1	191846.0	2.219308	[5.77841004000488, 11.55682008000976)	NaN
2	1	144414.0	2.322037	[11.55682008000976, 17.335230120014643)	NaN
3	1	147537.0	2.554825	[17.335230120014643, 23.11364016001952)	NaN
4	1	168619.0	2.741778	[23.11364016001952, 28.8920502000244)	NaN
...	...	...	...	...	...
95	2	9.0	9.577810	[52.00569036004392, 57.78410040004881)	[0.7980177770456596, 3.4235908346649175)
96	2	14.0	15.795427	[52.00569036004392, 57.78410040004881)	[3.4235908346649175, 6.0491638922841755)
97	2	5.0	6.308063	[52.00569036004392, 57.78410040004881)	[6.0491638922841755, 8.674736949903435)
98	2	0.0	NaN	[52.00569036004392, 57.78410040004881)	[8.674736949903435, 11.300310007522695)
99	2	0.0	NaN	[52.00569036004392, 57.78410040004881)	[11.300310007522695, 13.925883065141957)

120 rows × 5 columns

The third output is the 2D binning interpolant, i.e. an error function with the slope and maximum curvature (Note: below we divide the maximum curvature by 100 to convert it in m^-1 ):

for slope, maxc in [(0, 0), (40, 0), (0, 5), (40, 5)]:
    print(
        "Error for a slope of {:.0f} degrees and"
        " {:.2f} m-1 max. curvature: {:.1f} m".format(slope, maxc / 100, error_function((slope, maxc)))
    )

Error for a slope of 0 degrees and 0.00 m-1 max. curvature: 2.2 m
Error for a slope of 40 degrees and 0.00 m-1 max. curvature: 3.9 m
Error for a slope of 0 degrees and 0.05 m-1 max. curvature: 5.2 m
Error for a slope of 40 degrees and 0.05 m-1 max. curvature: 5.6 m

This pipeline will not always work optimally with default parameters: spread estimates can be affected by skewed distributions, the binning by extreme range of values, some DEMs do not have any error variability with terrain (e.g., terrestrial photogrammetry). To learn how to tune more parameters and use the subfunctions, see the gallery example: Estimation and modelling of heteroscedasticity!