{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# DEM subtraction\n\nSubtracting one DEM with another should be easy!\nThis is why ``xdem`` (with functionality from `geoutils <https://geoutils.readthedocs.io>`_) allows directly using the ``-`` or ``+`` operators on :class:`xdem.DEM` objects.\n\nBefore DEMs can be compared, they need to be reprojected/resampled/cropped to fit the same grid.\nThe :func:`xdem.DEM.reproject` method takes care of this.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport xdem"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dem_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)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can print the information about the DEMs for a \"sanity check\"\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(dem_2009)\nprint(dem_1990)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In this particular case, the two DEMs are already on the same grid (they have the same bounds, resolution and coordinate system).\nIf they don't, we need to reproject one DEM to fit the other.\n:func:`xdem.DEM.reproject` is a multi-purpose method that ensures a fit each time:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "_ = dem_1990.reproject(dem_2009)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Oops!\n``xdem`` just warned us that ``dem_1990`` did not need reprojection, but we asked it to anyway.\nTo hide this prompt, add ``.reproject(..., silent=True)``.\nBy default, :func:`xdem.DEM.reproject` uses \"nearest neighbour\" resampling (assuming resampling is needed).\nThis approach is great for speed, but it may often be erroneous when resampling.\nOther more accurate methods are \"bilinear\", \"cubic_spline\" and others.\nSee `geoutils.Raster.reproject() <https://geoutils.readthedocs.io/en/latest/api.html#geoutils.georaster.Raster.reproject>`_ and `rasterio.enums.Resampling <https://rasterio.readthedocs.io/en/latest/api/rasterio.enums.html#rasterio.enums.Resampling>`_ for more information about reprojection.\n\nNow, we are ready to generate the dDEM:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ddem = dem_2009 - dem_1990\n\nprint(ddem)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "It is a new :class:`xdem.DEM` instance, loaded in memory.\nLet's visualize it:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=(8, 5))\nplt.imshow(ddem.data.squeeze(), cmap=\"coolwarm_r\", vmin=-20, vmax=20)\n\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "... and that's it!\nFor missing values, ``xdem`` provides a number of interpolation methods which are shown in the other examples.\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
}