workshop_part_2.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",
    "\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",
    " * 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": [
    "# Import Python standard library and IPython packages we need.\n",
    "import os\n",
    "import subprocess\n",
    "import sys\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(os.path.expanduser(\"~/grassdata\"), \"dix_park\", \"PERMANENT\")"
   ]
  },
  {
   "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": [
    "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",
    "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:"
   ]
  },
  {
   "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 = 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)\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)\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",
    "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,\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": {},
   "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_grouped(type=\"raster\", pattern=\"viewshed_*\")[\"PERMANENT\"]"
   ]
  },
  {
   "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 = \"\"\"0% 70:70:70\n",
    "100% yellow\"\"\"\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(elevation_map=\"dsm\", resolution_fine=1, color_map=\"cumulative_viewshed\",\n",
    "           vline=\"umstead_drive\", vline_width=3, vline_color=\"white\",\n",
    "           position=[0.5, 0.8], height=2500, perspective=15)\n",
    "img.overlay.d_legend(raster=\"cumulative_viewshed\", at=(60, 97, 87, 92), color=\"white\", flags=\"f\")\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data reprojection and resampling\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",
    "1. We first reproject R and NIR Landsat bands from NCSPM sample dataset we already have available. Module 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",
    "2. We resample each band to our current resolution (1m) using bilinear interpolation. Note that we don't have to do that if nearest neighbor resampling is sufficient.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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)\n",
    "    gs.run_command(\"r.resamp.interp\", input=f\"lsat7_2002_{band}\", output=f\"lsat7_2002_{band}_1m\", method=\"bilinear\")"
   ]
  },
  {
   "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_1m\", nir=\"lsat7_2002_40_1m\", output=\"ndvi\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = gj.GrassRenderer()\n",
    "img.d_rast(map=\"ndvi\")\n",
    "img.d_legend(raster=\"ndvi\")\n",
    "img.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",
    "\n",
    "![Mask](img/mask.gif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.mask\", raster=\"cumulative_viewshed\", maskcats=\"1 thru 6\")\n",
    "data = gs.parse_command(\"r.univar\", map=\"ndvi\", flags=\"g\")\n",
    "print(f\"Average NDVI of visible cells: {data['mean']} ± {data['stddev']}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = gj.GrassRenderer()\n",
    "img.d_rast(map=\"ndvi\")\n",
    "img.d_legend(raster=\"ndvi\", flags=\"d\")\n",
    "img.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": [
    "img = gj.GrassRenderer()\n",
    "img.d_histogram(map=\"ndvi\", bgcolor=\"grey\")\n",
    "img.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",
   "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}