.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "basic_examples/plot_nuth_kaab.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_basic_examples_plot_nuth_kaab.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_basic_examples_plot_nuth_kaab.py:


Nuth and Kääb coregistration
============================

The Nuth and Kääb coregistration corrects horizontal and vertical shifts, and is especially performant for precise
sub-pixel alignment in areas with varying slope.
In xDEM, this approach is implemented through the :class:`xdem.coreg.NuthKaab` class.

See also the :ref:`nuthkaab` section in feature pages.

**Reference:** `Nuth and Kääb (2011) <https:https://doi.org/10.5194/tc-5-271-2011>`_.

.. GENERATED FROM PYTHON SOURCE LINES 13-19

.. code-block:: Python


    import geoutils as gu
    import numpy as np

    import xdem








.. GENERATED FROM PYTHON SOURCE LINES 20-21

We open example files.

.. GENERATED FROM PYTHON SOURCE LINES 21-28

.. 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"))

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








.. GENERATED FROM PYTHON SOURCE LINES 29-31

The DEM to be aligned (a 1990 photogrammetry-derived DEM) has some vertical and horizontal biases that we want to reduce.
These can be visualized by plotting a change map:

.. GENERATED FROM PYTHON SOURCE LINES 31-35

.. code-block:: Python


    diff_before = reference_dem - dem_to_be_aligned
    diff_before.plot(cmap="RdYlBu", vmin=-10, vmax=10, cbar_title="Elevation change (m)")




.. image-sg:: /basic_examples/images/sphx_glr_plot_nuth_kaab_001.png
   :alt: plot nuth kaab
   :srcset: /basic_examples/images/sphx_glr_plot_nuth_kaab_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 36-38

Horizontal and vertical shifts can be estimated using :class:`~xdem.coreg.NuthKaab`.
The shifts are estimated then applied to the to-be-aligned elevation data:

.. GENERATED FROM PYTHON SOURCE LINES 38-42

.. code-block:: Python


    nuth_kaab = xdem.coreg.NuthKaab()
    aligned_dem = nuth_kaab.fit_and_apply(reference_dem, dem_to_be_aligned, inlier_mask)








.. GENERATED FROM PYTHON SOURCE LINES 43-44

The shifts are stored in the affine metadata output

.. GENERATED FROM PYTHON SOURCE LINES 44-47

.. code-block:: Python


    print([nuth_kaab.meta["outputs"]["affine"][s] for s in ["shift_x", "shift_y", "shift_z"]])





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    [np.float64(9.180569669232375), np.float64(2.830744657104157), -1.9912575001893416]




.. GENERATED FROM PYTHON SOURCE LINES 48-49

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

.. GENERATED FROM PYTHON SOURCE LINES 49-53

.. code-block:: Python


    diff_after = reference_dem - aligned_dem
    diff_after.plot(cmap="RdYlBu", vmin=-10, vmax=10, cbar_title="Elevation change (m)")




.. image-sg:: /basic_examples/images/sphx_glr_plot_nuth_kaab_002.png
   :alt: plot nuth kaab
   :srcset: /basic_examples/images/sphx_glr_plot_nuth_kaab_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 54-55

We compare the median and NMAD to validate numerically that there was an improvement (see :ref:`robuststats-meanstd`):

.. GENERATED FROM PYTHON SOURCE LINES 55-64

.. code-block:: Python

    inliers_before = diff_before[inlier_mask]
    med_before, nmad_before = np.ma.median(inliers_before), gu.stats.nmad(inliers_before)

    inliers_after = diff_after[inlier_mask]
    med_after, nmad_after = np.ma.median(inliers_after), gu.stats.nmad(inliers_after)

    print(f"Error before: median = {med_before:.2f} - NMAD = {nmad_before:.2f} m")
    print(f"Error after: median = {med_after:.2f} - NMAD = {nmad_after:.2f} m")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Error before: median = -2.33 - NMAD = 3.42 m
    Error after: median = 0.00 - NMAD = 2.51 m




.. GENERATED FROM PYTHON SOURCE LINES 65-68

In the plot above, one may notice a positive (blue) tendency toward the east.
The 1990 DEM is a mosaic, and likely has a "seam" near there.
:ref:`sphx_glr_advanced_examples_plot_blockwise_coreg.py` tackles this issue, using a nonlinear coregistration approach.


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 2.820 seconds)


.. _sphx_glr_download_basic_examples_plot_nuth_kaab.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_nuth_kaab.ipynb <plot_nuth_kaab.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_nuth_kaab.py <plot_nuth_kaab.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_nuth_kaab.zip <plot_nuth_kaab.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_