In flee-group/watershed: Tools for Watershed Delineation

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning=FALSE, message=FALSE
)
knitr::opts_chunk$set(fig.width=5.5, fig.height=5.5, collapse = TRUE, comment = "##", dev="png")
par(mar = c(0,0,0,0))

0. Preparing the installation

Be sure you have followed the installation guide and have a working installation of watershed, rgrass7, and Grass GIS 7.4 or 7.6 (other version may work, but are untested).

0.5. Preparing the DEM

Many workflows advise either using a "burned-in" DEM or a "filled" DEM to avoid topological problems. The tools used by watershed should in theory not require either of those steps, so in general it is recommended you skip this step unless you encounter problems later on. Burning in a DEM is not yet supported, so if this step is required it is best to use another tool and then bring the already burned DEM into R. Filling the DEM is supported however; see ?fill_dem.

1. Delineating the stream

The pared-down workflow offered by watershed provides a relatively quick way to go from a dem to a stream map. We will make a first-draft map using the kamp_dem dataset provided by this package, for the Kamp river in Austria. The DEM is a raster, so we load this package as well.

The delineate() function is the main workhorse here. The resulting raster stack contains 3 layers: accum is the flow accumulation, drainage is the drainage direction, and stream is the stream map. We also produce a vector layer of the stream map for plotting. Some notes about our use of delineate here:

We do not set the threshold parameter, which is the minimum catchment area before something is marked as a stream. The default of 1e6 $m^2$ (10 $km^2$) works well; increasing the threshold will result in more streams, decreasing it will take longer and make more streams.
We restrict the output to a single river basin with outlet = NA (this chooses the largest basin). To select a different basin, we could instead use outlet = c(4704588, 2847762); the coordinates must be exactly on the stream for this to work.
Because outlet is provided, the output rasters also have their extent reduced to the catchment of our stream.

library(watershed)
library(raster)
library(sf)
data(kamp_dem)
kamp = delineate(kamp_dem, outlet = NA)
kv = vectorise_stream(kamp[["stream"]], Tp=kamp_Tp)
plot(kamp_dem, col=terrain.colors(20), axes = FALSE)
plot(st_geometry(kv), col='blue', add = TRUE)

2. Generating topologies

Many later features require a topology; this simply a square matrix T, with one row/column for each pixel or reach in the river network. T[i,j] == 1 indicates that pixel/reach i is directly upstream of pixel/reach j (in other words, j gets water from i without that water flowing through any other intermediate pixels/reaches).

Generating the topology can take some time for large and/or complex river networks, so it is best to do it once you have a reasonable first pass.

Our first step is to crop the rasters from the previous steps to only include our catchment of interest. Then we load the Matrix package to handle the sparseMatrix returned by pixel_topology. Finally, we create two topologies; one for the pixels, and one for reaches.

library(Matrix)
kamp_Tp = pixel_topology(kamp)
kamp_Tr = reach_topology(kamp, kamp_Tp)

Note the warning, indicating that there is a confluence where more than 2 streams flow into the same pixel. This is not a serious problem, so we will ignore it for now.

3. Producing smaller reaches

The default output from delineate is to assign each pixel to a reach, where a reach is a stretch of river between confluences. For many applications, it is useful to split these into smaller sections. watershed includes a function to split reaches into fixed (approximately) equal-length sections.

Note that resizing the reaches invalidates the old topologies, and a new ones must be made.
Alternatively, we could define the reach lengths directly in delineate by setting the reach_len parameter.
Setting trim=TRUE will remove small (by default, length < target_length/2) headwater reaches from the network

new_reaches = resize_reaches(kamp$stream, kamp_Tp, 2000, trim=TRUE)
new_v = vectorise_stream(new_reaches)
new_Tp = pixel_topology(drainage = kamp$drainage, stream = new_reaches)
new_Tr = reach_topology(Tp = new_Tp, stream = new_reaches)

## how many reaches, and what range of lengths?
c(nrow(kamp_Tr), nrow(new_Tr))
library(ggplot2)
pl = ggplot(data.frame(len = c(kamp_Tr@x, new_Tr@x), 
                  reach = c(rep("old", length(kamp_Tr@x)), rep("new", length(new_Tr@x)))), aes(x=len, fill=reach))
pl + geom_histogram(aes(y=..density..), position="dodge", binwidth = 200, colour="gray") + theme_minimal() + xlab("Reach length (m)")
## update the raster stack to use the new streams
kamp[["stream"]] = new_reaches
kv = new_v
kamp_Tp = new_Tp
kamp_Tr = new_Tr
rm(list=c("new_reaches", "new_v", "new_Tp", "new_Tr"))

4. Computing catchment area

watershed can compute catchment area for you at any set of points in the watershed; this will yield the total area draining into a given point. Note that this can take a very long time for large areas. Here, we will compute the catchment area per reach.

CA = catchment(kamp, type = 'reach', Tp = kamp_Tp)
rids = unique(kamp$stream)
## not all reaches survive the transition to vector, and the orders are not the same, so we match by rid
## We also convert to log scale to make the map more readable
kv$catch_area = CA[kv$reach_id]
kv$log_catchment_area = log(kv$catch_area)
plot(kv['log_catchment_area'])

Discharge and hydraulic geometry

watershed includes a simple function for estimating discharge and hyrdaulic geometry at all points/reaches in the stream given. See ?hydraulic_scaling for details and references on the scaling used.

We load the dataset kamp_q for observations of discharge. We first need to make sure the points are exactly on the stream, then compute the catchment area for each point. We set units using the units package to make sure unit conversions in the scaling equations are handled correctly. Finally, we can then use the reach-level catchment areas to compute discharge and geometry for every reach in the watershed.

library(units)
data(kamp_q)
kamp_q = st_snap(kamp_q, kv, tolerance = 200)
kamp_q$catch_area = catchment(kamp, type = "points", y = st_coordinates(kamp_q))
h <- hydraulic_geometry(set_units(kamp_q$catch_area, "m^2"), kamp_q$discharge, set_units(kv$catch_area, "m^2"))
head(h)