workshop_part_1_intro.ipynb

{
 "cells": [
  {
   "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 project (location) 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 project and `e` will _exit_ the command after creating the project. 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 project. 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": [
    "<div class=\"alert alert-info\">\n",
    "If you run this workshop online, data is already prepared in the environment. Otherwise download the data <a href=\"http://fatra.cnr.ncsu.edu/foss4g2021/\">here</a>.\n",
    "</div>\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": [
    "Schema of GRASS project dix_park's content:\n",
    "\n",
    "<img src=\"img/data_structure1.png\" alt=\"GRASS project dix_park\" width=\"200\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data visualization\n",
    "\n",
    "We will call d.rast/d.vect modules using the [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. It will work only locally.\n",
    "\n",
    "```\n",
    "g.manual r.univar\n",
    "```"
   ]
  },
  {
   "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, that allows you to fully control the extent and resolution of your raster computations.\n",
    "All raster computations will be performed in the specified extent and with the given resolution to ensure consistency.\n",
    "Among other things, computational region 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 (subprojects)\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",
    "\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"
   ]
  },
  {
   "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 projection, 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 alignment the resolution to 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=umstead_drive_segments align=dsm save=road_region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Saved region can be used later on, for example to set the rendering extent:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "elevation_map = gj.Map(saved_region=\"road_region\")\n",
    "elevation_map.d_rast(map=\"dsm\")\n",
    "elevation_map.d_vect(map=\"umstead_drive_segments\")\n",
    "elevation_map.show()"
   ]
  },
  {
   "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 practice!\n",
    "Compute viewshed in Python.\n",
    "1. Set computational region to match raster `dsm`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Compute binary viewshed (visible - 1, not visible - 0) using [r.viewshed](https://grass.osgeo.org/grass-stable/manuals/r.viewshed.html) using `dsm` raster from X=711260 and Y=3960860:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Compute size of the visible area, using [r.univar](https://grass.osgeo.org/grass-stable/manuals/r.univar.html) computing univariate statistics. Use `gs.parse_command()` together with r.univar's `-g` flag to get a number of visible cells. Use `gs.region()` to get the cell size. Compute percentage of area that is visible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Display the viewshed overlaid over the orthophoto `ortho`. Use [d.rast](https://grass.osgeo.org/grass-stable/manuals/d.rast.html)'s `values` to display only visible cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "Find answers when you edit this cell.\n",
    "<!--\n",
    "# 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\")\n",
    "\n",
    "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()\n",
    "-->"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}