{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Iterative Closest Point (ICP) coregistration.\nSome DEMs may for one or more reason be erroneously rotated in the X, Y or Z directions.\nEstablished coregistration approaches like `coregistration_nuthkaab` work great for X, Y and Z *translations*, but rotation is not accounted for at all.\n\nICP is one method that takes both rotation and translation into account.\nIt is however not as good as `coregistration_nuthkaab` when it comes to sub-pixel accuracy.\nFortunately, ``xdem`` provides the best of two worlds by allowing a combination of the two.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# sphinx_gallery_thumbnail_number = 2\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nimport xdem"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's load a DEM and crop it to a single mountain on Svalbard, called Battfjellet.\nIts aspects vary in every direction, and is therefore a good candidate for coregistration exercises.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dem = xdem.DEM(xdem.examples.get_path(\"longyearbyen_ref_dem\"), silent=True)\n\nsubset_extent = [523000, 8660000, 529000, 8665000]\ndem.crop(subset_extent)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's plot a hillshade of the mountain for context.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt_extent = [subset_extent[0], subset_extent[2], subset_extent[1], subset_extent[3]]\n\nplt.imshow(xdem.spatial_tools.hillshade(dem.data, resolution=dem.res), cmap=\"Greys_r\", extent=plt_extent)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To try the effects of rotation, we can artificially rotate the DEM using a transformation matrix.\nHere, a rotation of just one degree is attempted.\nBut 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!\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rotation = np.deg2rad(1)\nrotation_matrix = np.array(\n    [\n        [np.cos(rotation), 0, np.sin(rotation), 0],\n        [0, 1, 0, 0],\n        [-np.sin(rotation), 0, np.cos(rotation), 0],\n        [0, 0, 0, 1],\n    ]\n)\n\n# This will apply the matrix along the center of the DEM\nrotated_dem = xdem.coreg.apply_matrix(dem.data.squeeze(), transform=dem.transform, matrix=rotation_matrix)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can plot the difference between the original and rotated DEM.\nIt is now artificially tilting from east down to the west.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "diff_before = dem.data - rotated_dem\nplt.imshow(diff_before.squeeze(), cmap=\"coolwarm_r\", vmin=-20, vmax=20, extent=plt_extent)\n\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As previously mentioned, ``NuthKaab`` works well on sub-pixel scale but does not handle rotation.\n``ICP`` works with rotation but lacks the sub-pixel accuracy.\nLuckily, these can be combined!\nAny :class:`xdem.coreg.Coreg` subclass can be added with another, making a :class:`xdem.coreg.CoregPipeline`.\nWith a pipeline, each step is run sequentially, potentially leading to a better result.\nLet's try all three approaches: ``ICP``, ``NuthKaab`` and ``ICP + NuthKaab``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "approaches = [\n    (xdem.coreg.ICP(), \"ICP\"),\n    (xdem.coreg.NuthKaab(), \"NuthKaab\"),\n    (xdem.coreg.ICP() + xdem.coreg.NuthKaab(), \"ICP + NuthKaab\"),\n]\n\n\nplt.figure(figsize=(6, 12))\n\nfor i, (approach, name) in enumerate(approaches):\n    approach.fit(\n        reference_dem=dem.data,\n        dem_to_be_aligned=rotated_dem,\n        transform=dem.transform,\n    )\n\n    corrected_dem_data = approach.apply(dem=rotated_dem, transform=dem.transform)\n\n    diff = dem.data - corrected_dem_data\n\n    plt.subplot(3, 1, i + 1)\n    plt.title(name)\n    plt.imshow(diff.squeeze(), cmap=\"coolwarm_r\", vmin=-20, vmax=20, extent=plt_extent)\n\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The results show what we expected:\n\n* ``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).\n* ``NuthKaab`` barely helped at all, since the offset is purely rotational.\n* ``ICP + NuthKaab`` first handled the rotation, then fit the reference with sub-pixel accuracy.\n\nThe last result is an almost identical raster that was offset but then corrected back to its original position!\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
}