workshop_part_5_inundation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 3: Estimating inundation extent using HAND methodology\n",
    "\n",
    "\n",
    "In this example we will use some of GRASS GIS hydrology tools, namely:\n",
    "\n",
    "* [r.watershed](https://grass.osgeo.org/grass-stable/manuals/r.watershed.html): for computing flow accumulation, drainage direction, the location of streams and watershed basins; it does not need sink filling because of using the least-cost-path to route flow out of sinks\n",
    "* [r.lake](https://grass.osgeo.org/grass-stable/manuals/r.lake.html): to fill a lake to a target water level from a given start point or seed raster\n",
    "* [r.lake.series](https://grass.osgeo.org/grass-stable/manuals/r.lake.series.html): addon which runs r.lake for different water levels\n",
    "* [r.stream.distance](https://grass.osgeo.org/grass-stable/manuals/r.stream.distance.html): for computing the distance to streams or outlet, the relative elevation above streams; the distance and the elevation are calculated along watercourses\n",
    "\n",
    "First, let's create a new mapset *flooding* in nc_spm_08_grass7 sample dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "grass -c -e ~/grassdata/nc_spm_08_grass7/flooding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize GRASS session:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Python standard library and IPython packages we need.\n",
    "import subprocess\n",
    "import sys\n",
    "\n",
    "# Ask GRASS GIS where its Python packages are.\n",
    "sys.path.append(\n",
    "    subprocess.check_output([\"grass\", \"--config\", \"python_path\"], text=True).strip()\n",
    ")\n",
    "\n",
    "# Import the GRASS GIS packages we need.\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "# Start GRASS Session\n",
    "gj.init(\"~/grassdata\", \"nc_spm_08_grass7\", \"flooding\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Modules r.stream.distance and r.lake.series are addons and we need to install them first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.extension\", extension=\"r.stream.distance\")\n",
    "gs.run_command(\"g.extension\", extension=\"r.lake.series\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute HAND raster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will estimate inundation extent using the Height Above Nearest Drainage methodology ([A.D. Nobre, 2011](https://doi.org/10.1016/j.jhydrol.2011.03.051)). We will compute the HAND terrain model representing the differences in elevation between each grid cell and the elevations of the flowpath-connected downslope grid cells where the flow enters the channel.\n",
    "\n",
    "First we compute the flow accumulation, drainage and streams (with a threshold value of 100000). We convert the streams to vector for better visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", raster=\"elevation\")\n",
    "gs.run_command(\"r.watershed\", elevation=\"elevation\", accumulation=\"flowacc\",\n",
    "               drainage=\"drainage\", stream=\"rwatershed_streams\", threshold=100000)\n",
    "gs.run_command(\"r.to.vect\", input=\"rwatershed_streams\", output=\"rwatershed_streams\", type=\"line\")\n",
    "\n",
    "fllowacc_map = gj.Map()\n",
    "fllowacc_map.d_rast(map=\"flowacc\")\n",
    "fllowacc_map.d_vect(map=\"rwatershed_streams\", width=2, color=\"blue\")\n",
    "fllowacc_map.d_legend(raster=\"flowacc\", range=\"0,1000\")\n",
    "fllowacc_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's zoom in to see the flow accumulation raster better:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flowacc_map = gj.InteractiveMap()\n",
    "flowacc_map.add_raster(\"flowacc\")\n",
    "flowacc_map.add_vector(\"rwatershed_streams\")\n",
    "flowacc_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we use r.stream.distance with output parameter difference to compute a new raster map where each cell is the elevation difference between the cell and the the cell on the stream where the cell drains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.stream.distance\", stream_rast=\"rwatershed_streams\", direction=\"drainage\",\n",
    "               elevation=\"elevation\", method=\"downstream\", difference=\"above_stream\")\n",
    "gs.run_command(\"r.colors\", map=\"above_stream\", color=\"elevation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare original elevation and HAND raster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map3d = gj.Map3D()\n",
    "map3d.render(elevation_map=\"elevation\", color_map=\"elevation\",\n",
    "             position=(0.5, 1), height=3000, perspective=12, zexag=5)\n",
    "map3d.overlay.d_legend(raster=\"elevation\", at=(5, 80, 87, 92), flags=\"b\", border_color=\"none\")\n",
    "map3d.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map3d = gj.Map3D()\n",
    "map3d.render(elevation_map=\"above_stream\", color_map=\"above_stream\",\n",
    "             position=(0.5, 1), height=3000, perspective=12, zexag=5)\n",
    "map3d.overlay.d_legend(raster=\"above_stream\", at=(5, 80, 90, 95), flags=\"b\", border_color=\"none\")\n",
    "map3d.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inundation\n",
    "Before we compute the inundation, we will look at how r.lake works. We compute a lake from a specified coordinate pair and water level:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.lake\", elevation=\"elevation\", water_level=90, lake=\"lake\", coordinates=[637877, 218475])\n",
    "\n",
    "lake_map = gj.Map()\n",
    "lake_map.d_rast(map=\"elevation\")\n",
    "lake_map.d_rast(map=\"lake\")\n",
    "lake_map.d_legend(raster=\"lake\", label_values=\"0.1,5,10,15\", digits=0, at=(1, 40, 90, 95))\n",
    "lake_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now instead of the elevation raster we use the HAND raster to simulate 5-meter inundation and, as the seed we specify the entire stream."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.lake\", elevation=\"above_stream\", water_level=5, lake=\"flood\", seed=\"rwatershed_streams\")\n",
    "\n",
    "hand_map = gj.Map()\n",
    "hand_map.d_rast(map=\"above_stream\")\n",
    "hand_map.d_rast(map=\"flood\")\n",
    "hand_map.d_legend(raster=\"flood\", label_values=\"0.1,4,8\", digits=0, at=(1, 40, 90, 95))\n",
    "hand_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the r.lake.series addon we can create a series of inundation maps with rising water levels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from io import StringIO\n",
    "import pandas as pd\n",
    "\n",
    "gs.run_command(\"r.lake.series\", elevation=\"above_stream\", start_water_level=0, end_water_level=5, \n",
    "               water_level_step=0.5, output=\"inundation\", seed_raster=\"rwatershed_streams\", quiet=True)\n",
    "gs.run_command(\"t.rast.colors\", input=\"inundation\", color=\"water\")\n",
    "pd.read_csv(StringIO(gs.read_command(\"t.rast.list\", input=\"inundation\", separator=\"comma\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "r.lake.series creates a space-time dataset. We can use the [temporal modules](https://grass.osgeo.org/grass-stable/manuals/temporal.html) to further work with the data. For example, we could further compute the volume and extent of flood water using t.rast.univar:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.read_csv(StringIO(gs.read_command(\"t.rast.univar\", input=\"inundation\", separator=\"comma\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timemap = gj.TimeSeriesMap()\n",
    "timemap.d_rast(map=\"elevation_shade\")\n",
    "timemap.d_vect(map=\"streets_wake\", color=\"#4D4D4D\")\n",
    "timemap.add_raster_series(\"inundation\")\n",
    "timemap.d_legend(at=(1, 40, 90, 95))\n",
    "timemap.show()"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}