xdem.coreg¶
DEM coregistration classes and functions.
Contents
Functions¶
apply_matrix¶
- xdem.coreg.apply_matrix(dem, transform, matrix, invert=False, centroid=None, resampling='bilinear', dilate_mask=False)[source]¶
Apply a 3D transformation matrix to a 2.5D DEM.
The transformation is applied as a value correction using linear deramping, and 2D image warping.
Convert the DEM into a point cloud (not for gridding; for estimating the DEM shifts).
Transform the point cloud in 3D using the 4x4 matrix.
Measure the difference in elevation between the original and transformed points.
Estimate a linear deramp from the elevation difference, and apply the correction to the DEM values.
Convert the horizontal coordinates of the transformed points to pixel index coordinates.
Apply the pixel-wise displacement in 2D using the new pixel coordinates.
Apply the same displacement to a nodata-mask to exclude previous and/or new nans.
- Parameters
dem – The DEM to transform.
transform – The Affine transform object (georeferencing) of the DEM.
matrix – A 4x4 transformation matrix to apply to the DEM.
invert – Invert the transformation matrix.
centroid – The X/Y/Z transformation centroid. Irrelevant for pure translations. Defaults to the midpoint (Z=0)
resampling – The resampling method to use. Can be nearest, bilinear, cubic or an integer from 0-5.
dilate_mask – Dilate the nan mask to exclude edge pixels that could be wrong.
- Returns
The transformed DEM with NaNs as nodata values (replaces a potential mask of the input dem).
apply_xy_shift¶
- xdem.coreg.apply_xy_shift(ds, dx, dy)[source]¶
Apply horizontal shift to rio dataset using Transform affine matrix :type ds:
DatasetReader:param ds: DEM :type dx:float:param dx: dx shift value :type dy:float:param dy: dy shift valueReturns: Rio Dataset with updated transform
- Return type
apply_z_shift¶
- xdem.coreg.apply_z_shift(ds, dz)[source]¶
Apply vertical shift to rio dataset using Transform affine matrix :type ds:
DatasetReader:param ds: DEM :param dx: dz shift value
deramping¶
- xdem.coreg.deramping(elevation_difference, x_coordinates, y_coordinates, degree, verbose=False, metadata=None)[source]¶
Calculate a deramping function to account for rotational and non-rigid components of the elevation difference.
- Parameters
elevation_difference – The elevation difference array to analyse.
x_coordinates – x-coordinates of the above array (must have the same shape as elevation_difference)
y_coordinates – y-coordinates of the above array (must have the same shape as elevation_difference)
degree – The polynomial degree to estimate the ramp.
verbose – Print the least squares optimization progress.
metadata – Optional. A metadata dictionary that will be updated with the key “deramp”.
- Returns
A callable function to estimate the ramp.
get_horizontal_shift¶
- xdem.coreg.get_horizontal_shift(elevation_difference, slope, aspect, min_count=20)[source]¶
Calculate the horizontal shift between two DEMs using the method presented in Nuth and Kääb (2011).
- Parameters
elevation_difference – The elevation difference (reference_dem - aligned_dem).
slope – A slope map with the same shape as elevation_difference (units = pixels?).
aspect – An aspect map with the same shape as elevation_difference (units = radians).
min_count – The minimum allowed bin size to consider valid.
- Raises
ValueError – If very few finite values exist to analyse.
- Returns
The pixel offsets in easting, northing, and the c_parameter (altitude?).
get_terrainattr¶
- xdem.coreg.get_terrainattr(ds, attrib='slope_degrees')[source]¶
Derive terrain attribute for DEM opened with rasterio. One of “slope_degrees”, “slope_percentage”, “aspect”, “profile_curvature”, “planform_curvature”, “curvature” and others (see richDEM documentation) :type ds: rio.DatasetReader :param ds: DEM :param attrib: terrain attribute :rtype: rd.rdarray :return:
mask_as_array¶
warp_dem¶
- xdem.coreg.warp_dem(dem, transform, source_coords, destination_coords, resampling='cubic', trim_border=True)[source]¶
Warp a DEM using a set of source-destination 2D or 3D coordinates.
- Parameters
dem (
ndarray) – The DEM to warp. Allowed shapes are (1, row, col) or (row, col)transform (
Affine) – The Affine transform of the DEM.source_coords (
ndarray) – The source 2D or 3D points. must be X/Y/(Z) coords of shape (N, 2) or (N, 3).destination_coords (
ndarray) – The destination 2D or 3D points. Must have the exact same shape as ‘source_coords’resampling (
str) – The resampling order to use. Choices: [‘nearest’, ‘linear’, ‘cubic’].trim_border (
bool) – Remove values outside of the interpolation regime (True) or leave them unmodified (False).
- Raises
ValueError – If the inputs are poorly formatted.
AssertionError – For unexpected outputs.
- Return type
- Returns
A warped DEM with the same shape as the input.
Classes¶
BiasCorr¶
- class xdem.coreg.BiasCorr(bias_func=<function average>)[source]¶
Bases:
xdem.coreg.CoregDEM bias correction.
Estimates the mean (or median, weighted avg., etc.) offset between two DEMs.
BlockwiseCoreg¶
- class xdem.coreg.BlockwiseCoreg(coreg, subdivision, success_threshold=0.8, n_threads=None)[source]¶
Bases:
xdem.coreg.CoregBlock-wise coreg class for nonlinear estimations.
A coreg class of choice is run on an arbitrary subdivision of the raster. When later applying the coregistration, the optimal warping is interpolated based on X/Y/Z shifts from the coreg algorithm at the grid points.
E.g. a subdivision of 4 means to divide the DEM in four equally sized parts. These parts are then coregistered separately, creating four Coreg.fit results. If the subdivision is not divisible by the raster shape, subdivision is made as best as possible to have approximately equal pixel counts.
- __init__(coreg, subdivision, success_threshold=0.8, n_threads=None)[source]¶
Instantiate a blockwise coreg object.
- Parameters
coreg – An instantiated coreg object to fit in the subdivided DEMs.
subdivision – The number of chunks to divide the DEMs in. E.g. 4 means four different transforms.
success_threshold – Raise an error if fewer chunks than the fraction failed for any reason.
n_threads – The maximum amount of threads to use. Default=auto
- stats()[source]¶
Return statistics for each chunk in the blockwise coregistration.
center_{x,y,z}: The center coordinate of the chunk in georeferenced units.
{x,y,z}_off: The calculated offset in georeferenced units.
inlier_count: The number of pixels that were inliers in the chunk.
nmad: The NMAD after coregistration.
median: The bias after coregistration.
- Raises
ValueError – If no coregistration results exist yet.
- Return type
DataFrame- Returns
A dataframe of statistics for each chunk.
- subdivide_array(shape)[source]¶
Return the grid subdivision for a given DEM shape.
- Parameters
shape – The shape of the input DEM.
- Returns
An array of shape ‘shape’ with ‘self.subdivision’ unique indices.
- to_points()[source]¶
Convert the blockwise coregistration matrices to 3D (source -> destination) points.
- The returned shape is (N, 3, 2) where the dimensions represent:
The point index where N is equal to the amount of subdivisions.
The X/Y/Z coordinate of the point.
The old/new position of the point.
To acquire the first point’s original position: points[0, :, 0] To acquire the first point’s new position: points[0, :, 1] To acquire the first point’s Z difference: points[0, 2, 1] - points[0, 2, 0]
- Return type
- Returns
An array of 3D source -> destionation points.
Coreg¶
- class xdem.coreg.Coreg(meta=None, matrix=None)[source]¶
Bases:
objectGeneric Coreg class.
Made to be subclassed.
- error(reference_dem, dem_to_be_aligned, error_type='nmad', inlier_mask=None, transform=None)[source]¶
Calculate the error of a coregistration approach.
- Choices:
“nmad”: Default. The Normalized Median Absolute Deviation of the residuals.
“median”: The median of the residuals.
“mean”: The mean/average of the residuals
“std”: The standard deviation of the residuals.
“rms”: The root mean square of the residuals.
“mae”: The mean absolute error of the residuals.
“count”: The residual count.
- Parameters
reference_dem – 2D array of elevation values acting reference.
dem_to_be_aligned – 2D array of elevation values to be aligned.
error_type – The type of error meaure to calculate. May be a list of error types.
inlier_mask – Optional. 2D boolean array of areas to include in the analysis (inliers=True).
transform – Optional. Transform of the reference_dem. Mandatory in some cases.
- Returns
The error measure of choice for the residuals.
- fit(reference_dem, dem_to_be_aligned, inlier_mask=None, transform=None, weights=None, subsample=1.0, verbose=False)[source]¶
Estimate the coregistration transform on the given DEMs.
- Parameters
reference_dem (
Union[ndarray,MaskedArray]) – 2D array of elevation values acting reference.dem_to_be_aligned (
Union[ndarray,MaskedArray]) – 2D array of elevation values to be aligned.inlier_mask (
Optional[ndarray]) – Optional. 2D boolean array of areas to include in the analysis (inliers=True).transform (
Optional[Affine]) – Optional. Transform of the reference_dem. Mandatory in some cases.weights (
Optional[ndarray]) – Optional. Per-pixel weights for the coregistration.subsample (
Union[float,int]) – Subsample the input to increase performance. <1 is parsed as a fraction. >1 is a pixel count.verbose (
bool) – Print progress messages to stdout.
- classmethod from_matrix(matrix)[source]¶
Instantiate a generic Coreg class from a transformation matrix.
- Parameters
matrix (
ndarray) – A 4x4 transformation matrix. Shape must be (4,4).- Raises
ValueError – If the matrix is incorrectly formatted.
- Returns
The instantiated generic Coreg class.
- classmethod from_translation(x_off=0.0, y_off=0.0, z_off=0.0)[source]¶
Instantiate a generic Coreg class from a X/Y/Z translation.
- Parameters
x_off (
float) – The offset to apply in the X (west-east) direction.y_off (
float) – The offset to apply in the Y (south-north) direction.z_off (
float) – The offset to apply in the Z (vertical) direction.
- Raises
ValueError – If the given translation contained invalid values.
- Returns
An instantiated generic Coreg class.
- property is_affine: bool¶
Check if the transform be explained by a 3D affine transform.
- Return type
bool
- residuals(reference_dem, dem_to_be_aligned, inlier_mask=None, transform=None)[source]¶
Calculate the residual offsets (the difference) between two DEMs after applying the transformation.
- Parameters
reference_dem (
Union[ndarray,MaskedArray]) – 2D array of elevation values acting reference.dem_to_be_aligned (
Union[ndarray,MaskedArray]) – 2D array of elevation values to be aligned.inlier_mask (
Optional[ndarray]) – Optional. 2D boolean array of areas to include in the analysis (inliers=True).transform (
Optional[Affine]) – Optional. Transform of the reference_dem. Mandatory in some cases.
- Return type
- Returns
A 1D array of finite residuals.
CoregPipeline¶
- class xdem.coreg.CoregPipeline(pipeline)[source]¶
Bases:
xdem.coreg.CoregA sequential set of coregistration steps.
Deramp¶
- class xdem.coreg.Deramp(degree=1, subsample=500000.0)[source]¶
Bases:
xdem.coreg.CoregPolynomial DEM deramping.
Estimates an n-D polynomial between the difference of two DEMs.
- __init__(degree=1, subsample=500000.0)[source]¶
Instantiate a deramping correction object.
- Parameters
degree (
int) – The polynomial degree to estimate. degree=0 is a simple bias correction.subsample (
Union[int,float]) – Factor for subsampling the input raster for speed-up.
If <= 1, will be considered a fraction of valid pixels to extract. If > 1 will be considered the number of pixels to extract.
ICP¶
- class xdem.coreg.ICP(max_iterations=100, tolerance=0.05, rejection_scale=2.5, num_levels=6)[source]¶
Bases:
xdem.coreg.CoregIterative Closest Point DEM coregistration.
Estimates a rigid transform (rotation + translation) between two DEMs.
Requires ‘opencv’ See opencv docs for more info: https://docs.opencv.org/master/dc/d9b/classcv_1_1ppf__match__3d_1_1ICP.html
- __init__(max_iterations=100, tolerance=0.05, rejection_scale=2.5, num_levels=6)[source]¶
Instantiate an ICP coregistration object.
- Parameters
max_iterations – The maximum allowed iterations before stopping.
tolerance – The residual change threshold after which to stop the iterations.
rejection_scale – The threshold (std * rejection_scale) to consider points as outliers.
num_levels – Number of octree levels to consider. A higher number is faster but may be more inaccurate.
NuthKaab¶
- class xdem.coreg.NuthKaab(max_iterations=10, offset_threshold=0.05)[source]¶
Bases:
xdem.coreg.CoregNuth and Kääb (2011) DEM coregistration.
Implemented after the paper: https://doi.org/10.5194/tc-5-271-2011
ZScaleCorr¶
- class xdem.coreg.ZScaleCorr(degree=1, bin_count=100)[source]¶
Bases:
xdem.coreg.CoregCorrect linear or nonlinear elevation scale errors.
Often useful for nadir image DEM correction, where the focal length is slightly miscalculated.
DISCLAIMER: This function may introduce error when correcting non-photogrammetric biases. See Gardelle et al. (2012) (Figure 2), http://dx.doi.org/10.3189/2012jog11j175, for curvature-related biases.