Skip to content

r.dem.nk

Co-register an SfM DSM to a LiDAR DSM using a Nuth & Kaab-style model.

Estimates horizontal (dx, dy) and vertical (dz) offsets on stable terrain and applies a sub-cell translation.

r.dem.nk [-k] sfm=name lidar=name stable_mask=name output=name [interp=string] [slope_min=float] [slope_max=float] [iters=integer] [sigma=float] [max_iter=integer] [tol=float] [transform_output=name] [apply_transform=name] [--overwrite] [--verbose] [--quiet] [--qq] [--ui]

Example:

r.dem.nk sfm=name lidar=name stable_mask=name output=name

grass.tools.Tools.r_dem_nk(sfm, lidar, stable_mask, output, interp="bilinear", slope_min=2.0, slope_max=85.0, iters=2, sigma=2.5, max_iter=20, tol=0.01, transform_output=None, apply_transform=None, flags=None, overwrite=None, verbose=None, quiet=None, superquiet=None)

Example:

tools = Tools()
tools.r_dem_nk(sfm="name", lidar="name", stable_mask="name", output="name")

This grass.tools API is experimental in version 8.5 and expected to be stable in version 8.6.

grass.script.run_command("r.dem.nk", sfm, lidar, stable_mask, output, interp="bilinear", slope_min=2.0, slope_max=85.0, iters=2, sigma=2.5, max_iter=20, tol=0.01, transform_output=None, apply_transform=None, flags=None, overwrite=None, verbose=None, quiet=None, superquiet=None)

Example:

gs.run_command("r.dem.nk", sfm="name", lidar="name", stable_mask="name", output="name")

Parameters

sfm=name [required]
    Input SfM DSM raster
lidar=name [required]
    Input LiDAR DSM raster (reference)
stable_mask=name [required]
    Mask raster (1 for stable terrain, NULL elsewhere)
output=name [required]
    Output co-registered SfM DSM
interp=string
    Interpolation for sub-cell translation
    Allowed values: nearest, bilinear, bicubic
    Default: bilinear
slope_min=float
    Minimum slope (degrees) used for regression (avoid near-flat)
    Default: 2.0
slope_max=float
    Maximum slope (degrees) used for regression (avoid near-vertical)
    Default: 85.0
iters=integer
    Sigma-clipping iterations per co-registration pass (0 disables clip)
    Default: 2
sigma=float
    Sigma threshold for residual clipping (|resid| <= sigma * stddev)
    Default: 2.5
max_iter=integer
    Maximum outer co-registration passes (re-warp and re-solve)
    Default: 20
tol=float
    Convergence tolerance in map units for the outer passes
    Default: 0.01
transform_output=name
    Write the solved transform (dz, dx, dy) to a file
apply_transform=name
    Apply a saved transform (dz, dx, dy) instead of solving
-k
    Keep intermediate rasters
--overwrite
    Allow output files to overwrite existing files
--help
    Print usage summary
--verbose
    Verbose module output
--quiet
    Quiet module output
--qq
    Very quiet module output
--ui
    Force launching GUI dialog

sfm : str | np.ndarray, required
    Input SfM DSM raster
    Used as: input, raster, name
lidar : str | np.ndarray, required
    Input LiDAR DSM raster (reference)
    Used as: input, raster, name
stable_mask : str | np.ndarray, required
    Mask raster (1 for stable terrain, NULL elsewhere)
    Used as: input, raster, name
output : str | type(np.ndarray) | type(np.array) | type(gs.array.array), required
    Output co-registered SfM DSM
    Used as: output, raster, name
interp : str, optional
    Interpolation for sub-cell translation
    Allowed values: nearest, bilinear, bicubic
    Default: bilinear
slope_min : float, optional
    Minimum slope (degrees) used for regression (avoid near-flat)
    Default: 2.0
slope_max : float, optional
    Maximum slope (degrees) used for regression (avoid near-vertical)
    Default: 85.0
iters : int, optional
    Sigma-clipping iterations per co-registration pass (0 disables clip)
    Default: 2
sigma : float, optional
    Sigma threshold for residual clipping (|resid| <= sigma * stddev)
    Default: 2.5
max_iter : int, optional
    Maximum outer co-registration passes (re-warp and re-solve)
    Default: 20
tol : float, optional
    Convergence tolerance in map units for the outer passes
    Default: 0.01
transform_output : str, optional
    Write the solved transform (dz, dx, dy) to a file
    Used as: output, file, name
apply_transform : str | io.StringIO, optional
    Apply a saved transform (dz, dx, dy) instead of solving
    Used as: input, file, name
flags : str, optional
    Allowed values: k
    k
        Keep intermediate rasters
overwrite : bool, optional
    Allow output files to overwrite existing files
    Default: None
verbose : bool, optional
    Verbose module output
    Default: None
quiet : bool, optional
    Quiet module output
    Default: None
superquiet : bool, optional
    Very quiet module output
    Default: None

Returns:

result : grass.tools.support.ToolResult | np.ndarray | tuple[np.ndarray] | None
If the tool produces text as standard output, a ToolResult object will be returned. Otherwise, None will be returned. If an array type (e.g., np.ndarray) is used for one of the raster outputs, the result will be an array and will have the shape corresponding to the computational region. If an array type is used for more than one raster output, the result will be a tuple of arrays.

Raises:

grass.tools.ToolError: When the tool ended with an error.

sfm : str, required
    Input SfM DSM raster
    Used as: input, raster, name
lidar : str, required
    Input LiDAR DSM raster (reference)
    Used as: input, raster, name
stable_mask : str, required
    Mask raster (1 for stable terrain, NULL elsewhere)
    Used as: input, raster, name
output : str, required
    Output co-registered SfM DSM
    Used as: output, raster, name
interp : str, optional
    Interpolation for sub-cell translation
    Allowed values: nearest, bilinear, bicubic
    Default: bilinear
slope_min : float, optional
    Minimum slope (degrees) used for regression (avoid near-flat)
    Default: 2.0
slope_max : float, optional
    Maximum slope (degrees) used for regression (avoid near-vertical)
    Default: 85.0
iters : int, optional
    Sigma-clipping iterations per co-registration pass (0 disables clip)
    Default: 2
sigma : float, optional
    Sigma threshold for residual clipping (|resid| <= sigma * stddev)
    Default: 2.5
max_iter : int, optional
    Maximum outer co-registration passes (re-warp and re-solve)
    Default: 20
tol : float, optional
    Convergence tolerance in map units for the outer passes
    Default: 0.01
transform_output : str, optional
    Write the solved transform (dz, dx, dy) to a file
    Used as: output, file, name
apply_transform : str, optional
    Apply a saved transform (dz, dx, dy) instead of solving
    Used as: input, file, name
flags : str, optional
    Allowed values: k
    k
        Keep intermediate rasters
overwrite : bool, optional
    Allow output files to overwrite existing files
    Default: None
verbose : bool, optional
    Verbose module output
    Default: None
quiet : bool, optional
    Quiet module output
    Default: None
superquiet : bool, optional
    Very quiet module output
    Default: None

DESCRIPTION

r.dem.nk implements the Nuth and Kääb (2011) algorithm for co-registering two Digital Elevation Models (DEMs). The algorithm estimates and corrects for vertical and horizontal offsets between the DEMs by minimizing elevation differences on stable terrain.

Model

For stable terrain pixels:

dh = SfM - LiDAR
   approx. delta_x * tan(slope) * cos(aspect)
         + delta_y * tan(slope) * sin(aspect)
         + delta_z

Where delta_x (east), delta_y (north), delta_z (vertical) are solved by ordinary least squares using raster-wide sums computed internally. Only stable-terrain cells whose slope lies between slope_min and slope_max (degrees; defaults 2 and 85) enter the regression, excluding near-flat and near-vertical cells.

Iterative solve

The linear model above is only first order, so a single pass under-estimates shifts larger than one cell. The module therefore solves the increment on a working surface, accumulates it into the running transform, re-warps the working surface from the original SfM by the accumulated transform, and repeats until the increment falls below tol (in map units) or max_iter outer passes are reached. If the passes do not converge, a warning is issued and the last estimate is used.

Horizontal application (native)

The solved offsets are applied as a single inverse warp: each output cell at map coordinate (x, y) samples the original SfM at (x + delta_x, y + delta_y) (using interp) and subtracts delta_z. No region shifting or external resampling tool is involved.

Robustness

Optional iterative sigma-clipping on residuals (iters per outer pass) reduces outlier influence during each solve.

The module always writes a residual raster named output_resid which contains output - lidar on the stable-terrain mask used for regression. When -k is provided, additional intermediate rasters are written: output_slope, output_aspect, and output_mask.

Saving and reusing a transform

transform_output writes the solved offsets (dz, dx, dy) to a small text file. apply_transform reads such a file and applies it directly, skipping the regression. This lets a transform solved on one surface (for example a clean bare-earth DTM) be replayed onto another surface from the same acquisition (for example its DSM) so both share the same horizontal alignment. In apply mode the stable_mask is used only to define the reported residual raster, so a near-flat mask no longer triggers the "not enough valid pixels" error.

NOTES

The model is first order in the elevation difference, so a smooth long-wavelength vertical bias (e.g., photogrammetric doming) is partly degenerate with a horizontal shift: over sloped terrain a gentle tilt and a translation produce a similar dh pattern. When both are present the solve splits the signal between them and the reported dx and dy absorb part of the doming.

Estimate the alignment on a surface where the long-wavelength component is small, or remove it first, and treat the horizontal offsets with suspicion when the residual after co-registration still shows a broad, smoothly varying pattern. r.dem.bias method=spline is the tool for the long-wavelength part, and it operates on the difference rather than on the DEM pair, so it runs after this step.

The stable_mask must cover broad, sloped, unchanged terrain. Flat features (e.g., roads and parking lots) are filtered out by slope_min and carry no aspect information, so a mask built from them alone leaves the horizontal offsets poorly constrained. Those features belong in the PGCP vertical step of r.dem.coregister instead.

EXAMPLES

The commands below use the example scene built in the r.dem toolset manual, which is derived from the North Carolina sample dataset. Build it there first.

Solve the offset of the misregistered surface against the lidar reference. The stable mask must be broad, sloped terrain, not the flat roads used for the PGCP step:

g.region raster=elev_lid792_1m

r.dem.nk sfm=dsm_offset lidar=elev_lid792_1m \
    stable_mask=stable_terrain output=dsm_nk \
    transform_output=nk_transform.txt

The applied offset was 0.4596 m east, 0.4596 m north, and 1.32 m up, and the solve returns it:

Converged transform: dz=1.320551 dx=0.449873 dy=0.457072

Replay the saved transform onto another surface from the same acquisition, so a DSM and a DTM end up sharing one horizontal alignment:

r.dem.nk sfm=dsm_offset lidar=elev_lid792_1m stable_mask=stable_terrain \
    output=dsm_nk_replay apply_transform=nk_transform.txt

r.dem.nk example
Figure: Stable-terrain residual before and after r.dem.nk.

REFERENCES

  • Nuth, C., and A. Kääb. 2011. "Co-Registration and Bias Corrections of Satellite Elevation Data Sets for Quantifying Glacier Thickness Change." The Cryosphere 5 (1): 271-90. https://doi.org/10.5194/tc-5-271-2011

SEE ALSO

r.dem, r.dem.coregister, r.dem.icp

AUTHORS

Corey T. White, Center for Geospatial Analytics, NC State University

SOURCE CODE

Available at: r.dem.nk source code (history)
Latest change: Thursday Aug 20 19:32:37 2026 in commit d08d385