|
@@ -0,0 +1,382 @@
|
|
|
|
+{
|
|
|
|
+ "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`, a *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",
|
|
|
|
+ "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(\"../../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 GrassRenderer map\n",
|
|
|
|
+ "# GrassRenderer makes non-interactive maps using a PNG image\n",
|
|
|
|
+ "img = gj.GrassRenderer()\n",
|
|
|
|
+ "\n",
|
|
|
|
+ "# Add a raster, vector and legend to the map\n",
|
|
|
|
+ "img.d_rast(map=\"elevation\")\n",
|
|
|
|
+ "img.d_legend(raster=\"elevation\", at=(65, 90, 85, 88), fontsize=12, flags=\"b\", title=\"DTM\")\n",
|
|
|
|
+ "\n",
|
|
|
|
+ "# Display map\n",
|
|
|
|
+ "img.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": [
|
|
|
|
+ "fig = 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",
|
|
|
|
+ "fig.add_raster(\"elevation\")\n",
|
|
|
|
+ "fig.add_raster(\"drainage\", opacity=1.0)\n",
|
|
|
|
+ "fig.add_raster(\"flowacc\", opacity=1.0)\n",
|
|
|
|
+ "fig.add_raster(\"watersheds\", opacity=1.0)\n",
|
|
|
|
+ "fig.add_vector(\"streams\")\n",
|
|
|
|
+ "\n",
|
|
|
|
+ "fig.add_layer_control()\n",
|
|
|
|
+ "\n",
|
|
|
|
+ "fig.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. We'll use GrassRenderer here\n",
|
|
|
|
+ "gs.run_command(\"r.colors\", map=\"watershed_area\", color=\"plasma\")\n",
|
|
|
|
+ " \n",
|
|
|
|
+ "watershed_map = gj.GrassRenderer()\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.GrassRenderer()\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_map2 = gj.GrassRenderer()\n",
|
|
|
|
+ "watershed_map2.d_rast(map=\"elevation\")\n",
|
|
|
|
+ "watershed_map2.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_map2.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
|
|
|
|
+}
|