Cheatsheet: How to correct… ?#

In elevation data analysis, the problem generally starts with identifying what correction method to apply when observing a specific pattern of error in your own data.

Below, we summarize a cheatsheet that links what method is likely to correct a pattern of error you can visually identify on a map of elevation differences with another elevation dataset (looking at static surfaces)!

Cheatsheet#

The patterns of errors categories listed in this spreadsheet are linked to visual examples further below that you can use to compare to your own elevation differences.

Pattern

Description

Cause and correction

Notes

Sharp landforms

Positive and negative errors that are larger near high slopes and symmetric with opposite slope orientation, making landforms appear visually.

Likely horizontal shift due to geopositioning errors, use a Coregistration such as NuthKaab.

Even a tiny horizontal misalignment can be visually identified! To not confuse with Peak cuts and cavity fills.

Peak cuts and cavity fills

Positive and negative errors, with one sign located exclusively near peaks and the other exclusively near cavities.

Likely resolution-type errors, use a Bias correction such as TerrainBias.

Can be over-corrected, sometimes better to simply ignore during analysis. Or avoid by downsampling all elevation data to the lowest resolution, rather than upsampling to the highest.

Smooth-field offset

Smooth offsets varying at scale of 10 km+, often same sign (either positive or negative).

Likely wrong Vertical referencing, can set and transform with set_vcrs() and to_vcrs().

Vertical references often only exists in a user guide, they are not coded in the raster CRS and need to be set manually.

Ramp or dome

Ramping errors, often near the edge of the data extent, sometimes with a center dome.

Likely ramp/rotations due to camera errors, use either a Coregistration such as ICP or a Bias correction such as Deramp.

Can sometimes be more rigorously fixed ahead of DEM generation with bundle adjustment.

Undulations

Positive and negative errors undulating patterns at one or several frequencies well larger than pixel size.

Likely jitter-type errors, use a Bias correction such as DirectionalBias.

Can sometimes be more rigorously fixed ahead of DEM generation with jitter correction.

Point alternating

Point data errors that alternate between negative and positive, higher on steeper slopes.

Likely wrong point-raster comparison, use point interpolation or reduction on the raster instead such as interp_points().

Rasterizing point data introduces spatially correlated random errors, instead it is recommended to interpolate raster data at the point coordinates.

Visual patterns of errors#

Important

The patterns of errors below are created synthetically to examplify the effect of a single source of error. In your own elevation differences, those will be mixed together and with random errors inherent to the data.

For examples on real data, see the Basic and Advanced gallery examples!

It is often crucial to relate the location of your errors on static surfaces to the terrain distribution (in particular, its slope and curvature), which can usually be inferred visually from a hillshade.

Hide the code for opening example file
import xdem
import geoutils as gu
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt

# Open an example DEM
dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
hs = dem.hillshade()
hs.plot(cmap="Greys_r", cbar_title="Hillshade")
_images/98df991705e0cc7359bee917d0ca46bab74c7df2ece037455d5839430d109b70.png

Sharp landforms#

Example of sharp landforms appearing with a horizontal shift due to geolocation errors. We here translate the DEM horizontally by 1/10th of a pixel, for a pixel resolution of 20 m.

Hide code to simulate horizontal shift errors
# Simulate a translation of 1/10th of a pixel
x_shift = 0.1
y_shift = 0.1
dem_shift = dem.translate(x_shift, y_shift, distance_unit="pixel")

# Resample and plot the elevation differences of the horizontal shift
dh = dem - dem_shift.reproject(dem)
dh.plot(cmap='RdYlBu', vmin=-3, vmax=3, cbar_title="Elevation differences of\nhorizontal shift (m)")
_images/e5353779e32a52ce30eb2141f7dcab92161ac0f3392821f769e0a59e367a1725.png

Peak cuts and cavity fills#

Example of peak cutting and cavity filling errors. We here downsampled our DEM from 20 m to 100 m to simulate a lower native resolution, then upsample it again to 20 m, to show the errors affect areas near high curvatures such as peaks and cavities.

Hide code to simulate resolution errors
# Downsample DEM (bilinear)
dem_100m = dem.reproject(res=100)

# Upsample (bilinear again) and compare
dh = dem - dem_100m.reproject(dem)
dh.plot(cmap='RdYlBu', vmin=-40, vmax=40, cbar_title="Elevation differences of\nresolution change (m)")
_images/4f09d48d9f891cb779fa1b5d9f4c31322c5b1bfced51b412ef83616e7d752bf2.png

Smooth-field offset#

Example of smooth large offset field created by a wrong vertical CRS. We here show the difference due to the EGM96 geoid added on top of the ellipsoid.

Hide code to simulate vertical referencing errors
# Set current vertical CRS as ellipsoid
dem.set_vcrs("EGM96")
# Transform vertical reference to geoid
trans_dem = dem.to_vcrs("Ellipsoid")

# Plot the elevation differences of the vertical transformation
dh = dem - trans_dem
dh.plot(cmap='RdYlBu', cbar_title="Elevation differences of\nvertical transform (m)")
_images/d964c0efaf632fac095aab19a35ac09d454e41f2e7e238108689c2f941b8fa39.png

Ramp or dome#

Example of ramp created by a rotation due to camera errors. We here show just a slight rotation of 0.02 degrees.

Hide code to simulate rotation errors
# Apply a rotation of 0.02 degrees
rotation = np.deg2rad(0.02)
# Affine matrix for 3D transformation
matrix = np.array(
    [
        [1, 0, 0, 0],
        [0, np.cos(rotation), -np.sin(rotation), 0],
        [0, np.sin(rotation), np.cos(rotation), 0],
        [0, 0, 0, 1],
    ]
)
# Select a centroid
centroid = [dem.bounds.left + 5000, dem.bounds.top - 2000, np.median(dem) + 100]
# We rotate the elevation data
dem_rotated = xdem.coreg.apply_matrix(dem, matrix, centroid=centroid)

# Plot the elevation differences of the rotation
dh = dem - dem_rotated
dh.plot(cmap='RdYlBu', cbar_title="Elevation differences of\nrotation (m)")
_images/b3675ace2fc561dd8357b286b390295b2c55ca3bd8a7d4b36eb4a52826393292.png

Undulations#

Example of undulations resembling jitter errors. We here artificially create a sinusoidal signal at a certain angle.

Hide code to simulate undulation errors
# Get rotated coordinates along an angle
angle = -20
xx = gu.raster.get_xy_rotated(dem, along_track_angle=angle)[0]

# One sinusoid: amplitude, phases and frequencies
params = np.array([(2, 3000, np.pi)]).flatten()

# Create a sinusoidal bias and add to the DEM
from xdem.fit import sumsin_1d
synthetic_bias_arr = sumsin_1d(xx.flatten(), *params).reshape(np.shape(dem.data))

# Plot the elevation differences of the undulations
synthetic_bias = dem.copy(new_array=synthetic_bias_arr)
synthetic_bias.plot(cmap='RdYlBu', vmin=-3, vmax=3, cbar_title="Elevation differences of\nundulations (m)")
_images/63dbaaf91aaae300cadc0375c3558d9dd6fc8a89b56471f4b4e40b9cb111b519.png

Point alternating#

An example of alternating point errors created by wrong point-raster comparison by rasterization of the points, which are especially large around steep slopes.

Hide code to simulate point-raster comparison errors
# Simulate swath coordinates of an elevation point cloud
x = np.linspace(dem.bounds.left, dem.bounds.right, 100)
y = np.linspace(dem.bounds.top - 5000, dem.bounds.bottom + 5000, 100)

# Interpolate DEM at these coordinates to build the point cloud
# (to approximate the real elevation at these coordinates,
# which has negligible impact compared to rasterization)
z = dem.interp_points((x,y))
epc = gu.Vector(gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=x, y=y, crs=dem.crs), data={"z": z}))

# Rasterize point cloud back on the DEM grid
epc_rast = epc.rasterize(dem, in_value=z, out_value=np.nan)

# For easier plotting, convert the valid dh values to points
dh = dem - epc_rast
dh_pc = dh.to_pointcloud(data_column_name="dh")

# Plot the elevation differences of the rasterization on top of a hillshade
hs.plot(cmap="Greys_r", add_cbar=False)
dh_pc.plot(column="dh", cmap='RdYlBu', vmin=-10, vmax=10, legend=True, cbar_title="Elevation differences of\npoint-raster differencing (m)")
/home/docs/checkouts/readthedocs.org/user_builds/xdem/conda/stable/lib/python3.12/site-packages/geoutils/raster/raster.py:1719: UserWarning: Setting default nodata -99999 to mask non-finite values found in the array, as no nodata value was defined.
  warnings.warn(
_images/c0d2ea9436bc81ce04a8f42af4169c13578c0e1fc4fa5f823ab370a46e767d42.png