{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hydrology with GRASS GIS\n",
    "\n",
    "This is a short introduction to common hydrologic workflows in *GRASS GIS* in *Jupyter Notebook*. In addition to common *Python* packages, it demonstrates the usage of `grass.script`, the *Python* API for GRASS GIS, and `grass.jupyter`, an experimental *Jupyter Notebook* specific package that helps with the launch of *GRASS GIS* and with displaying maps.\n",
    "\n",
    "This interactive notebook is available online thanks to the [https://mybinder.org](Binder) service. To run the select part (called a *cell*), hit `Shift + Enter`.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Starting GRASS in Jupyter Notebooks"
   ]
  },
  {
   "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",
    "import csv\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import defaultdict\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(\"../../data/grassdata\", \"nc_basic_spm_grass7\", \"user1\")\n",
    "\n",
    "# Set computational region to elevation raster\n",
    "gs.run_command(\"g.region\", raster=\"elevation\", flags=\"pg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's view the elevation raster to get an overview of the area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Start a Map\n",
    "elev_map = gj.Map()\n",
    "\n",
    "# Add a raster, vector and legend to the map\n",
    "elev_map.d_rast(map=\"elevation\")\n",
    "elev_map.d_legend(raster=\"elevation\", at=(65, 90, 85, 88), fontsize=12, flags=\"b\", title=\"DTM\")\n",
    "\n",
    "# Display map\n",
    "elev_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Depression Filling\n",
    "\n",
    "Depression filling is often necessary for certain flow routing algorithms. In this section, we'll find out how extensive the depressions are in our DEM using `r.fill.dir`. Note that r.watershed doesn't need any depression filling thanks to its underlying algorithm which uses least cost path to get over depressions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.fill.dir\", input=\"elevation\", output=\"elev_fill1\", direction=\"dir1\", areas=\"area1\")\n",
    "gs.run_command(\"r.fill.dir\", input=\"elev_fill1\", output=\"elev_fill2\", direction=\"dir2\", areas=\"area2\")\n",
    "gs.run_command(\"r.fill.dir\", input=\"elev_fill2\", output=\"elev_fill3\", direction=\"dir3\", areas=\"area3\")\n",
    "gs.mapcalc(\"depr_bin = if((elevation-elev_fill3) < 0., 1, null())\")\n",
    "gs.run_command(\"r.colors\", map=\"depr_bin\", color=\"blues\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Display the depressions with InteractiveMap to see how they compare to existing waterbodies\n",
    "depr_map = gj.InteractiveMap()\n",
    "depr_map.add_raster(\"depr_bin\")\n",
    "depr_map.add_layer_control()\n",
    "depr_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing Watersheds, Drainage Direction, Flow Accumulation, and Streams\n",
    "\n",
    "From the elevation raster, we compute the watersheds, drainage direction and flow accumulation and display the results. Since `r.watershed` uses a least cost algorithm, we don't need to use the depression-filled raster; instead, we'll use the original elevation raster.\n",
    "\n",
    "It may take a minute for this cell to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.watershed\", \n",
    "               elevation=\"elevation@PERMANENT\",\n",
    "               drainage=\"drainage\",    # Drainage Direction\n",
    "               accumulation=\"flowacc\", # Flow Accumulation\n",
    "               basin=\"watersheds\",\n",
    "               stream=\"streams\",\n",
    "               threshold=80000)\n",
    "\n",
    "# Convert streams raster to vector\n",
    "gs.run_command(\"r.to.vect\", input=\"streams\", output=\"streams\", type=\"line\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, to view and compare the ouputs of `r.watersheds`, we'll use `grass.jupyter`'s `InteractiveMap` class which allows us to toggle between layers and zoom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "hydro_map = gj.InteractiveMap(height=400, width=600)\n",
    "\n",
    "# We can modify with color table for rasters with `r.colors`.\n",
    "# Note that if the raster is located in a different mapset (for example,\n",
    "# elevation is in PERMANENT, not user1), the `r.colors` will not change \n",
    "# the color in InteractiveMap.\n",
    "gs.run_command(\"r.colors\", map=\"drainage\", color=\"aspect\")\n",
    "\n",
    "# Add elements to map\n",
    "# We set opacity to 1.0 (default is 0.8) so layers won't interfere with eachother.\n",
    "hydro_map.add_raster(\"elevation\")\n",
    "hydro_map.add_raster(\"drainage\", opacity=1.0)\n",
    "hydro_map.add_raster(\"flowacc\", opacity=1.0)\n",
    "hydro_map.add_raster(\"watersheds\", opacity=1.0)\n",
    "hydro_map.add_vector(\"streams\")\n",
    "\n",
    "hydro_map.add_layer_control()\n",
    "\n",
    "hydro_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Watershed Area\n",
    "\n",
    "With our watersheds, we can compute some zonal statistics. In this section, we use the `count` method in `r.stats.zonal` to make a map of watershed area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Count cells in each watershed\n",
    "gs.run_command(\"r.stats.zonal\", base=\"watersheds\", cover=\"elevation\", method=\"count\", output=\"watersheds_count\")\n",
    "\n",
    "# Get projection resolution\n",
    "proj=gs.parse_command(\"g.region\", flags=\"m\")\n",
    "\n",
    "# Multiply N-S resollution by E-W resolution to get cell area\n",
    "cell_area = float(proj[\"nsres\"])*float(proj[\"ewres\"])\n",
    "\n",
    "# Calculate watersheds areas and convert from m2 to km2\n",
    "gs.mapcalc(\"'watershed_area' = float('watersheds_count'*{})/1000000\".format(cell_area))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create choropleth map of watershed area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display a map of watershed areas.\n",
    "gs.run_command(\"r.colors\", map=\"watershed_area\", color=\"plasma\")\n",
    "               \n",
    "watershed_map = gj.Map()\n",
    "watershed_map.d_rast(map=\"watershed_area\")\n",
    "watershed_map.d_legend(raster=\"watershed_area\",\n",
    "                       bgcolor=\"none\",\n",
    "                       color=\"black\",\n",
    "                       border_color=\"none\",\n",
    "                       at=(3, 40, 84, 88),\n",
    "                       lines=2,\n",
    "                       fontsize=15,\n",
    "                       title=\"Area\",\n",
    "                       units=\" km2\")\n",
    "watershed_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zonal Statistics: Average Slope by Watershed\n",
    "\n",
    "In this section, we compute average slope and standard deviation in each watershed then make a bar plot to compare them. Each watershed is a zone. We use `r.univar` to find compute a table of univariate statistics. An alternative approach would be to use `r.stats.zonal` which returns a raster. \n",
    "\n",
    "We start by computing the slope."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute Slope\n",
    "gs.run_command(\"r.slope.aspect\", elevation=\"elevation\", slope=\"slope\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display slope map\n",
    "slope_map = gj.Map()\n",
    "slope_map.d_rast(map=\"slope\")\n",
    "slope_map.d_legend(raster=\"slope\", at=(65, 90, 85, 90), fontsize=15, flags=\"b\", title=\"Slope\", units=\"°\")\n",
    "slope_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we use `r.univar` to calculate the average slope in each watershed and return a csv."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "separator = \"|\"\n",
    "\n",
    "columns = defaultdict(list) # each value in each column is appended to a list\n",
    "\n",
    "text = gs.read_command(\"r.univar\", map=\"elevation\", zones=\"watersheds\", separator=separator, flags=\"t\")\n",
    "reader = csv.DictReader(text.splitlines(), delimiter=separator)\n",
    "for row in reader: # read a row as {column1: value1, column2: value2,...}\n",
    "    for (k,v) in row.items(): # go over each column name and value \n",
    "        columns[k].append(v) # append the value into the appropriate list\n",
    "                             # based on column name k\n",
    "\n",
    "watersheds = columns[\"zone\"]\n",
    "means = np.array(columns[\"mean\"], dtype=np.float32)\n",
    "stddevs = np.array(columns[\"stddev\"], dtype=np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a bar plot of average slope by watershed\n",
    "bar_positions = np.arange(len(watersheds))\n",
    "plt.style.use(\"ggplot\")\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_title(\"Average Slope\", fontsize=16)\n",
    "ax.set_xlabel(\"Watershed\")\n",
    "ax.set_ylabel(\"Slope [degrees]\")\n",
    "ax.bar(bar_positions, means)\n",
    "ax.set_xticks(bar_positions)\n",
    "ax.set_xticklabels(watersheds)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Converting to Vectors\n",
    "\n",
    "Convert watersheds from raster to vector to make a nice map. Label the watersheds so we can compare to the bar chart above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert to vector\n",
    "gs.run_command(\"r.to.vect\", flags=\"s\", input=\"watersheds\", output=\"watersheds_vector\", type=\"area\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display\n",
    "watershed_vect_map = gj.Map()\n",
    "watershed_vect_map.d_rast(map=\"elevation\")\n",
    "watershed_vect_map.d_vect(map=\"watersheds_vector\",\n",
    "                      fill_color=\"none\",\n",
    "                      width=1.5,\n",
    "                      color=\"black\",\n",
    "                      attribute_column=\"value\",\n",
    "                      label_bgcolor=\"black\",\n",
    "                      label_color=\"white\",\n",
    "                      label_size=10)\n",
    "watershed_vect_map.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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}