.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "advanced_examples/plot_blockwise_coreg.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_advanced_examples_plot_blockwise_coreg.py: Blockwise coregistration ======================== Often, biases are spatially variable, and a "global" shift may not be enough to coregister a DEM properly. In the :ref:`sphx_glr_basic_examples_plot_nuth_kaab.py` 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 :class:`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. .. GENERATED FROM PYTHON SOURCE LINES 19-29 .. code-block:: Python 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 .. GENERATED FROM PYTHON SOURCE LINES 31-32 We open example files. .. GENERATED FROM PYTHON SOURCE LINES 32-40 .. code-block:: Python 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) .. GENERATED FROM PYTHON SOURCE LINES 41-44 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: .. GENERATED FROM PYTHON SOURCE LINES 44-50 .. code-block:: Python diff_before = reference_dem - dem_to_be_aligned diff_before.plot(cmap="RdYlBu", vmin=-10, vmax=10) plt.show() .. image-sg:: /advanced_examples/images/sphx_glr_plot_blockwise_coreg_001.png :alt: plot blockwise coreg :srcset: /advanced_examples/images/sphx_glr_plot_blockwise_coreg_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 51-54 Horizontal and vertical shifts can be estimated using :class:`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. .. GENERATED FROM PYTHON SOURCE LINES 54-62 .. code-block:: Python # 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 ) .. GENERATED FROM PYTHON SOURCE LINES 63-64 Coregistration is performed with the ``.fit()`` method. .. GENERATED FROM PYTHON SOURCE LINES 64-70 .. code-block:: Python blockwise.fit(reference_dem, dem_to_be_aligned, inlier_mask) blockwise.apply(dem_to_be_aligned) aligned_dem = xdem.DEM("aligned_dem.tif") .. GENERATED FROM PYTHON SOURCE LINES 71-74 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. .. GENERATED FROM PYTHON SOURCE LINES 74-107 .. code-block:: Python 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() .. image-sg:: /advanced_examples/images/sphx_glr_plot_blockwise_coreg_002.png :alt: shifts in X, shifts in Y, shifts in Z :srcset: /advanced_examples/images/sphx_glr_plot_blockwise_coreg_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 108-109 Then, the new difference can be plotted to validate that it improved. .. GENERATED FROM PYTHON SOURCE LINES 109-115 .. code-block:: Python diff_after = reference_dem - aligned_dem diff_after.plot(cmap="RdYlBu", vmin=-10, vmax=10) plt.show() .. image-sg:: /advanced_examples/images/sphx_glr_plot_blockwise_coreg_003.png :alt: plot blockwise coreg :srcset: /advanced_examples/images/sphx_glr_plot_blockwise_coreg_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 116-117 We can compare the NMAD to validate numerically that there was an improvement: .. GENERATED FROM PYTHON SOURCE LINES 117-120 .. code-block:: Python 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") .. rst-class:: sphx-glr-script-out .. code-block:: none Error before: 3.42 m Error after: 2.29 m .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 26.327 seconds) .. _sphx_glr_download_advanced_examples_plot_blockwise_coreg.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_blockwise_coreg.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_blockwise_coreg.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_blockwise_coreg.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_