.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "basic_examples/plot_icp_coregistration.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_basic_examples_plot_icp_coregistration.py: Iterative Closest Point coregistration ====================================== Some DEMs may for one or more reason be erroneously rotated in the X, Y or Z directions. Established coregistration approaches like :ref:`coregistration-nuthkaab` work great for X, Y and Z *translations*, but rotation is not accounted for at all. Iterative Closest Point (ICP) is one method that takes both rotation and translation into account. It is however not as good as :ref:`coregistration-nuthkaab` when it comes to sub-pixel accuracy. Fortunately, ``xdem`` provides the best of two worlds by allowing a combination of the two. **Reference**: `Besl and McKay (1992) `_. .. GENERATED FROM PYTHON SOURCE LINES 13-18 .. code-block:: default import matplotlib.pyplot as plt import numpy as np import xdem .. GENERATED FROM PYTHON SOURCE LINES 20-22 Let's load a DEM and crop it to a single mountain on Svalbard, called Battfjellet. Its aspects vary in every direction, and is therefore a good candidate for coregistration exercises. .. GENERATED FROM PYTHON SOURCE LINES 22-27 .. code-block:: default dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem")) subset_extent = [523000, 8660000, 529000, 8665000] dem.crop(subset_extent) .. GENERATED FROM PYTHON SOURCE LINES 28-29 Let's plot a hillshade of the mountain for context. .. GENERATED FROM PYTHON SOURCE LINES 29-31 .. code-block:: default xdem.terrain.hillshade(dem).show(cmap="gray") .. image-sg:: /basic_examples/images/sphx_glr_plot_icp_coregistration_001.png :alt: plot icp coregistration :srcset: /basic_examples/images/sphx_glr_plot_icp_coregistration_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 32-35 To try the effects of rotation, we can artificially rotate the DEM using a transformation matrix. Here, a rotation of just one degree is attempted. But keep in mind: the window is 6 km wide; 1 degree of rotation at the center equals to a 52 m vertical difference at the edges! .. GENERATED FROM PYTHON SOURCE LINES 35-50 .. code-block:: default rotation = np.deg2rad(1) rotation_matrix = np.array( [ [np.cos(rotation), 0, np.sin(rotation), 0], [0, 1, 0, 0], [-np.sin(rotation), 0, np.cos(rotation), 0], [0, 0, 0, 1], ] ) # This will apply the matrix along the center of the DEM rotated_dem_data = xdem.coreg.apply_matrix(dem.data.squeeze(), transform=dem.transform, matrix=rotation_matrix) rotated_dem = xdem.DEM.from_array(rotated_dem_data, transform=dem.transform, crs=dem.crs, nodata=-9999) .. GENERATED FROM PYTHON SOURCE LINES 51-53 We can plot the difference between the original and rotated DEM. It is now artificially tilting from east down to the west. .. GENERATED FROM PYTHON SOURCE LINES 53-57 .. code-block:: default diff_before = dem - rotated_dem diff_before.show(cmap="coolwarm_r", vmin=-20, vmax=20) plt.show() .. image-sg:: /basic_examples/images/sphx_glr_plot_icp_coregistration_002.png :alt: plot icp coregistration :srcset: /basic_examples/images/sphx_glr_plot_icp_coregistration_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 58-64 As previously mentioned, ``NuthKaab`` works well on sub-pixel scale but does not handle rotation. ``ICP`` works with rotation but lacks the sub-pixel accuracy. Luckily, these can be combined! Any :class:`xdem.coreg.Coreg` subclass can be added with another, making a :class:`xdem.coreg.CoregPipeline`. With a pipeline, each step is run sequentially, potentially leading to a better result. Let's try all three approaches: ``ICP``, ``NuthKaab`` and ``ICP + NuthKaab``. .. GENERATED FROM PYTHON SOURCE LINES 64-92 .. code-block:: default approaches = [ (xdem.coreg.ICP(), "ICP"), (xdem.coreg.NuthKaab(), "NuthKaab"), (xdem.coreg.ICP() + xdem.coreg.NuthKaab(), "ICP + NuthKaab"), ] plt.figure(figsize=(6, 12)) for i, (approach, name) in enumerate(approaches): approach.fit( reference_dem=dem, dem_to_be_aligned=rotated_dem, ) corrected_dem = approach.apply(dem=rotated_dem) diff = dem - corrected_dem ax = plt.subplot(3, 1, i + 1) plt.title(name) diff.show(cmap="coolwarm_r", vmin=-20, vmax=20, ax=ax) plt.tight_layout() plt.show() .. image-sg:: /basic_examples/images/sphx_glr_plot_icp_coregistration_003.png :alt: ICP, NuthKaab, ICP + NuthKaab :srcset: /basic_examples/images/sphx_glr_plot_icp_coregistration_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 93-100 The results show what we expected: * ``ICP`` alone handled the rotational offset, but left a horizontal offset as it is not sub-pixel accurate (in this case, the resolution is 20x20m). * ``NuthKaab`` barely helped at all, since the offset is purely rotational. * ``ICP + NuthKaab`` first handled the rotation, then fit the reference with sub-pixel accuracy. The last result is an almost identical raster that was offset but then corrected back to its original position! .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 4.505 seconds) .. _sphx_glr_download_basic_examples_plot_icp_coregistration.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_icp_coregistration.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_icp_coregistration.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_