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.
This interactive notebook is available online thanks to the https://mybinder.org service. To run the select part (called a cell), hit Shift + Enter
.
# Import Python standard library and IPython packages we need.
import os
import subprocess
import sys
import csv
import numpy as np
import matplotlib.pyplot as plt
from collections import defaultdict
# Ask GRASS GIS where its Python packages are.
sys.path.append(
subprocess.check_output(["grass", "--config", "python_path"], text=True).strip()
)
# Import the GRASS GIS packages we need.
import grass.script as gs
import grass.jupyter as gj
# Start GRASS Session
session = gj.init("../../data/grassdata", "nc_basic_spm_grass7", "user1")
# Set computational region to elevation raster
gs.run_command("g.region", raster="elevation", flags="pg")
First, let's view the elevation raster to get an overview of the area
# Start a Map
elev_map = gj.Map()
# Add a raster, vector and legend to the map
elev_map.d_rast(map="elevation")
elev_map.d_legend(raster="elevation", at=(65, 90, 85, 88), fontsize=12, flags="b", title="DTM")
# Display map
elev_map.show()
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.
gs.run_command("r.fill.dir", input="elevation", output="elev_fill1", direction="dir1", areas="area1")
gs.run_command("r.fill.dir", input="elev_fill1", output="elev_fill2", direction="dir2", areas="area2")
gs.run_command("r.fill.dir", input="elev_fill2", output="elev_fill3", direction="dir3", areas="area3")
gs.mapcalc("depr_bin = if((elevation-elev_fill3) < 0., 1, null())")
gs.run_command("r.colors", map="depr_bin", color="blues")
# Display the depressions with InteractiveMap to see how they compare to existing waterbodies
depr_map = gj.InteractiveMap()
depr_map.add_raster("depr_bin")
depr_map.add_layer_control()
depr_map.show()
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.
It may take a minute for this cell to run.
gs.run_command("r.watershed",
elevation="elevation@PERMANENT",
drainage="drainage", # Drainage Direction
accumulation="flowacc", # Flow Accumulation
basin="watersheds",
stream="streams",
threshold=80000)
# Convert streams raster to vector
gs.run_command("r.to.vect", input="streams", output="streams", type="line")
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.
hydro_map = gj.InteractiveMap(height=400, width=600)
# We can modify with color table for rasters with `r.colors`.
# Note that if the raster is located in a different mapset (for example,
# elevation is in PERMANENT, not user1), the `r.colors` will not change
# the color in InteractiveMap.
gs.run_command("r.colors", map="drainage", color="aspect")
# Add elements to map
# We set opacity to 1.0 (default is 0.8) so layers won't interfere with eachother.
hydro_map.add_raster("elevation")
hydro_map.add_raster("drainage", opacity=1.0)
hydro_map.add_raster("flowacc", opacity=1.0)
hydro_map.add_raster("watersheds", opacity=1.0)
hydro_map.add_vector("streams")
hydro_map.add_layer_control()
hydro_map.show()
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.
# Count cells in each watershed
gs.run_command("r.stats.zonal", base="watersheds", cover="elevation", method="count", output="watersheds_count")
# Get projection resolution
proj=gs.parse_command("g.region", flags="m")
# Multiply N-S resollution by E-W resolution to get cell area
cell_area = float(proj["nsres"])*float(proj["ewres"])
# Calculate watersheds areas and convert from m2 to km2
gs.mapcalc("'watershed_area' = float('watersheds_count'*{})/1000000".format(cell_area))
Create choropleth map of watershed area.
# Display a map of watershed areas.
gs.run_command("r.colors", map="watershed_area", color="plasma")
watershed_map = gj.Map()
watershed_map.d_rast(map="watershed_area")
watershed_map.d_legend(raster="watershed_area",
bgcolor="none",
color="black",
border_color="none",
at=(3, 40, 84, 88),
lines=2,
fontsize=15,
title="Area",
units=" km2")
watershed_map.show()
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.
We start by computing the slope.
# Compute Slope
gs.run_command("r.slope.aspect", elevation="elevation", slope="slope")
# Display slope map
slope_map = gj.Map()
slope_map.d_rast(map="slope")
slope_map.d_legend(raster="slope", at=(65, 90, 85, 90), fontsize=15, flags="b", title="Slope", units="°")
slope_map.show()
Now, we use r.univar
to calculate the average slope in each watershed and return a csv.
separator = "|"
columns = defaultdict(list) # each value in each column is appended to a list
text = gs.read_command("r.univar", map="elevation", zones="watersheds", separator=separator, flags="t")
reader = csv.DictReader(text.splitlines(), delimiter=separator)
for row in reader: # read a row as {column1: value1, column2: value2,...}
for (k,v) in row.items(): # go over each column name and value
columns[k].append(v) # append the value into the appropriate list
# based on column name k
watersheds = columns["zone"]
means = np.array(columns["mean"], dtype=np.float32)
stddevs = np.array(columns["stddev"], dtype=np.float32)
# Make a bar plot of average slope by watershed
bar_positions = np.arange(len(watersheds))
plt.style.use("ggplot")
fig, ax = plt.subplots()
ax.set_title("Average Slope", fontsize=16)
ax.set_xlabel("Watershed")
ax.set_ylabel("Slope [degrees]")
ax.bar(bar_positions, means)
ax.set_xticks(bar_positions)
ax.set_xticklabels(watersheds)
plt.show()
Convert watersheds from raster to vector to make a nice map. Label the watersheds so we can compare to the bar chart above.
# Convert to vector
gs.run_command("r.to.vect", flags="s", input="watersheds", output="watersheds_vector", type="area")
# Display
watershed_vect_map = gj.Map()
watershed_vect_map.d_rast(map="elevation")
watershed_vect_map.d_vect(map="watersheds_vector",
fill_color="none",
width=1.5,
color="black",
attribute_column="value",
label_bgcolor="black",
label_color="white",
label_size=10)
watershed_vect_map.show()