GRIPPA_al_2017_AnOpenSourceSemiAutomatedProcessingChain_OBIA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This processing chain is available under [Creative Common Licence (CC-BY)](https://creativecommons.org/licenses/by/4.0/), so feel free to re-use, adapt or enhance it to match your own needs. \n",
    "\n",
    "![alt text](https://i.creativecommons.org/l/by/4.0/88x31.png)\n",
    "\n",
    "This processing chain is available on [GitHub.com](https://github.com/tgrippa/Opensource_OBIA_processing_chain). Don't hesitate to create \"Pull requests\" to propose corrections, modifications or enhancements.\n",
    "\n",
    "This processing chain is linked to a publication in [MDPI - Remote Sensing](http://www.mdpi.com/2072-4292/9/4/358/htm). If you need to refer to this processing chain, you can refer directly to this publication."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Table of Contents**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell is a Javascript section of code for building the Jupyter notebook's table of content."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"toc\"></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Instructions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to use the semi-automated processing chain, by calling up GRASS GIS functions [without starting grass explicitly](https://grasswiki.osgeo.org/wiki/Working_with_GRASS_without_starting_it_explicitly) in this [Jupyter notebook](http://jupyter.org/), please:\n",
    " \n",
    "1. Make sure that [GRASS GIS](https://grasswiki.osgeo.org/wiki/Installation_Guide) is installed and fully functional on your computer.\n",
    "2. Make sure that [Anaconda with Python 2.7](https://www.continuum.io/downloads) is installed and fully functional on your computer.\n",
    "3. Make sure that [R software](https://www.r-project.org/) is installed and fully functional on your computer.\n",
    "4. Adjust the **\"Define working environment\"** part to match your system configuration, for both [GRASS GIS' environment variables](https://grass.osgeo.org/grass64/manuals/variables.html) and [R' environment variables](https://stat.ethz.ch/R-manual/R-devel/library/base/html/EnvVar.html).\n",
    "5. Run this notebook (cell-by-cell running is higly recommanded on first time to control the process and adapt steps, variables and parameters to your own needs).\n",
    "    \n",
    "*Note: This script was developed using Windows7 x64, Anaconda 2 with Python 2.7, GRASS 7.3, R 3.3.0 (Caret package v. 6.0-70).*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some useful resources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- [Book: Learning IPython for interactive computing and data visualisation](http://it-ebooks.info/book/7021/) (previous name of Jupyter notebook)\n",
    "- [Wiki: GRASS and Python](https://grasswiki.osgeo.org/wiki/GRASS_and_Python)\n",
    "- [Wiki: GRASS Python scripting library](https://grasswiki.osgeo.org/wiki/GRASS_Python_Scripting_Library)\n",
    "\n",
    "\n",
    "- For a nice and easy first view of possibilities of GRASS scripting usin Python, see this [workshop video](http://www.youtube.com/watch?feature=player_embedded&v=PX2UpMhp2hc) on Youtube \n",
    "<a href=\"http://www.youtube.com/watch?feature=player_embedded&v=PX2UpMhp2hc\n",
    "\" target=\"_blank\"><img src=\"http://img.youtube.com/vi/PX2UpMhp2hc/0.jpg\" \n",
    "alt=\"IMAGE ALT TEXT HERE\" width=\"240\" height=\"180\" border=\"10\" /></a>\n",
    "\n",
    "\n",
    "- For more information about the GRASS GIS, please refer to  [Neteler and Mitasova, 2008](http://link.springer.com/book/10.1007%2F978-0-387-68574-8)\n",
    "![alt text](https://static-content.springer.com/cover/book/978-0-387-68574-8.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define the working environment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cells are used to: \n",
    "- Import required libraries\n",
    "- Set the environment variables for Python, Anaconda, GRASS GIS and R \n",
    "- Define the [\"GRASSDATA\" folder](https://grass.osgeo.org/grass73/manuals/helptext.html), along with the name of the \"location\" and the \"mapset\" in which you will work."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Import libraries**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Import libraries needed for setting parameters of operating system \n",
    "import os\n",
    "import sys\n",
    "\n",
    "## Import library for temporary files creation \n",
    "import tempfile \n",
    "\n",
    "## Import Pandas library\n",
    "import pandas as pd\n",
    "\n",
    "## Import Numpy library\n",
    "import numpy as np\n",
    "\n",
    "## Import subprocess\n",
    "import subprocess\n",
    "\n",
    "## Import multiprocessing\n",
    "import multiprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Set 'Python' and 'GRASS GIS' environment variables**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we set [the environment variables allowing to use of GRASS GIS](https://grass.osgeo.org/grass64/manuals/variables.html) inside this Jupyter notebook. Please modify the directory paths, so that they match your own system configuration. \n",
    "\n",
    "If you are working on Windows, with the GRASS GIS [stand-alone installation](https://grass.osgeo.org/download/software/ms-windows/), the paths displayed below should be similar. \n",
    "\n",
    "The setting of environmental variables could be improved like proposed on [this GRASS wiki page](https://grasswiki.osgeo.org/wiki/Working_with_GRASS_without_starting_it_explicitly#Python:_GRASS_GIS_7_without_existing_location_using_metadata_only)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Define path to the 'grass7x.bat' file\n",
    "grass7bin_win = 'C:\\\\Program Files\\\\GRASS GIS 7.3.svn\\\\grass73svn.bat'\n",
    "\n",
    "## Define GRASS GIS environment variables\n",
    "os.environ['GISBASE'] = 'C:\\\\Program Files\\\\GRASS GIS 7.3.svn'\n",
    "os.environ['PATH'] = 'C:\\\\Program Files\\\\GRASS GIS 7.3.svn\\\\lib;C:\\\\Program Files\\\\GRASS GIS 7.3.svn\\\\bin;C:\\\\Program Files\\\\GRASS GIS 7.3.svn\\\\extrabin' + os.pathsep + os.environ['PATH']\n",
    "os.environ['PATH'] = 'C:\\\\Program Files\\\\GRASS GIS 7.3.svn\\\\etc;C:\\\\Program Files\\\\GRASS GIS 7.3.svn\\\\etc\\\\python;C:\\\\Python27' + os.pathsep + os.environ['PATH']\n",
    "os.environ['PATH'] = 'C:\\\\Program Files\\\\GRASS GIS 7.3.svn\\\\Python27;C:\\\\Users\\\\Admin_ULB\\\\AppData\\\\Roaming\\\\GRASS7\\\\addons\\\\scripts' + os.pathsep + os.environ['PATH']\n",
    "os.environ['PATH'] = 'C:\\\\Program Files\\\\Anaconda2\\\\lib\\\\site-packages' + os.pathsep + os.environ['PATH']\n",
    "os.environ['PYTHONLIB'] = 'C:\\\\Python27'\n",
    "os.environ['PYTHONPATH'] = 'C:\\\\Program Files\\\\GRASS GIS 7.3.svn\\\\etc\\\\python'\n",
    "os.environ['GIS_LOCK'] = '$$'\n",
    "os.environ['GISRC'] = 'C:\\\\Users\\\\Admin_ULB\\\\AppData\\\\Roaming\\\\GRASS7\\\\rc'\n",
    "os.environ['GDAL_DATA'] = 'C:\\\\Program Files\\\\GRASS GIS 7.3.svn\\\\share\\\\gdal'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Define GRASS-Python environment\n",
    "sys.path.append(os.path.join(os.environ['GISBASE'],'etc','python'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please notice that paths will differ if you installed GRASS through the 'OSGeo4W package'. Here are some identified environment variables to use with a OSGeo4W installation: \n",
    "- grass7bin_win = 'C:\\\\OSGeo4W64\\\\bin\\\\grass73svn.bat'\n",
    "- os.environ['GISBASE'] = 'C:\\\\OSGeo4W64\\\\apps\\\\grass\\\\grass-7.3.svn'\n",
    "- os.environ['PATH'] = 'C:\\\\OSGeo4W64\\\\bin' + os.pathsep + os.environ['PATH']\n",
    "- os.environ['PYTHONLIB'] = 'C:\\\\OSGeo4W64\\\\apps\\\\Python27'\n",
    "- os.environ['GDAL_DATA'] = 'C:\\\\OSGeo4W64\\\\share\\\\gdal'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Set 'R statistical computing software' environment variables**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we set [the environment variables allowing to use the R statistical computing software](https://stat.ethz.ch/R-manual/R-devel/library/base/html/EnvVar.html) inside this Jupyter notebook. Please change the directory path to match your system configuration. If you are working on Windows, the paths below should be similar. \n",
    "\n",
    "Please notice that you will probably have to set the path of R_LIBS_USER also directly in R interface. For that, open R software (or [Rstudio software](https://www.rstudio.com/)) and enter the following command in the command prompt (you should adapt this path to match your own configuration: **.libPaths('C:\\\\R_LIBS_USER\\\\win-library\\\\3.3')**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Add the R software directory to the general PATH\n",
    "os.environ['PATH'] = 'C:\\\\Program Files\\\\R\\\\R-3.3.0\\\\bin' + os.pathsep + os.environ['PATH']\n",
    "\n",
    "## Set R software specific environment variables\n",
    "os.environ['R_HOME'] = 'C:\\Program Files\\R\\R-3.3.0'\n",
    "os.environ['R_ENVIRON'] = 'C:\\Program Files\\R\\R-3.3.0\\etc\\x64'\n",
    "os.environ['R_DOC_DIR'] = 'C:\\Program Files\\R\\R-3.3.0\\doc'\n",
    "os.environ['R_LIBS'] = 'C:\\Program Files\\R\\R-3.3.0\\library'\n",
    "os.environ['R_LIBS_USER'] = 'C:\\R_LIBS_USER\\win-library\\\\3.3'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Display current environment variables of your computer**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Display the current defined environment variables\n",
    "for key in os.environ.keys():\n",
    "    print \"%s = %s \\t\" % (key,os.environ[key])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# User inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Define a empty dictionnary for saving user inputs\n",
    "user={}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here after:\n",
    "- Enter the path to the directory you want to use as \"[GRASSDATA](https://grass.osgeo.org/programming7/loc_struct.png)\". \n",
    "- Enter the name of the location in which you want to work and its projection information in [EPSG code](http://spatialreference.org/ref/epsg/) format. Please note that the GRASSDATA folder and locations will be automatically created if they do not yet exist. If the location name already exists, the projection information will not be used.  \n",
    "- Enter the name you want for the mapsets which will be used later for Unsupervised Segmentation Parameter Optimization (USPO), Segmentation and Classification steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Enter the path to GRASSDATA folder\n",
    "user[\"gisdb\"] = \"F:\\\\myusername\\\\GRASSDATA\"\n",
    "\n",
    "## Enter the name of the location (existing or for a new one)\n",
    "user[\"location\"] = \"mycityname_32630\"\n",
    "\n",
    "## Enter the EPSG code for this location \n",
    "user[\"locationepsg\"] = \"32630\"\n",
    "\n",
    "## Enter the name of the mapset to use for Unsupervised Segmentation Parameter Optimization (USPO) step\n",
    "user[\"uspo_mapsetname\"] = \"TEST_USPO\"\n",
    "\n",
    "## Enter the name of the mapset to use for segmentation step\n",
    "user[\"segmentation_mapsetname\"] = \"TEST_SEGMENT\"\n",
    "\n",
    "## Enter the name of the mapset to use for classification step\n",
    "user[\"classification_mapsetname\"] = \"TEST_CLASSIF\"\n",
    "\n",
    "## Enter the maximum number of processes to run in parallel\n",
    "user[\"nb_proc\"] = 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if user[\"nb_proc\"] > multiprocessing.cpu_count():\n",
    "    print \"The requiered number of cores is higher than the amount available. Please fix it\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define the GRASSDATA folder and create GRASS location and mapsets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here after, the python script will check if the GRASSDATA folder, locations and mapsets already exist. If not, they will be automatically created."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Import GRASS Python packages**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Import libraries needed to launch GRASS GIS in the jupyter notebook\n",
    "import grass.script.setup as gsetup\n",
    "\n",
    "## Import libraries needed to call GRASS using Python\n",
    "import grass.script as grass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Define GRASSDATA folder and create location and mapsets**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Automatic creation of GRASSDATA folder\n",
    "if os.path.exists(user[\"gisdb\"]):\n",
    "    print \"GRASSDATA folder already exists\" \n",
    "else: \n",
    "    os.makedirs(user[\"gisdb\"]) \n",
    "    print \"GRASSDATA folder created in \"+user[\"gisdb\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Automatic creation of GRASS location if it doesn't exist\n",
    "if os.path.exists(os.path.join(user[\"gisdb\"],user[\"location\"])):\n",
    "    print \"Location \"+user[\"location\"]+\" already exists\" \n",
    "else : \n",
    "    if sys.platform.startswith('win'):\n",
    "        grass7bin = grass7bin_win\n",
    "        startcmd = grass7bin + ' -c epsg:' + user[\"locationepsg\"] + ' -e ' + os.path.join(user[\"gisdb\"],user[\"location\"])\n",
    "        p = subprocess.Popen(startcmd, shell=True, \n",
    "                             stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "        out, err = p.communicate()\n",
    "        if p.returncode != 0:\n",
    "            print >>sys.stderr, 'ERROR: %s' % err\n",
    "            print >>sys.stderr, 'ERROR: Cannot generate location (%s)' % startcmd\n",
    "            sys.exit(-1)\n",
    "        else:\n",
    "            print 'Created location %s' % os.path.join(user[\"gisdb\"],user[\"location\"])\n",
    "    else:\n",
    "        print 'This notebook was developed for use with Windows. It seems you are using another OS.'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Automatic creation of GRASS GIS mapsets\n",
    "\n",
    "## Import library for file copying \n",
    "import shutil\n",
    "\n",
    "## USPO mapset\n",
    "mapsetname=user[\"uspo_mapsetname\"]\n",
    "if os.path.exists(os.path.join(user[\"gisdb\"],user[\"location\"],mapsetname)):\n",
    "    print \"'\"+mapsetname+\"' mapset already exists\" \n",
    "else: \n",
    "    os.makedirs(os.path.join(user[\"gisdb\"],user[\"location\"],mapsetname))\n",
    "    shutil.copy(os.path.join(user[\"gisdb\"],user[\"location\"],'PERMANENT','WIND'),\n",
    "                os.path.join(user[\"gisdb\"],user[\"location\"],mapsetname,'WIND'))\n",
    "    print \"'\"+mapsetname+\"' mapset created in location \"+user[\"gisdb\"]\n",
    "\n",
    "## SEGMENTATION mapset\n",
    "mapsetname=user[\"segmentation_mapsetname\"]\n",
    "if os.path.exists(os.path.join(user[\"gisdb\"],user[\"location\"],mapsetname)):\n",
    "    print \"'\"+mapsetname+\"' mapset already exists\" \n",
    "else: \n",
    "    os.makedirs(os.path.join(user[\"gisdb\"],user[\"location\"],mapsetname))\n",
    "    shutil.copy(os.path.join(user[\"gisdb\"],user[\"location\"],'PERMANENT','WIND'),\n",
    "                os.path.join(user[\"gisdb\"],user[\"location\"],mapsetname,'WIND'))\n",
    "    print \"'\"+mapsetname+\"' mapset created in location \"+user[\"gisdb\"]\n",
    "\n",
    "## CLASSIFICATION mapset\n",
    "mapsetname=user[\"classification_mapsetname\"]\n",
    "if os.path.exists(os.path.join(user[\"gisdb\"],user[\"location\"],mapsetname)):\n",
    "    print \"'\"+mapsetname+\"' mapset already exists\" \n",
    "else: \n",
    "    os.makedirs(os.path.join(user[\"gisdb\"],user[\"location\"],mapsetname))\n",
    "    shutil.copy(os.path.join(user[\"gisdb\"],user[\"location\"],'PERMANENT','WIND'),\n",
    "                os.path.join(user[\"gisdb\"],user[\"location\"],mapsetname,'WIND'))\n",
    "    print \"'\"+mapsetname+\"' mapset created in location \"+user[\"gisdb\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section of the notebook is dedicated to defining functions which will then be called later in the script. If you want to create your own functions, define them here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Function for computing processing time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"print_processing_time\" function is used to calculate and display the processing time for various stages of the processing chain. At the beginning of each major step, the current time is stored in a new variable, using [time.time() function](https://docs.python.org/2/library/time.html). At the end of the stage in question, the \"print_processing_time\" function is called and takes as an argument, the name of this new variable containing the recorded time at the beginning of the stage, and an output message."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Import library for managing time in python\n",
    "import time  \n",
    "\n",
    "## Function \"print_processing_time()\" compute processing time and print it.\n",
    "# The argument \"begintime\" wait for a variable containing the begintime (result of time.time()) of the process for which to compute processing time.\n",
    "# The argument \"printmessage\" wait for a string format with information about the process. \n",
    "def print_processing_time(begintime, printmessage):    \n",
    "    endtime=time.time()           \n",
    "    processtime=endtime-begintime\n",
    "    remainingtime=processtime\n",
    "\n",
    "    days=int((remainingtime)/86400)\n",
    "    remainingtime-=(days*86400)\n",
    "    hours=int((remainingtime)/3600)\n",
    "    remainingtime-=(hours*3600)\n",
    "    minutes=int((remainingtime)/60)\n",
    "    remainingtime-=(minutes*60)\n",
    "    seconds=round((remainingtime)%60,1)\n",
    "\n",
    "    if processtime<60:\n",
    "        finalprintmessage=str(printmessage)+str(seconds)+\" seconds\"\n",
    "    elif processtime<3600:\n",
    "        finalprintmessage=str(printmessage)+str(minutes)+\" minutes and \"+str(seconds)+\" seconds\"\n",
    "    elif processtime<86400:\n",
    "        finalprintmessage=str(printmessage)+str(hours)+\" hours and \"+str(minutes)+\" minutes and \"+str(seconds)+\" seconds\"\n",
    "    elif processtime>=86400:\n",
    "        finalprintmessage=str(printmessage)+str(days)+\" days, \"+str(hours)+\" hours and \"+str(minutes)+\" minutes and \"+str(seconds)+\" seconds\"\n",
    "    \n",
    "    return finalprintmessage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "begintime_full=time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1 - Importing data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usually, original data are imported and stored in the \"PERMANENT\" mapset (automatically created when creating a new location)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Launch GRASS GIS working session**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Save the name of the mapset in which to import the data\n",
    "mapsetname='PERMANENT'\n",
    "\n",
    "## Launch GRASS GIS working session in the mapset\n",
    "if os.path.exists(os.path.join(user[\"gisdb\"],user[\"location\"],mapsetname)):\n",
    "    gsetup.init(os.environ['GISBASE'], user[\"gisdb\"], user[\"location\"], mapsetname)\n",
    "    print \"You are now working in mapset '\"+mapsetname+\"'\" \n",
    "else: \n",
    "    print \"'\"+mapsetname+\"' mapset doesn't exists in \"+user[\"gisdb\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "begintime_importingdata=time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import raw data in PERMANENT mapset "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For optical and nDSM import, please adapt the input of the ['r.in.gdal' commands](https://grass.osgeo.org/grass73/manuals/r.in.gdal.html) to match your own data location."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import optical raster imagery "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please adapt the number of lines in the loop to match the number of layers stacked in your imagery file (ensuring the number of ['g.rename' commands](https://grass.osgeo.org/grass73/manuals/g.rename.html) equal the number of layers stacked). Ensure the names of the layers match their position in the stack (e.g. \"opt_blue\" for the first layer).\n",
    "\n",
    "Please note that it is assumed that your data has at least a red band layer, called \"opt_red\". If not, you will have to change several parameters through this notebook, notably when defining [computation region](https://grasswiki.osgeo.org/wiki/Computational_region). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "begintime_optical=time.time()\n",
    "\n",
    "## Import optical imagery and rename band with color name\n",
    "print (\"Importing optical raster imagery at \" + time.ctime())\n",
    "\n",
    "grass.run_command('r.in.gdal', input=\"F:\\\\....\\\\.....\\\\mosaique_georef_ordre2.tif\", output=\"optical\", overwrite=True)\n",
    "for rast in grass.list_strings(\"rast\"):\n",
    "    if rast.find(\"1\")!=-1: grass.run_command(\"g.rename\", overwrite=True, rast=(rast,\"opt_blue\"))\n",
    "    elif rast.find(\"2\")!=-1: grass.run_command(\"g.rename\", overwrite=True, rast=(rast,\"opt_green\"))\n",
    "    elif rast.find(\"3\")!=-1: grass.run_command(\"g.rename\", overwrite=True, rast=(rast,\"opt_red\"))\n",
    "    elif rast.find(\"4\")!=-1: grass.run_command(\"g.rename\", overwrite=True, rast=(rast,\"opt_nir\"))\n",
    "        \n",
    "print_processing_time(begintime_optical ,\"Optical imagery has been imported in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import nDSM raster imagery "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have null value in your nDSM raster, please be careful to define the \"setnull\" parameter of ['r.null' command](https://grass.osgeo.org/grass73/manuals/r.null.html) well, according to your own data. If you didn't have any null values in your nDSM raster, simply comment the r.null command line with an '#' as first character (to put it in [comment](http://www.pythonforbeginners.com/comments/comments-in-python)). Notice that you can display the line's number by pressing the L key when cell edge is in blue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "begintime_ndsm=time.time()\n",
    "\n",
    "## Import nDSM imagery \n",
    "print (\"Importing nDSM raster imagery at \" + time.ctime())\n",
    "grass.run_command('r.in.gdal', input=\"F:\\\\MAUPP\\\\.....\\\\Orthorectified\\\\mosaique_georef\\\\nDSM\\\\nDSM_mosaik_georef_ordre2.tif\", output=\"ndsm\", overwrite=True)\n",
    "\n",
    "## Define null value for specific value in nDSM raster. Adapt the value to your own data. \n",
    "# If there is no null value in your data, comment the next line\n",
    "grass.run_command('r.null', map=\"ndsm\", setnull=\"-999\")\n",
    "\n",
    "print_processing_time(begintime_ndsm, \"nDSM imagery has been imported in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Update imagery group \"optical\" with optical rasters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In GRASS GIS several operations, mainly segmentation when dealing with OBIA, require an 'imagery group'.\n",
    "\n",
    "In the next cell, a new imagery group called 'optical' containing the optical raster layers is created with ['i.group command'](https://grass.osgeo.org/grass73/manuals/i.group.html). It is impossible to create a new imagery group, if the name already exists, therefore existing imagery groups with the same name are removed (with ['g.remove command'](https://grass.osgeo.org/grass73/manuals/g.remove.html))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print (\"Updating imagery group 'optical' with optical rasters at \" + time.ctime())\n",
    "\n",
    "## Remove existing imagery group nammed \"optical\". This group was created when importing multilayer raster data\n",
    "grass.run_command(\"g.remove\", type=\"group\", name=\"optical\", flags=\"f\")\n",
    "\n",
    "## Add each raster which begin with the prefix \"opt\" into a new imagery group \"optical\"\n",
    "for rast in grass.list_strings(\"rast\", pattern=\"opt\", flag=\"r\"):\n",
    "    grass.run_command(\"i.group\", group=\"optical\", input=rast)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save default GRASS GIS' computational region for the whole extent of optical imagery"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In GRASS GIS, the concept of computational region is fundamental. We highly recommend reading [information about the computation region in GRASS GIS](https://grasswiki.osgeo.org/wiki/Computational_region) to be sure to understand the concept. \n",
    "\n",
    "Here after, the 'default' computational region is defined as corresponding to the red band image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Save default computational region to match the full extend of optical imagery\n",
    "grass.run_command('g.region', flags=\"s\", raster=\"opt_red@PERMANENT\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute pseudo-band raster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set computational region and mask layer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ['r.mask' command](https://grass.osgeo.org/grass73/manuals/r.mask.html) is used not to perform further processing on 'nodata' pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Set computational region to default\n",
    "grass.run_command('g.region', flags=\"d\")\n",
    "\n",
    "## Apply mask to not compute out-of-AOI pixels \n",
    "grass.run_command('r.mask', overwrite=True, raster=\"opt_red\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute NDVI (Normalized difference vegetation index) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here after, we use the [r.mapcalc command](https://grass.osgeo.org/grass73/manuals/r.mapcalc.html) to compute NDVI (Normalized difference vegetation index) indice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "print (\"Begin compute NDVI on \"+time.ctime())\n",
    "begintime_ndvi=time.time()\n",
    "\n",
    "## Compute NDVI\n",
    "formula=\"NDVI=(float(opt_nir@PERMANENT)-float(opt_red@PERMANENT))/(float(opt_nir@PERMANENT)+float(opt_red@PERMANENT))\"\n",
    "grass.mapcalc(formula, overwrite=True)\n",
    "    \n",
    "print_processing_time(begintime_ndvi, \"NDVI has been computed in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute NDWI (Normalized difference water index) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here after, we use the [r.mapcalc command](https://grass.osgeo.org/grass73/manuals/r.mapcalc.html) to compute the normalized difference water index (NDWI) as indice. The formula used was proposed by [McFeeters in 1996](http://www.tandfonline.com/doi/abs/10.1080/01431169608948714)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Saving current time for processing time management\n",
    "print (\"Begin compute NDWI on \"+time.ctime())\n",
    "begintime_ndvi=time.time()\n",
    "\n",
    "## Compute NDVI\n",
    "formula=\"NDWI=(float(opt_green@PERMANENT)-float(opt_nir@PERMANENT))/(float(opt_green@PERMANENT)+float(opt_nir@PERMANENT))\"\n",
    "grass.mapcalc(formula, overwrite=True)\n",
    "    \n",
    "print_processing_time(begintime_ndvi, \"NDWI has been computed in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute brightness "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here after, we use the [r.mapcalc command](https://grass.osgeo.org/grass73/manuals/r.mapcalc.html) to compute the brightness indice defined as the sum of visible bands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Saving current time for processing time management\n",
    "print (\"Begin compute brightness on \"+time.ctime())\n",
    "begintime_brightness=time.time()\n",
    "\n",
    "## Compute Brightness\n",
    "formula=\"Brightness=opt_blue@PERMANENT+opt_green@PERMANENT+opt_red@PERMANENT\"\n",
    "grass.mapcalc(formula, overwrite=True)\n",
    "\n",
    "print_processing_time(begintime_brightness, \"Brightness has been computed in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute texture "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we use ['r.texture' command](https://grass.osgeo.org/grass73/manuals/r.texture.html) to compute angular second moment with 5x5 moving window. This texture layer is not currently being used in the process, but you can adapt the script if you want to use it. Other textures can also be computed using the 'r.texture' command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Saving current time for processing time management\n",
    "print (\"Computing a 5x5 window angular second moment texture on \"+time.ctime())\n",
    "begintime_texture=time.time()\n",
    "\n",
    "## Compute Angular second moment texture\n",
    "grass.run_command('r.texture', overwrite=True, input=\"Brightness\", output=\"texture\", method=\"asm\", size=\"5\")\n",
    "\n",
    "print_processing_time(begintime_texture, \"Angular second moment texture has been computed in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remove current mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Check if there is a raster layer named \"MASK\"\n",
    "if not grass.list_strings(\"rast\", pattern=\"MASK\", flag='r'):\n",
    "    print 'There is currently no MASK'\n",
    "else:\n",
    "    ## Remove the current MASK layer\n",
    "    grass.run_command('r.mask',flags='r')\n",
    "    print 'The current MASK has been removed'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### End of part 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(\"Importation of data ends at \"+ time.ctime())\n",
    "print_processing_time(begintime_importingdata, \"Importation of data has been achieved in :\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-* </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-* </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2 - Unsupervised segmentation parameter optimization (USPO)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Launch GRASS GIS working session**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Set the name of the mapset in which to work\n",
    "mapsetname=user[\"uspo_mapsetname\"]\n",
    "\n",
    "## Launch GRASS GIS working session in the mapset\n",
    "if os.path.exists(os.path.join(user[\"gisdb\"],user[\"location\"],mapsetname)):\n",
    "    gsetup.init(os.environ['GISBASE'], user[\"gisdb\"], user[\"location\"], mapsetname)\n",
    "    print \"You are now working in mapset '\"+mapsetname+\"'\" \n",
    "else: \n",
    "    print \"'\"+mapsetname+\"' mapset doesn't exists in \"+user[\"gisdb\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "begintime_USPO_full=time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create USPO's working regions based on several image subsets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The ['i.segment.uspo' add-on](https://grass.osgeo.org/grass70/manuals/addons/i.segment.uspo.html) could use the whole area of your data for segmentation parameters optimization, but it could take days depending on the size of your data and the number of parameter combinations you ask for testing. For this reason, we use small rectangular polygons (recommended size: fewer than 4 million pixels), representing the diversity of landscapes in your scene, for which i.segment.uspo will find optimized segmentation parameters. Please refer to the \"region\" parameter of the [i.segment.uspo help](https://grass.osgeo.org/grass70/manuals/addons/i.segment.uspo.html) for more explanations. Please create this polygon layer (shapefile) focusing on subsets representing the diversity of landscapes in your scene. You could use [QuantumGIS](http://www.qgis.org/en/site/) to perform this step.\n",
    "\n",
    "- Please adapt the \"input\" parameter of the 'v.import' command to match the path to your own data. \n",
    "\n",
    "- We call here \"USPO's regions\" the GRASS's computational regions where i.segment.uspo will perform segmentation and compute optimization function in order to find optimized segmentation parameter(s). \n",
    "\n",
    "- The ['v.import' command](https://grass.osgeo.org/grass72/manuals/v.import.html) is used to import vector data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Import shapefile with polygons corresponding to computational region's extension for USPO \n",
    "print (\"Importing vector data with USPO's regions at \" + time.ctime())\n",
    "\n",
    "grass.run_command('g.region', flags=\"d\")\n",
    "\n",
    "grass.run_command('v.import', overwrite=True, \n",
    "                  input=\"F:/...../region_USPO/Ouaga_region_USPO.shp\", \n",
    "                  output=\"region_uspo\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the \"region_uspo\" layer contains polygons corresponding to the computational regions to be used in i.segment.uspo, this layer is used here to define (and save) a GRASS computational region for each polygon. We use here [v.extract command](https://grass.osgeo.org/grass72/manuals/v.extract.html) to extract each polygon temporarily from the \"region uspo\" layer,  [g.region command](https://grass.osgeo.org/grass72/manuals/g.region.html) to create computational region corresponding to the polygon and save it with a specific name and [g.remove command](https://grass.osgeo.org/grass72/manuals/g.remove.html) to remove the temporarily created polygon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Create a computional region for each polygon in the 'region_uspo' layer\n",
    "print (\"Defining a GRASS region for each polygon at \" + time.ctime())\n",
    "\n",
    "for cat in grass.parse_command('v.db.select', map='region_uspo',  columns='cat', flags='c'):\n",
    "    condition=\"cat=\"+cat\n",
    "    outputname=\"region_uspo_\"+cat\n",
    "    regionname=\"subset_uspo_\"+cat\n",
    "    grass.run_command('v.extract', overwrite=True, quiet=True, \n",
    "                      input=\"region_uspo\", type=\"area\", where=condition, output=outputname)\n",
    "    grass.run_command('g.region', overwrite=True, vector=outputname, save=regionname, align=\"opt_red@PERMANENT\", flags=\"u\")\n",
    "    grass.run_command('g.remove', type=\"vector\", name=outputname, flags=\"f\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unsupervised Segmentation Parameter Optimization with i.segment.uspo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [i.segment.uspo add-on](https://grass.osgeo.org/grass70/manuals/addons/i.segment.uspo.html) is used to automatically find the optimized segmentation parameter(s) for specific computational regions. We highly recommend reading of [the add-on help](https://grass.osgeo.org/grass70/manuals/addons/i.segment.uspo.html) to be sure to understand the different fonctionalities well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define a imagery group for i.segment.uspo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to create an imagery group for i.segment.uspo. In GRASS GIS, these imagery groups contain raster layers to use to a speficic task. Here, for i.segment.uspo, we define and use the same imagery group of layers for both segmentation and parameter optimization. You could use different imagery groups for segmentation and parameter optimization (e.g. you could segment using only optical data and optimize the segmentation parameter based on a specific indice like NDVI i.e.). \n",
    "\n",
    "Here, we use the 4 optical layers and the NDVI layer for the segmentation and the USPO step. If you want to use other layers, please adapt the script accordingly. If you make changes, please be careful to use the same layer for the segmentation in i.segment.uspo and for further segmentation step (with i.segment). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "begintime_USPO_FULL=time.time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print (\"Defining a imagery group with raster used for i.segment.uspo at \" + time.ctime())\n",
    "\n",
    "## Remove existing imagery group named \"group\"\n",
    "grass.run_command('g.remove', flags=\"rf\", type=\"group\", name=\"group\")\n",
    " \n",
    "## Add all optical imagery in the imagery group\n",
    "for rast in grass.list_strings(\"rast\", pattern=\"opt\", flag='r'):\n",
    "    grass.run_command('i.group', group=\"group\", input=rast)\n",
    "\n",
    "## Add NDVI imagery in the imagery group\n",
    "grass.run_command('i.group', group=\"group\", input=\"NDVI@PERMANENT\")\n",
    "\n",
    "## list files in the group\n",
    "print grass.read_command('i.group', group=\"group\", flags=\"l\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instal GRASS extensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GRASS GIS have both a core part (the one installed by default on your computer) and add-ons (which have to be installed using the extension manager ['g.extension'](https://grass.osgeo.org/grass72/manuals/g.extension.html)).\n",
    "\n",
    "In the next cell, 'i.segment.uspo' will be installed (if not yet) and also other add-ons ['r.neighborhoodmatrix'](https://grass.osgeo.org/grass70/manuals/addons/r.neighborhoodmatrix.html) and ['i.segment.hierarchical'](https://grass.osgeo.org/grass70/manuals/addons/i.segment.hierarchical.html) required for running i.segment.uspo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Instal r.neighborhoodmatrix if not yet installed\n",
    "if \"r.neighborhoodmatrix\" not in grass.parse_command('g.extension', flags=\"a\"):\n",
    "    grass.run_command('g.extension', extension=\"r.neighborhoodmatrix\")\n",
    "    print \"r.neighborhoodmatrix have been installed on your computer\"\n",
    "else: print \"r.neighborhoodmatrix is already installed on your computer\" \n",
    "\n",
    "## Instal i.segment.hierarchical if not yet installed\n",
    "if \"i.segment.hierarchical\" not in grass.parse_command('g.extension', flags=\"a\"):\n",
    "    grass.run_command('g.extension', extension=\"i.segment.hierarchical\")\n",
    "    print \"i.segment.hierarchical have been installed on your computer\"\n",
    "else: print \"i.segment.hierarchical is already installed on your computer\" \n",
    "\n",
    "## Instal i.segment.uspo if not yet installed\n",
    "if \"i.segment.uspo\" not in grass.parse_command('g.extension', flags=\"a\"):\n",
    "    grass.run_command('g.extension', extension=\"i.segment.uspo\")\n",
    "    print \"i.segment.uspo have been installed on your computer\"\n",
    "else: print \"i.segment.uspo is already installed on your computer\" "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unsupervised segmentation parameter optimisation with i.segment.uspo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Please read the [i.segment.uspo](https://grass.osgeo.org/grass70/manuals/addons/i.segment.uspo.html) help to ensure you understand the extension well and you choose correctly the different parameters for your case. \n",
    "- It is recommended to quickly identify (by visual check) thresholds resulting in clearly oversegmentedand undersegmented segments, and use those thresholds as 'threshold_start' and 'threshold_start', respectively. \n",
    "- If the number of threshold/minsize combination to test is too big, running i.segment.uspo could be very time-consuming or even failed. In this case, please reduce the range and/or steps of thresholds and/or minsize to be tested. \n",
    "- Choose between the following optimization function: \"sum\" of [Espindola (2006)](http://www.tandfonline.com/doi/abs/10.1080/01431160600617194) or \"F function\" of [Johnson (2015)](http://www.mdpi.com/2220-9964/4/4/2292). \n",
    "- If you choose the \"F function\", please use an adapted alpha parameter. A value of alpha>1 is *\"(...) appropriate for selecting the parameters for finer segmentation levels (to ensure that smaller objects of interest or objects spectrally-similar to their surroundings are not undersegmented at these levels), while values of a[alpha] < 1 may be more appropriate for selecting parameters for coarser segmentation levels (to ensure that larger/more heterogeneous objects of interest are not oversegmented at these levels)\"* (Johnson et al., 2015, pp. 2295).\n",
    "- Please notice that, for our requirements, the minsize parameter was set according to our MMU (Minimum Mapping Unit) and was not optimized with i.segment.uspo. You can change the script if you want but notice that you should then add a step further to select the optimized minsize as it is done currently for the threshold. \n",
    "- Please notice that the RAM allowed to the segmentation has been set to 2Gb but can be modified according to the capacity of your own system. The same is true for the number of processors to use. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Define the optimization function name (\"sum\" or \"f\")\n",
    "opti_f=\"f\"\n",
    "\n",
    "## Define the alpha, only if selected optimization function is \"f\"\n",
    "if opti_f==\"f\":\n",
    "    alpha=1.25\n",
    "else : alpha=\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next cell, please adapt the path to the directory where you want to save the .csv output of i.segment.uspo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Define the csv output file name, according to the optimization function selected\n",
    "outputcsv=\"F:\\\\.....\\\\Segmentation_param\\\\ouaga_uspo_\"+str(opti_f)+str(alpha)+\".csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Defining a list of GRASS GIS' computational regions where i.segment.uspo will optimize the segmentation parameters\n",
    "regions_uspo=grass.list_strings(\"region\", pattern=\"subset_uspo_\", flag='r')[0]\n",
    "for region in grass.list_strings(\"region\", pattern=\"subset_uspo_\", flag='r')[1:]:\n",
    "    regions_uspo+=\",\"+region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Running i.segment.uspo\n",
    "print (\"Runing i.segment.uspo at \" + time.ctime())\n",
    "begintime_USPO=time.time()\n",
    "\n",
    "grass.run_command('i.segment.uspo', overwrite=True, group='group', \n",
    "                  output=outputcsv, segment_map=\"best\", \n",
    "                  regions=regions_uspo, threshold_start=\"0.001\", threshold_stop=\"0.03\", threshold_step=\"0.001\", minsizes=\"8\", \n",
    "                  optimization_function=opti_f, f_function_alpha=alpha, memory=\"2000\", processes=str(user[\"nb_proc\"]))\n",
    "\n",
    "## Create a .csvt file containing each colomn type of i.segment.uspo' csv output. Required for further import of .csv file\n",
    "model_output_desc = outputcsv + \"t\"\n",
    "f = open(model_output_desc, 'w')\n",
    "header_string = '\"String\",\"Real\",\"Integer\",\"Real\",\"Real\",\"Real\"'\n",
    "f.write(header_string)\n",
    "f.close()\n",
    "\n",
    "## Print\n",
    "print_processing_time(begintime_USPO, \"USPO process achieved in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Export of i.segment.uspo results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section is optional and can be used to export the segmentation results of i.segment.uspo in shapefile format. ['r.to.vect' command](https://grass.osgeo.org/grass72/manuals/r.to.vect.html) is used to vectorize the \"best\" segmentation results coming from i.segment.uspo. If you don't want to export those vectors and just visualize it in GRASS GIS, you can use the r.to.vect command. The ['v.out.ogr' command](https://grass.osgeo.org/grass72/manuals/v.out.ogr.html) is used to perform the export of vector layers in conventional vector formats like shapefile."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert i.segment.uspo raster outputs in vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please adapt the path to the directory where you want to save the segmentation results of i.segment.uspo, in shapefile format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Define the output folder name, according to the optimization function selected\n",
    "outputfolder=\"F:\\\\.....\\\\Segmentation_param\\\\USPO_bestsegment_\"+str(opti_f)+str(alpha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "print (\"Begin to export i.segment.uspo results at \" + time.ctime())\n",
    "begintime_exportUSPO=time.time()\n",
    "\n",
    "count=0\n",
    "for rast in grass.list_strings(\"rast\", pattern=\"best_\", flag='r'):\n",
    "    count+=1\n",
    "    print \"Working on raster '\"+str(rast)+\"' - \"+str(count)+\"/\"+str(len(grass.list_strings(\"rast\", pattern=\"best_\", flag='r')))\n",
    "                                                 \n",
    "    strindex=rast.find(\"subset_uspo_\")\n",
    "    subregion=rast[strindex: strindex+14]\n",
    "    vectname=\"temp_bestsegment_\"+subregion\n",
    "    \n",
    "    print (\"Converting raster layer into vector\")\n",
    "    grass.run_command('r.to.vect', overwrite=True, input=rast, output=vectname, type='area')\n",
    "    \n",
    "    print (\"Exporting shapefile\")\n",
    "    grass.run_command('v.out.ogr', overwrite=True, input=vectname, type='area', \n",
    "                      output=outputfolder, format='ESRI_Shapefile')\n",
    "    \n",
    "    print (\"Remove newly created vector layer\")                                                      \n",
    "    grass.run_command(\"g.remove\", type=\"vector\", pattern=\"temp_bestsegment_\", flags=\"rf\")\n",
    "    \n",
    "## Print\n",
    "print_processing_time(begintime_exportUSPO, \"Export of i.segment.uspo results done in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### End of part 2 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(\"The script ends at \"+ time.ctime())\n",
    "print_processing_time(begintime_USPO_FULL, \"Entire process has been achieved in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-* </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-* </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3 - Segmentation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Launch GRASS GIS working session**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Set the name of the mapset in which to work\n",
    "mapsetname=user[\"segmentation_mapsetname\"]\n",
    "\n",
    "## Launch GRASS GIS working session in the mapset\n",
    "if os.path.exists(os.path.join(user[\"gisdb\"],user[\"location\"],mapsetname)):\n",
    "    gsetup.init(os.environ['GISBASE'], user[\"gisdb\"], user[\"location\"], mapsetname)\n",
    "    print \"You are now working in mapset '\"+mapsetname+\"'\" \n",
    "else: \n",
    "    print \"'\"+mapsetname+\"' mapset doesn't exists in \"+user[\"gisdb\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "begintime_segmentation_full=time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Import the csv output file from i.segment.uspo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Here, the results of i.segment.uspo are imported and manipulated to select, for each 'USPO's region', the segmentation parameter achieving the highest optimization score. Then, the threshold to be used to segment the whole scene is selected as the lowest threshold among the different \"best\" thresholds (notice that you can easily change it by using ['median'](https://docs.scipy.org/doc/numpy/reference/generated/numpy.median.html) function or  instead of ['amin'](https://docs.scipy.org/doc/numpy/reference/generated/numpy.amin.html) function of numpy library). \n",
    "- Be careful to control this step well, as selection of the lowest threshold could result in over-segmented results if the optimized thresholds are very different among the 'USPO's region'. In that case, median threshold could be preferred. \n",
    "- Please notice that the minsize have been fixed and not optimized with i.segment.uspo\n",
    "- If you made several tests, please be sure to import the .csv file corresponding to the wanted optimization function and alpha parameter !\n",
    "- Python's [Pandas](http://pandas.pydata.org/) library for managing .csv table in dataframe. We specifically used [.read_csv](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.read_csv.html), [.to_csv](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.to_csv.html), [.merge](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.merge.html), [.loc](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.loc.html), [.head](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.head.html), [.sort_values](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.sort_values.html), [.groupby](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.groupby.html), [.max](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.max.html) Pandas' dataframe functions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Import Pandas library\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Import the optimization results of i.segment.uspo in a dataframe\n",
    "print (\"Import .csv file with results of i.segment.uspo on \" + time.ctime())\n",
    "ouaga_uspo=pd.read_csv(outputcsv, sep=',',header=0)\n",
    "ouaga_uspo.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Create temporary dataframe with maximum value of optimization criteria for \"USPO's region\"\n",
    "temp=ouaga_uspo.loc[:,['region','optimization_criteria']].groupby('region').max()\n",
    "temp.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Merge between dataframes for identification of threshold corresponding to the maximum optimizaion criteria of each \"USPO's region\"\n",
    "uspo_parameters = pd.merge(ouaga_uspo, temp, on='optimization_criteria').loc[:,['region','threshold','optimization_criteria']].sort_values(by='region')\n",
    "uspo_parameters.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Save the optimized threshold of each \"USPO's region\" in a list\n",
    "uspo_parameters_list=uspo_parameters['threshold'].tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Save the minimum of optimized threshold in a new variable called \"optimized_threshold\"\n",
    "optimized_threshold=round(np.amin(uspo_parameters_list),3)\n",
    "print \"The lowest of the 'USPOs region' optimized threshold is \"+str(optimized_threshold)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Copy imagery group created for i.segment.uspo to the current mapset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the imagery group to use for segmentation (which is going to be performed with ['i.segment' module](https://grass.osgeo.org/grass72/manuals/i.segment.html)) should be in the working mapset, we copy the one created for i.segment.uspo in the previous step. We use ['g.copy' command](https://grass.osgeo.org/grass72/manuals/g.copy.html) for this purpose. ['i.group' command](https://grass.osgeo.org/grass72/manuals/i.group.html) is used to print, as a reminder, the list of raster in the imagery group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Copy of imagery group to the current mapset\n",
    "grass.run_command('g.copy', overwrite=True, group='group@USPO,group')\n",
    "print grass.read_command('i.group', group=\"group\", flags=\"l\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import processing tiles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to allow processing large areas, this processing chain is designed to work with polygons tiles. Please, create a polygon layer (shapefile) containing at least 2 tiles of your area of interest (it could have decades if you're dealing with very large dataset). Please add a column in the attribute table, called \"area_km2\" and containing the area of the tiles (in km2), which will be used to inform about the progress of the segmentation process. You can use [QuantumGIS](http://www.qgis.org/en/site/) to perform this step.\n",
    "\n",
    "Please adapt the \"input\" parameter of the 'v.import' command to match the path to your own data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print ('Importing processing tiles at '+ time.ctime())\n",
    "\n",
    "## Import vectorial tiles zones layers\n",
    "grass.run_command('v.import', overwrite=True, input=\"F:\\\\.....\\\\processing_tiles.shp\", output=\"processing_tiles\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Segmentation with optimized segmentation parameter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- For each processing tile, a segmentation is made (using ['i.segment'](https://grass.osgeo.org/grass72/manuals/i.segment.html)).\n",
    "\n",
    "- In the sequence of processes, each processing tile is extracted (with ['v.extract'](https://grass.osgeo.org/grass72/manuals/v.extract.html)), then is use to define the computational region (with ['g.region'](https://grass.osgeo.org/grass72/manuals/g.region.html)) and finally used as mask (with ['r.mask'](https://grass.osgeo.org/grass72/manuals/r.mask.html)). The processing tile is then segmented according to the optimized parameter, previously defined by i.segment.uspo. \n",
    "\n",
    "- Please notice that the \"minsize\" parameter of 'i.segment' command is here fixed, but you can adapt the script to your own need.\n",
    "\n",
    "- Please notice that the RAM allowed to the segmentation has been set to 2Gb but can be modified according to the capacity of your own system.\n",
    "\n",
    "- Note that it is necessary to use the *align* parameter when setting the \"computational region\" to match the polygon extention, as otherwise it could result in misalignment of segmentation results regardind to raster imagery.\n",
    "\n",
    "- We use also ['v.db.select'](https://grass.osgeo.org/grass72/manuals/v.db.select.html) and ['g.remove'](https://grass.osgeo.org/grass72/manuals/g.remove.html) modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Compute total area to be segmented for process progression information\n",
    "total_area=0\n",
    "processed_area=0\n",
    "\n",
    "for id in grass.parse_command('v.db.select', map='processing_tiles@SEGMENTATION', columns='cat', flags='c'):\n",
    "    condition='cat='+id\n",
    "    size=float(grass.read_command('v.db.select', map=\"processing_tiles@SEGMENTATION\", columns=\"area_km2\", where=condition,flags=\"c\"))\n",
    "    total_area+=size\n",
    "    \n",
    "print total_area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print (\"Begin segmentation process on \" + time.ctime())\n",
    "## Saving current time for processing time management\n",
    "begintime_segmentation=time.time()\n",
    "\n",
    "## Initialize a variable for process progression purpose \n",
    "processed_area=0.0\n",
    "    \n",
    "## Segmentation step. We use here a loop trhough the different polygons\n",
    "for id in grass.parse_command('v.db.select', map='processing_tiles@SEGMENTATION', columns='cat', flags='c'):\n",
    "    ## Save current time at loop' start. \n",
    "    begintime_current_id=time.time()\n",
    "    \n",
    "    ## Build condition for selection in attributes of vector layer\n",
    "    condition='cat='+id\n",
    "    \n",
    "    ## Save size of the current polygon \n",
    "    size=float(grass.read_command('v.db.select', map=\"processing_tiles@SEGMENTATION\", columns=\"area_km2\", where=condition,flags=\"c\"))\n",
    "\n",
    "    ## Build name for temporary vector layer used for computational region and mask definition\n",
    "    tempvector=\"temp_polygon_\"+id\n",
    "    \n",
    "    ## Extract the current polygon in a new vector layer\n",
    "    grass.run_command('v.extract', overwrite=True, input=\"processing_tiles@SEGMENTATION\", type=\"area\", where=condition, output=tempvector)\n",
    "    ## Define computational region to match the current polygon vector layer and align the computational reigion with optical imagery\n",
    "    grass.run_command('g.region', overwrite=True, vector=tempvector, align=\"opt_red@PERMANENT\")\n",
    "    ## Apply mask using the current polygon \n",
    "    grass.run_command('r.mask', overwrite=True, vector=tempvector)\n",
    "\n",
    "    \n",
    "    ## Segmentation of current polygon with i.segment\n",
    "    print (\"Segmenting tile number \"+str(id)+\" corresponding to \"+str(size)+\" km2\" )\n",
    "    outputsegment=\"segmentation_tile_\"+id\n",
    "    grass.run_command('i.segment', overwrite=True, group=\"group\", output=outputsegment, threshold=optimized_threshold, minsize=\"8\", memory=\"2000\")\n",
    "\n",
    "    ## Delete the current polygon vector layer\n",
    "    grass.run_command('g.remove', type=\"vector\", name=tempvector, flags=\"f\")\n",
    "    ## Remove current mask\n",
    "    grass.run_command('r.mask', flags=\"r\")\n",
    "    \n",
    "    ## Add size of the current polygon to the already processed area \n",
    "    processed_area+=size\n",
    "    \n",
    "    ## Print of what happened \n",
    "    print(\"Tile \"+str(id)+\" processed.\")\n",
    "    print_processing_time(begintime_current_id, \" Process achieved in \")\n",
    "    print (\"Progress = \"+str((processed_area/total_area)*100)+\" percent of the total area segmented\")   \n",
    "\n",
    "## Compute processing time and print it\n",
    "print_processing_time(begintime_segmentation, \"Segmentation process on all tiles achieved in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Merging each individual segmentation raster one \"patched\" raster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The different segmentation results of 'i.segment' (raster layers) for each processing tiles will be \"patched\" (merged) in one resulting raster. \n",
    "['r.mapcalc' command](https://grass.osgeo.org/grass71/manuals/r.mapcalc.html) is used to combine all the segmentation rasters together. The 'nmax' expression is used to keep the maximum value of input rasters, excluding the NULL values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "## Setting the computational region extend to all rasters to be merged \n",
    "groupraster=grass.list_strings(\"rast\", pattern=\"segmentation_tile_\", mapset=\"SEGMENTATION\", flag='r')[0]\n",
    "count=1\n",
    "for rast in grass.list_strings(\"rast\", pattern=\"segmentation_tile_\", mapset=\"SEGMENTATION\", flag='r')[1:]:\n",
    "    groupraster+=\",\"+rast\n",
    "    count+=1\n",
    "\n",
    "## Define computational region \n",
    "grass.run_command('g.region', overwrite=True, raster=groupraster)\n",
    "\n",
    "## Print and saving current time for processing time management\n",
    "print (\"Begin to merge \"+str(count)+\" individual segmentation maps on \" + time.ctime())\n",
    "begintime_merge=time.time()\n",
    "\n",
    "## Defining the formula for r.mapcalc\n",
    "formula=\"unclumped_raster= nmax(\"+grass.list_strings(\"rast\", pattern=\"segmentation_tile_\", mapset=\"SEGMENTATION\", flag='r')[0]\n",
    "for rast in grass.list_strings(\"rast\", pattern=\"segmentation_tile_\", mapset=\"SEGMENTATION\", flag='r')[1:]:\n",
    "    formula+=\",\"+rast\n",
    "formula+=\")\"\n",
    "\n",
    "## Running r.mapcalc to merge all raster together\n",
    "grass.mapcalc(formula, overwrite=True)\n",
    "\n",
    "## Compute processing time and print it\n",
    "print(str(count)+\" individual segment maps have been merge with 'r.mapcalc'\")\n",
    "print_processing_time(begintime_merge, \" Merging process achieved in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clump patched raster  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use here the ['r.clump' command](https://grass.osgeo.org/grass71/manuals/r.clump.html) to allow a new (unique) ID for each group of pixels with different values from their neighbors (because segments resulting from 'i.segment' on different processing tiles could have the same ID after being patched in the precedent step)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Print and saving current time for processing time management\n",
    "print (\"Begin clump of raster on \" + time.ctime())\n",
    "begintime_clump=time.time()\n",
    "\n",
    "## Generate new individual values for group of pixels\n",
    "grass.run_command('r.clump', overwrite=True, input=\"unclumped_raster@SEGMENTATION\", output=\"segmentation_raster@SEGMENTATION\")\n",
    "\n",
    "## Compute processing time and print it\n",
    "print_processing_time(begintime_clump, \"Segmentation raster have been clumped in \")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Compute basic statistics about the clumped raster. The maximum value correspond to the number of objets (patchs)\n",
    "nbrobject=grass.raster_info(\"segmentation_raster\")\n",
    "print \"The segmentation raster contain \"+str(int(nbrobject.max))+\" objects\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Erase intermediate maps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell will remove all the temporary files needed for the different previous steps. Be careful to be sure that your 'segmentation_raster' has correctly been processed before running this part, otherwise you should start all the segmentation part of this script again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Print \n",
    "print (\"Begin deleting temporary maps on \" + time.ctime())\n",
    "\n",
    "## Delete individual segmentation rasters \n",
    "grass.run_command('g.remove', flags=\"rf\", type=\"raster\", pattern=\"segmentation_tile_\")\n",
    "\n",
    "## Delete unclumped segmentation rasters \n",
    "grass.run_command('g.remove', flags=\"f\", type=\"raster\", name=\"unclumped_raster\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### End of part 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(\"The script ends at \"+ time.ctime())\n",
    "print_processing_time(begintime_segmentation_full, \"Entire process has been achieved in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-* </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-* </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4 - Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Launch GRASS GIS working session**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Set the name of the mapset in which to work\n",
    "mapsetname=user[\"classification_mapsetname\"]\n",
    "\n",
    "## Launch GRASS GIS working session in the mapset\n",
    "if os.path.exists(os.path.join(user[\"gisdb\"],user[\"location\"],mapsetname)):\n",
    "    gsetup.init(os.environ['GISBASE'], user[\"gisdb\"], user[\"location\"], mapsetname)\n",
    "    print \"You are now working in mapset '\"+mapsetname+\"'\" \n",
    "else: \n",
    "    print \"'\"+mapsetname+\"' mapset doesn't exists in \"+user[\"gisdb\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "begintime_classif_full=time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Copy data from other mapset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some data need to be copied from other mapsets into the current mapset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Copy segmentation raster in the current mapset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Copy segmentation raster layer from SEGMENTATION mapset to current mapset\n",
    "grass.run_command('g.copy', overwrite=True, raster=\"segmentation_raster@SEGMENTATION,segments\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Copy nDSM raster in the current mapset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-  Specifically to our data, our nDSM have *null values* which have been defined during the importation of the raw data, those *null values* (which correspond, in our case, most of the time to missed pixels in the stereo-photogrammetry process) have to be set to zero values of elevation (with the ['r.null' command](https://grass.osgeo.org/grass72/manuals/r.null.html)). Those missed pixels are for almost water surfaces (0 elevation on nDSM) or hidden side of buildings.\n",
    "\n",
    "-  If you are working with a nDSM data which didn't have any null values, please skip the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Copy nDSM raster layer from PERMANENT mapset to current mapset\n",
    "grass.run_command('g.copy', overwrite=True, raster=\"nDSM@PERMANENT,nDSM\")\n",
    "\n",
    "## Define computational region to match the extention of GEOBIA Subset\n",
    "grass.run_command('g.region', overwrite=True, raster=\"nDSM@CLASSIFICATION\")\n",
    "\n",
    "## Replace null values of nDSM with zero values \n",
    "grass.run_command('r.null', map=\"nDSM@CLASSIFICATION\", null=\"0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Display list of raster available in the current mapset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## List of raster available in CLASSIFICATION mapset\n",
    "print grass.list_strings(\"raster\", mapset=\"CLASSIFICATION\", flag='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creation of training/validation and test sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, you can choose between different procedures for building/creation of your training and test sets : \n",
    "\n",
    "1. Method A: Pre-defined training and test sets.\n",
    "2. Method B: Stratified random split of training and test sets.\n",
    "3. Method C: Spatial split of training and test sets.\n",
    "\n",
    "**Please run only the cells corresponding to the procedure you choose!**\n",
    "\n",
    "The creation of training/test sets is the most human-labour intensive work of the processing. Be careful that the quality of the training set and test set are important to achieve satisfying results. \n",
    "\n",
    "Create a shapefile of points with an attribute column called **'Class_num'** (INT type) and containing the class of the objet. Use numbers as class categories (1,10,2,20,25...) instead of text. \n",
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**\n",
    "Please notice that we use here [the following definitions](https://en.wikipedia.org/wiki/Test_set) when speaking about \"training set\", \"validation set\" and \"test set\":\n",
    "\n",
    "- Training set: A set of examples used for learning, that is to fit the parameters [i.e., weights] of the classifier.\n",
    "- Validation set: A set of examples used to tune the hyperparameters [i.e., architecture, not weights] of a classifier, for example to choose the number of hidden units in a neural network.\n",
    "- Test set: A set of examples used only to assess the performance [generalization] of a fully-specified classifier.\n",
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**\n",
    "Here, we build only a training set and an independent test set which will be used later for model's performance evaluation. The validation set to use for tunning machine learning model's parameters will be automatically generated from the training set provided to the GRASS GIS' classification add-on ['v.class.mlR'](https://grass.osgeo.org/grass70/manuals/addons/v.class.mlR.html) which implement cross-validation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Method A: Pre-defined training and test sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This procedure is the simplest. Training and test sets are should be created in distinct shapefile and are imported as separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Import vector shapefile with the training set\n",
    "grass.run_command('v.in.ogr', overwrite=True, input='F:\\\\.....\\\\Training_test\\\\training.shp', output='training_set')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Import vector shapefile with the test set\n",
    "grass.run_command('v.in.ogr', overwrite=True, input='F:\\\\.....\\\\Training_test\\\\test_set.shp', output='test_set')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Method B: Stratified random split of training and test sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this procedure the training set is the result of a (nonspatial) random selection from all available sample points, with stratification based on LULC classes. The user has to choose the ratio of available points to be used for the training set (e.g. 0.5 to split in two equal parts ; 0.75 to select 3/4 of available points for training). Then, the test set is defined as the opposite of the training set (the points still available after selection of training points). For this part, the following commands are used: ['g.remove'](https://grass.osgeo.org/grass72/manuals/g.remove.html), ['v.db.select'](https://grass.osgeo.org/grass72/manuals/v.db.select.html), ['v.extract'](https://grass.osgeo.org/grass72/manuals/v.extract.html), ['v.patch'](https://grass.osgeo.org/grass72/manuals/v.patch.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Import sample of points to be divided into training and test set**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Set computational region to match the default region\n",
    "grass.run_command('g.region', flags=\"d\")\n",
    "\n",
    "## Import sample data (points)\n",
    "grass.run_command('v.in.ogr', overwrite=True, input='F:\\\\.....\\\\Training_test\\\\sample_point.shp', output='samples')\n",
    "\n",
    "## Print\n",
    "print \"Point sample imported on \"+time.ctime()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Creation of training set as a ramdom stratified selection from available samples**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, you can modify the 'ratio' variable to change the percentage of available points which are going to be randomly selected in the training set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "print (\"Start building training set on \" + time.ctime())\n",
    "begintime_trainingset=time.time()\n",
    "\n",
    "## Set the ratio of available sample to select for training (between 0 and 1) \n",
    "ratio=0.5\n",
    "\n",
    "## Erase potential existing vector \n",
    "grass.run_command('g.remove', flags=\"rf\", type=\"vector\", pattern=\"temp_sample_\")\n",
    "grass.run_command('g.remove', flags=\"rf\", type=\"vector\", pattern=\"training_\")\n",
    "grass.run_command('g.remove', flags=\"rf\", type=\"vector\", pattern=\"training_set\")\n",
    "\n",
    "## Loop through all class label\n",
    "for classnum in grass.parse_command('v.db.select', map='samples', columns='Class_num', flags='c'):\n",
    "    ## Extract one vector layer per class\n",
    "    condition=\"Class_num='\"+str(classnum)+\"'\"\n",
    "    tempvectorname=\"temp_sample_\"+str(classnum)   \n",
    "    grass.run_command('v.extract', overwrite=True, input=\"samples\", type=\"point\", where=condition, output=tempvectorname)\n",
    "    \n",
    "    ## Extract class-based training sample (one for each class layer)\n",
    "    nbravailable=grass.vector_info(tempvectorname).points\n",
    "    nbrextract=int(nbravailable*ratio)\n",
    "    outputname=\"training_\"+classnum\n",
    "    grass.run_command('v.extract', overwrite=True, input=tempvectorname, output=outputname, type=\"point\", random=nbrextract)\n",
    "    print str(nbrextract)+\" training samples extracted from the \"+str(nbravailable)+\" available for class '\"+str(classnum)+\"'\"\n",
    "\n",
    "    grass.run_command('g.remove', flags=\"f\", type=\"vector\", name=tempvectorname)\n",
    "    \n",
    "    \n",
    "## Setting the list of vector to be patched \n",
    "inputlayers=grass.list_strings(\"vector\", pattern=\"training_\", mapset=\"CLASSIFICATION\", flag='r')[0]\n",
    "count=1\n",
    "for vect in grass.list_strings(\"vector\", pattern=\"training_\", mapset=\"CLASSIFICATION\", flag='r')[1:]:\n",
    "    inputlayers+=\",\"+vect\n",
    "    count+=1\n",
    "    \n",
    "## Patch of class-based trainings samples in one unique training set\n",
    "grass.run_command('g.remove', flags=\"f\", type=\"vector\", name=\"training_set\")\n",
    "grass.run_command('v.patch', flags=\"ne\", overwrite=True, input=inputlayers, output=\"training_set\")\n",
    "print str(count)+\" vector layers patched in one unique training set\"\n",
    "\n",
    "## Erase individual class-based training sample\n",
    "for vect in grass.list_strings(\"vector\", pattern=\"training_\", mapset=\"CLASSIFICATION\", flag='r'):\n",
    "    grass.run_command('g.remove', flags=\"f\", type=\"vector\", name=vect)\n",
    "    \n",
    "## Save number of records in the training set\n",
    "nbtraining=len(grass.parse_command('v.db.select', map='training_set', columns='Id', flags='c'))\n",
    "\n",
    "## Print number of records in the training set and processing time\n",
    "print(str(nbtraining)+\" points in the training set\")\n",
    "print_processing_time(begintime_trainingset, \"Training set build in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Creation of test set as the opposite of training set**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "print (\"Start building test set on \" + time.ctime())\n",
    "begintime_testset=time.time()\n",
    "\n",
    "## Erase existing vector \n",
    "grass.run_command('g.remove', flags=\"rf\", type=\"vector\", pattern=\"test_set\")\n",
    "\n",
    "## Save the id of training point \n",
    "list_id=[]\n",
    "for point_id in grass.parse_command('v.db.select', map='training_set', columns='Id', flags='c'):\n",
    "    list_id.append(str(point_id))\n",
    "\n",
    "## Build SQL statement for 'v.extract' command \n",
    "condition=\"Id not in (\"+str(list_id[0])\n",
    "for point_id in list_id[1:]:\n",
    "    condition+=\",\"+str(point_id)\n",
    "condition+=\")\"\n",
    "\n",
    "## From sample point, extract point not yet selected in training set\n",
    "grass.run_command('g.remove', flags=\"f\", type=\"vector\", name=\"test_set\")\n",
    "grass.run_command('v.extract', overwrite=True, input=\"samples\", type=\"point\", where=condition, output=\"test_set\")\n",
    "\n",
    "## Save number of records in the test set\n",
    "nbvalidation=len(grass.parse_command('v.db.select', map='test_set', columns='Id', flags='c'))\n",
    "\n",
    "## Print number of records in the test set and processing time\n",
    "print(str(nbvalidation)+\" points in the test set\")\n",
    "print_processing_time(begintime_testset, \"Test set build in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Method C: Spatial split of training and test sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this procedure, the training and test sets are split in order not to have training points inside a specified image subset (called here after 'GEOBIA_subset'). This subset is determined by a polygon vector layer to be imported (shapefile). All training points will be outside, while available points inside the image subset will be used for the test set. This approach could avoid spatial autocorrelation between training and test points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Import polygon to be used as image subset**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Import vector shapefile of the extend of image subset\n",
    "grass.run_command('v.in.ogr', overwrite=True, input='F:\\\\.....\\\\Training_test\\\\image_subset_polygon.shp', output='GEOBIA_subset')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Import sample of points to be divided into training/validation and test set**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Set computational region to match the default region\n",
    "grass.run_command('g.region', flags=\"d\")\n",
    "\n",
    "## Import sample data (points)\n",
    "grass.run_command('v.in.ogr', overwrite=True, input='F:\\\\.....\\\\Training_test\\\\sample_point.shp', output='samples')\n",
    "\n",
    "## Print\n",
    "print \"Point sample imported on \"+time.ctime()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Creation of training set as all available samples outside the image subset**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "print (\"Start building training set on \" + time.ctime())\n",
    "begintime_trainingset=time.time()\n",
    "\n",
    "## Erase existing vector \n",
    "grass.run_command('g.remove', flags=\"rf\", type=\"vector\", pattern=\"training_set\")\n",
    "\n",
    "## Select points inside the image subset\n",
    "grass.run_command('v.select', overwrite=True, ainput=\"samples\", atype=\"point\", binput=\"GEOBIA_subset\", btype=\"area\", \n",
    "                  output=\"training_set\", operator=\"overlap\", flags=\"r\")\n",
    "\n",
    "## Save number of records in the training set\n",
    "nbtraining=len(grass.parse_command('v.db.select', map='training_set', columns='Id', flags='c'))\n",
    "\n",
    "## Print number of records in the training set and processing time\n",
    "print(str(nbtraining)+\" points in the training set\")\n",
    "print_processing_time(begintime_trainingset, \"Training set build in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Creation of test set as the opposite of training set**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "print (\"Start building test set on \" + time.ctime())\n",
    "begintime_testset=time.time()\n",
    "\n",
    "## Erase existing vector \n",
    "grass.run_command('g.remove', flags=\"rf\", type=\"vector\", pattern=\"test_set\")\n",
    "\n",
    "## Save the id of training point \n",
    "list_id=[]\n",
    "for point_id in grass.parse_command('v.db.select', map='training_set', columns='Id', flags='c'):\n",
    "    list_id.append(str(point_id))\n",
    "\n",
    "## Build SQL statement for v.extract \n",
    "condition=\"Id not in (\"+str(list_id[0])\n",
    "for point_id in list_id[1:]:\n",
    "    condition+=\",\"+str(point_id)\n",
    "condition+=\")\"\n",
    "\n",
    "## Extract point not in training from all sample points\n",
    "grass.run_command('v.extract', overwrite=True, input=\"samples\", type=\"point\", where=condition, output=\"test_set\")\n",
    "\n",
    "## Save number of records in the test set\n",
    "nbvalidation=len(grass.parse_command('v.db.select', map='test_set', columns='Id', flags='c'))\n",
    "\n",
    "## Print number of records in the test set and processing time\n",
    "print(str(nbvalidation)+\" points in the test set\")\n",
    "print_processing_time(begintime_testset, \"Test set build in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Display by-class sample points distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following part is used to count up the number of points in training and test set. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Create temporary .csv file with attribute table (all columns) of \"training_set\" vector layer\n",
    "grass.run_command('v.db.select', overwrite=True, map=\"training_set@CLASSIFICATION\",\n",
    "                  file=os.path.join(tempfile.gettempdir(),\"tempfile.csv\"),separator=\"comma\")\n",
    "\n",
    "## Import .csv file into Jupyter notebook (with panda)\n",
    "dataframe=pd.read_csv(os.path.join(tempfile.gettempdir(),\"tempfile.csv\"), sep=',',header=0)\n",
    "print str(len(dataframe))+\" points in training_set layer\\n\"\n",
    "\n",
    "## Delete temporary .csv file\n",
    "os.remove(os.path.join(tempfile.gettempdir(),\"tempfile.csv\"))\n",
    "\n",
    "## Display the number of points per class in sample\n",
    "print \"Number of points per class in training_set\"\n",
    "print dataframe.groupby(\"Class_num\").size()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Create temporary .csv file with attribute table (all columns) of \"test_set\" vector layer\n",
    "grass.run_command('v.db.select', overwrite=True, map=\"test_set@CLASSIFICATION\",\n",
    "                  file=os.path.join(tempfile.gettempdir(),\"tempfile.csv\"),separator=\"comma\")\n",
    "\n",
    "## Import .csv file into Jupyter notebook (with panda)\n",
    "dataframe=pd.read_csv(os.path.join(tempfile.gettempdir(),\"tempfile.csv\"), sep=',',header=0)\n",
    "print str(len(dataframe))+\" points in test_set layer\\n\"\n",
    "\n",
    "## Delete temporary .csv file\n",
    "os.remove(os.path.join(tempfile.gettempdir(),\"tempfile.csv\"))\n",
    "\n",
    "## Display the number of points per class in sample\n",
    "print \"Number of points per class in test_set\"\n",
    "print dataframe.groupby(\"Class_num\").size()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Export of training and test sets in shapefile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This part is optional. If you built the training and test sets from a simple set of points (Method B or C), you can run the following cells to exports those sets as shapefiles in desired folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Export training set \n",
    "grass.run_command('v.out.ogr', flags=\"sc\", overwrite=True, input=\"training_set\", \n",
    "                  output=\"F:\\\\.....\\\\sample_shapefiles\\\\new_training_set.shp\",\n",
    "                  format=\"ESRI_Shapefile\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Export test set \n",
    "grass.run_command('v.out.ogr', flags=\"sc\", overwrite=True, input=\"test_set\", \n",
    "                  output=\"F:\\\\.....\\\\sample_shapefiles\\\\new_test_set.shp\",\n",
    "                  format=\"ESRI_Shapefile\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Compute statistics of training objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Identify which segment correspond to each training point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our processing chain, the training and test sets are formed of points. However, objects are needed to train the supervised classification in OBIA context. In this section, each point in the training set is used to identify the underlying object in the segmentation layer, and save its unique ID. We use the ['v.db.addcolumn' command](https://grass.osgeo.org/grass72/manuals/v.db.addcolumn.html), ['v.what.rast' command](https://grass.osgeo.org/grass72/manuals/v.what.rast.html) for this purpose."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "begintime_whatrast=time.time()\n",
    "\n",
    "## Add a column \"seg_id\" in training_set layer\n",
    "grass.run_command('v.db.addcolumn', map=\"training_set\", columns=\"seg_id int\")\n",
    "\n",
    "## Set computational region to the default region\n",
    "grass.run_command('g.region', flags=\"d\")\n",
    "\n",
    "## For each training point, add the value of the underlying segmentation raster pixel in column \"seg_id\"\n",
    "grass.run_command('v.what.rast', map=\"training_set\", raster=\"segments\", column=\"seg_id\")\n",
    "\n",
    "## Compute processing time and print it\n",
    "print_processing_time(begintime_whatrast, \"Segment iD added in attribute table of the 'training_set' vector layer in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create dataframe with \"seg_id\" and \"class\" of segments in training set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we create a [Pandas' dataframe](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.html) containing only two columns, the segment ID (named 'cat') and class. This dataframe will be used further for joint with the computed statistics of each segment. Please notice that the number of (distinct) segments to be used for training could be different of the number of points in initial training sample, as some points could refer to the same segment depending of the segmentation results.\n",
    "\n",
    "**In the 'columns' parameter, please set only the segment iD and the class to be used in the classification process (only two columns).**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Create a temporary .csv file containing segment iD and class of each training point\n",
    "grass.run_command('v.db.select', overwrite=True, map=\"training_set@CLASSIFICATION\", columns=\"seg_id,Class_num\",\n",
    "                  file=os.path.join(tempfile.gettempdir(),\"temp_train_segid_class.csv\"),separator=\"comma\")\n",
    "\n",
    "## Import .csv file in a temporary Pandas' dataframe\n",
    "temp=pd.read_csv(os.path.join(tempfile.gettempdir(),\"temp_train_segid_class.csv\"), sep=',',header=0)\n",
    "\n",
    "## Erase the temporary .csv file\n",
    "os.remove(os.path.join(tempfile.gettempdir(),\"temp_train_segid_class.csv\"))\n",
    "\n",
    "## Rename columns \"seg_id\" in \"cat\" for joint further \n",
    "temp.rename(columns={'seg_id': 'cat'}, inplace=True)\n",
    "\n",
    "## Keep only distinct value of column \"cat\" \n",
    "seg_id_class=temp.drop_duplicates(subset='cat', keep=False)\n",
    "\n",
    "## Print\n",
    "print \"Dataframe created with \"+str(len(seg_id_class))+\" distinct segments' ID for training set from the \"+str(len(temp))+\" point provided in the initial training sample\"\n",
    "\n",
    "## Display table\n",
    "seg_id_class.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a new raster layer with segments to be used for training"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we build a new raster layer containing only segments to be used for training. It will be used after to compute statistics of training objects. \n",
    "\n",
    "This raster is created by reclassifying the original segmentation layer (with segments for the whole area). For that, segments not included in the training set will be replaced with *NULL* values. The ['r.reclass' command](https://grass.osgeo.org/grass72/manuals/r.reclass.html) is used for this purpose. Before reclassification, a 'reclass rule file' containing instructions for reclassification is created.\n",
    "\n",
    "In GRASS GIS, a reclassified raster is only a specific rule assigned to another existing raster. When dealing with very large dataset, display a reclassified raster could be very long. If you want to ensure a faster display of a reclassified raster, you can write a new raster based on the reclassified one. Please note that writing a new raster will use more disk space. The last part of the following cell is dedicated to this purpose. It is optional."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Compute reclassification rules and build a raster of training segments**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "print (\"Bulding a raster map with training segments on \" + time.ctime())\n",
    "begintime_reclassify=time.time()\n",
    "\n",
    "## Define reclass rule\n",
    "rule=\"\"\n",
    "for seg_id in grass.parse_command('v.db.select', map='training_set', columns='seg_id', flags='c'):  #note that parse_command provide a list of DISTINCT values\n",
    "    rule+=str(seg_id)\n",
    "    rule+=\"=\"\n",
    "    rule+=str(seg_id)\n",
    "    rule+=\"\\n\"\n",
    "rule+=\"*\"\n",
    "rule+=\"=\"\n",
    "rule+=\"NULL\"\n",
    "\n",
    "## Create a temporary 'reclass_rule.csv' file\n",
    "outputcsv=os.path.join(tempfile.gettempdir(),\"reclass_rules.csv\") # Define the csv output file name\n",
    "f = open(outputcsv, 'w')\n",
    "f.write(rule)\n",
    "f.close()\n",
    "\n",
    "## Set computational region to the default region\n",
    "grass.run_command('g.region', flags=\"d\")\n",
    "\n",
    "## Reclass segments raster layer to keep only training segments, using the reclas_rule.csv file\n",
    "grass.run_command('r.reclass', overwrite=True, input=\"segments\", output=\"segments_training\", rules=outputcsv)\n",
    "\n",
    "## Erase the temporary 'reclass_rule.csv' file\n",
    "os.remove(outputcsv)\n",
    "\n",
    "## Create the same raster with r.mapcalc (to ensure fast display) \n",
    "##### Comment the following lines if you want to save disk space instead of fast display\n",
    "formula=\"segments_training_temp=segments_training\"\n",
    "grass.mapcalc(formula, overwrite=True)\n",
    "## Rename the new raster with the name of the original one (will be overwrited)\n",
    "grass.run_command('g.rename', overwrite=True, raster=\"segments_training_temp,segments_training\")\n",
    "# Remove the existing GRASS colortable (for faster display in GRASS map display)\n",
    "grass.run_command('r.colors', flags=\"r\", map=\"segments_training\", color=\"random\")\n",
    "\n",
    "## Compute processing time and print it\n",
    "print_processing_time(begintime_reclassify, \"Raster map with training segments builted in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute statistics on training segments with i.segment.stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use here the ['i.segment.stats' add-on](https://grass.osgeo.org/grass70/manuals/addons/i.segment.stats.html) to compute statistics for each object. As this add-on is not by-default installed, the first cell is there to install it with ['g.extension' command](https://grass.osgeo.org/grass72/manuals/g.extension.html). Another add-on, ['r.object.geometry'](https://grass.osgeo.org/grass70/manuals/addons/r.object.geometry.html) is also installed and is required for computing morphological statistics by i.segment.stats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Instal i.segment.stats if not yet installed\n",
    "if \"i.segment.stats\" not in grass.parse_command('g.extension', flags=\"a\"):\n",
    "    grass.run_command('g.extension', extension=\"i.segment.stats\")\n",
    "    print \"i.segment.stats have been installed on your computer\"\n",
    "else: print \"i.segment.stats is already installed on your computer\" \n",
    "    \n",
    "## Instal r.object.geometry if not yet installed\n",
    "if \"r.object.geometry\" not in grass.parse_command('g.extension', flags=\"a\"):\n",
    "    grass.run_command('g.extension', extension=\"r.object.geometry\")\n",
    "    print \"r.object.geometry have been installed on your computer\"\n",
    "else: print \"r.object.geometry is already installed on your computer\" "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Set list of raster from which to compute statistics with i.segment.stats**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here after, a list of raster layer on which to compute statistics is saved. Please adapt those layers according to the raster you want to use for object statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Display the name of rasters available in PERMANENT and CLASSIFICATION mapset\n",
    "print grass.list_strings(\"raster\", mapset=\"PERMANENT\", flag='r')\n",
    "print grass.list_strings(\"raster\", mapset=\"CLASSIFICATION\", flag='r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Define the list of raster layers for which statistics will be computed\n",
    "inputstats=\"opt_blue@PERMANENT\"\n",
    "inputstats+=\",opt_green@PERMANENT\"\n",
    "inputstats+=\",opt_nir@PERMANENT\"\n",
    "inputstats+=\",opt_red@PERMANENT\"\n",
    "inputstats+=\",NDVI@PERMANENT\"\n",
    "inputstats+=\",Brightness@PERMANENT\"\n",
    "inputstats+=\",nDSM@CLASSIFICATION\"\n",
    "print inputstats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Compute statistics of segments with i.segment.stats**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following section, ['i.segment.stats' add-on](https://grass.osgeo.org/grass70/manuals/addons/i.segment.stats.html) is used to compute object statistics. Please refer to the official help if you want to modify the parameters. Other raster statistics and morphological features could be used according to your needs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Define computational region to match the extention of segmentation raster\n",
    "grass.run_command('g.region', overwrite=True, raster=\"segments@CLASSIFICATION\")\n",
    "\n",
    "## Saving current time for processing time management\n",
    "print (\"Start computing statistics for training segments, using i.segment.stats on \" + time.ctime())\n",
    "begintime_isegmentstats=time.time()\n",
    "\n",
    "## Compute statistics of objets using i.segment.stats only with .csv output (no vectormap output)\n",
    "grass.run_command('i.segment.stats', overwrite=True, map=\"segments_training@CLASSIFICATION\", \n",
    "                  rasters=inputstats, \n",
    "                  raster_statistics=\"min,max,range,mean,stddev,sum,coeff_var,first_quart,median,third_quart,perc_90\", \n",
    "                  area_measures=\"area,perimeter,compact_circle\",\n",
    "                  csvfile=\"F:\\\\.....\\\\Classification\\\\i.segment.stats\\\\stats_training_sample.csv\", \n",
    "                  processes=str(user[\"nb_proc\"]))\n",
    "\n",
    "## Compute processing time and print it\n",
    "print_processing_time(begintime_isegmentstats, \"Segment statistics computed in :\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Remove temporary raster layer**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Remove \"segment_training\" raster layer\n",
    "grass.run_command('g.remove', flags=\"f\", type=\"raster\", name=\"segments_training@CLASSIFICATION\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check for unwanted values (Null/Inf values) in data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The purpose of the following section is to check presence of unwanted values in the statistics previously computed, like *null values* or *infinite values*. CSV file with object statistics just created with i.segment.stats is imported into a Pandas' dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Import .csv file\n",
    "temp_stat_train=pd.read_csv(\"F:\\\\.....\\\\Classification\\\\i.segment.stats\\\\stats_training_sample.csv\", sep=',',header=0)\n",
    "print \"The .csv file with results of i.segment.stats for \"+str(len(temp_stat_train))+\" training segments imported in a new dataframe\"\n",
    "    \n",
    "## Check and count for NaN values by column in the table\n",
    "if temp_stat_train.isnull().any().any():\n",
    "    for colomn in list(temp_stat_train.columns.values):\n",
    "        if temp_stat_train[colomn].isnull().any():\n",
    "            print \"Column '\"+str(colomn)+\"' have \"+str(temp_stat_train[colomn].isnull().sum())+\" NULL values\"\n",
    "else: print \"No missing values in dataframe\" \n",
    "        \n",
    "## Check and count for Inf values by column in the table\n",
    "if np.isinf(temp_stat_train).any().any():\n",
    "    for colomn in list(temp_stat_train.columns.values):\n",
    "        if np.isinf(temp_stat_train[colomn]).any():\n",
    "            print \"Column '\"+str(colomn)+\"' have \"+str(np.isinf(temp_stat_train[colomn]).sum())+\" Infinite values\"\n",
    "else: print \"No infinite values in dataframe\" \n",
    "    \n",
    "## Display table\n",
    "temp_stat_train.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<left> <font size=4> <b> Replace Null values in data with zero values (PLEASE USE CAREFULLY) </b> </font> </left> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "FOR EXPERIENCED USERS ONLY! \n",
    "\n",
    "The following section is dedicated to replace Null values in data with zero values (or other values, according to your needs). Use it only if you are sure that the missing values in your data could be replaced by another value!\n",
    "\n",
    "If you want to use the following cell, change its type in \"code\" instead of \"Markdown\" in the \"Jupyter notebook\" interface. Also, you will have to remove the first and the last line."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "```python\n",
    "## Fill NaN values with Zero value \n",
    "if temp_stat_train.isnull().any().any():\n",
    "    nbnan=temp_stat_train.isnull().sum().sum()\n",
    "    temp_stat_train.fillna(0, inplace=True)\n",
    "    print str(nbnan)+\" NaN value have been filled with Zero values\"\n",
    "else: print \"No missing values in dataframe\" \n",
    "    \n",
    "## Check and count for NaN values by column in the table\n",
    "if temp_stat_train.isnull().any().any():\n",
    "    for colomn in list(temp_stat_train.columns.values):\n",
    "        if temp_stat_train[colomn].isnull().any():\n",
    "            print \"Column '\"+str(colomn)+\"' still have \"+str(temp_stat_train[colomn].isnull().sum())+\" NULL values\"\n",
    "else: print \"No more missing values in dataframe\" \n",
    "        \n",
    "## Display table\n",
    "temp_stat_train.head()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<left> <font size=4> <b> Inf values in data </b> </font> </left> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have infinite values in your data, please find and solve this problem. The dataset could'nt have any Null or infinite values for classification process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Check and count for Inf values by column in the table\n",
    "if np.isinf(temp_stat_train).any().any():\n",
    "    for colomn in list(temp_stat_train.columns.values):\n",
    "        if np.isinf(temp_stat_train[colomn]).any():\n",
    "            print \"Column '\"+str(colomn)+\"' still have \"+str(np.isinf(temp_stat_train[colomn]).sum())+\" Infinite values\"\n",
    "else: print \"No more infinite values in dataframe\" "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Building final training set table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here after, dataframe of training segments' classes with dataframe of training segments' statistics are merged together and saved into a .csv file. This one will be used further in the machine learning classification add-on 'v.class.mlR'. [The merge function of Pandas](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.merge.html) is used to perform the joint between dataframes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Join between tables (pandas dataframe) on column 'cat'\n",
    "training_sample=pd.merge(seg_id_class, temp_stat_train, on='cat')\n",
    "\n",
    "## Check if there are NaN values in the table and print basic information\n",
    "if training_sample.isnull().any().any():\n",
    "    print \"WARNING: Some values are missing in the dataset\"\n",
    "else: \n",
    "    # Write dataframe in a .csv file\n",
    "    training_sample.to_csv(path_or_buf=\"F:\\\\.....\\\\Classification\\\\i.segment.stat\\\\stats_training_set.csv\", \n",
    "                       sep=',', header=True,  quoting=None, decimal='.', index=False)\n",
    "    print \"A new csv table called 'stats_training_set', to be used for training, have been created with \"+str(len(training_sample))+\" rows.\"\n",
    "    \n",
    "## Display table\n",
    "training_sample.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### Compute number of points per class in training sample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell could be used to see the distribution of training segments by LULC classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Number of points per class in training sample\n",
    "print \"Number of segments per class in training sample\\n\"\n",
    "print training_sample.groupby(\"Class_num\").size()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exclude some specific classes from training sample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section is optional and has to be used only if some specific classes have not to be used in the classification process.\n",
    "\n",
    "If you want to use the following cells, change their type in \"code\" instead of \"Markdown\" in the \"Jupyter notebook\" interface. Also, you will have to remove the first and the last line."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "```python\n",
    "## Import \n",
    "samples_attributes=pd.read_csv(\"F:\\\\MAUPP\\\\.....\\\\Classification\\\\i.segment.stat\\\\stats_training_set.csv\", sep=',',header=0)\n",
    "\n",
    "## Exclude specific row corresponding to some classes.\n",
    "samples_attributes.drop(samples_attributes[samples_attributes.Class_num==12].index, inplace=True)\n",
    "    \n",
    "## Write .csv file\n",
    "samples_attributes.to_csv(path_or_buf=\"F:\\\\MAUPP\\\\.....\\\\Classification\\\\i.segment.stat\\\\stats_training_set.csv\", \n",
    "                          sep=',', header=True,  quoting=None, decimal='.', index=False)\n",
    "print \"A new csv table called 'sample_training', with samples to be used for training, have been created with \"+str(len(samples_attributes))+\" rows.\"\n",
    "\n",
    "## Display table\n",
    "samples_attributes.head()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "## Number of points per class in training sample\n",
    "print \"Number of segments per class in training sample\\n\"\n",
    "print samples_attributes.groupby(\"Class_num\").size()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Compute statistics for segments to be classified"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section uses the ['i.segment.stats' add-on](https://grass.osgeo.org/grass70/manuals/addons/i.segment.stats.html) to compute statistics for each object to be classified. In the context of the GEOBIA conference, just an image subset (GEOBIA_subset) has been classified. You can adapt this part of script according to your own needs.\n",
    "\n",
    "**Please be careful that the statistic you will compute for objects to be classified should be the same as those computed previously for the training set!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Create raster of segments to be classified which are in the image subset**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define computational region to match the extention of image Subset\n",
    "grass.run_command('g.region', overwrite=True, vector=\"GEOBIA_subset@CLASSIFICATION\", align=\"segments@CLASSIFICATION\")\n",
    "\n",
    "# Create a new raster layer with segments inside the current computational region, using r.map.calc\n",
    "formula=\"GEOBIA_segments=segments@CLASSIFICATION\"\n",
    "grass.mapcalc(formula, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Set list of raster from which to compute statistics with i.segment.stats**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Display the name of rasters available in PERMANENT and CLASSIFICATION mapset\n",
    "print grass.read_command('g.list',type=\"raster\", mapset=\"PERMANENT\", flags='rp')\n",
    "print grass.read_command('g.list',type=\"raster\", mapset=\"CLASSIFICATION\", flags='rp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Define the list of raster layers for which statistics will be computed\n",
    "inputstats=\"opt_blue@PERMANENT\"\n",
    "inputstats+=\",opt_green@PERMANENT\"\n",
    "inputstats+=\",opt_nir@PERMANENT\"\n",
    "inputstats+=\",opt_red@PERMANENT\"\n",
    "inputstats+=\",NDVI@PERMANENT\"\n",
    "inputstats+=\",Brightness@PERMANENT\"\n",
    "inputstats+=\",nDSM@CLASSIFICATION\"\n",
    "print inputstats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Compute statistics of segments to be classified (with i.segment.stats)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Define computational region to match the extention of segmentation raster\n",
    "grass.run_command('g.region', overwrite=True, raster=\"GEOBIA_segments@CLASSIFICATION\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "print (\"Start computing statistics for segments to be classified, using i.segment.stats on \" + time.ctime())\n",
    "begintime_isegmentstats=time.time()\n",
    "\n",
    "## Compute statistics of objets using i.segment.stats only with .csv output (no vectormap output).\n",
    "grass.run_command('i.segment.stats', overwrite=True, map=\"GEOBIA_segments@CLASSIFICATION\", \n",
    "                  rasters=inputstats, \n",
    "                  raster_statistics=\"min,max,range,mean,stddev,sum,coeff_var,first_quart,median,third_quart,perc_90\", \n",
    "                  area_measures=\"area,perimeter,compact_circle\",\n",
    "                  csvfile=\"F:\\\\.....\\\\Classification\\\\i.segment.stat\\\\stats_segments.csv\", \n",
    "                  processes=str(user[\"nb_proc\"]))\n",
    "\n",
    "## Compute processing time and print it\n",
    "print_processing_time(begintime_isegmentstats, \"Segment statistics computed in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check for unwanted values (Null/NaN/Inf values) in data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The purpose of the following section is to check presence of unwanted values in the statistics previously computed, like *null values* or *infinite values*. CSV file with object statistics just created with i.segment.stats is imported into a Pandas' dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Import .csv file\n",
    "stats_segments=pd.read_csv(\"F:\\\\.....\\\\Classification\\\\i.segment.stat\\\\stats_segments.csv\", sep=',',header=0)\n",
    "print \"The .csv file with results of i.segment.stats for the \"+str(len(stats_segments))+\" segments to be classified imported in a new dataframe\"\n",
    "\n",
    "## Check and count for NaN values by column in the table\n",
    "if stats_segments.isnull().any().any():\n",
    "    for colomn in list(stats_segments.columns.values):\n",
    "        if stats_segments[colomn].isnull().any():\n",
    "            print \"Column '\"+str(colomn)+\"' have \"+str(stats_segments[colomn].isnull().sum())+\" NULL values\"\n",
    "else: print \"No missing values in dataframe\" \n",
    "        \n",
    "## Check and count for Inf values by column in the table\n",
    "if np.isinf(stats_segments).any().any():\n",
    "    for colomn in list(stats_segments.columns.values):\n",
    "        if np.isinf(stats_segments[colomn]).any():\n",
    "            print \"Column '\"+str(colomn)+\"' have \"+str(np.isinf(stats_segments[colomn]).sum())+\" Infinite values\"\n",
    "else: print \"No infinite values in dataframe\" \n",
    "\n",
    "## Display table\n",
    "stats_segments.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<left> <font size=4> <b> Replace Null/NaN values in data with zero values (PLEASE USE CAREFULLY) </b> </font> </left> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "FOR EXPERIENCED USERS ONLY! \n",
    "\n",
    "The following section is dedicated to replace Null values in data with zero values (or other values, according to your needs). Use it only if you are sure that the missing values in your data could be replaced by another value!\n",
    "\n",
    "If you want to use the following cell, change its type in \"code\" instead of \"Markdown\" in the \"Jupyter notebook\" interface. Also, you will have to remove the first and the last line."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "```python\n",
    "## Fill NaN values with Zero value \n",
    "if stats_segments.isnull().any().any():\n",
    "    nbnan=stats_segments.isnull().sum().sum()\n",
    "    stats_segments.fillna(0, inplace=True)\n",
    "    print str(nbnan)+\" NaN value have been filled with Zero values\"\n",
    "else: print \"No missing values in dataframe\" \n",
    "    \n",
    "## Check and count for NaN values by column in the table\n",
    "if stats_segments.isnull().any().any():\n",
    "    for colomn in list(stats_segments.columns.values):\n",
    "        if stats_segments[colomn].isnull().any():\n",
    "            print \"Column '\"+str(colomn)+\"' still have \"+str(stats_segments[colomn].isnull().sum())+\" NULL values\"\n",
    "else: print \"No more missing values in dataframe\" \n",
    "        \n",
    "## Display table\n",
    "stats_segments.head()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<left> <font size=4> <b> Inf values in data </b> </font> </left> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have infinite values in your data, please find and solve this problem. The dataset could'nt have any Null or infinite values for classification process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Check and count for Inf values by column in the table\n",
    "if np.isinf(stats_segments).any().any():\n",
    "    for colomn in list(stats_segments.columns.values):\n",
    "        if np.isinf(stats_segments[colomn]).any():\n",
    "            print \"Column '\"+str(colomn)+\"' still have \"+str(np.isinf(stats_segments[colomn]).sum())+\" Infinite values\"\n",
    "else: print \"No infinite values in dataframe\" "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Classification using multiple machine learning classifiers system "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Classification with v.class.mlR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ['v.class.mlR' add-on](https://grass.osgeo.org/grass70/manuals/addons/v.class.mlR.html) is used here in order to classify the segments using the training data. You can choose between several machine learning classifiers and several majority-voting systems. Please read the [add-on's help](https://grass.osgeo.org/grass70/manuals/addons/v.class.mlR.html) for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Instal v.class.mlR if not yet installed\n",
    "if \"v.class.mlR\" not in grass.parse_command('g.extension', flags=\"a\"):\n",
    "    grass.run_command('g.extension', extension=\"v.class.mlR\")\n",
    "    print \"v.class.mlR have been installed on your computer\"\n",
    "else: print \"v.class.mlR is already installed on your computer\" "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Classification**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please notice that if the classification process failed, it could be due to : \n",
    "- R software is not installed on your computer.\n",
    "- the \"R_LIBS_USER\" environment variable defined in this notebook and in R is not the same. Please return on the begin of this notebook and read the instructions.\n",
    "- You didn't respect the syntax of the folder path in the following cell (/,//) resulting in failure when running the R script.\n",
    "- The .csv files containing the statistics of objects to be classified and of training objects present some unaccepted values like *Null* or *infinite*. It couldn't have any 'hole' in the dataset.\n",
    "\n",
    "Please read the [official help](https://grass.osgeo.org/grass70/manuals/addons/v.class.mlR.html) to know which parameter adapt or not according to your needs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "print (\"Start classification process, using v.class.mlR on \" + time.ctime())\n",
    "begintime_vclassmlr=time.time()\n",
    "\n",
    "## Classification using v.class.mlR\n",
    "grass.run_command('v.class.mlR', flags=\"fi\", overwrite=True, \n",
    "separator=\"comma\",\n",
    "segments_file=\"F:/...../Classification/i.segment.stat/stats_segments.csv\", \n",
    "training_file=\"F:/...../Classification/i.segment.stat/training_sample.csv\", \n",
    "raster_segments_map=\"GEOBIA_segments@CLASSIFICATION\",\n",
    "classified_map=\"indiv_classification\", \n",
    "train_class_column=\"Class_num\",\n",
    "output_class_column=\"vote\", \n",
    "output_prob_column=\"prob\", \n",
    "classifiers=\"svmRadial,rf,rpart,knn\", \n",
    "folds=\"5\", \n",
    "partitions=\"10\", \n",
    "tunelength=\"10\", \n",
    "weighting_modes=\"smv,swv,bwwv,qbwwv\", \n",
    "weighting_metric=\"accuracy\", \n",
    "classification_results=\"F://.....//Classification//all_results.csv\", \n",
    "accuracy_file=\"F://.....//Classification//accuracy.csv\", \n",
    "model_details=\"F://.....//Classification//classifier_runs.txt\", \n",
    "bw_plot_file=\"F://.....//Classification//box_whisker\",\n",
    "r_script_file=\"F://.....//Classification//Rscript_mlR.R\", \n",
    "processes=\"2\")\n",
    "\n",
    "## Compute processing time and print it\n",
    "print_processing_time(begintime_vclassmlr, \"Classification process achieved in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import results of v.class.mlR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Import accuracy results of individual classifiers, resulting from cross-validation of tuning**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Import .csv file\n",
    "accuracy=pd.read_csv(\"F:\\\\.....\\\\Classification\\\\accuracy.csv\", sep=',',header=0)\n",
    "    \n",
    "## Display table\n",
    "accuracy.head(15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Import classifiers tuning parameters and individual classifier confusion matrix**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Open file\n",
    "classifier_runs = open('F:\\\\.....\\\\Classification\\\\classifier_runs.txt', 'r')\n",
    "    \n",
    "## Read file\n",
    "print classifier_runs.read()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Copy classified raster as 'real raster' in the current mapset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As classified maps from v.class.mlR are reclassed map of the original segmented map, display in GRASS GIS can be too slow. If you want, you can copy this classified maps as \"real raster\" in the current mapset. Please notice that it will use more disk space !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Display the list of raster available in the current mapset\n",
    "print grass.read_command('g.list', type=\"raster\", mapset=\"CLASSIFICATION\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When copying the classified raster, we change the color table in the same time, using [r.colors](https://grass.osgeo.org/grass72/manuals/r.colors.html). You can adapt the R:G:B values for the color table according to the colors you want each class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Make a copy of the classified maps of faster display in GRASS GIS\n",
    "\n",
    "## Saving current time for processing time management\n",
    "print (\"Making a copy of classified maps in current mapset on \" + time.ctime())\n",
    "begintime_copyraster=time.time()\n",
    "\n",
    "for classif in grass.list_strings(\"rast\", pattern=\"indiv_classification_\", flag='r'):\n",
    "    ## Create the same raster with r.mapcalc\n",
    "    formula=str(classif[:-15])+\"_temp=\"+str(classif[:-15])\n",
    "    grass.mapcalc(formula, overwrite=True)\n",
    "       \n",
    "    ## Rename the new raster with the name of the original one (will be overwrited)\n",
    "    renameformula=str(classif[:-15])+\"_temp,\"+str(classif[:-15])\n",
    "    grass.run_command('g.rename', overwrite=True, raster=renameformula)\n",
    "    \n",
    "    ## Define color table. Replace with the RGB values of wanted colors of each class\n",
    "    color_table=\"11  227:26:28\"+\"\\n\"\n",
    "    color_table+=\"12  255:141:1\"+\"\\n\"\n",
    "    color_table+=\"13  94:221:227\"+\"\\n\"\n",
    "    color_table+=\"14  102:102:102\"+\"\\n\"\n",
    "    color_table+=\"21  246:194:142\"+\"\\n\"\n",
    "    color_table+=\"22  211:217:173\"+\"\\n\"\n",
    "    color_table+=\"31  0:128:0\"+\"\\n\"\n",
    "    color_table+=\"32  189:255:185\"+\"\\n\"\n",
    "    color_table+=\"33  88:190:141\"+\"\\n\"\n",
    "    color_table+=\"34  29:220:0\"+\"\\n\"\n",
    "    color_table+=\"41  30:30:192\"+\"\\n\"\n",
    "    color_table+=\"51  0:0:0\"+\"\\n\"\n",
    "    \n",
    "    ## Create a temporary 'color_table.txt' file\n",
    "    outputcsv=\"F:\\\\.....\\\\Classification\\\\Results_maps\\\\temp_color_table.txt\" # Define the csv output file name\n",
    "    f = open(outputcsv, 'w')\n",
    "    f.write(color_table)\n",
    "    f.close()\n",
    "    \n",
    "    ## Apply new color the existing GRASS colortable (for faster display in GRASS map display)\n",
    "    grass.run_command('r.colors', map=classif, rules=outputcsv)\n",
    "    \n",
    "    ## Erase the temporary 'color_table.txt' file\n",
    "    os.remove(\"F:\\\\.....\\\\Classification\\\\Results_maps\\\\temp_color_table.txt\")\n",
    "    \n",
    "## Compute processing time and print it\n",
    "print_processing_time(begintime_copyraster, \"Classified raster maps have been copied in current mapset in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Export of classification raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "print (\"Export classified raster maps on \" + time.ctime())\n",
    "begintime_exportraster=time.time()\n",
    "\n",
    "for classif in grass.list_strings(\"rast\", pattern=\"indiv_classification_\", flag='r'):\n",
    "    outputname=\"F:\\\\.....\\\\Classification\\\\classified_raster\\\\\"+str(classif[21:-15])+\".tif\"\n",
    "    grass.run_command('r.out.gdal', overwrite=True, input=classif, output=outputname, format='GTiff')\n",
    "    \n",
    "## Compute processing time and print it\n",
    "print_processing_time(begintime_exportraster, \"Classified raster maps exported in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### End of part 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(\"The script ends at \"+ time.ctime())\n",
    "print_processing_time(begintime_classif_full, \"Entire process has been achieved in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-* </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> *-*-*-*-*-*-*-*-*-*-*-* </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5 - Performance evaluation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Launch GRASS GIS working session**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Set the name of the mapset in which to work\n",
    "mapsetname=user[\"classification_mapsetname\"]\n",
    "\n",
    "## Launch GRASS GIS working session in the mapset\n",
    "if os.path.exists(os.path.join(user[\"gisdb\"],user[\"location\"],mapsetname)):\n",
    "    gsetup.init(os.environ['GISBASE'], user[\"gisdb\"], user[\"location\"], mapsetname)\n",
    "    print \"You are now working in mapset '\"+mapsetname+\"'\" \n",
    "else: \n",
    "    print \"'\"+mapsetname+\"' mapset doesn't exists in \"+user[\"gisdb\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "begintime_perform=time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Import validation sample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following section is dedicated to the importation of test set points. Please adapt the path to you own data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Set computational region\n",
    "grass.run_command('g.region', overwrite=True, raster=\"segments\")\n",
    "\n",
    "## Import points sample\n",
    "grass.run_command('v.in.ogr', overwrite=True, \n",
    "                  input='F:\\\\.....\\\\Training_test\\\\test_set.shp', output='test_set')\n",
    "\n",
    "## Print\n",
    "print \"Validation sample imported on \"+time.ctime()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can run the next cell if you want to see the attribute table of the \"test_set\" vector layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Create temporary .csv file with columns of \"test_set\" vector layer\n",
    "grass.run_command('v.db.select', overwrite=True, map=\"test_set@CLASSIFICATION\",\n",
    "                  file=\"F:\\\\.....\\\\Training_validation\\\\test_set.csv\",separator=\"comma\")\n",
    "\n",
    "## Import .csv file into Jupyter notebook (with panda)\n",
    "validation_samples_attributes=pd.read_csv(\"F:\\\\.....\\\\Training_validation\\\\test_set.csv\", sep=',',header=0)\n",
    "print str(len(validation_samples_attributes))+\" points in sample layer imported\"\n",
    "\n",
    "## Delete temporary .csv file\n",
    "os.remove(\"F:\\\\.....\\\\Training_validation\\\\test_set.csv\")\n",
    "\n",
    "## Display table\n",
    "validation_samples_attributes.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## Number of points per class in validation sample\n",
    "print \"Number of points per class in validation sample\\n\"\n",
    "print validation_samples_attributes.groupby(\"Class_num\").size()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Build dataframe with predicted class for each validation point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Add prediction of each classifier for validations points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, predictions of each classifier is saved in the attribute table of 'test_set' vector layer. The ['v.db.addcolumn' command](https://grass.osgeo.org/grass72/manuals/v.db.addcolumn.html) and the ['v.what.rast' command](https://grass.osgeo.org/grass72/manuals/v.what.rast.html) are used for this purpose. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "begintime_whatrast=time.time()\n",
    "\n",
    "## Initialize a empty list\n",
    "allclassif=[]\n",
    "\n",
    "## Loop through all individual classification results\n",
    "for classif in grass.list_strings(\"rast\", pattern=\"indiv_classification_\", flag='r'):\n",
    "    nameclassif=str(classif[21:-15])  # Save the name of classifier\n",
    "    allclassif.append(nameclassif)  # Add the name of classifier in the list\n",
    "    \n",
    "    ## Add a \"int\" column in test_set layer, for each classification result\n",
    "    grass.run_command('v.db.addcolumn', map=\"test_set\", columns=nameclassif+\" int\")\n",
    "    \n",
    "    ## For each validation point, add the value of the underlying classifier raster pixel in column \"seg_id\"\n",
    "    grass.run_command('v.what.rast', map=\"test_set\", raster=\"indiv_classification_\"+nameclassif+\"@CLASSIFICATION\", column=nameclassif)\n",
    "    \n",
    "## Compute processing time and print it\n",
    "print(\"Predicted classes for '\"+', '.join(allclassif)+\"' added in the 'test_set' layer\")\n",
    "print_processing_time(begintime_whatrast, \"Prossess achieved in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next cell, the 'test_set' vector layer attribute table is exported in .csv file. Please replace the \"columnstoexport\" variable according to your own data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Export 'test_set' vector layer attribute table in .csv file.\n",
    "columnstoexport=\"Class_num\"+\",\"\n",
    "columnstoexport+=', '.join(allclassif)\n",
    "grass.run_command('v.db.select', overwrite=True, map=\"test_set@CLASSIFICATION\", columns=columnstoexport,\n",
    "                  file=\"F:\\\\.....\\\\Classification\\\\Validation\\\\predicted_gtruth.csv\",separator=\"comma\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exclude some specific classes from test set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section is optional and has to be used only if some specific classes not have to be used in the performance evaluation process.\n",
    "\n",
    "Please notice that you need to have the same classes in your training set and test set. \n",
    "\n",
    "If you want to use the following cell, change its type in \"code\" instead of \"Markdown\" in the \"Jupyter notebook\" interface. Also, you will have to remove the first and the last line."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "```python\n",
    "## Extract only test set points which have to be used for performance evaluation.PLEASE REPLACE THE \"WHERE\" CONDITION ACCORDING TO YOUR OWN DATA\n",
    "grass.run_command('v.extract', overwrite=True, input=\"test_set\", where=\"Class_num is not 12\", output=\"test_set_filter\")\n",
    "\n",
    "## Export 'test_set' vector layer attribute table.PLEASE REPLACE THE \"COLUMNS\" PARAMETER ACCORDING TO YOUR OWN DATA\n",
    "columnstoexport=\"Class_num\"+\",\"\n",
    "columnstoexport+=', '.join(allclassif)\n",
    "grass.run_command('v.db.select', overwrite=True, map=\"test_set_filter\", columns=columnstoexport,\n",
    "                  file=\"F:\\\\.....\\\\Classification\\\\Validation\\\\predicted_gtruth.csv\",separator=\"comma\")\n",
    "\n",
    "## Import data in dataframe\n",
    "validation_samples_attributes=pd.read_csv(\"F:\\\\.....\\\\Classification\\\\Validation\\\\predicted_gtruth.csv\", sep=',',header=0)\n",
    "\n",
    "## Number of points per class in training sample\n",
    "print \"Number of points per class in validation sample: \"+str(len(validation_samples_attributes))+\"\\n\"\n",
    "print validation_samples_attributes.groupby(\"Class_num\").size()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the end of this section, performance evaluation of classifications is performed. As mentioned in our article, our legend scheme is designed in two hierarchical levels. The second level is the most detailed (11 classes). The first level classes are derived from classification results on the first level classes. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Create inputs for classification perfomance evaluation (Level-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Define computational region to match the extention of segmentation raster\n",
    "grass.run_command('g.region', overwrite=True, raster=\"segments@CLASSIFICATION\")\n",
    "\n",
    "## Create raster layers with one pixel corresponding to each object. Pixels values representing either the ground thruth or the prediction of a specific classifier\n",
    "grass.run_command('v.to.rast', overwrite=True, input='test_set_filter', output='PE_L2_Class_num', use='attr', attribute_column='Class_num')    \n",
    "for result in allclassif:\n",
    "    outputname=\"PE_L2_\"+str(result)\n",
    "    grass.run_command('v.to.rast', overwrite=True, input='test_set_filter', output=outputname, use='attr', attribute_column=result)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Create inputs of classification perfomance evaluation (Level-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Loop through all_raster used for PE at level2\n",
    "for result in grass.list_strings('rast', pattern=\"PE_L2\", flag=\"r\"):\n",
    "    ## Reclass the pixels of inputs of level2 to match the level1 classes\n",
    "    rule=\"\"\n",
    "    for pixel_value in grass.parse_command('v.db.select', map='test_set_filter', columns=result[6:-15], flags='c'):  #note that parse_command provide a list of distinct values\n",
    "        rule+=str(pixel_value)\n",
    "        rule+=\"=\"\n",
    "        rule+=str(pixel_value[:-1])\n",
    "        rule+=\"\\n\"\n",
    "    rule+=\"*\"\n",
    "    rule+=\"=\"\n",
    "    rule+=\"NULL\"\n",
    "\n",
    "    ## Create a temporary 'reclass_rule.csv' file\n",
    "    outputcsv=\"F:\\\\.....\\\\Classification\\\\Validation\\\\reclass_rules.csv\" # Define the csv output file name\n",
    "    f = open(outputcsv, 'w')\n",
    "    f.write(rule)\n",
    "    f.close()\n",
    "\n",
    "    #### Reclass level2 classes to match level1 classes\n",
    "    outputname=\"PE_L1_\"+str(result[6:-15])\n",
    "    grass.run_command('r.reclass', overwrite=True, input=result, output=outputname, rules=outputcsv)\n",
    "    ## Erase the temporary 'reclass_rule.csv' file\n",
    "    os.remove(\"F:\\\\.....\\\\Classification\\\\Validation\\\\reclass_rules.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Classification performance evaluation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Performance evaluation is conducted here with the ['r.kappa' module](https://grass.osgeo.org/grass73/manuals/r.kappa.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Level 2**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "begintime_kappa_L2=time.time()\n",
    "\n",
    "## Classification perfomance evalutation using r.kappa (compute per-class kappa)\n",
    "for result in grass.list_strings('rast', pattern=\"PE_L2\", flag=\"r\", exclude=\"PE_L2_Class_num\"):\n",
    "    outputfile=\"F:\\\\.....\\\\Classification\\\\Validation\\\\rkappa_\"+str(result[3:-15])+\".txt\"\n",
    "    grass.run_command('r.kappa', flags=\"w\", overwrite=True, classification=result, reference=\"PE_L2_Class_num\", output=outputfile)\n",
    "    \n",
    "## Compute processing time and print it\n",
    "print_processing_time(begintime_kappa_L2, \"Performance evaluation for Level 2 achieved in :\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Level 1**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Saving current time for processing time management\n",
    "begintime_kappa_L1=time.time()\n",
    "\n",
    "## Classification perfomance evalutation using r.kappa (compute per-class kappa)\n",
    "for result in grass.list_strings('rast', pattern=\"PE_L1\", flag=\"r\", exclude=\"PE_L1_Class_num\"):\n",
    "    outputfile=\"F:\\\\.....\\\\Classification\\\\Validation\\\\rkappa_\"+str(result[3:-15])+\".txt\"\n",
    "    grass.run_command('r.kappa', flags=\"w\", overwrite=True, classification=result, reference=\"PE_L1_Class_num\", output=outputfile)\n",
    "\n",
    "## Compute processing time and print it\n",
    "print_processing_time(begintime_kappa_L1, \"Performance evaluation for Level 1 achieved in :\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Clean mapset "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Erase temporary files no more needed\n",
    "grass.run_command('g.remove', flags=\"rf\", type=\"raster\", pattern=\"PE_\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### End of part 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(\"The script ends at \"+ time.ctime())\n",
    "print_processing_time(begintime_perform, \"Entire process has been achieved in \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# End of this Jupyter notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(\"The script ends at \"+ time.ctime())\n",
    "print_processing_time(begintime_full,\"Entire process has been achieved in \")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}