colab_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# install dependencies\n",
    "apt install -y autoconf2.13 autotools-dev bison flex g++ gettext git imagemagick libblas-dev libbz2-dev libcairo2-dev libfftw3-dev libfreetype6-dev libgdal-dev libgeos-dev libglu1-mesa-dev libjpeg-dev liblapack-dev liblas-c-dev libncurses5-dev libnetcdf-dev libpng-dev libpq-dev libproj-dev libreadline-dev libsqlite3-dev libtiff-dev libxmu-dev libzstd-dev make netcdf-bin p7zip proj-bin sqlite3 unixodbc-dev xvfb zlib1g-dev\n",
    "\n",
    "pip install Pillow\n",
    "pip install ply\n",
    "pip install PyVirtualDisplay\n",
    "\n",
    "git clone --depth 1 --branch Jupyter-Timeseries https://github.com/chaedri/grass.git grass-src\n",
    "\n",
    "# enter the directory with source code\n",
    "cd grass-src\n",
    "\n",
    "# compile\n",
    "./configure \\\n",
    "    --enable-largefile=yes \\\n",
    "    --with-nls \\\n",
    "    --with-cxx \\\n",
    "    --with-readline \\\n",
    "    --with-bzlib \\\n",
    "    --with-pthread \\\n",
    "    --with-proj-share=/usr/share/proj \\\n",
    "    --with-geos=/usr/bin/geos-config \\\n",
    "    --with-cairo \\\n",
    "    --with-opengl-libs=/usr/include/GL \\\n",
    "    --with-freetype=yes --with-freetype-includes=\"/usr/include/freetype2/\" \\\n",
    "    --with-sqlite=yes\n",
    "make -j2\n",
    "make -j2 install\n",
    "\n",
    "# leave the directory with source code\n",
    "cd ..\n",
    "\n",
    "# download sample data\n",
    "mkdir -p grassdata\n",
    "wget http://fatra.cnr.ncsu.edu/foss4g2021/dsm.tif\n",
    "wget http://fatra.cnr.ncsu.edu/foss4g2021/ground.tif\n",
    "wget http://fatra.cnr.ncsu.edu/foss4g2021/ortho.tif\n",
    "wget http://fatra.cnr.ncsu.edu/foss4g2021/ortho.tif.aux.xml\n",
    "wget http://fatra.cnr.ncsu.edu/foss4g2021/roads_backup.geojson -O roads.geojson"
   ]
  },
  {
   "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.\n",
    "We will first use the GUI and then we will replicate it in the notebook environment running on mybinder.com."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GUI intro"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GRASS GIS in a notebook"
   ]
  },
  {
   "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",
    "\n",
    "v = sys.version_info\n",
    "print(f\"We are using Python {v.major}.{v.minor}.{v.micro}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, we will also show how to use GRASS from Bash shell, for that we will use IPython magic:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "pwd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Start GRASS GIS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create new empty location (project) that uses projection [UTM zone 17 N](https://epsg.io/6346). Here, we place it directory called _grassdata_ under the home directory which already exists in this Binder environment (create it if you are running the notebook locally and change the path here accordingly)."
   ]
  },
  {
   "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 os\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/dix_park/PERMANENT\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import prepared DSM, DTM and ortho maps. If you run this notebook locally (not in Binder), download the data [here](http://fatra.cnr.ncsu.edu/foss4g2021/). In this case no reprojection is done since the data CRS matches the location's CRS.\n",
    "\n",
    "We also set the computational region based on the imported raster (more on that later)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.import\", input=\"dsm.tif\", output=\"dsm\", resample=\"bilinear\")\n",
    "gs.run_command(\"r.import\", input=\"ground.tif\", output=\"ground\")\n",
    "gs.run_command(\"r.import\", input=\"ortho.tif\", output=\"ortho\")\n",
    "\n",
    "# set the region to match ortho raster\n",
    "gs.run_command(\"g.region\", raster=\"ortho\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`roads.geojson` We extract roads for our study area using Download OSM data using Overpass Turbo.\n",
    "The data was <https://overpass-turbo.eu/s/1aGu> The query should run upon running the cell below, so\n",
    "it should be enough to *press Export - download as GeoJSON* and save it as `roads.geojson`.\n",
    "This will download the file on your local machine, so if you are running this in Jupyter Lab on Binder, upload it there with *Upload files* in the top left corner."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we import the downloaded file. This time the data will be automatically reprojected because the CRS doesn't match."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "v.import input=roads.geojson output=roads --quiet"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will call d.rast/d.vect modules using the new GRASS Jupyter API:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start a GrassRenderer\n",
    "img = gj.GrassRenderer()\n",
    "# Add a raster and vector to the map\n",
    "img.d_rast(map=\"ground\")\n",
    "img.d_vect(map=\"roads\")\n",
    "img.d_legend(raster=\"ground\", flags=\"b\")\n",
    "# Display map\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = gj.GrassRenderer()\n",
    "img.d_rast(map=\"ortho\")\n",
    "img.d_vect(map=\"roads\", width=\"2\", color=\"yellow\", where=\"name = 'Umstead Drive'\")\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is how we can visualize data with folium:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = gj.InteractiveMap(width=600)\n",
    "# Add vector and layer control to map\n",
    "fig.add_vector(\"roads\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also visualize data in 3D. Here we drape the ortho over the DSM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = gj.Grass3dRenderer()\n",
    "img.render(\n",
    "    elevation_map=\"dsm\",\n",
    "    color_map=\"ortho\",\n",
    "    position=(0.5, 1),\n",
    "    height=3000,\n",
    "    perspective=12,\n",
    ")\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = gj.Grass3dRenderer()\n",
    "img.render(\n",
    "    elevation_map=\"dsm\",\n",
    "    resolution_fine=1,\n",
    "    color_map=\"ortho\",\n",
    "    light_position=(1, 0, 0.5),\n",
    "    position=(0.75, 0.35),\n",
    "    height=1500,\n",
    "    perspective=10,\n",
    ")\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GRASS GIS modules"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GRASS functionality is available through modules (tools, functions). There are over 500 different modules in the core distribution and over 300 addon modules that can be used to prepare and analyze data.\n",
    "\n",
    "Modules 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 modules have multiple dots in their names. For example, modules staring with v.net.* deal with vector network analysis and r.in.* modules import raster data into GRASS GIS spatial database.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finding and running a module\n",
    "\n",
    "![Find modules](img/search_module.gif)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find a module for your analysis, type the term into the **search box within the Modules tab** in the Layer Manager or just browse the module tree.\n",
    "Most modules are also available from the **main menu**. For example, to find information about a raster map, use: *Raster → Reports and statistics → Basic raster metadata*.\n",
    "\n",
    "If you already know the name of the module, you can just use it in the command line. The GUI offers a **Console tab with command line** specifically built for running GRASS GIS modules. If you type the module name there, you will get suggestions for automatic completion of the name. After pressing Enter, you will get the GUI dialog for the module.\n",
    "\n",
    "Last but not least, there is a module for finding modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.search.modules keyword=zonal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "r.univar --help"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GUI  vs command line vs Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GRASS modules can be executed either through the GUI, command line or Python interfaces. The GUI offers a user-friendly approach to execute modules where the user can navigate to data layers that they would like to analyze and modify processing options with simple check boxes.\n",
    "The GUI also offers an easily accessible manual on how to execute a model. The command line interface allows users to execute a module using command prompts specific to that module. This is handy when you are running similar analyses with minor modification or are familiar with the module commands for quick efficient processing.\n",
    "\n",
    "In this example we will show how the same module can be executed with the GUI, in shell and in Python:\n",
    "\n",
    "![v.extract](img/v.extract.png) &nbsp;&nbsp;&nbsp;&nbsp;\n",
    "![v.extract](img/v.extract.gif)"
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\n",
    "    \"v.extract\",\n",
    "    input=\"roads\",\n",
    "    where=\"name = 'Umstead Drive'\",\n",
    "    output=\"umstead_drive_segments\",\n",
    ")"
   ]
  },
  {
   "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 must 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` or in menu *Settings - Region - Display region* to see current region settings\n",
    " \n",
    " ![Computational region](img/region.gif)"
   ]
  },
  {
   "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. If the raster is added to Layer Manager, right click on it and select *Set computational region from selected map*. Or run the g.region module (we include -p flag to always see the resulting region):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", 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 GUI, this can be done in the same way as for the raster map. 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 -n 10"
   ]
  },
  {
   "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": [
    "print(gs.read_command(\"g.region\", vector=\"roads\", res=2, flags=\"ap\"))"
   ]
  },
  {
   "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": [
    "gs.parse_command(\"g.region\", vector=\"roads\", align=\"dsm\", flags=\"g\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python API"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two Python APIs for accessing a module's functionality - [GRASS GIS Python Scripting Library](https://grass.osgeo.org/grass-stable/manuals/libpython/script_intro.html) (_grass.script_) and [PyGRASS](https://grass.osgeo.org/grass78/manuals/libpython/pygrass_index.html) (_grass.pygrass_).\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",
    "We are also using an API for Jupyter (`import grass.jupyter as gj`) which makes our workflows in notebooks more convenient.\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",
    " * [script.core.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",
    " * [script.core.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",
    " * [script.core.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",
    " * [script.core.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 common modules to simplify basic tasks. The list of convenient wrapper functions with examples includes:\n",
    "\n",
    " * Raster metadata using [script.raster.raster_info()](https://grass.osgeo.org/grass-stable/manuals/libpython/script.html#script.raster.raster_info): `gs.raster_info('dsm')`\n",
    " * Vector metadata using [script.vector.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 [script.core.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 [script.core.region()](https://grass.osgeo.org/grass-stable/manuals/libpython/script.html#script.core.region): `gs.region()`\n",
    "\n",
    "**Note for GUI:** The simplest way, especially on Windows, to execute the Python code which uses GRASS GIS packages is to use the _Simple Python Editor_ integrated in GRASS GIS (accessible from the toolbar or the _Python_ tab in the _Layer Manager_). Another option is to write the Python code in your favorite plain text editor like _Notepad++_ (note that Python editors are plain text editors). Then run the script in GRASS GIS using the main menu *File -> Launch script*.\n",
    "\n",
    "![Run Python in GUI](img/python.gif)"
   ]
  },
  {
   "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",
    "\n",
    "Topics covered:\n",
    " * viewshed computation ([r.viewshed](https://grass.osgeo.org/grass-stable/manuals/r.viewshed.html))\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",
    " * 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",
    " * reprojecting ([r.proj](https://grass.osgeo.org/grass-stable/manuals/r.proj.html))\n",
    " * resampling ([r.resample.interp](https://grass.osgeo.org/grass-stable/manuals/r.resample.interp.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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "grass -e -c ~/grassdata/dix_park/viewsheds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Python standard library and IPython packages we need.\n",
    "import os\n",
    "import subprocess\n",
    "import sys\n",
    "\n",
    "# Import package for progress bars, see binder/requirements.txt.\n",
    "from tqdm import tqdm\n",
    "\n",
    "# Ask GRASS GIS where its Python packages are.\n",
    "gisbase = subprocess.check_output([\"grass\", \"--config\", \"path\"], text=True).strip()\n",
    "os.environ[\"GISBASE\"] = gisbase\n",
    "sys.path.append(os.path.join(gisbase, \"etc\", \"python\"))\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/viewsheds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute viewshed and visible area\n",
    "\n",
    "Set computational region based on DSM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", raster=\"dsm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute viewshed, flag 'b' is for binary (0 and 1) output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\n",
    "    \"r.viewshed\",\n",
    "    input=\"dsm\",\n",
    "    output=\"viewshed\",\n",
    "    flags=\"b\",\n",
    "    coordinates=(711260, 3960860),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute basic univariate statistics, flag 'g' is for parsable output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "univar = gs.parse_command(\"r.univar\", map=\"viewshed\", flags=\"g\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get current region settings to get cell size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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} ha\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = gj.GrassRenderer()\n",
    "img.d_rast(map=\"ortho\")\n",
    "# select only cells with value 1 to visualize\n",
    "img.d_rast(map=\"viewshed\", values=1)\n",
    "img.show()"
   ]
  },
  {
   "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(\n",
    "    \"v.build.polylines\",\n",
    "    input=\"umstead_drive_segments\",\n",
    "    output=\"umstead_drive\",\n",
    "    cats=\"first\",\n",
    ")\n",
    "gs.run_command(\n",
    "    \"v.to.points\", input=\"umstead_drive\", type=\"line\", output=\"viewpoints\", dmax=50\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", raster=\"dsm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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",
    "3. Compute height above ground (DSM - DTM).\n",
    "4. Extract height above ground for the viewpoint locations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.mapcalc(\"diff = dsm - ground\")\n",
    "gs.run_command(\"v.what.rast\", map=\"viewpoints\", layer=2, raster=\"diff\", column=\"height\")\n",
    "gs.run_command(\"r.colors\", map=\"diff\", color=\"differences\")"
   ]
  },
  {
   "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",
    "\n",
    "pd.DataFrame(\n",
    "    json.loads(\n",
    "        gs.read_command(\n",
    "            \"v.db.select\",\n",
    "            map=\"viewpoints\",\n",
    "            columns=\"cat,height\",\n",
    "            layer=2,\n",
    "            format=\"json\",\n",
    "        )\n",
    "    )[\"records\"]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize the viewpoints with the height-above-ground raster. You can filter the points based on the height:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = gj.GrassRenderer()\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=10, 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 = (\n",
    "    gs.read_command(\n",
    "        \"v.out.ascii\",\n",
    "        input=\"viewpoints\",\n",
    "        separator=\"comma\",\n",
    "        layer=2,\n",
    "        where=\"height < 2\",\n",
    "    )\n",
    "    .strip()\n",
    "    .splitlines()\n",
    ")\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(\n",
    "        \"r.viewshed\", input=\"dsm\", output=name, coordinates=(x, y), max_distance=300\n",
    "    )\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": "markdown",
   "metadata": {},
   "source": [
    "List of map names is ordered based on the points on the road (this will be important later on):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "for point in tqdm(viewpoints):\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(\n",
    "        align=\"dsm\",\n",
    "        e=x + max_distance,\n",
    "        w=x - max_distance,\n",
    "        n=y + max_distance,\n",
    "        s=y - max_distance,\n",
    "    )\n",
    "    name = f\"viewshed_{cat}\"\n",
    "    gs.run_command(\n",
    "        \"r.viewshed\",\n",
    "        input=\"dsm\",\n",
    "        output=name,\n",
    "        coordinates=(x, y),\n",
    "        max_distance=max_distance,\n",
    "        env=env,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See the difference in number of cells between one of the viewsheds and the DSM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Viewshed num cells: {gs.raster_info(maps[0])['cells']}\")\n",
    "print(f\"DSM num cells: {gs.raster_info('dsm')['cells']}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cumulative viewshed\n",
    "Finally, 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",
    "\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_*\", mapset=\".\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since by default viewshed rasters have no data in areas that are not visible, we will use r.series method *count*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cumulative viewshed\n",
    "gs.run_command(\"r.series\", input=maps, output=\"cumulative_viewshed\", method=\"count\")\n",
    "# set color of cumulative viewshed from grey to yellow\n",
    "color_rules = \"\"\"\\\n",
    "0% 70:70:70\n",
    "100% yellow\n",
    "\"\"\"\n",
    "gs.write_command(\"r.colors\", map=\"cumulative_viewshed\", rules=\"-\", stdin=color_rules)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the results in an interactive map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = gj.InteractiveMap(width=600)\n",
    "# Add vector and layer control to map\n",
    "fig.add_raster(\"cumulative_viewshed\", opacity=0.6)\n",
    "fig.add_vector(\"umstead_drive\")\n",
    "fig.add_layer_control(position=\"bottomright\")\n",
    "fig.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": [
    "img = gj.Grass3dRenderer()\n",
    "img.render(\n",
    "    elevation_map=\"dsm\",\n",
    "    resolution_fine=1,\n",
    "    color_map=\"cumulative_viewshed\",\n",
    "    vline=\"umstead_drive\",\n",
    "    vline_width=3,\n",
    "    vline_color=\"white\",\n",
    "    position=[0.5, 0.8],\n",
    "    height=2500,\n",
    "    perspective=15,\n",
    ")\n",
    "img.overlay.d_legend(\n",
    "    raster=\"cumulative_viewshed\", at=(60, 97, 87, 92), color=\"white\", flags=\"f\"\n",
    ")\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Timeseries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\n",
    "    \"t.create\",\n",
    "    output=\"viewshed_walk\",\n",
    "    type=\"strds\",\n",
    "    temporaltype=\"absolute\",\n",
    "    title=\"Average temperature\",\n",
    "    description=\"Monthly temperature average in NC [deg C]\",\n",
    ")\n",
    "gs.run_command(\n",
    "    \"t.register\",\n",
    "    input=\"viewshed_walk\",\n",
    "    type=\"raster\",\n",
    "    start=\"2000-01-01 14:00\",\n",
    "    increment=\"5 minutes\",\n",
    "    maps=maps,\n",
    "    flags=\"i\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_rules = \"\"\"\\\n",
    "0% red\n",
    "100% yellow\n",
    "\"\"\"\n",
    "gs.write_command(\"t.rast.colors\", input=\"viewshed_walk\", rules=\"-\", stdin=color_rules)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = gj.TimeSeries(\"viewshed_walk\")\n",
    "img.baselayer.d_rast(map=\"ortho\")\n",
    "img.time_slider()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img.animate(duration=500, label=True, text_size=16, text_color=\"gray\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}