.. 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 ====================================== Iterative Closest Point (ICP) is a registration methods accounting for both rotation and translation. It is used primarily to correct rotations, as it performs worse than :ref:`coregistration-nuthkaab` for sub-pixel shifts. Fortunately, xDEM provides the best of two worlds by allowing a combination of the two methods in a pipeline, demonstrated below! **Reference**: `Besl and McKay (1992) `_. .. GENERATED FROM PYTHON SOURCE LINES 13-18 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np import xdem .. GENERATED FROM PYTHON SOURCE LINES 20-22 We 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:: Python dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem")) subset_extent = [523000, 8660000, 529000, 8665000] dem = 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:: Python xdem.terrain.hillshade(dem).plot(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-49 .. code-block:: Python 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 = xdem.coreg.apply_matrix(dem, matrix=rotation_matrix) .. rst-class:: sphx-glr-script-out .. code-block:: none Residual check at iteration number 1: Mean diff x: 0.26929990302364193 Mean diff y: 0.0 Points not within tolerance: 73498 for X; 0 for Y Residual check at iteration number 5: Mean diff x: 1.4429699573081453e-09 Mean diff y: 0.0 Points not within tolerance: 0 for X; 0 for Y .. GENERATED FROM PYTHON SOURCE LINES 50-52 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 52-56 .. code-block:: Python diff_before = dem - rotated_dem diff_before.plot(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 57-63 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 63-91 .. code-block:: Python 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_elev=dem, to_be_aligned_elev=rotated_dem, ) corrected_dem = approach.apply(elev=rotated_dem) diff = dem - corrected_dem ax = plt.subplot(3, 1, i + 1) plt.title(name) diff.plot(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 .. rst-class:: sphx-glr-script-out .. code-block:: none Residual check at iteration number 1: Mean diff x: 0.270728172614538 Mean diff y: 6.042520672455082e-05 Points not within tolerance: 72045 for X; 0 for Y Residual check at iteration number 5: Mean diff x: 1.3032642358797546e-09 Mean diff y: 1.551234769294988e-13 Points not within tolerance: 0 for X; 0 for Y Residual check at iteration number 1: Mean diff x: 0.270728172614538 Mean diff y: 6.042520672455082e-05 Points not within tolerance: 72045 for X; 0 for Y Residual check at iteration number 5: Mean diff x: 1.3032642358797546e-09 Mean diff y: 1.551234769294988e-13 Points not within tolerance: 0 for X; 0 for Y Residual check at iteration number 1: Mean diff x: 0.270728172614538 Mean diff y: 6.042520672455082e-05 Points not within tolerance: 72045 for X; 0 for Y Residual check at iteration number 5: Mean diff x: 1.3032642358797546e-09 Mean diff y: 1.551234769294988e-13 Points not within tolerance: 0 for X; 0 for Y .. GENERATED FROM PYTHON SOURCE LINES 92-99 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 8.153 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-jupyter :download:`Download Jupyter notebook: plot_icp_coregistration.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_icp_coregistration.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_