Blockwise coregistration

Often, biases are spatially variable, and a “global” shift may not be enough to coregister a DEM properly. In the Nuth and Kääb coregistration example, we saw that the method improved the alignment significantly, but there were still possibly nonlinear artefacts in the result. Clearly, nonlinear coregistration approaches are needed. One solution is xdem.coreg.BlockwiseCoreg, a helper to run any Coreg class over an arbitrary grid. Thanks to this tool, local errors are detected more effectively and memory usage is reduced. Indeed, the entire DEM does not need to be loaded into memory, the processes run for each block, also enabling multiprocessing.

The BlockwiseCoreg class runs in four steps:

1. Generate a subdivision grid to divide the DEM in N blocks.

2. Run the requested coregistration approach in each block and save the results.

3. Based on these results, calculation of the overall offset.

4. Apply to the entire DEM also by block, but this time with an overlap, this overlap corresponds to the maximum offset calculated in step 2.

import geoutils as gu

import matplotlib.pyplot as plt
import numpy as np
from geoutils.raster import ClusterGenerator
from geoutils.raster.distributed_computing import MultiprocConfig

import xdem

We open example files.

reference_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
dem_to_be_aligned = xdem.DEM(xdem.examples.get_path("longyearbyen_tba_dem"))
glacier_outlines = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))

# Create a stable ground mask (not glacierized) to mark "inlier data"
inlier_mask = ~glacier_outlines.create_mask(reference_dem)

The DEM to be aligned (a 1990 photogrammetry-derived DEM) has some vertical and horizontal biases that we want to avoid, as well as possible nonlinear distortions. The product is a mosaic of multiple DEMs, so “seams” may exist in the data. These can be visualized by plotting a change map:

diff_before = reference_dem - dem_to_be_aligned

diff_before.plot(cmap="RdYlBu", vmin=-10, vmax=10)
plt.show()

Horizontal and vertical shifts can be estimated using xdem.coreg.NuthKaab. Let’s prepare a coregistration class with a tiling configuration BlockwiseCoreg is also available without mp_config by specifying the block size and the save directory.

# Create a configuration with two processing nodes
mp_config = MultiprocConfig(chunk_size=500, outfile="aligned_dem.tif", cluster=ClusterGenerator("multi", nb_workers=2))
# Currently, with mp_config, block_size_fit must be equal to chunk_size
blockwise = xdem.coreg.BlockwiseCoreg(
    xdem.coreg.NuthKaab(), mp_config=mp_config, block_size_fit=500, block_size_apply=1000
)

Coregistration is performed with the .fit() method.

blockwise.fit(reference_dem, dem_to_be_aligned, inlier_mask)
blockwise.apply(dem_to_be_aligned)

aligned_dem = xdem.DEM("aligned_dem.tif")

The estimated shifts can be visualized by applying the coregistration to a completely flat surface. This shows the estimated shifts that would be applied in elevation; additional horizontal shifts will also be applied if the method supports it.

rows, cols, _ = blockwise.shape_tiling_grid

matrix_x = np.full((rows, cols), np.nan)
matrix_y = np.full((rows, cols), np.nan)
matrix_z = np.full((rows, cols), np.nan)

for key, value in blockwise.meta["outputs"].items():
    row, col = map(int, key.split("_"))
    matrix_x[row, col] = value["shift_x"]
    matrix_y[row, col] = value["shift_y"]
    matrix_z[row, col] = value["shift_z"]


def plot_heatmap(matrix, title, cmap, ax):
    im = ax.imshow(matrix, cmap=cmap)
    for (i, j), val in np.ndenumerate(matrix):
        ax.text(j, i, f"{val:.2f}", ha="center", va="center", color="black")
    ax.set_title(title)
    ax.set_xticks(np.arange(cols))
    ax.set_yticks(np.arange(rows))
    ax.invert_yaxis()
    plt.colorbar(im, ax=ax)


fig, axes = plt.subplots(1, 3, figsize=(18, 6))
plot_heatmap(matrix_x, "shifts in X", "Reds", axes[0])
plot_heatmap(matrix_y, "shifts in Y", "Greens", axes[1])
plot_heatmap(matrix_z, "shifts in Z", "Blues", axes[2])

plt.tight_layout()
plt.show()

Then, the new difference can be plotted to validate that it improved.

diff_after = reference_dem - aligned_dem

diff_after.plot(cmap="RdYlBu", vmin=-10, vmax=10)
plt.show()

We can compare the NMAD to validate numerically that there was an improvement:

print(f"Error before: {gu.stats.nmad(diff_before[inlier_mask]):.2f} m")
print(f"Error after: {gu.stats.nmad(diff_after[inlier_mask]):.2f} m")

Error before: 3.42 m
Error after: 2.22 m