.. 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 arbitrarily small grid, and then "puppet warp" the DEM to fit the reference best. The ``BlockwiseCoreg`` class runs in five steps: 1. Generate a subdivision grid to divide the DEM in N blocks. 2. Run the requested coregistration approach in each block. 3. Extract each result as a source and destination X/Y/Z point. 4. Interpolate the X/Y/Z point-shifts into three shift-rasters. 5. Warp the DEM to apply the X/Y/Z shifts. .. GENERATED FROM PYTHON SOURCE LINES 19-26 .. code-block:: Python import geoutils as gu import matplotlib.pyplot as plt import numpy as np import xdem .. GENERATED FROM PYTHON SOURCE LINES 28-29 **Example files** .. GENERATED FROM PYTHON SOURCE LINES 29-44 .. 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) plt_extent = [ reference_dem.bounds.left, reference_dem.bounds.right, reference_dem.bounds.bottom, reference_dem.bounds.top, ] .. GENERATED FROM PYTHON SOURCE LINES 45-48 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 48-54 .. code-block:: Python diff_before = reference_dem - dem_to_be_aligned diff_before.plot(cmap="coolwarm_r", 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 55-57 Horizontal and vertical shifts can be estimated using :class:`xdem.coreg.NuthKaab`. Let's prepare a coregistration class that calculates 64 offsets, evenly spread over the DEM. .. GENERATED FROM PYTHON SOURCE LINES 57-61 .. code-block:: Python blockwise = xdem.coreg.BlockwiseCoreg(xdem.coreg.NuthKaab(), subdivision=64) .. GENERATED FROM PYTHON SOURCE LINES 62-64 The grid that will be used can be visualized with a helper function. Coregistration will be performed in each block separately. .. GENERATED FROM PYTHON SOURCE LINES 64-69 .. code-block:: Python plt.title("Subdivision grid") plt.imshow(blockwise.subdivide_array(dem_to_be_aligned.shape), cmap="gist_ncar") plt.show() .. image-sg:: /advanced_examples/images/sphx_glr_plot_blockwise_coreg_002.png :alt: Subdivision grid :srcset: /advanced_examples/images/sphx_glr_plot_blockwise_coreg_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 70-72 Coregistration is performed with the ``.fit()`` method. This runs in multiple threads by default, so more CPU cores are preferable here. .. GENERATED FROM PYTHON SOURCE LINES 72-77 .. code-block:: Python blockwise.fit(reference_dem, dem_to_be_aligned, inlier_mask=inlier_mask) aligned_dem = blockwise.apply(dem_to_be_aligned) .. GENERATED FROM PYTHON SOURCE LINES 78-81 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. The :func:`xdem.coreg.BlockwiseCoreg.stats` method can be used to annotate each block with its associated Z shift. .. GENERATED FROM PYTHON SOURCE LINES 81-90 .. code-block:: Python z_correction = blockwise.apply( np.zeros_like(dem_to_be_aligned.data), transform=dem_to_be_aligned.transform, crs=dem_to_be_aligned.crs )[0] plt.title("Vertical correction") plt.imshow(z_correction, cmap="coolwarm_r", vmin=-10, vmax=10, extent=plt_extent) for _, row in blockwise.stats().iterrows(): plt.annotate(round(row["z_off"], 1), (row["center_x"], row["center_y"]), ha="center") .. image-sg:: /advanced_examples/images/sphx_glr_plot_blockwise_coreg_003.png :alt: Vertical correction :srcset: /advanced_examples/images/sphx_glr_plot_blockwise_coreg_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 91-92 Then, the new difference can be plotted to validate that it improved. .. GENERATED FROM PYTHON SOURCE LINES 92-98 .. code-block:: Python diff_after = reference_dem - aligned_dem diff_after.plot(cmap="coolwarm_r", vmin=-10, vmax=10) plt.show() .. image-sg:: /advanced_examples/images/sphx_glr_plot_blockwise_coreg_004.png :alt: plot blockwise coreg :srcset: /advanced_examples/images/sphx_glr_plot_blockwise_coreg_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 99-100 We can compare the NMAD to validate numerically that there was an improvment: .. GENERATED FROM PYTHON SOURCE LINES 100-103 .. code-block:: Python print(f"Error before: {xdem.spatialstats.nmad(diff_before):.2f} m") print(f"Error after: {xdem.spatialstats.nmad(diff_after):.2f} m") .. rst-class:: sphx-glr-script-out .. code-block:: none Error before: 3.81 m Error after: 2.24 m .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 9.922 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 ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_