{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Normalized regional hypsometric interpolation\n\nThere are many ways of interpolating gaps in a dDEM.\nIn the case of glaciers, one very useful fact is that elevation change is generally varies with elevation.\nThis means that if valid pixels exist in a certain elevation bin, their values can be used to fill other pixels in the same approximate elevation.\nFilling gaps by elevation is the main basis of \"hypsometric interpolation approaches\", of which there are many variations of.\n\nOne problem with simple hypsometric approaches is that they may not work glaciers with different elevation ranges and scales.\nLet's say we have two glaciers: one gigantic reaching from 0-1000 m, and one small from 900-1100 m.\nUsually in the 2000s, glaciers thin rapidly at the bottom, while they may be neutral or only thin slightly in the top.\nIf we extrapolate the hypsometric signal of the gigantic glacier to use on the small one, it may seem like the smaller glacier has almost no change whatsoever.\nThis may be right, or it may be catastrophically wrong!\n\nNormalized regional hypsometric interpolation solves the scale and elevation range problems in one go. It:\n\n    1. Calculates a regional signal using the weighted average of each glacier's normalized signal:\n\n        a. The glacier's elevation range is scaled from 0-1 to be elevation-independent.\n        b. The glaciers elevation change is scaled from 0-1 to be magnitude-independent.\n        c. A weight is assigned by the amount of valid pixels (well-covered large glaciers gain a higher weight)\n\n    2. Re-scales that signal to fit each glacier once determined.\n\nThe consequence is a much more accurate interpolation approach that can be used in a multitude of glacierized settings.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# sphinx_gallery_thumbnail_number = 2\nimport matplotlib.pyplot as plt\n\nimport numpy as np\nimport xdem\nimport xdem.misc\nimport geoutils as gu"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Example files**\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "xdem.examples.download_longyearbyen_examples()\n\ndem_2009 = xdem.DEM(xdem.examples.get_path(\"longyearbyen_ref_dem\"), silent=True)\ndem_1990 = xdem.DEM(xdem.examples.get_path(\"longyearbyen_tba_dem_coreg\"), silent=True)\n\nglacier_outlines = gu.Vector(xdem.examples.get_path(\"longyearbyen_glacier_outlines\"))\n\n# Rasterize the glacier outlines to create an index map.\n# Stable ground is 0, the first glacier is 1, the second is 2, etc.\nglacier_index_map = glacier_outlines.rasterize(dem_2009)\n\nplt_extent = [\n    dem_2009.bounds.left,\n    dem_2009.bounds.right,\n    dem_2009.bounds.bottom,\n    dem_2009.bounds.top,\n]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To test the method, we can generate a semi-random mask to assign nans to glacierized areas.\nLet's remove 30% of the data.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "np.random.seed(42)\nrandom_nans = (xdem.misc.generate_random_field(dem_1990.shape, corr_size=5) > 0.7) & (glacier_index_map > 0)\n\nplt.imshow(random_nans)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The normalized hypsometric signal shows the tendency for elevation change as a function of elevation.\nThe magnitude may vary between glaciers, but the shape is generally similar.\nNormalizing by both elevation and elevation change, and then re-scaling the signal to every glacier, ensures that it is as accurate as possible.\n**NOTE**: The hypsometric signal does not need to be generated separately; it will be created by :func:`xdem.volume.norm_regional_hypsometric_interpolation`.\nGenerating it first, however, allows us to visualize and validate it.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ddem = (dem_2009 - dem_1990).data\nddem_voided = np.where(random_nans, np.nan, ddem)\n\nsignal = xdem.volume.get_regional_hypsometric_signal(\n    ddem=ddem_voided,\n    ref_dem=dem_2009.data,\n    glacier_index_map=glacier_index_map,\n)\n\nplt.fill_between(signal.index.mid, signal[\"sigma-1-lower\"], signal[\"sigma-1-upper\"], label=\"Spread (+- 1 sigma)\")\nplt.plot(signal.index.mid, signal[\"w_mean\"], color=\"black\", label=\"Weighted mean\")\nplt.ylabel(\"Normalized elevation change\")\nplt.xlabel(\"Normalized elevation\")\nplt.legend()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The signal can now be used (or simply estimated again if not provided) to interpolate the DEM.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ddem_filled = xdem.volume.norm_regional_hypsometric_interpolation(\n    voided_ddem=ddem_voided,\n    ref_dem=dem_2009.data,\n    glacier_index_map=glacier_index_map,\n    regional_signal=signal\n)\n\n\nplt.figure(figsize=(8, 5))\nplt.imshow(ddem_filled, cmap=\"coolwarm_r\", vmin=-10, vmax=10, extent=plt_extent)\nplt.colorbar()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can plot the difference between the actual and the interpolated values, to validate the method.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "difference = (ddem_filled - ddem).squeeze()[random_nans].filled(np.nan)\nmedian = np.nanmedian(difference)\nnmad = xdem.spatial_tools.nmad(difference)\n\nplt.title(f\"Median: {median:.2f} m, NMAD: {nmad:.2f} m\")\nplt.hist(difference, bins=np.linspace(-15, 15, 100))\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As we see, the median is close to zero, while the `spatial_stats_nmad` varies slightly more.\nThis is expected, as the regional signal is good for multiple glaciers at once, but it cannot account for difficult local topography and meteorological conditions.\nIt is therefore highly recommended for large regions; just don't zoom in too close!\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
}