colab_workshop.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "qibM9kYhQ1uw"
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "# install dependencies\n",
    "DEBIAN_FRONTEND=noninteractive sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable && apt update && apt install grass subversion grass-dev parallel libglu1-mesa-dev xvfb && apt remove libproj22\n",
    "\n",
    "pip install PyVirtualDisplay\n",
    "\n",
    "cd ~\n",
    "\n",
    "# download sample data\n",
    "mkdir -p grassdata\n",
    "curl -SL https://grass.osgeo.org/sampledata/north_carolina/nc_spm_08_grass7.zip > nc_spm_08_grass7.zip\n",
    "unzip -qq nc_spm_08_grass7.zip\n",
    "mv nc_spm_08_grass7 grassdata\n",
    "rm nc_spm_08_grass7.zip\n",
    "wget https://zenodo.org/record/6967699/files/FOSS4G_2022_GRASS_GIS_workshop.zip\n",
    "unzip FOSS4G_2022_GRASS_GIS_workshop.zip\n",
    "rm FOSS4G_2022_GRASS_GIS_workshop.zip"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "grt_FB6epyVk"
   },
   "outputs": [],
   "source": [
    "import os\n",
    "os.chdir(os.path.expanduser(\"~\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 1: GRASS GIS intro\n",
    "\n",
    "In this first part, we will demonstrate starting GRASS GIS, creating new project and basic data visualization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notebook environment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default all cells are running Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "v = sys.version_info\n",
    "print(f\"We are using Python {v.major}.{v.minor}.{v.micro}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cells can also run bash using IPython magic. Later, we will use this to show how GRASS can be used from shell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "pwd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or we can use ! to run individual lines in bash:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pwd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Start GRASS GIS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we create new empty location (project) called *dix_park* that uses projection [UTM zone 17 N](https://epsg.io/6346) with EPSG:6346.\n",
    "\n",
    "Flag `c` stands for _creating_ new location and `e` will _exit_ the command after creating the location. See [manual](https://grass.osgeo.org/grass-stable/manuals/grass.html) for more examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "grass -c EPSG:6346 -e ~/grassdata/dix_park/ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start GRASS GIS session in the newly created location. Load Python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Python standard library and IPython packages we need.\n",
    "import subprocess\n",
    "import sys\n",
    "\n",
    "# Ask GRASS GIS where its Python packages are.\n",
    "sys.path.append(\n",
    "    subprocess.check_output([\"grass\", \"--config\", \"python_path\"], text=True).strip()\n",
    ")\n",
    "\n",
    "# Import the GRASS GIS packages we need.\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "# Start GRASS Session\n",
    "session = gj.init(\"~/grassdata\", \"dix_park\", \"PERMANENT\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_If you run this workshop online, data is already prepared in the environment. Otherwise download the data [here](http://fatra.cnr.ncsu.edu/foss4g2021/)._\n",
    "\n",
    "We will import prepared digitial surface model (DSM), bare ground (digital terrain model, DTM) and ortho maps. \n",
    "The data CRS matches the CRS of the *dix_park* location, so we don't need to reproject it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "r.import input=dsm.tif output=dsm resample=bilinear\n",
    "r.import input=ground.tif output=ground\n",
    "r.import input=ortho.tif output=ortho"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will import pre-downloaded OSM data of roads restricted to our study area. We obtained the roads using Overpass Turbo with this [query](https://overpass-turbo.eu/s/1aGu) and exported to GeoJSON. \n",
    "\n",
    "The data comes in EPSG:4326, so it will be automatically reprojected to UTM during the import."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "v.import input=roads.geojson output=roads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the available data in our location:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.list type=raster,vector -m -t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data visualization\n",
    "\n",
    "We will call d.rast/d.vect modules using the new [GRASS Jupyter API](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.jupyter.html). The `Map` class creates and displays GRASS maps as static PNG images.\n",
    "\n",
    "First let's display bare ground elevation with the roads dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Map instance\n",
    "elevation_map = gj.Map()\n",
    "# Add a raster and vector to the map\n",
    "elevation_map.d_rast(map=\"ground\")\n",
    "elevation_map.d_vect(map=\"roads\")\n",
    "elevation_map.d_legend(raster=\"ground\")\n",
    "# Display map\n",
    "elevation_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, display road _Umstead Drive_ in yellow on top of ortho. Method `d_vect` can be called with options of [d.vect](https://grass.osgeo.org/grass-stable/manuals/d.vect.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ortho_map = gj.Map()\n",
    "ortho_map.d_rast(map=\"ortho\")\n",
    "ortho_map.d_vect(map=\"roads\", width=\"2\", color=\"yellow\", where=\"name = 'Umstead Drive'\")\n",
    "ortho_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is how we can visualize data interactively with folium:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roads_map = gj.InteractiveMap()\n",
    "roads_map.add_raster(\"dsm\", opacity=0.5)\n",
    "roads_map.add_vector(\"roads\")\n",
    "roads_map.add_layer_control(position=\"bottomright\")\n",
    "roads_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also visualize data in 3D. Here we drape the ortho over the DSM.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map3d = gj.Map3D()\n",
    "map3d.render(elevation_map=\"dsm\", color_map=\"ortho\",\n",
    "             position=(0.5, 1), height=3000, perspective=12)\n",
    "map3d.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map3d = gj.Map3D()\n",
    "map3d.render(elevation_map=\"dsm\", resolution_fine=1, color_map=\"ortho\",\n",
    "             light_position=(1, 0, 0.5),\n",
    "             position=(0.75, 0.35), height=1500, perspective=10)\n",
    "map3d.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GRASS GIS tools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GRASS functionality is available through tools (also called modules). There are over 500 different tools in the core distribution and over 300 addon tools that can be used to prepare and analyze data.\n",
    "\n",
    "Tools respect the following naming conventions:\n",
    "\n",
    "Prefix | Function | Example\n",
    "------ | -------- | -------\n",
    "r.* | raster processing | r.mapcalc: map algebra\n",
    "v.*\t| vector processing\t| v.clean: topological cleaning\n",
    "i.*\t| imagery processing | i.segment: object recognition\n",
    "db.* | database management | db.select: select values from table\n",
    "r3.* | 3D raster processing | r3.stats: 3D raster statistics\n",
    "t.* | temporal data processing | t.rast.aggregate: temporal aggregation\n",
    "g.* | general data management | g.rename: renames map\n",
    "d.* | display | d.rast: display raster map\n",
    "\n",
    "Note also that some tools have multiple dots in their names. For example, tools staring with v.net.* deal with vector network analysis and r.in.* tools import raster data into GRASS GIS spatial database.\n",
    "\n",
    "There is a tool for finding other tools:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.search.modules keyword=zonal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is how to get all options and flags of a GRASS tool:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "r.univar --help"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will open the tool's manual page. Remove -m flag to open HTML page in your browser (will work only locally)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.manual r.univar -m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GRASS modules can be executed either through the GUI, command line or Python interfaces. This is an example how to execute a tool in command line. Specifically, it will extract road _Umstead Drive_ into a new vector `umstead_drive_segments`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "v.extract input=roads where=\"name = 'Umstead Drive'\" output=umstead_drive_segments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the same in Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"v.extract\", input=\"roads\", where=\"name = 'Umstead Drive'\", output=\"umstead_drive_segments\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computational region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computational region is an important raster concept in GRASS GIS.\n",
    "Before we use a module to compute a new raster map, we need to properly set the computational region. All raster computations will be performed in the specified extent and with the given resolution.\n",
    "Among other things, this allows us to easily subset larger extent data for quicker testing of analysis, or to run an analysis of specific regions given by e.g. administrative units.\n",
    "\n",
    "A few points to keep in mind:\n",
    "\n",
    " * computational region is defined by region extent and raster resolution\n",
    " * applies to all raster operations and vector to raster operations\n",
    " * persists between GRASS sessions, can be different for different mapsets\n",
    " * advantages: keeps your results consistent, avoids clipping, for computationally demanding tasks set region to smaller extent, check that your result is good and then set the computational region to the entire study area and rerun analysis\n",
    " * run `g.region -p` to see current region settings\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.region -p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most common way to set region is **based on a raster map** - both extent and resolution. Run again g.region (we include -p flag to always see the resulting region):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.region -p raster=dsm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computational region can be set also **using a vector map**. In that case, only extent is set (as vector maps do not have any resolution - at least not in the way raster maps do). In the command line, it looks like this:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.region -p vector=roads | tail"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However now the resolution was adjusted based on the extent of the vector map, it is no longer a nice rounded number. If that's not desired, we can set it explicitly using -a flag and parameter res. Now the resolution is aligned to even multiples of 2 (the units are the units of the current location, in our case meters):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.region -p vector=roads res=2 -a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often we need to set the computational extent based on a vector map, but take the resolution and alignment from a raster map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.region -p vector=roads align=dsm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can save a specific region for later. This won't actually modify the current region:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.region vector=roads align=dsm save=roads_region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to later set that named region, we can do it simply with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.region region=roads_region -p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python API"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two Python APIs for accessing a tool's functionality - [GRASS GIS Python Scripting Library](https://grass.osgeo.org/grass-stable/manuals/libpython/script_intro.html) and [PyGRASS](https://grass.osgeo.org/grass-stable/manuals/libpython/pygrass_index.html).\n",
    "PyGRASS is advantageous for more advanced workflows. Here we will be using Python Scripting Library (`import grass.script as gs`)\n",
    "as it is simple and straightforward to use.\n",
    " \n",
    "\n",
    "GRASS GIS Python Scripting Library provides functions to call GRASS modules within scripts as subprocesses. The most often used functions include:\n",
    "\n",
    " * [run_command()](https://grass.osgeo.org/grass-stable/manuals/libpython/script.html#script.core.run_command): used with modules which output raster/vector data where text output is not expected\n",
    " * [read_command()](https://grass.osgeo.org/grass-stable/manuals/libpython/script.html#script.core.read_command): used when we are interested in text output\n",
    " * [parse_command()](https://grass.osgeo.org/grass-stable/manuals/libpython/script.html#script.core.parse_command): used with modules producing text output as key=value pair\n",
    " * [write_command()](https://grass.osgeo.org/grass-stable/manuals/libpython/script.html#script.core.write_command): for modules expecting text input from either standard input or file\n",
    "\n",
    "\n",
    "It also provides several wrapper functions for often called modules. The list of convenient wrapper functions with examples includes:\n",
    "\n",
    " * Raster metadata using [raster_info()](https://grass.osgeo.org/grass-stable/manuals/libpython/script.html#script.raster.raster_info): `gs.raster_info('dsm')`\n",
    " * Vector metadata using [vector_info()](https://grass.osgeo.org/grass-stable/manuals/libpython/script.html#script.vector.vector_info): `gs.vector_info('roads')`\n",
    " * List raster data in current location using [list_grouped()](https://grass.osgeo.org/grass-stable/manuals/libpython/script.html#script.core.list_grouped): `gs.list_grouped(type=['raster'])`\n",
    " * Get current computational region using [region()](https://grass.osgeo.org/grass-stable/manuals/libpython/script.html#script.core.region): `gs.region()`\n",
    " * Run raster algebra using [mapcalc()](https://grass.osgeo.org/grass-stable/manuals/libpython/script.html#script.raster.mapcalc): `gs.mapcalc()`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use Python to compute viewshed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1. set computational region based on DSM\n",
    "gs.run_command(\"g.region\", raster=\"dsm\")\n",
    "# 2. Compute viewshed, flag 'b' is for binary (0 and 1) output \n",
    "gs.run_command(\"r.viewshed\", input=\"dsm\", output=\"viewshed\", flags=\"b\", coordinates=(711260, 3960860))\n",
    "# 3. Compute basic univariate statistics, flag 'g' is for parsable output\n",
    "univar = gs.parse_command(\"r.univar\", map=\"viewshed\", flags='g')\n",
    "# 4. Get current region settings to get cell size\n",
    "region = gs.region()\n",
    "cell_size = region[\"nsres\"] * region[\"ewres\"]\n",
    "percentage = 100 * float(univar['sum']) / float(univar['n'])\n",
    "area = cell_size * float(univar['sum'])\n",
    "print(f\"Percentage of visible area is {percentage:.2f}%, which is {area / 10000:.2f} ha\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewshed_map = gj.Map()\n",
    "viewshed_map.d_rast(map=\"ortho\")\n",
    "# select only cells with value 1 to visualize\n",
    "viewshed_map.d_rast(map=\"viewshed\", values=1)\n",
    "viewshed_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2: Viewshed case study\n",
    "\n",
    "In the second part, we will demonstrate the use of GRASS for a small viewshed case study.\n",
    "The goal is to **compute and analyze the area a driver would see from a road**.\n",
    "This notebook can be run only after notebook Part 1 was executed.\n",
    "\n",
    "Topics covered:\n",
    " * Python scripting\n",
    " * manipulating vector data ([v.build.polylines](https://grass.osgeo.org/grass-stable/manuals/v.build.polylines.html), [v.to.points](https://grass.osgeo.org/grass-stable/manuals/v.to.points.html))\n",
    " * vector attributes ([v.db.select](https://grass.osgeo.org/grass-stable/manuals/v.db.select.html))\n",
    " * viewshed computation ([r.viewshed](https://grass.osgeo.org/grass-stable/manuals/r.viewshed.html))\n",
    " * simple parallelization ([multiprocessing.Pool](https://docs.python.org/3/library/multiprocessing.html))\n",
    " * region handling ([grass.script.region_env](https://grass.osgeo.org/grass-stable/manuals/libpython/script.html#script.core.region_env))\n",
    " * raster algebra ([r.mapcalc](https://grass.osgeo.org/grass-stable/manuals/r.mapcalc.html))\n",
    " * temporal data handling\n",
    " * reprojecting ([r.proj](https://grass.osgeo.org/grass-stable/manuals/r.proj.html))\n",
    " * raster mask ([r.mask](https://grass.osgeo.org/grass-stable/manuals/r.mask.html))\n",
    " * raster as numpy array ([grass.script.array](https://grass.osgeo.org/grass-stable/manuals/libpython/script.html#module-script.array))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous notebook (Part 1) we created new location (project) *dix_park*. This automatically created new mapset (subproject) _PERMANENT_ where we then imported our base data. Now it's time to create a new mapset for our viewshed analysis, we will name it _viewshed_:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "grass -c -e ~/grassdata/dix_park/viewshed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Python standard library and IPython packages we need.\n",
    "import subprocess\n",
    "import sys\n",
    "from tqdm import tqdm\n",
    "\n",
    "# Ask GRASS GIS where its Python packages are.\n",
    "sys.path.append(\n",
    "    subprocess.check_output([\"grass\", \"--config\", \"python_path\"], text=True).strip()\n",
    ")\n",
    "\n",
    "# Import the GRASS GIS packages we need.\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "# Start GRASS Session\n",
    "gj.init(\"~/grassdata\", \"dix_park\", \"viewshed\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data preparation\n",
    "We will first derive viewpoints along the road *Umstead Drive* (vector `umstead_drive_segments`) that we extracted in the first part of the workshop.\n",
    "\n",
    "1. Because the road consists of several segments, we will merge them into one.\n",
    "2. Create new vector of points along the line with distance 50 m."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"v.build.polylines\", input=\"umstead_drive_segments\", output=\"umstead_drive\", cats=\"first\")\n",
    "gs.run_command(\"v.to.points\", input=\"umstead_drive\", type=\"line\", output=\"viewpoints\", dmax=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize the points with InteractiveMap with OSM tiles (see [other tile options](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.jupyter.html#module-grass.jupyter.interactivemap)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "road_map = gj.InteractiveMap(tiles=\"OpenStreetMap\")\n",
    "road_map.add_vector(\"umstead_drive\")\n",
    "road_map.add_vector(\"viewpoints\")\n",
    "road_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next part of analysis is raster-based, so we need to make sure we set computational region as we need. Specifically, we set it to match the DSM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", raster=\"dsm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to compute the visibility using DSM, however some points may fall on top of a tree, so we need to filter those out.\n",
    "\n",
    "First compute height above ground (DSM - DTM):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.mapcalc(\"diff = dsm - ground\")\n",
    "gs.run_command(\"r.colors\", map=\"diff\", color=\"differences\")\n",
    "\n",
    "diff_map = gj.Map()\n",
    "diff_map.d_rast(map=\"diff\")\n",
    "diff_map.d_vect(map=\"umstead_drive\")\n",
    "diff_map.d_legend(raster=\"diff\")\n",
    "diff_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract height above ground for the viewpoint locations to identify points that fall on top of a tree growing next to the road:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"v.what.rast\", map=\"viewpoints\", layer=2, raster=\"diff\", column=\"height\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See the newly computed attribute data. This example shows how the attribute data can be loaded into pandas:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "import pandas as pd\n",
    "pd.DataFrame(json.loads(gs.read_command(\"v.db.select\", map=\"viewpoints\", columns=\"cat,height\", layer=2, format=\"json\"))[\"records\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize the viewpoints with the height-above-ground raster. You can filter the points based on the height above ground, we won't display points with height > 2.\n",
    "Additionally, we will render the result larger (`width=1000`) and we will render the map zoomed in to the area with the points\n",
    "by saving a region and using it in Map (`saved_region=\"umstead_drive_region\"`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", vector=\"umstead_drive\", align=\"dsm\", save=\"umstead_drive_region\")\n",
    "\n",
    "img = gj.Map(width=1000, saved_region=\"umstead_drive_region\")\n",
    "img.d_rast(map=\"diff\")\n",
    "img.d_vect(map=\"umstead_drive\")\n",
    "img.d_vect(map=\"viewpoints\", layer=2, where=\"height >= 2\", size=15, icon=\"basic/pin\", fill_color=\"red\")\n",
    "img.d_vect(map=\"viewpoints\", layer=2, where=\"height < 2\", size=15, icon=\"basic/pin\")\n",
    "img.d_legend(raster=\"diff\")\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Viewshed computation\n",
    "To get the cumulative viewshed, we will compute viewsheds from all the viewpoints we generated earlier.\n",
    "First, we get the list coordinates of the viewpoints that are likely lying on the ground:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewpoints = gs.read_command('v.out.ascii', input='viewpoints',\n",
    "                             separator='comma', layer=2, where=\"height < 2\").strip().splitlines()\n",
    "viewpoints = [p.split(\",\") for p in viewpoints]\n",
    "viewpoints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now compute the viewshed from each viewpoint in a loop. We set max distance of 300 m. Each viewshed will be named `viewshed_cat`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "maps = []\n",
    "for x, y, cat in tqdm(viewpoints):\n",
    "    name = f\"viewshed_{cat}\"\n",
    "    gs.run_command(\"r.viewshed\", input=\"dsm\", output=name,\n",
    "                   coordinates=(x, y), max_distance=300, flags=\"b\")\n",
    "    maps.append(name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since these are independent runs, we can easily parallelize the r.viewshed calls using Python multiprocessing.\n",
    "We define a function that computes the viewshed and returns the name of the output or None in case of error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "from grass.exceptions import CalledModuleError\n",
    "from multiprocessing import Pool, cpu_count\n",
    "\n",
    "\n",
    "def viewshed(point):\n",
    "    x, y, cat = point\n",
    "    x, y = float(x), float(y)\n",
    "    name = f\"viewshed_{cat}\"\n",
    "    try:\n",
    "        gs.run_command(\"r.viewshed\", input=\"dsm\", output=name,\n",
    "                       coordinates=(x, y), max_distance=300, flags=\"b\")\n",
    "        return f\"viewshed_{cat}\"\n",
    "    except CalledModuleError:\n",
    "        return None\n",
    "\n",
    "# run with the number of CPUs available\n",
    "# proc = cpu_count()\n",
    "proc = 1\n",
    "with Pool(processes=proc) as pool:\n",
    "    maps = pool.map_async(viewshed, viewpoints).get()\n",
    "print(maps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One trick to speedup viewshed computation is to limit the computation only\n",
    "to the actual area given by the maxdistance option. To do that we will locally modify the computational region\n",
    "and pass the environment to the module directly. The current computational region won't be affected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "import os\n",
    "from grass.exceptions import CalledModuleError\n",
    "from multiprocessing import Pool, cpu_count\n",
    "\n",
    "\n",
    "def viewshed(point):\n",
    "    x, y, cat = point\n",
    "    x, y = float(x), float(y)\n",
    "    max_distance = 300\n",
    "    # copy current environment\n",
    "    env = os.environ.copy()\n",
    "    # set GRASS_REGION variable using region_env function\n",
    "    env[\"GRASS_REGION\"] = gs.region_env(align=\"dsm\",\n",
    "                                        e=x + max_distance,\n",
    "                                        w=x - max_distance,\n",
    "                                        n=y + max_distance,\n",
    "                                        s=y - max_distance)\n",
    "    name = f\"viewshed_{cat}\"\n",
    "    try:\n",
    "        gs.run_command(\"r.viewshed\", input=\"dsm\", output=name, flags=\"b\",\n",
    "                      coordinates=(x, y), max_distance=max_distance, env=env)\n",
    "        return f\"viewshed_{cat}\"\n",
    "    except CalledModuleError:\n",
    "        return None\n",
    "\n",
    "# run with the number of CPUs available\n",
    "# proc = cpu_count()\n",
    "proc = 1\n",
    "with Pool(processes=proc) as pool:\n",
    "    maps = pool.map_async(viewshed, viewpoints).get()\n",
    "print(maps)\n",
    "print(f\"Viewshed num cells: {gs.raster_info(maps[0])['cells']}\")\n",
    "print(f\"DSM num cells: {gs.raster_info('dsm')['cells']}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Temporal dataset of viewsheds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this part we will create, analyze and visualize a temporal dataset of viewsheds using [temporal tools](https://grass.osgeo.org/grass-stable/manuals/temporal.html). \n",
    "\n",
    "First, let's check we have the viewshed rasters ready:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.list_strings(type=\"raster\", pattern=\"viewshed_*\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will create an empty space-time raster dataset called _viewsheds_ with relative temporal type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"t.create\", output=\"viewsheds\", type=\"strds\", temporaltype=\"relative\",\n",
    "              title=\"Viewshed series\", description=\"Series of viewsheds along a road\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we register the viewshed rasters with start time 1 and 1-minute increment to simulate a change of view of a car driving slowly along the road:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"t.register\", input=\"viewsheds\", maps=\",\".join(maps), start=1, unit=\"minutes\", increment=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's print basic dataset info. We will use this info later on to set computational region covering the entire dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pprint import pprint\n",
    "info = gs.parse_command(\"t.info\", input=\"viewsheds\", flags=\"g\")\n",
    "pprint(info)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To list the individual rasters, we will use t.rast.list. Notice there is no end time, because we don't use interval data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from io import StringIO\n",
    "\n",
    "pd.read_csv(StringIO(gs.read_command(\"t.rast.list\", input=\"viewsheds\", separator=\"comma\", columns=\"name,start_time,end_time\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can quickly get basic statistics such as the size of the viewsheds (see _sum_ column for the number of visible cells):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(StringIO(gs.read_command(\"t.rast.univar\", input=\"viewsheds\", separator=\"comma\")))\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's find and visualize largest and smallest viewshed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "largest = df.iloc[df[['sum']].idxmax()[\"sum\"]].id\n",
    "smallest = df.iloc[df[['sum']].idxmin()[\"sum\"]].id\n",
    "\n",
    "gs.run_command(\"g.region\", raster=[largest, smallest], save=\"zoom_region\")\n",
    "viewshed_map = gj.Map(saved_region=\"zoom_region\")\n",
    "viewshed_map.d_rast(map=\"ortho\")\n",
    "viewshed_map.d_rast(map=largest, values=1)\n",
    "viewshed_map.d_rast(map=smallest, values=1)\n",
    "viewshed_map.d_vect(map=\"umstead_drive\", color=\"white\")\n",
    "viewshed_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute a temporal dataset where values of each viewshed will represent the registered start time.\n",
    "\n",
    "Before we do that, let's set computational region to match the bounding box of the entire dataset, here we use previous t.info output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", n=info[\"north\"], s=info[\"south\"], e=info[\"east\"], w=info[\"west\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use temporal raster algebra. Here we compute a new temporal dataset _viewsheds_start_ so that for example viewshed with start time 5 has value 5 for visible area and no data for invisible area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"t.rast.mapcalc\", inputs=\"viewsheds\", output=\"viewsheds_start\", basename=\"viewshed_start\",\n",
    "               expression=\"if (viewsheds == 0, null(), start_time())\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set color of the newly computed time series:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"t.rast.colors\", input=\"viewsheds_start\", color=\"plasma\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With TimeSeriesMap, we can interactively visualize the time series:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "timemap = gj.TimeSeriesMap(width=800)\n",
    "timemap.d_rast(map=\"ortho\")\n",
    "timemap.d_vect(map=\"umstead_drive\", )\n",
    "timemap.add_raster_series(\"viewsheds_start\")\n",
    "timemap.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can export an animated GIF:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "\n",
    "Image(timemap.save(\"animation.gif\", duration=300))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cumulative viewshed\n",
    "We can compute the cumulative viewshed, which aggregates viewsheds from multiple viewpoints. In this way you can e.g., identify the most frequently visible areas from the road.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since our viewshed rasters are binary (0 invisible, 1 visible), we will use r.series method *sum*. Then we replace zeros with no data using r.null and set a new color ramp:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cumulative viewshed\n",
    "gs.run_command(\"r.series\", input=maps, output=\"cumulative_viewshed\", method=\"sum\")\n",
    "gs.run_command(\"r.null\", map=\"cumulative_viewshed\", setnull=0)\n",
    "gs.run_command(\"r.colors\", map=\"cumulative_viewshed\", color=\"plasma\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cumulative_map = gj.InteractiveMap()\n",
    "cumulative_map.add_raster(\"cumulative_viewshed\", opacity=0.8)\n",
    "cumulative_map.add_vector(\"umstead_drive\")\n",
    "cumulative_map.add_layer_control(position=\"bottomright\")\n",
    "cumulative_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And create a 3D rendering with draped cumulative viewshed over the DSM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map3d = gj.Map3D()\n",
    "map3d.render(elevation_map=\"dsm\", resolution_fine=1, color_map=\"cumulative_viewshed\",\n",
    "           vline=\"umstead_drive\", vline_width=3, vline_color=\"white\", light_brightness=50,\n",
    "           position=[0.4, 0.8], height=3000, perspective=10)\n",
    "map3d.overlay.d_legend(raster=\"cumulative_viewshed\", at=(0, 30, 1, 7), use=[1, 2, 3, 4, 5, 6], flags=\"fb\")\n",
    "map3d.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data reprojection\n",
    "Next, we will analyze the cumulative viewshed to see how much greenery a driver would see on the way. To do that we compute NDVI:\n",
    "\n",
    "We reproject R and NIR Landsat bands from NCSPM sample dataset we already have available. Tool r.proj respects the current region (extent and resolution), but you can set resolution to certain value, we use 28.5 m which is the original resolution.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", raster=\"dsm\")\n",
    "for band in [30, 40]:\n",
    "    gs.run_command(\"r.proj\", location=\"nc_spm_08_grass7\", mapset=\"PERMANENT\", input=f\"lsat7_2002_{band}\", method=\"nearest\", resolution=28.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute NDVI:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"i.vi\", viname=\"ndvi\", red=\"lsat7_2002_30\", nir=\"lsat7_2002_40\", output=\"ndvi\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi_map = gj.Map()\n",
    "ndvi_map.d_rast(map=\"ndvi\")\n",
    "ndvi_map.d_legend(raster=\"ndvi\")\n",
    "ndvi_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mask\n",
    "Now let's analyze what is the distribution of NDVI within the visible area. We will mask the data by the visible area:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.mask\", raster=\"cumulative_viewshed\")\n",
    "data = gs.parse_command(\"r.univar\", map=\"ndvi\", flags=\"g\")\n",
    "print(f\"Average NDVI of visible cells: {float(data['mean']):.2f} ± {float(data['stddev']):.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi_map = gj.Map()\n",
    "ndvi_map.d_rast(map=\"ndvi\")\n",
    "ndvi_map.d_legend(raster=\"ndvi\", flags=\"d\")\n",
    "ndvi_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see the histogram of visible NDVI using d.histogram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "histo = gj.Map(width=800, height=400)\n",
    "histo.d_histogram(map=\"ndvi\", bgcolor=\"grey\")\n",
    "histo.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read as numpy array\n",
    "It is also easy to use the results as a numpy array and then use other Python libraries to analyze the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from grass.script import array as garray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi = garray.array(mapname=\"ndvi\", null='nan')\n",
    "ndvi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('darkgrid')\n",
    "sns.histplot(ndvi.ravel(), kde=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, remove the mask:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.mask\", flags=\"r\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 3: Estimating inundation extent using HAND methodology\n",
    "\n",
    "\n",
    "In this example we will use some of GRASS GIS hydrology tools, namely:\n",
    "\n",
    "* [r.watershed](https://grass.osgeo.org/grass-stable/manuals/r.watershed.html): for computing flow accumulation, drainage direction, the location of streams and watershed basins; it does not need sink filling because of using the least-cost-path to route flow out of sinks\n",
    "* [r.lake](https://grass.osgeo.org/grass-stable/manuals/r.lake.html): to fill a lake to a target water level from a given start point or seed raster\n",
    "* [r.lake.series](https://grass.osgeo.org/grass-stable/manuals/r.lake.series.html): addon which runs r.lake for different water levels\n",
    "* [r.stream.distance](https://grass.osgeo.org/grass-stable/manuals/r.stream.distance.html): for computing the distance to streams or outlet, the relative elevation above streams; the distance and the elevation are calculated along watercourses\n",
    "\n",
    "First, let's create a new mapset *flooding* in nc_spm_08_grass7 sample dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "grass -c -e ~/grassdata/nc_spm_08_grass7/flooding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize GRASS session:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Python standard library and IPython packages we need.\n",
    "import subprocess\n",
    "import sys\n",
    "\n",
    "# Ask GRASS GIS where its Python packages are.\n",
    "sys.path.append(\n",
    "    subprocess.check_output([\"grass\", \"--config\", \"python_path\"], text=True).strip()\n",
    ")\n",
    "\n",
    "# Import the GRASS GIS packages we need.\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "# Start GRASS Session\n",
    "gj.init(\"~/grassdata\", \"nc_spm_08_grass7\", \"flooding\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Modules r.stream.distance and r.lake.series are addons and we need to install them first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.extension\", extension=\"r.stream.distance\")\n",
    "gs.run_command(\"g.extension\", extension=\"r.lake.series\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute HAND raster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will estimate inundation extent using the Height Above Nearest Drainage methodology ([A.D. Nobre, 2011](https://doi.org/10.1016/j.jhydrol.2011.03.051)). We will compute the HAND terrain model representing the differences in elevation between each grid cell and the elevations of the flowpath-connected downslope grid cells where the flow enters the channel.\n",
    "\n",
    "First we compute the flow accumulation, drainage and streams (with a threshold value of 100000). We convert the streams to vector for better visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", raster=\"elevation\")\n",
    "gs.run_command(\"r.watershed\", elevation=\"elevation\", accumulation=\"flowacc\",\n",
    "               drainage=\"drainage\", stream=\"rwatershed_streams\", threshold=100000)\n",
    "gs.run_command(\"r.to.vect\", input=\"rwatershed_streams\", output=\"rwatershed_streams\", type=\"line\")\n",
    "\n",
    "fllowacc_map = gj.Map()\n",
    "fllowacc_map.d_rast(map=\"flowacc\")\n",
    "fllowacc_map.d_vect(map=\"rwatershed_streams\", width=2, color=\"blue\")\n",
    "fllowacc_map.d_legend(raster=\"flowacc\", range=\"0,1000\")\n",
    "fllowacc_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's zoom in to see the flow accumulation raster better:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flowacc_map = gj.InteractiveMap()\n",
    "flowacc_map.add_raster(\"flowacc\")\n",
    "flowacc_map.add_vector(\"rwatershed_streams\")\n",
    "flowacc_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we use r.stream.distance with output parameter difference to compute a new raster map where each cell is the elevation difference between the cell and the the cell on the stream where the cell drains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.stream.distance\", stream_rast=\"rwatershed_streams\", direction=\"drainage\",\n",
    "               elevation=\"elevation\", method=\"downstream\", difference=\"above_stream\")\n",
    "gs.run_command(\"r.colors\", map=\"above_stream\", color=\"elevation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare original elevation and HAND raster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map3d = gj.Map3D()\n",
    "map3d.render(elevation_map=\"elevation\", color_map=\"elevation\",\n",
    "             position=(0.5, 1), height=3000, perspective=12, zexag=5)\n",
    "map3d.overlay.d_legend(raster=\"elevation\", at=(5, 80, 87, 92), flags=\"b\", border_color=\"none\")\n",
    "map3d.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map3d = gj.Map3D()\n",
    "map3d.render(elevation_map=\"above_stream\", color_map=\"above_stream\",\n",
    "             position=(0.5, 1), height=3000, perspective=12, zexag=5)\n",
    "map3d.overlay.d_legend(raster=\"above_stream\", at=(5, 80, 90, 95), flags=\"b\", border_color=\"none\")\n",
    "map3d.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inundation\n",
    "Before we compute the inundation, we will look at how r.lake works. We compute a lake from a specified coordinate pair and water level:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.lake\", elevation=\"elevation\", water_level=90, lake=\"lake\", coordinates=[637877, 218475])\n",
    "\n",
    "lake_map = gj.Map()\n",
    "lake_map.d_rast(map=\"elevation\")\n",
    "lake_map.d_rast(map=\"lake\")\n",
    "lake_map.d_legend(raster=\"lake\", label_values=\"0.1,5,10,15\", digits=0, at=(1, 40, 90, 95))\n",
    "lake_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now instead of the elevation raster we use the HAND raster to simulate 5-meter inundation and, as the seed we specify the entire stream."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.lake\", elevation=\"above_stream\", water_level=5, lake=\"flood\", seed=\"rwatershed_streams\")\n",
    "\n",
    "hand_map = gj.Map()\n",
    "hand_map.d_rast(map=\"above_stream\")\n",
    "hand_map.d_rast(map=\"flood\")\n",
    "hand_map.d_legend(raster=\"flood\", label_values=\"0.1,4,8\", digits=0, at=(1, 40, 90, 95))\n",
    "hand_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the r.lake.series addon we can create a series of inundation maps with rising water levels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from io import StringIO\n",
    "import pandas as pd\n",
    "\n",
    "gs.run_command(\"r.lake.series\", elevation=\"above_stream\", start_water_level=0, end_water_level=5, \n",
    "               water_level_step=0.5, output=\"inundation\", seed_raster=\"rwatershed_streams\", quiet=True)\n",
    "gs.run_command(\"t.rast.colors\", input=\"inundation\", color=\"water\")\n",
    "pd.read_csv(StringIO(gs.read_command(\"t.rast.list\", input=\"inundation\", separator=\"comma\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "r.lake.series creates a space-time dataset. We can use the [temporal modules](https://grass.osgeo.org/grass-stable/manuals/temporal.html) to further work with the data. For example, we could further compute the volume and extent of flood water using t.rast.univar:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.read_csv(StringIO(gs.read_command(\"t.rast.univar\", input=\"inundation\", separator=\"comma\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timemap = gj.TimeSeriesMap()\n",
    "timemap.d_rast(map=\"elevation_shade\")\n",
    "timemap.d_vect(map=\"streets_wake\", color=\"#4D4D4D\")\n",
    "timemap.add_raster_series(\"inundation\")\n",
    "timemap.d_legend(at=(1, 40, 90, 95))\n",
    "timemap.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 4: Image segmentation\n",
    "\n",
    "\n",
    "In this example we will show segmentation of a Landsat scene.\n",
    "We show two segmentation modules: [i.segment](https://grass.osgeo.org/grass-stable/manuals/i.segment.html) and the addon [i.superpixels.slic](https://grass.osgeo.org/grass-stable/manuals/addons/i.superpixels.slic.html).\n",
    "Note that each segmentation algorithm is designed for different purpose, so we can't directly compare them.\n",
    "\n",
    "First, let's create a new mapset *segmentation* in nc_spm_08_grass7 location:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "grass -c -e ~/grassdata/nc_spm_08_grass7/segmentation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize GRASS session:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Python standard library and IPython packages we need.\n",
    "import subprocess\n",
    "import sys\n",
    "\n",
    "# Ask GRASS GIS where its Python packages are.\n",
    "sys.path.append(\n",
    "    subprocess.check_output([\"grass\", \"--config\", \"python_path\"], text=True).strip()\n",
    ")\n",
    "\n",
    "# Import the GRASS GIS packages we need.\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "# Start GRASS Session\n",
    "gj.init(\"~/grassdata\", \"nc_spm_08_grass7\", \"segmentation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Install the addon:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.extension\", extension=\"i.superpixels.slic\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data preparation\n",
    "Imagery modules typically work with *imagery groups*. We first list the landsat raster data and then create an imagery group:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maps = gs.list_grouped(type=\"raster\", pattern=\"lsat*\")[\"PERMANENT\"]\n",
    "print(maps)\n",
    "gs.run_command(\"i.group\", group=\"landsat\", subgroup=\"landsat\", input=maps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we derive NDVI to see some of the effects of segmentation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", raster=\"lsat7_2002_30\")\n",
    "gs.run_command(\"i.vi\", red=\"lsat7_2002_30\",  nir=\"lsat7_2002_40\", output=\"ndvi\", viname=\"ndvi\")\n",
    "\n",
    "ndvi_map = gj.Map()\n",
    "ndvi_map.d_rast(map=\"ndvi\")\n",
    "ndvi_map.d_legend(raster=\"ndvi\", at=(2, 50, 1, 5))\n",
    "ndvi_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Segmentation with i.superpixels.slic\n",
    "Superpixels can be defined as a group of pixels that share common characteristics and are useful in computer vision and image processing.\n",
    "Here we run [i.superpixels.slic](https://grass.osgeo.org/grass-stable/manuals/addons/i.superpixels.slic.html) and convert the resulting raster to vector for better visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"i.superpixels.slic\", input=\"landsat\", output=\"superpixels\", num_pixels=1000, compactness=0.5)\n",
    "gs.run_command(\"r.to.vect\", input=\"superpixels\", output=\"superpixels\", type=\"area\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can play with *compactness* and *num_pixels* parameters and see how the resulting segmentation changes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "superpixels_map = gj.Map()\n",
    "superpixels_map.d_rast(map=\"ndvi\")\n",
    "superpixels_map.d_vect(map=\"superpixels\", width=1, color=\"black\", fill_color=\"none\")\n",
    "superpixels_map.d_legend(raster=\"ndvi\", at=(2, 50, 1, 5))\n",
    "superpixels_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For fun, let's do zonal statistics on the results. We compute the median NDVI value within each segment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.stats.quantile\", base=\"superpixels\", cover=\"ndvi\", output=\"superpixels_ndvi\")\n",
    "\n",
    "superpixels_map = gj.Map()\n",
    "superpixels_map.d_rast(map=\"superpixels_ndvi\")\n",
    "superpixels_map.d_legend(raster=\"superpixels_ndvi\", at=(2, 50, 1, 5))\n",
    "superpixels_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Segmentation with i.segment\n",
    "\n",
    "Next, we do the same, but with i.segment to see the different behavior. Note that i.segment uses *region growing* algorithm by default, but *mean shift* is also available:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"i.segment\", group=\"landsat\", output=\"segments\", threshold=0.5, minsize=50)\n",
    "gs.run_command(\"r.to.vect\", input=\"segments\", output=\"segments\", type=\"area\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "segments_map = gj.Map()\n",
    "segments_map.d_rast(map=\"ndvi\")\n",
    "segments_map.d_vect(map=\"segments\", width=1, color=\"black\", fill_color=\"none\")\n",
    "segments_map.d_legend(raster=\"ndvi\", at=(2, 50, 1, 5))\n",
    "segments_map.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.stats.quantile\", base=\"segments\", cover=\"ndvi\", output=\"segments_ndvi\")\n",
    "\n",
    "segments_map = gj.Map()\n",
    "segments_map.d_rast(map=\"segments_ndvi\")\n",
    "segments_map.d_legend(raster=\"superpixels_ndvi\", at=(2, 50, 1, 5))\n",
    "segments_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 5: Parallelization examples\n",
    "\n",
    "This part will briefly cover how to run GRASS computations in parallel. First, create a new mapset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "grass -c -e ~/grassdata/nc_spm_08_grass7/parallelization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import subprocess\n",
    "import sys\n",
    "\n",
    "# Ask GRASS GIS where its Python packages are.\n",
    "sys.path.append(\n",
    "    subprocess.check_output([\"grass\", \"--config\", \"python_path\"], text=True).strip()\n",
    ")\n",
    "\n",
    "# Import the GRASS GIS packages we need.\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "# Start GRASS Session\n",
    "gj.init(\"~/grassdata\", \"nc_spm_08_grass7\", \"parallelization\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tool-level parallelization\n",
    "There are several [internally parallelized tools](https://grass.osgeo.org/grass-stable/manuals/keywords.html#parallel), either using OpenMP or Python multiprocessing library. We can use `nprocs` option to set the number of cores to be used for processing.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set computational region to match _elevation_ raster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", raster=\"elevation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute moving window analysis and measure time first with one core, than with 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "gs.run_command(\"r.neighbors\", input=\"elevation\", output=\"elevation_smoothed\", method=\"average\", size=25, nprocs=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "gs.run_command(\"r.neighbors\", input=\"elevation\", output=\"elevation_smoothed\", method=\"average\", size=25, nprocs=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize original and smoothed raster (turn layers on and off):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neighbors_map = gj.InteractiveMap()\n",
    "neighbors_map.add_raster(\"elevation\")\n",
    "neighbors_map.add_raster(\"elevation_smoothed\")\n",
    "neighbors_map.add_layer_control(position=\"bottomright\")\n",
    "neighbors_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GridModule for tiling\n",
    "Some compute-intensive tasks can benefit from spatially splitting the task into tiles, and then running the task in parallel. [GridModule](https://grass.osgeo.org/grass-stable/manuals/libpython/pygrass.modules.grid.html) can automate this splitting-computing-merging procedure and execute the computation in parallel.\n",
    "\n",
    "In this example, we will interpolate an elevation surface from vector points using IDW interpolation. First, set computational region to match the extent of vector points and set the resolution to 1 meter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", vector=\"elev_lid792_bepts\", res=1, flags=\"a\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Measure the time without using GridModule:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "gs.run_command(\"v.surf.idw\", input=\"elev_lid792_bepts\", output=\"elev_lid792_interp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now with GridModule:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile interpolation.py\n",
    "from grass.pygrass.modules.grid import GridModule\n",
    "\n",
    "\n",
    "grid = GridModule(\n",
    "    \"v.surf.idw\",\n",
    "    input=\"elev_lid792_bepts\",\n",
    "    output=\"elev_lid792_interp\",\n",
    "    processes=3,\n",
    "    overlap=12,\n",
    "    quiet=True,\n",
    ")\n",
    "grid.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "!python interpolation.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.colors\", map=\"elev_lid792_interp\", color=\"elevation\")\n",
    "neighbors_map = gj.Map()\n",
    "neighbors_map.d_rast(map=\"elev_lid792_interp\")\n",
    "neighbors_map.d_vect(map=\"elev_lid792_bepts\", size=1, color=\"black\")\n",
    "neighbors_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running multiple independent computations\n",
    "In this example, our goal is to compute multiple viewsheds and export them as picture. Since these are independent computations, we can run them in parallel.\n",
    "The first part implements this task in Python using _multiprocessing_ library (also explained in Part 2)\n",
    "and the second part will run each computation using `grass --exec` interface in separate mapsets that allows us to potentially distribute the computation across multiple nodes on an HPC."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First compute a shaded relief raster for visualization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", raster=\"elevation\")\n",
    "gs.run_command(\"r.relief\", input=\"elevation\", output=\"relief\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will compute viewsheds from vector points _firestations_:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewpoints = gs.read_command('v.out.ascii', input='firestations', separator='comma', flags=\"r\").strip().splitlines()\n",
    "viewpoints = [p.split(\",\") for p in viewpoints]\n",
    "viewpoints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will extend script from Part 2 of this workshop to include rendering to file. Set `nprocs` to more than 1 when possible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from grass.exceptions import CalledModuleError\n",
    "from multiprocessing import Pool, cpu_count\n",
    "\n",
    "\n",
    "def viewshed(point):\n",
    "    x, y, cat = point\n",
    "    x, y = float(x), float(y)\n",
    "    max_distance = 2000\n",
    "    # copy current environment\n",
    "    env = os.environ.copy()\n",
    "    # set GRASS_REGION variable using region_env function\n",
    "    env[\"GRASS_REGION\"] = gs.region_env(align=\"elevation\",\n",
    "                                        e=x + max_distance,\n",
    "                                        w=x - max_distance,\n",
    "                                        n=y + max_distance,\n",
    "                                        s=y - max_distance)\n",
    "    name = f\"viewshed_{cat}\"\n",
    "    try:\n",
    "        gs.run_command(\"r.viewshed\", input=\"elevation\", output=name, flags=\"b\",\n",
    "                      coordinates=(x, y), max_distance=max_distance, env=env)\n",
    "        # create visualization\n",
    "        viewshed_map = gj.Map(use_region=True, env=env)\n",
    "        viewshed_map.d_rast(map=\"relief\")\n",
    "        viewshed_map.d_rast(map=f\"viewshed_{cat}\", values=1)\n",
    "        viewshed_map.d_vect(map=\"firestations\", cat=cat, size=15, icon=\"basic/pin\")\n",
    "        viewshed_map.save(f\"viewshed_{cat}.png\")\n",
    "        return f\"viewshed_{cat}\"\n",
    "    except CalledModuleError:\n",
    "        return None\n",
    "\n",
    "# run with the number of CPUs available\n",
    "# proc = cpu_count()\n",
    "nprocs = 1\n",
    "with Pool(processes=nprocs) as pool:\n",
    "    maps = pool.map_async(viewshed, viewpoints).get()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at one of the computed and rendered viewsheds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "\n",
    "Image(\"viewshed_22.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this way, we can't distribute the computation across multiple nodes (hundreds of cores).\n",
    "We will do the same thing differently, using `grass --exec` [interface](https://grass.osgeo.org/grass-stable/manuals/grass.html), running each task in a separate mapset. This way, the tasks could be distributed across multiple nodes.\n",
    "\n",
    "`--exec` interface allows GRASS tools and user scripts to be executed in a GRASS GIS non-interactive session. For example, here is a simple call to list all available vectors in PERMANENT mapset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "grass ~/grassdata/nc_spm_08_grass7/PERMANENT --exec g.list type=vector mapset=PERMANENT -t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will create a Python script `myscript.py` computing and rendering viewsheds similarly as in the previous example. The script requires 3 parameters (x and y coordinate, and category). Note that we can set computational region in a standard way, because each script will run in separate mapset, so the different regions won't interfere with each other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile myscript.py\n",
    "import sys\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "\n",
    "def main(x, y, cat):\n",
    "    max_distance = 2000\n",
    "    x, y = float(x), float(y)\n",
    "    name = f\"viewshed_{cat}\"\n",
    "    gs.run_command(\"g.region\", align=\"elevation\", e=x + max_distance,\n",
    "                   w=x - max_distance, n=y + max_distance, s=y - max_distance)\n",
    "    gs.run_command(\"r.viewshed\", input=\"elevation\", output=name, coordinates=(x, y),\n",
    "                   observer_elevation=3, max_distance=max_distance, flags=\"b\")\n",
    "    # create visualization\n",
    "    viewshed_map = gj.Map(use_region=True)\n",
    "    viewshed_map.d_rast(map=\"relief@parallelization\")\n",
    "    viewshed_map.d_rast(map=f\"viewshed_{cat}\", values=1)\n",
    "    viewshed_map.d_vect(map=\"firestations\", cat=cat, size=15, icon=\"basic/pin\")\n",
    "    viewshed_map.save(f\"viewshed_{cat}.png\")\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    args = sys.argv[1:]\n",
    "    main(*args)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will generate a file `jobs.sh` with one command per line. We will run each task in a temporary mapset so all computed data will be deleted afterwards. That is fine for our example where we need only the final PNG files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"jobs.sh\", \"w\") as f:\n",
    "    for viewpoint in viewpoints:\n",
    "        f.write(f\"grass --tmp-mapset ~/grassdata/nc_spm_08_grass7 --exec python myscript.py {viewpoint[0]} {viewpoint[1]} {viewpoint[2]}\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the content of the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat jobs.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To execute these commands in parallel, we can use e.g. [GNU Parallel](https://www.gnu.org/software/parallel/):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "parallel -j 2 < jobs.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check one of the resulting PNG files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image(\"viewshed_22.png\")"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "Untitled3.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}