sentinel2_grass_gis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [Part 2] Accessing Sentinel EO data\n",
    "\n",
    "## Introduction\n",
    "\n",
    "In this notebook we will demonstrate the use of GRASS GIS 8.2+ in combination with Python within a Jupyter Notebook environment. We will use some GRASS GIS modules and related Python libraries which facilitate scripting (`grass.script`) and which let us connect and interact with Jupyter Notebooks (`grass.jupyter`). \n",
    "\n",
    "In this notebook we'll process Sentinel-2 data: NDVI computation and some time series processing.\n",
    "\n",
    "### Table of contents\n",
    "\n",
    "1. Why Jupyter Notebooks and how to use them?\n",
    "2. GRASS GIS & Python\n",
    "3. Paths and variables, connecting to GRASS GIS backend\n",
    "4. Data upload to notebook session\n",
    "5. Initialization of GRASS GIS in the Jupyter notebook session\n",
    "6. Creating an area of interest map\n",
    "7. Importing geodata into GRASS GIS\n",
    "8. Sentinel-2 processing overview\n",
    "9. Computing NDVI\n",
    "10. Time series data processing\n",
    "11. Creating an image stack (imagery group)\n",
    "12. Object recognition with image segmentation\n",
    "13. Supervised Classification: RandomForest\n",
    "14. What's next?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Why Jupyter Notebooks and how to use them?\n",
    "\n",
    "Jupyter Notebooks are server-client applications that allow code written in a notebook document to be edited and executed through a web browser. They can be run on a local computer that does not require Internet access, as well as used to control computations on a remote server accessed via the Internet ([documentation](https://jupyter-notebook-beginner-guide.readthedocs.io/en/latest/what_is_jupyter.html)).\n",
    "\n",
    "Jupyter notebooks can be interactive and are run through a web browser. They provide the ability to combine live code, explanatory text, and computational results into a single document. Jupyter Notebooks can be easily shared as documents.\n",
    "\n",
    "Jupyter Notebooks are \n",
    "\n",
    "* convenient for initial code development (prototyping)\n",
    "* ideal for code segmentation with the ability to re-run segments\n",
    "* able to store values of variables from already executed segments.\n",
    "\n",
    "The notebook can be saved as an executable Python script in addition to the native `.ipynb` format, or exported to various documentation formats such as PDF or Sphinx RST with nice styling.\n",
    "\n",
    "#### Editing and interactive use\n",
    "\n",
    "Editing a Jupyter notebook is very easy: in the web browser, you can navigate between text or code blocks (\"cells\") using the mouse or keyboard shortcuts (see Menu > Help > Keyboard Shortcuts). You can execute small code segments cell by cell, save the notebook in its current state, or modify and recalculate cells or return them to their previous state. In addition to executable code cells, you can use Markdown in documentation cells to make them presentable to others.\n",
    "\n",
    "And: there are now many Jupyter notebooks on the Internet for inspiration!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. GRASS GIS & Python\n",
    "\n",
    "<!-- This cell has been written by Veronica Andreo, https://veroandreo.gitlab.io/ -->\n",
    "\n",
    "### Python library \"grass.script\"\n",
    "\n",
    "The **grass.script** or GRASS GIS Python Scripting Library provides functions for calling GRASS modules within Python scripts as sub-processes. The most commonly used functions include:\n",
    "\n",
    "- `run_command()`: used when the output of the modules is a raster or vector, no text type output is expected\n",
    "- `read_command()`: used when the output of the modules is a text\n",
    "- `parse_command()`: used with modules whose output can be converted to `key=value` pairs\n",
    "- `write_command()`: used with modules that expect text input, either in the form of a file or from stdin\n",
    "\n",
    "It also provides several wrapper functions for frequently used modules:\n",
    "\n",
    "- To get info from a raster, script.raster.raster_info() is used: `gs.raster_info('dsm')`\n",
    "- To get info of a vector, script.vector.vector_info() is used: `gs.vector_info('roads')`\n",
    "- To list the raster in a location, script.core.list_grouped() is used: `gs.list_grouped(type=['raster'])`\n",
    "- To obtain the computational region, script.core.region() is used: `gs.region()`\n",
    "\n",
    "More info: https://grass.osgeo.org/grass-stable/manuals/libpython/script_intro.html\n",
    "\n",
    "### Python library \"grass.jupyter\"\n",
    "\n",
    "The **grass.jupyter** library improves the integration of GRASS and Jupyter, and provides different classes to visualize maps.\n",
    "\n",
    "- `init()`: starts a GRASS session and sets up necessary environment variables\n",
    "- `Map()`: 2D rendering\n",
    "- `Map3D()`: 3D rendering\n",
    "- `InteractiveMap()`: interactive visualization with folium\n",
    "- `TimeSeriesMap()`: visualization for spatio-temporal data\n",
    "\n",
    "More info: https://grass.osgeo.org/grass-stable/manuals/libpython/grass.jupyter.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once set, we will test if the installation works."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Paths and variables, connecting to GRASS GIS backend\n",
    "\n",
    "For the ease of working in this notebook, we define some session variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define some needed path variables\n",
    "\n",
    "import os\n",
    "\n",
    "# data directory\n",
    "homedir=os.path.join(os.path.expanduser('~'), \"odse_workshop2022\")\n",
    "\n",
    "# GRASS GIS related variables, adapt path to your's\n",
    "grassdata=os.path.join(homedir, \"grassdata\")\n",
    "location=\"odse_workshop2022\"\n",
    "mapset=\"PERMANENT\"\n",
    "\n",
    "# be sure to overwrite existing data when re-running a Jupyter notebook cell\n",
    "GRASS_OVERWRITE=1\n",
    "\n",
    "# Sentinel-2 related directories, adapt path to your's\n",
    "s2_data=os.path.join(homedir, \"sentinel\")\n",
    "s2_timestamps=os.path.join(homedir, \"sentinel-timestamps.txt\")\n",
    "\n",
    "# create directories if not already existing\n",
    "!mkdir -p $grassdata\n",
    "!mkdir -p $s2_data\n",
    "\n",
    "# the variables are also accessible via Python\n",
    "print(homedir)\n",
    "\n",
    "# list content\n",
    "!ls -la $homedir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Next we check the GRASS GIS installation on the updated system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Next we test the installation: you will see version and copyright statements\n",
    "!grass --version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create an empty GRASS GIS location and mapset\n",
    "\n",
    "The GRASS GIS database comes with a structure of folders and subfolders. We create an initial empty folder structure, for now using the GRASS GIS command line (in future also possible with Python). To match the projection of the Sentinel-2 data of Greece used later on, we use [EPSG:32634](http://epsg.io/32634)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create the GRASS GIS location with CRS set to EPSG:32634, then exit\n",
    "!grass -c epsg:32634 -e $homedir/grassdata/$location"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Data upload to notebook session\n",
    "\n",
    "Next we will get the following data set: https://geo.fsv.cvut.cz/geoharmonizer/odse_workshop_2022/grass/geodata.zip\n",
    "\n",
    "Here we transfer the ZIP file to our notebook instance with `wget`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download workshop data into target directory homedir\n",
    "!wget -c https://geo.fsv.cvut.cz/geoharmonizer/odse_workshop_2022/grass/geodata.zip -O $homedir/geodata.zip\n",
    "\n",
    "# unpack workshop data into target directory homedir\n",
    "!unzip -o -q -d $homedir $homedir/geodata.zip\n",
    "\n",
    "# get extra OSM map\n",
    "!wget -c https://data.neteler.org/tmp/osm_greece_landuse.gpkg -O $homedir/geodata/osm_greece_landuse.gpkg\n",
    "\n",
    "print(\"List uploaded file(s) in target directory homedir:\")\n",
    "!ls $homedir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The data structure is as follows:\n",
    "\n",
    "- geodata/highways.gpkg - Greece\n",
    "- geodata/odse_tiles.gpkg - Geo-harmonizer tiles, The Netherlands and Greece\n",
    "- geodata/dtm_5606.tif - elevation model, Greece (30 km x 30 km)\n",
    "- geodata/dtm_14580.tif - elevation model, The Netherlands (30 km x 30 km)\n",
    "- geodata/osm_greece_landuse.gpkg - OpenStreetMap \"landuse\" extract, Greece\n",
    "- geodata/t34sgh_20_60m/* - three Sentinel-2 scenes, Greece (reduced in size)\n",
    "\n",
    "We will later on look at the data here in the Jupyter notebook session.\n",
    "\n",
    "We are now settled and ready to start the GRASS GIS analysis part."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Initialization of GRASS GIS in the Jupyter notebook session\n",
    "\n",
    "While we have already an installation of GRASS GIS, we still need to connect to it to this session."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import subprocess\n",
    "import sys\n",
    "\n",
    "# We ask GRASS GIS where its Python packages are to be able to start them from the notebook\n",
    "sys.path.append(\n",
    "    subprocess.check_output([\"grass\", \"--config\", \"python_path\"], text=True).strip()\n",
    ")\n",
    "\n",
    "# Import GRASS packages\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "# Start GRASS Session\n",
    "session = gj.init(grassdata, location, mapset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point we are in a running GRASS GIS session, let's verify."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show current GRASS GIS settings\n",
    "print(gs.read_command(\"g.gisenv\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Creating an area of interest map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a vector map defined by a 1km buffer around a tile. We first import the tiles vector map `odse_tiles.gpkg`. Note that the CRS of the tiles map (`\"ETRS89 Lambert Azimutal Equal Area\"`) differs from the current location's CRS (EPSG `32634`). It means that GRASS GIS will perform the reprojection of the data into the current location CRS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import\n",
    "gs.run_command(\"v.import\", input=homedir+\"/geodata/odse_tiles.gpkg\", output=\"odse_tiles\")\n",
    "\n",
    "# extract Greece tile from vector map\n",
    "gs.run_command(\"v.extract\", input=\"odse_tiles\", where=\"tile_id = 5606\", output=\"odse_tile_greece\")\n",
    "\n",
    "# show attributes\n",
    "gs.vector_db_select('odse_tile_greece')['values']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we create the 1km buffer around a selected tile, followed by a [Python API](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html?highlight=info%20vector#grass.script.vector.vector_info) metadata query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# buffer single tile\n",
    "gs.run_command(\"v.buffer\", input=\"odse_tile_greece\", where=\"tile_id = 5606\", distance=\"1000\",\n",
    "               output=\"tile_5606_1km\")\n",
    "\n",
    "# show metadata\n",
    "gs.vector_info(\"tile_5606_1km\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Importing geodata into GRASS GIS\n",
    "\n",
    "### Importing the highways vector map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we import the highways vector map `highways.gpkg`. Note that the CRS of highways (EPSG `4326`) differs from the current location's CRS (EPSG `32634`). It means that GRASS GIS will perform the reprojection of the data into the current location CRS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"v.import\", input=homedir+\"/geodata/highways.gpkg\", output=\"highways\")\n",
    "gs.vector_info(\"highways\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing the landuse vector map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we import the landuse vector map `osm_greece_landuse.gpkg` which is based on selected OpenStreetMap tags. Note that the CRS of highways (EPSG `4326`) differs from the current location's CRS (EPSG `32634`). It means that GRASS GIS will perform the reprojection of the data into the current location CRS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"v.import\", input=homedir+\"/geodata/osm_greece_landuse.gpkg\", output=\"osm_greece_landuse\")\n",
    "\n",
    "# show topological metadata\n",
    "gs.vector_info_topo(\"osm_greece_landuse\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Display of imported maps\n",
    "\n",
    "We make use of the `InteractiveMap` tool - interactive visualization with folium."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "greecemap = gj.InteractiveMap(width = 400, tiles=\"OpenStreetMap\")\n",
    "# greecemap.add_vector(\"odse_tile_greece\")\n",
    "# greecemap.add_vector(\"tile_5606_1km\")\n",
    "greecemap.add_vector(\"highways\")\n",
    "greecemap.add_vector(\"osm_greece_landuse\")\n",
    "greecemap.add_layer_control(position = \"bottomright\")\n",
    "greecemap.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Sentinel-2 processing overview\n",
    "\n",
    "There are plenty of libraries or tools which allow downloading\n",
    "Sentinel products from [Copernicus Open Access Hub](https://scihub.copernicus.eu/).\n",
    "\n",
    "For GRASS GIS there is the [i.sentinel](https://grass.osgeo.org/grass-stable/manuals/addons/i.sentinel.html) toolbox. It consists of six GRASS addon modules:\n",
    "\n",
    "* [i.sentinel.download](https://grass.osgeo.org/grass-stable/manuals/addons/i.sentinel.download.html)\n",
    "* [i.sentinel.import](https://grass.osgeo.org/grass-stable/manuals/addons/i.sentinel.import.html)\n",
    "* [i.sentinel.preproc](https://grass.osgeo.org/grass-stable/manuals/addons/i.sentinel.preproc.html)\n",
    "* [i.sentinel.mask](https://grass.osgeo.org/grass-stable/manuals/addons/i.sentinel.mask.html)\n",
    "* [i.sentinel.coverage](https://grass.osgeo.org/grass-stable/manuals/addons/i.sentinel.coverage.html)\n",
    "* [i.sentinel.parallel.download](https://grass.osgeo.org/grass-stable/manuals/addons/i.sentinel.parallel.download.html)\n",
    "\n",
    "Addons modules can be easily installed via [g.extension](https://grass.osgeo.org/grass-stable/manuals/addons/g.extension.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.extension\", extension=\"i.sentinel\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that [i.sentinel.download](https://grass.osgeo.org/grass-stable/manuals/addons/i.sentinel.download.html) requires also the [sentinelsat library](https://pypi.python.org/pypi/sentinelsat) to be installed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip3 install sentinelsat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select and download of Sentinel-2 data\n",
    "\n",
    "Pre-downloaded Sentinel-2 scenes are available in the [sample dataset](https://geo.fsv.cvut.cz/geoharmonizer/odse_workshop_2022/grass/geodata.zip). Workshop participants can skip the download and continue with section \"Import Sentinel-2 data\", see below.\n",
    "\n",
    "#### Steps to select and download of S2 data (please use the pre-downloaded Sentinel-2 scenes instead!)\n",
    "\n",
    "Let’s download suitable Sentinel products for our area of interest (AOI) and perform a NDVI calculation. AOI region is defined by `tile_5606_1km` created above.\n",
    "\n",
    "[Sentinel-2 L2A products](https://www.sentinel-hub.com/blog/sentinel-2-l2a-products-available-sentinel-hub)\n",
    "will be used to avoid computing atmospheric corrections. Let’s\n",
    "search for the latest available product by means of [i.sentinel.download](https://grass.osgeo.org/grass-stable/manuals/addons/i.sentinel.download.html). Setting the `-l` flag, the result will only\n",
    "be printed. The download procedure will be performed later. In order to\n",
    "search and download Sentinel products from the Copernicus Open Access Hub,\n",
    "you have to create an account first. See the manual page of [i.sentinel.download](https://grass.osgeo.org/grass-stable/manuals/addons/i.sentinel.download.html) module for details. Upload or create a new text\n",
    "file in the data directory (homedir) named `esa_credentials.txt` containing two lines (username and password).\n",
    "\n",
    "#### Note\n",
    "\n",
    "To get username and password you need to register at the [Copernicus Open Access Hub](https://scihub.copernicus.eu/), see the [Register new account](https://scihub.copernicus.eu/dhus/#/self-registration) page for signing up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# list available Sentinel-2 L2A scenes for AOI\n",
    "# note that we use parse_command() in order to intercept the output for display in the browser\n",
    "gs.parse_command(\"i.sentinel.download\", flags=\"l\", map=\"tile_5606_1km\", producttype=\"S2MSI2A\",\n",
    "               settings=homedir+\"/esa_credentials.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This should result into something like:\n",
    "\n",
    "```bash\n",
    "c48f04f2-81d1-4fcc-bb32-97405925d6d4 ... 2022-05-25T09:05:59Z  0% S2MSI2A 1.02 GB\n",
    "4720c028-0163-4aba-92a9-3e53fc4047c5 ... 2022-05-15T09:05:59Z  1% S2MSI2A 985.18 MB\n",
    "0eb10137-025e-402b-a050-7e2919602e61 ... 2022-05-15T09:05:59Z  1% S2MSI2A 1.01 GB\n",
    "...\n",
    "```\n",
    "\n",
    "By default the module returns products for the last 60 days. Let’s change\n",
    "the search period setting `start` and `end` options. We will also\n",
    "limit products by `clouds` coverage percentage threshold and `sort` products\n",
    "by ingestion date.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.parse_command(\"i.sentinel.download\", flags=\"l\", map=\"tile_5606_1km\",\n",
    "               producttype=\"S2MSI2A\", settings=homedir+\"/esa_credentials.txt\",\n",
    "               start=\"2022-02-01\", end=\"2022-05-31\", clouds=\"5\",\n",
    "               sort=\"ingestiondate\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This should result in:\n",
    "\n",
    "```bash\n",
    "a119f7b0-0020-4229-b0fb-efd0f5bdfb37 ... 2022-02-19T09:10:21Z  0% S2MSI2A 1.02 GB\n",
    "d6e4253d-80e3-49ff-899e-cee19a0ba935 ... 2022-02-19T09:10:21Z  0% S2MSI2A 998.70 MB\n",
    "cdb37429-1c5d-4c67-8ceb-b6c68a33eef6 ... 2022-03-26T09:06:09Z  0% S2MSI2A 1.02 GB\n",
    "...\n",
    "```\n",
    "\n",
    "#### Tip: If more products have been found, you can limit the amount with the `limit` option.\n",
    "\n",
    "\n",
    "Let’s download the desired product(s). Just remove the `-l` flag and add the `output` option in order to define the path to the output directory where data should be saved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# note: s2_data directory and esa_credentials.txt have been defined above\n",
    "\n",
    "## TO BE SKIPPED, we use pre-downloaded data below!\n",
    "#print(\"Storing S2 data in <\"+s2_data+\"> - this takes a few minutes...\")\n",
    "#gs.run_command(\"i.sentinel.download\", map=\"tile_5606_1km\", producttype=\"S2MSI2A\",\n",
    "#               settings=homedir+\"/esa_credentials.txt\", start=\"2022-02-01\", end=\"2022-05-31\",\n",
    "#               clouds=\"5\", output=s2_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing Sentinel-2 data\n",
    "\n",
    "Before importing or linking Sentinel-2 data we print a list of\n",
    "filtered raster files including projection match (second column, 1 for\n",
    "match, otherwise 0). If the CRS of input data differs from the current location\n",
    "consider reprojection (`-r`) or creating a new location for import.\n",
    "\n",
    "*Important*: Data will be imported into the new location by means of the [i.sentinel.import](https://grass.osgeo.org/grass-stable/manuals/addons/i.sentinel.import.html) tool. The command will import **all** Sentinel bands from `input` directory\n",
    "recursively. Before importing data let’s check content of the input directory by `-p` flag."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check list of pre-downloaded Sentinel-2 scenes, with i.sentinel.import (-p: print)\n",
    "gs.parse_command(\"i.sentinel.import\", flags=\"p\", input=homedir+\"/geodata/t34sgh_20_60m\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This should result in (projection match in second column: 1 for match, otherwise 0):\n",
    "\n",
    "```bash\n",
    ".../R20m/T34SGH_20210828T090549_B03_20m.jp2 1 (EPSG: 32634)\n",
    ".../R20m/T34SGH_20210828T090549_B04_20m.jp2 1 (EPSG: 32634)\n",
    ".../R20m/T34SGH_20210828T090549_B8A_20m.jp2 1 (EPSG: 32634)\n",
    "...\n",
    "```\n",
    "\n",
    "In the example below, we limit the S2 data import to the RGB and NIR bands (2, 3, 4, 8A) in 20m spatial resolution by the `pattern` option."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.parse_command(\"i.sentinel.import\", flags=\"p\", input=homedir+\"/geodata/t34sgh_20_60m\", pattern=\"B(02|03|04|8A)_20m\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This should result in:\n",
    "\n",
    "```bash\n",
    ".../R20m/T34SGH_20210828T090549_B8A_20m.jp2 1 (EPSG: 32634)\n",
    ".../R20m/T34SGH_20210828T090549_B04_20m.jp2 1 (EPSG: 32634)\n",
    ".../R20m/T34SGH_20210624T090601_B04_20m.jp2 1 (EPSG: 32634)\n",
    ".../R20m/T34SGH_20210624T090601_B8A_20m.jp2 1 (EPSG: 32634)\n",
    ".../R20m/T34SGH_20210729T090559_B04_20m.jp2 1 (EPSG: 32634)\n",
    ".../R20m/T34SGH_20210729T090559_B8A_20m.jp2 1 (EPSG: 32634)\n",
    "```\n",
    "\n",
    "By default, input data are imported into GRASS data format.\n",
    "Alternatively, data can be linked if `-l` is given. It is also\n",
    "useful to import cloud mask vector features by `-c` flag. We also use\n",
    "`register_output` option to produce a timestamp plain text file\n",
    "which will be used in [section 07](07.rst)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use 2GB of RAM for faster operations, s2_timestamps defined above\n",
    "# this takes up to a few minutes...\n",
    "gs.parse_command(\"i.sentinel.import\", flags=\"c\", input=homedir+\"/geodata/t34sgh_20_60m\",\n",
    "               pattern=\"B(02|03|04|8A)_20m\", memory=2000, register_output=s2_timestamps)\n",
    "\n",
    "# print timestamp file for inspection\n",
    "f = open(s2_timestamps, 'r')\n",
    "content = f. read()\n",
    "print(content)\n",
    "f. close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Semantic labels\n",
    "\n",
    "A fairly new concept within GRASS GIS is semantic labels. These are especially relevant for satellite imagery as they allow us to identify to which sensor and band a given raster corresponds. These labels are particularly relevant when working with satellite image collections and also when classifying different scenes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print semantic labels\n",
    "for i in gs.list_grouped(type=\"raster\")['PERMANENT']:\n",
    "    label = gs.raster_info(map=i)['semantic_label']\n",
    "    print('Map: {}, Semantic label: {}'.format(i,label))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9. Computing NDVI\n",
    "\n",
    "To compute [NDVI - Normalized difference vegetation index](http://en.wikipedia.org/wiki/NDVI),\n",
    "the red (VIS) and near-infrared (NIR) channels are required. In the case\n",
    "of Sentinel-2, these are the 4th and 8th band, respectively.\n",
    "\n",
    "<code>NDVI = (NIR - VIS) / (NIR  + VIS)</code>\n",
    "\n",
    "Before computing, let’s set up the computational region based on the tile map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# note that we use parse_command() to immediately see the printed computational region\n",
    "gs.parse_command(\"g.region\", vector=\"tile_5606_1km\", flags=\"p\")\n",
    "gs.run_command(\"v.in.region\", output=\"greece_box\")\n",
    "gs.parse_command(\"g.region\", vector=\"greece_box\", align=\"T34SGH_20210624T090601_B04_20m\", flags=\"p\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NDVI can be computed via the [i.vi](https://grass.osgeo.org/grass-stable/manuals/i.vi.html)\n",
    "tool which allows computing various vegetation indices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"i.vi\", red=\"T34SGH_20210624T090601_B04_20m\", viname=\"ndvi\",\n",
    "               nir=\"T34SGH_20210624T090601_B8A_20m\", output=\"ndvi\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add semantic label \"NDVI\" to the resulting `ndvi` map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.support\", map=\"ndvi\", semantic_label=\"S2_NDVI\")\n",
    "gs.raster_info(map=\"ndvi\")[\"semantic_label\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Displaying the NDVI map in Folium (Leaflet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display newly created NDVI map\n",
    "ndvimap = gj.InteractiveMap(width = 400, tiles=\"OpenStreetMap\")\n",
    "ndvimap.add_raster(\"ndvi\", opacity=0.8)\n",
    "ndvimap.add_vector(\"tile_5606_1km\")\n",
    "# ndvimap.add_vector(\"greece_box\")\n",
    "ndvimap.add_layer_control(position = \"bottomright\")\n",
    "ndvimap.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 10. NDVI time series data processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A few concepts of time series data processing in GRASS GIS\n",
    "\n",
    "GRASS GIS offers specialized tools for spatio-temporal data\n",
    "processing, see GRASS documentation [temporalintro](https://grass.osgeo.org/grass-stable/manuals/temporalintro.html) for details.\n",
    "\n",
    "GRASS introduces three special data types that are designed to handle\n",
    "time-series data:\n",
    "\n",
    "* *Space-time raster datasets* (`strds`) for managing raster map\n",
    "  time series.\n",
    "\n",
    "* *Space-time 3D raster datasets* (`str3ds`) for managing 3D raster\n",
    "  map time series.\n",
    "\n",
    "* *Space-time vector datasets* (`stvds`) for managing vector map time\n",
    "  series."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create space-time dataset\n",
    "\n",
    "At this moment a new space-time dataset can be created by means of [t.create](https://grass.osgeo.org/grass-stable/manuals/t.create.html) and all imported Sentinel bands registered with [t.register](https://grass.osgeo.org/grass-stable/manuals/t.register.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"t.create\", output=\"s2_tile_5606\", title=\"Sentinel L2A 2021\", desc=\"Tile 5606\")\n",
    "gs.run_command(\"t.register\", input=\"s2_tile_5606\", file=s2_timestamps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s check basic metadata (see [t.info](https://grass.osgeo.org/grass-stable/manuals/t.info.html)) and list the registered maps ([t.rast.list](https://grass.osgeo.org/grass-stable/manuals/t.rast.list.html))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.parse_command(\"t.info\", input=\"s2_tile_5606\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This should result in:\n",
    "```bash\n",
    "...\n",
    "| ... ------ Space Time Raster Dataset ------- ...\n",
    "...\n",
    "| Id: ...... s2_tile_5606@PERMANENT\n",
    "...\n",
    "| ... ------ Absolute time ------- ...\n",
    "...\n",
    "| Start time:................. 2021-06-24 09:19:52.607078\n",
    "| End time:................... 2021-08-28 09:19:49.080855\n",
    "...\n",
    "| ... ------ Metadata information ------- ...\n",
    "...\n",
    "| Number of registered maps:.. 12\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### List registered bands in space-time cube"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.parse_command(\"t.rast.list\", input=\"s2_tile_5606\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This should result in:\n",
    "\n",
    "```bash\n",
    "name|mapset|start_time|end_time\n",
    "T34SGH_20210624T090601_B04_20m|PERMANENT|2021-06-24 09:19:52.607078|None\n",
    "T34SGH_20210624T090601_B8A_20m|PERMANENT|2021-06-24 09:19:52.607078|None\n",
    "...\n",
    "T34SGH_20210729T090559_B04_20m|PERMANENT|2021-07-29 09:19:53.186492|None\n",
    "T34SGH_20210729T090559_B8A_20m|PERMANENT|2021-07-29 09:19:53.186492|None\n",
    "T34SGH_20210828T090549_B04_20m|PERMANENT|2021-08-28 09:19:49.080855|None\n",
    "T34SGH_20210828T090549_B8A_20m|PERMANENT|2021-08-28 09:19:49.080855|None\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NDVI Space-Time computation\n",
    "\n",
    "For NDVI computation 4th and 8th bands are required. Map algebra for spatio-temporal data is performed by [t.rast.algebra](https://grass.osgeo.org/grass-stable/manuals/t.rast.algebra.html) which requires bands separated into different spatio-temporal datasets. Such datasets can\n",
    "be prepared by [t.rast.extract](https://grass.osgeo.org/grass-stable/manuals/t.rast.extract.html).\n",
    "\n",
    "Note that data aren't replicated by this extraction as the space-time support in GRASS GIS is metadata based."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"t.rast.extract\", input=\"s2_tile_5606\", where=\"name like '%B04%'\", output=\"s2_b4\")\n",
    "gs.run_command(\"t.rast.extract\", input=\"s2_tile_5606\", where=\"name like '%B8A%'\", output=\"s2_b8\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s check content of the new datasets by [t.rast.list](https://grass.osgeo.org/grass-stable/manuals/t.rast.list.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.parse_command(\"t.rast.list\", input=\"s2_b4\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.parse_command(\"t.rast.list\", input=\"s2_b8\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the computational region by [g.region](https://grass.osgeo.org/grass-stable/manuals/g.region.html)\n",
    "including a mask for the area of interest by [r.mask](https://grass.osgeo.org/grass-stable/manuals/r.mask.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.parse_command(\"g.region\", vector=\"tile_5606_1km\", align=\"T34SGH_20210624T090601_B04_20m\", flags=\"p\")\n",
    "gs.run_command(\"r.mask\", vector=\"tile_5606_1km\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: NDVI (see above) computation on spatio-temporal datasets can be performed in parallel (`nproc` parameter)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"t.rast.algebra\", basename=\"ndvi\", nproc=\"2\",\n",
    "               expression=\"ndvi = float(s2_b8 - s2_b4) / ( s2_b8 + s2_b4 )\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When computation is finished *ndvi* color table can be set with [t.rast.colors](https://grass.osgeo.org/grass-stable/manuals/t.rast.colors.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"t.rast.colors\", input=\"ndvi\", color=\"ndvi\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Show metadata of NDVI time series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.parse_command(\"t.info\", input=\"ndvi\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Query time series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# query region center coordinates for query, in UTM34N\n",
    "gs.parse_command(\"g.region\", flags=\"c\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# query map at center coordinates\n",
    "gs.parse_command(\"t.rast.what\", strds=\"ndvi\", coordinates=\"748800,4226800\", layout=\"col\", flags=\"n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Display NDVI time series in Folium (Leaflet)\n",
    "\n",
    "\n",
    "Note: [TimeSeriesMap()](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.jupyter.html?highlight=timeseriesmap#module-grass.jupyter.timeseriesmap) of `grass.jupyter` is still experimental and under development.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### TO BE SKIPPED - in GRASS GIS 8.2.0 it takes \"forever\", bugfix pending.\n",
    "\n",
    "\n",
    "## reduce resolution for faster display of time series, save original first for later\n",
    "#gs.parse_command(\"g.region\", save=\"default_res\")\n",
    "#gs.parse_command(\"g.region\", flags=\"pa\", res=50)\n",
    "# \n",
    "## Display newly created NDVI time series map\n",
    "#ndviseries = gj.TimeSeriesMap(use_region=True)\n",
    "#ndviseries.add_raster_series(\"ndvi\", fill_gaps=False)\n",
    "#ndviseries.d_legend(color=\"black\", at=(10,40,2,6))\n",
    "#ndviseries.d_barscale()\n",
    "#ndviseries.show()  # Create TimeSlider\n",
    "\n",
    "# optionally, write out to animated GIF\n",
    "# ndviseries.save(\"image.gif\")\n",
    "\n",
    "## restore original region\n",
    "#gs.parse_command(\"g.region\", region=\"default_res\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 11. Creating an image stack (imagery group)\n",
    "\n",
    "**Stack of maps = imagery group**\n",
    "\n",
    "When you work with a stack of raster maps (e.g., R-G-B channels or more) in GRASS GIS, you can best handle this stack by creating a group with [i.group](https://grass.osgeo.org/grass-stable/manuals/i.group.html). It is just based on metadata, so it does not take up more disk space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate list of selected S2 maps\n",
    "s2_maps = gs.list_grouped(type=\"raster\", pattern=\"*20210624T090601*\")['PERMANENT']\n",
    "print(s2_maps)\n",
    "\n",
    "# since imagery groups can not be overwritten, we first delete any leftover \"s2\" group from previous runs\n",
    "gs.run_command(\"g.remove\", type=\"group\", name=\"s2\", flags=\"f\")\n",
    "\n",
    "# create group and subgroup with S2 bands\n",
    "gs.run_command(\"i.group\", group=\"s2\", subgroup=\"s2\", input=s2_maps)\n",
    "gs.parse_command(\"i.group\", group=\"s2\", flags=\"l\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12. Object recognition with image segmentation\n",
    "\n",
    "As a first step for the image segmentation of Sentinel-2 bands with [i.segment](https://grass.osgeo.org/grass-stable/manuals/i.segment.html), we set the computational region to our area of interest (here driven by the available landuse data), while matching the pixel geometry of the Sentinel-2 red band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define pixel geometry for new raster map: use extent from vector map,\n",
    "# and align pixel position and size from Sentinel-2 20m band\n",
    "gs.parse_command(\"g.region\", vector=\"osm_greece_landuse\",\n",
    "                 align=\"T34SGH_20210624T090601_B04_20m\", flags=\"p\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to perform an image segmentation (btw.: it works with any raster map, not only image bands!), we add the previously created NDVI map to the `s2` imagery group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add NDVI map to group and subgroup already populated with S2 bands\n",
    "gs.run_command(\"i.group\", group=\"s2\", subgroup=\"s2\", input=\"ndvi\")\n",
    "# list content of group\n",
    "gs.parse_command(\"i.group\", group=\"s2\", flags=\"l\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to perform the segmentation. In addition to the segmentation map we also create a goodness of fit estimate map from the data contained in the \"s2\" group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Threshold = 0 merges only identical segments; threshold = 1 merges all\n",
    "# (you may try different values)\n",
    "gs.run_command(\"i.segment\", group=\"s2\", threshold=\"0.05\", minsize=\"100\", output=\"sentinel_segments_min100\", goodness=\"sentinel_segments_goodness_min100\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# display newly created segments raster map\n",
    "segments = gj.InteractiveMap(width = 400, use_region=True)\n",
    "# choose low opacity to better compare segments with OpenStreetMap\n",
    "segments.add_raster(\"sentinel_segments_min100\", opacity=0.3)\n",
    "# segments.add_raster(\"ndvi\", opacity=0.3)\n",
    "segments.add_layer_control(position = \"bottomright\")\n",
    "segments.show()\n",
    "\n",
    "# use layer control to only see selected maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show univariate statistics of goodness-of-fit raster map, with extended statistics (quartiles)\n",
    "gs.parse_command(\"r.univar\", map=\"sentinel_segments_goodness_min100\", flags=\"ge\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add color table (low fit values: red; high fit values: green)\n",
    "gs.run_command(\"r.colors\", map=\"sentinel_segments_goodness_min100\", color=\"byg\", flags=\"e\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# display newly created goodness-of-fit raster map\n",
    "segments = gj.InteractiveMap(width = 400, use_region=True)\n",
    "segments.add_raster(\"sentinel_segments_goodness_min100\", opacity=0.8)\n",
    "segments.add_layer_control(position = \"bottomright\")\n",
    "segments.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 13. Supervised Classification: RandomForest\n",
    "\n",
    "We will now demonstrate a very much simplified workflow to perform a supervised [RandomForest classification](https://en.wikipedia.org/wiki/Random_forest).\n",
    "\n",
    "We will feed the following data into the model:\n",
    "\n",
    "- NDVI map (created above)\n",
    "- image segmentation (created above)\n",
    "- random points extracted from the training polygons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we inspect the available raster bands, just to recall their names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.parse_command(\"g.list\", type=\"raster\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show the R-G-B composite of our chosen Sentinel-2 scene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s2map = gj.Map(width=400, use_region=True)\n",
    "# show as a static map\n",
    "s2map.d_rgb(red=\"T34SGH_20210624T090601_B04_20m\", \n",
    "              green=\"T34SGH_20210624T090601_B03_20m\", \n",
    "              blue=\"T34SGH_20210624T090601_B02_20m\")\n",
    "s2map.show()\n",
    "\n",
    "# ... this will expectedly look very dark. Color balancing next!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# perform color auto-balancing for RGB bands (also other band combinations may be auto-balanced)\n",
    "gs.run_command(\"i.colors.enhance\", \n",
    "               red=\"T34SGH_20210624T090601_B04_20m\",\n",
    "               green=\"T34SGH_20210624T090601_B03_20m\", \n",
    "               blue=\"T34SGH_20210624T090601_B02_20m\",\n",
    "               strength=90)\n",
    "\n",
    "# display again\n",
    "s2map = gj.Map(width=400, use_region=True)\n",
    "s2map.d_rgb(red=\"T34SGH_20210624T090601_B04_20m\", \n",
    "              green=\"T34SGH_20210624T090601_B03_20m\", \n",
    "              blue=\"T34SGH_20210624T090601_B02_20m\")\n",
    "s2map.d_vect(map=\"osm_greece_landuse\")\n",
    "s2map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creation of a training map by sampling from existing data\n",
    "\n",
    "We will perform stratified sampling on the `osm_greece_landuse` map which is an incomplete extract from [OpenStreetMap](https://www.openstreetmap.org/#map=11/38.2050/23.8266) by using selected tags (aiming at matching the CORINE legend).\n",
    "\n",
    "As a first step we rasterize the map, matching the pixel geometry of the Sentinel-2 red band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (already done, but... :-) define pixel geometry for new raster map: use extent from vector map,\n",
    "# and align pixel position and size from Sentinel-2 20m band\n",
    "gs.parse_command(\"g.region\", vector=\"osm_greece_landuse\",\n",
    "                 align=\"T34SGH_20210624T090601_B04_20m\", flags=\"p\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# list column names of vector map\n",
    "gs.vector_columns(\"osm_greece_landuse\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now we convert the landuse map from vector to raster model to simplify the next steps\n",
    "# (names may be identical since the data type is different)\n",
    "gs.run_command(\"v.to.rast\", input=\"osm_greece_landuse\", output=\"osm_greece_landuse\", \n",
    "               use=\"attr\", attribute_column=\"CORINE_class\", label_column=\"landuse\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simplification of the training classes\n",
    "\n",
    "Since we have (in this example) only a few Sentinel-2 bands available, we will simplify the legend prior to using the training data for the training of the model:\n",
    "\n",
    "- 11 12 = 11 urban\n",
    "- 13    = 13 agriculture\n",
    "- 14 21 22 23 32 = 14 greenland\n",
    "- 31    = 31 forest\n",
    "- 51    = 51 water\n",
    "- all other  = NULL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to simplify teaching life, we just download the class redefinition file\n",
    "!wget -c https://data.neteler.org/tmp/lulc_rules.recl -O $homedir/lulc_rules.recl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Way 1: use Python\n",
    "#rules = [\"11 12 = 11 urban\", \"13 = 13 agriculture\", \"14 21 22 23 32 = 14 greenland\", \"31 = 31 forest\", \"* = NULL\"]\n",
    "#p = gs.feed_command(\"r.reclass\", input=\"landuse_train\", output=\"landuse_train_simplified2\",\n",
    "#                          rules=\"-\")\n",
    "#p.stdin.write(\"\\n\".join(rules).encode())\n",
    "#p.stdin.close()\n",
    "\n",
    "# Way 2: use the rules file which we downloaded\n",
    "# raster map: apply class simplification\n",
    "gs.run_command(\"r.reclass\", input=\"osm_greece_landuse\", output=\"osm_greece_landuse_simplified\",\n",
    "                 rules=homedir+\"/lulc_rules.recl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check raster categories of simplified training map\n",
    "gs.parse_command(\"r.category\", map=\"osm_greece_landuse_simplified\", separator=\"comma\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set colors\n",
    "colours = [\"11 195:20:0\", \"13 255,255,168\", \"14 100:240:100\", \"31 77:105:0\", \"51 0:204:242\"]\n",
    "colourise = gs.feed_command(\"r.colors\", map=\"osm_greece_landuse_simplified\", rules=\"-\", quiet=True)\n",
    "colourise.stdin.write(\"\\n\".join(colours).encode())\n",
    "colourise.stdin.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# display newly created raster map\n",
    "osmlulc = gj.InteractiveMap(width = 400, use_region=True, tiles=\"OpenStreetMap\")\n",
    "osmlulc.add_raster(\"osm_greece_landuse_simplified\", opacity=0.8)\n",
    "osmlulc.add_layer_control(position = \"bottomright\")\n",
    "osmlulc.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# show simple legend\n",
    "legend = gj.Map(width=400, use_region=True)\n",
    "# at=bottom,top,left,right, percentage of screen coordinates (0,0 is lower left)\n",
    "legend.d_legend(raster=\"osm_greece_landuse_simplified\", title=\"Classes\",\n",
    "                fontsize=20, at=(20, 80, 20, 80), flags=\"n\")\n",
    "legend.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random sampling from rasterized simplified landuse map\n",
    "\n",
    "We now perform stratified sampling, i.e. we extract for each land use class `n` sampling points, using the GRASS GIS addon [r.sample.category](https://grass.osgeo.org/grass-stable/manuals/addons/r.sample.category.html).\n",
    "\n",
    "First, we install this addon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.extension\", extension=\"r.sample.category\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# stratified random sampling, generated vector points\n",
    "gs.run_command(\"r.sample.category\", input=\"osm_greece_landuse_simplified\", output=\"landuse_train\", n=\"150\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# display newly created vector points map\n",
    "train = gj.InteractiveMap(width = 400, use_region=True)\n",
    "train.add_raster(\"osm_greece_landuse_simplified\", opacity=0.7)\n",
    "train.add_vector(\"landuse_train\")\n",
    "train.add_layer_control(position = \"bottomright\")\n",
    "train.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# list column names of vector points map\n",
    "gs.vector_columns(\"landuse_train\", getDict=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show vector attribute table\n",
    "gs.vector_db_select(\"landuse_train\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the machine learning classifier expects raster points, we convert the vector sampling points accordingly using  [v.to.rast](https://grass.osgeo.org/grass-stable/manuals/v.to.rast.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert points from vector to raster model\n",
    "# (names may remain identical since the data type is different)\n",
    "gs.run_command(\"v.to.rast\", input=\"landuse_train\", output=\"landuse_train\", \n",
    "               use=\"attr\", attribute_column=\"osm_greece_landuse_simplified\", label_column=\"label\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check raster categories\n",
    "gs.parse_command(\"r.category\", map=\"landuse_train\", separator=\"comma\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show univariate statistics\n",
    "gs.parse_command(\"r.univar\", map=\"landuse_train\", flags=\"g\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# display newly created raster map\n",
    "# NOTE: zoom in to spot the raster sampling points ;-)\n",
    "train = gj.InteractiveMap(width = 400, use_region=True)\n",
    "train.add_raster(\"landuse_train\", opacity=0.8)\n",
    "train.add_layer_control(position = \"bottomright\")\n",
    "train.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add segmentation map created above to group and subgroup already populated with S2 bands and NDVI\n",
    "gs.run_command(\"i.group\", group=\"s2\", subgroup=\"s2\", input=\"sentinel_segments_min100\")\n",
    "# list content of group\n",
    "gs.parse_command(\"i.group\", group=\"s2\", flags=\"l\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform machine learning model training (RandomForest)\n",
    "\n",
    "First we have to install the [r.learn.ml2](https://grass.osgeo.org/grass-stable/manuals/addons/r.learn.ml2.html) extention. It consists of two modules: `r.learn.train` and `r.learn.predict`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# install ML extension\n",
    "gs.run_command(\"g.extension\", extension=\"r.learn.ml2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now train the ML model using [r.learn.train](https://grass.osgeo.org/grass-stable/manuals/addons/r.learn.train.html), with model \"RandomForestClassifier\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# train a random forest classification model using r.learn.train\n",
    "gs.run_command(\"r.learn.train\", group=\"s2\", training_map=\"landuse_train\",\n",
    "               model_name=\"RandomForestClassifier\",\n",
    "               n_estimators=\"500\", save_model=homedir+\"/rf_model.gz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " The model has been stored in the file `rf_model.gz` for use in the prediction step of the supervised classification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls -la $homedir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform ML supervised classification\n",
    "\n",
    "The trained model will now be applied to the entire dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# perform prediction using r.learn.predict, takes a few minutes\n",
    "gs.run_command(\"r.learn.predict\", group=\"s2\", load_model=homedir+\"/rf_model.gz\", output=\"sentinel_rf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add some styling: category labels and colors (colors inspired by CORINE land cover map legend)\n",
    "\n",
    "class | R:G:B | label\n",
    "- 11 195:20:0 urban\n",
    "- 13 255,255,168 agriculture\n",
    "- 14 100:240:100 greenland\n",
    "- 31 77:105:0 forest\n",
    "- 51 0:204:242 water"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set colors\n",
    "colours = [\"11 195:20:0\", \"13 255,255,168\", \"14 100:240:100\", \"31 77:105:0\", \"51 0:204:242\"]\n",
    "colourise = gs.feed_command(\"r.colors\", map=\"sentinel_rf\", rules=\"-\", quiet=True)\n",
    "colourise.stdin.write(\"\\n\".join(colours).encode())\n",
    "colourise.stdin.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this the (oversimplified) supervised classification has been completed and we can display the result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Supervised classification: reporting and display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display newly created sentinel_rf map\n",
    "map = gj.InteractiveMap(width = 400, tiles=\"OpenStreetMap\")\n",
    "map.add_raster(\"sentinel_rf\", opacity=0.8)\n",
    "map.add_layer_control(position = \"bottomright\")\n",
    "map.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show simple legend\n",
    "legend = gj.Map(width=400, use_region=True)\n",
    "# at=bottom,top,left,right, percentage of screen coordinates (0,0 is lower left)\n",
    "legend.d_legend(raster=\"sentinel_rf\", title=\"Classes\",\n",
    "                fontsize=20, at=(10, 90, 10, 90), flags=\"n\")\n",
    "legend.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show class distribution in percent\n",
    "gs.parse_command(\"r.report\", map=\"sentinel_rf\", units=\"p\", flags=\"h\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# export map to COG\n",
    "gs.run_command(\"r.out.gdal\", flags=\"fmt\", input=\"sentinel_rf\", output=homedir+\"/greece_sentinel2_RF.tif\",\n",
    "               format=\"COG\", overviews=\"4\")\n",
    "\n",
    "# see the output file\n",
    "!ls -la $homedir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Keep in mind, this classification was just a simplified example to show how the procedure works."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 14. What's next?\n",
    "\n",
    "You may enjoy more Jupyter notebooks at\n",
    "\n",
    "https://github.com/OSGeo/grass/tree/main/doc/notebooks\n",
    "\n",
    "### Talk to us\n",
    "\n",
    "- Martin Landa, PhD, https://geo.fsv.cvut.cz/gwiki/Ing._Martin_Landa,_Ph.D.\n",
    "- Markus Neteler, PhD, https://www.mundialis.de/en/neteler/"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}