workshop_part_2_viewshed.ipynb

{
 "cells": [
  {
   "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 project *dix_park*. This automatically created new default 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": [
    "Schema of GRASS project dix_park's content:\n",
    "\n",
    "<img src=\"img/data_structure2.png\" alt=\"GRASS project dix_park\" width=\"400\"/>"
   ]
  },
  {
   "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",
    "Because the road consists of several segments, we will first merge them into one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"v.build.polylines\", input=\"umstead_drive_segments\", output=\"umstead_drive\", cats=\"first\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then 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.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) using raster algebra:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.mapcalc(\"diff = dsm - ground\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the color ramp of the raster to \"differences\", which will highlight in red areas with vegetation and buildings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "We will compute viewsheds from all the viewpoints we generated earlier and from those we compute a cumulative viewshed.\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": [
    "<div class=\"alert alert-info\">\n",
    "In part 3 of this workshop we will parallelize the viewshed computation.\n",
    "<div>"
   ]
  },
  {
   "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": [
    "info = gs.parse_command(\"t.info\", input=\"viewsheds\", flags=\"g\")\n",
    "pd.DataFrame(info.values(), index=info.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To list the individual rasters, we will use t.rast.list."
   ]
  },
  {
   "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\")))"
   ]
  },
  {
   "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 North Carolina sample dataset (State Plane meters, EPSG:3358) 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\")"
   ]
  }
 ],
 "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
}