{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Nuth & K\u00e4\u00e4b (2011) coregistration\n\nNuth and K\u00e4\u00e4b (`2011 <https:https://doi.org/10.5194/tc-5-271-2011>`_) coregistration allows for horizontal and vertical shifts to be estimated and corrected for.\nIn ``xdem``, this approach is implemented through the :class:`xdem.coreg.NuthKaab` class.\n\nFor more information about the approach, see `coregistration_nuthkaab`.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import geoutils as gu\nimport matplotlib.pyplot as plt\n\nimport xdem"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Example files**\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reference_dem = xdem.DEM(xdem.examples.get_path(\"longyearbyen_ref_dem\"), silent=True)\ndem_to_be_aligned = xdem.DEM(xdem.examples.get_path(\"longyearbyen_tba_dem\"), silent=True)\nglacier_outlines = gu.Vector(xdem.examples.get_path(\"longyearbyen_glacier_outlines\"))\n\n# Create a stable ground mask (not glacierized) to mark \"inlier data\"\ninlier_mask = ~glacier_outlines.create_mask(reference_dem)\n\nplt_extent = [\n    reference_dem.bounds.left,\n    reference_dem.bounds.right,\n    reference_dem.bounds.bottom,\n    reference_dem.bounds.top,\n]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The DEM to be aligned (a 1990 photogrammetry-derived DEM) has some vertical and horizontal biases that we want to avoid.\nThese can be visualized by plotting a change map:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "diff_before = (reference_dem - dem_to_be_aligned).data\n\nplt.figure(figsize=(8, 5))\nplt.imshow(diff_before.squeeze(), cmap=\"coolwarm_r\", vmin=-10, vmax=10, extent=plt_extent)\nplt.colorbar()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Horizontal and vertical shifts can be estimated using :class:`xdem.coreg.NuthKaab`.\nFirst, the shifts are estimated, and then they can be applied to the data:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nuth_kaab = xdem.coreg.NuthKaab()\n\nnuth_kaab.fit(reference_dem.data, dem_to_be_aligned.data, inlier_mask, transform=reference_dem.transform)\n\naligned_dem_data = nuth_kaab.apply(dem_to_be_aligned.data, transform=dem_to_be_aligned.transform)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then, the new difference can be plotted to validate that it improved.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "diff_after = reference_dem.data - aligned_dem_data\n\nplt.figure(figsize=(8, 5))\nplt.imshow(diff_after.squeeze(), cmap=\"coolwarm_r\", vmin=-10, vmax=10, extent=plt_extent)\nplt.colorbar()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can compare the `spatial_stats_nmad` to validate numerically that there was an improvment:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(f\"Error before: {xdem.spatial_tools.nmad(diff_before):.2f} m\")\nprint(f\"Error after: {xdem.spatial_tools.nmad(diff_after):.2f} m\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In the plot above, one may notice a positive (blue) tendency toward the east.\nThe 1990 DEM is a mosaic, and likely has a \"seam\" near there.\n`sphx_glr_auto_examples_plot_blockwise_coreg.py` tackles this issue, using a nonlinear coregistration approach.\n\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}