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 Ybbs river in Austria. The DEM is a raster, so we load this package as well.
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
thresholdparameter, 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 useoutlet = c(4704588, 2847762); the coordinates must be exactly on the stream for this to work.
+ - Because
outletis provided, the output rasters also have their extent reduced to the catchment of our stream.
+
library(raster)
-## Loading required package: sp
+library(sf)
data(kamp_dem)
-plot(kamp_dem, col=terrain.colors(20), axes = FALSE)
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. Note that the plotted stream looks “speckled” because it is a raster, and we don’t plot every pixel to save time; later you will see how to produce a vector map of streams.
kamp = delineate(kamp_dem)
-plot(kamp_dem, col=terrain.colors(20), axes = FALSE)
-plot(kamp$stream, col='blue', add = TRUE, legend = FALSE)
The defaults produce pretty good results in this case. Perhaps we know from on-the-ground work, or aerial photo inspection, that there are more streams that what is produced here. We can improve this with the threshold parameter, which sets the default minimum basin size. This defaults to 1e6 m^2, or 1 km^2. Here we reduce it by half; note that this gives us a warning about possible long computation times; it’s not such a problem here, because the Ybbs is a relatively small area, but with large DEMs this can make the computation very slow.
kamp2 = delineate(kamp_dem, threshold = 0.5e5)
-## Warning in delineate(kamp_dem, threshold = 50000): Small threshold;
-## excessive computation time and memory usage are possible if threshold not
-## increased
-plot(kamp_dem, col=terrain.colors(20), axes = FALSE)
-plot(kamp2$stream, col='blue', add = TRUE, legend = FALSE)
This produces a lot of streams! From our fieldwork in the Kamp, we know that the default produces a pretty good network, so we will keep the first one.
-
