{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nrasterMath with several rasters as inputs\n===============================================================\n\nCompute substract and addition between two raster bands.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Import librairies\n-------------------------------------------\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from museotoolbox.processing import RasterMath,image_mask_from_vector\nfrom museotoolbox import datasets\nimport numpy as np"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Load HistoricalMap dataset\n-------------------------------------------\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "raster,vector = datasets.load_historical_data()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Initialize rasterMath with raster\n------------------------------------\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "If invert is set to True, it means polygons will be set to nodata\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "image_mask_from_vector(vector,raster,'/tmp/mask.tif',invert=True)\nrM = RasterMath(in_image = raster,in_image_mask='/tmp/mask.tif',return_3d=True)\nrM.add_image(raster)\n\nprint('Number of rasters : '+str(len(rM.get_random_block())))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's suppose you want compute the substraction between the blue and green band of two inputs\nI suggest you to define type in numpy array to save space while creating the raster!\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "x = rM.get_block(0)\n\ndef sub(x):\n    firstBandOfFirstRaster = x[0][...,0]\n    thirdBandOfSecondRaster = x[1][...,2]\n    difference = np.subtract(firstBandOfFirstRaster,thirdBandOfSecondRaster)\n    return difference\n\nrM.add_function(sub,out_image='/tmp/sub_2inputs.tif')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run the script\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "#rM.run()\nrM.run()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot result\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from osgeo import gdal\nfrom matplotlib import pyplot as plt \n\nsrc = gdal.Open('/tmp/sub_2inputs.tif')\narr = src.ReadAsArray()\narr = np.where(arr==0,np.nan,arr)\nplt.imshow(arr)"
      ]
    }
  ],
  "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.7"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}