{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nModal class and number of agreements\n===============================================================\n\nCreate a raster with the modal class and the number of agreements.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Import librairies\n-------------------------------------------\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import museotoolbox as mtb\nfrom scipy.stats import mode\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 = mtb.datasets.load_historical_data(low_res=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Initialize rasterMath with raster\n-----------------------------------------\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "########\n# In case you want to add a mask\nmask = '/tmp/maskFromPolygons.tif'\n\nmtb.processing.image_mask_from_vector(vector,raster,out_image = mask)\n\nrM = mtb.processing.RasterMath(raster,in_image_mask=mask)\n\nprint(rM.get_random_block())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's suppose you want compute the modal classification between several predictions\nThe first band will be the most predicted class, and the second the number of times it has been predicted.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "x = rM.get_random_block()\n\ndef modal_class(x):    \n    tmpStack = np.column_stack(mode(x,axis=1)).astype(np.int16)\n    return tmpStack\n\nmodal_class(x)\n\nrM.add_function(modal_class,out_image='/tmp/modal.tif',out_nodata=0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run the script\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rM.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/modal.tif')\ndata = src.ReadAsArray()[0,:,:]\ndata = np.where(data== 0,np.nan,data)\nplt.imshow(data)"
      ]
    }
  ],
  "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
}