Full example: download DEM data and extract a river network
traudem can be used in combination with R-package elevatr to extract
a river network from elevation data retrieved from the web.
The following example shows how to derive the contour of a catchment and drainage area values for all cells within that contour by using the D8 flow direction algorithm.
The only required inputs are:
- coordinates of the catchment outlet;
- coordinates of the lower-left and top-right corners of a rectangular area containing the catchment of interest;
- EPSG code corresponding to the projection used to express coordinate data (see also https://epsg.org/home.html);
- zoom level at which the DEM data are extracted (see
elevatrdocumentation).
library(elevatr)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(terra)
#> terra 1.6.17
library(traudem)
library(shapefiles)
#> Loading required package: foreign
#>
#> Attaching package: 'shapefiles'
#> The following objects are masked from 'package:foreign':
#>
#> read.dbf, write.dbf
library(withr)
taudem_sitrep() # verify that TauDEM is correctly installed
#> ✔ Found GDAL version 3.0.4.
#> ✔ Found mpiexec (OpenRTE) 4.0.3 (MPI).
#> ✔ Found TauDEM path (/usr/local/taudem).
#> ✔ Found TauDEM executables directory (/usr/local/taudem).
#> ✔ Found all TauDEM executables.
#> ℹ Testing TauDEM on an example file (please wait a bit)...
#> ── TauDEM output ────────────────────────────────────────────────────────────────────
#> PitRemove version 5.3.9
#> Input file DEM.tif has projected coordinate system.
#> Nodata value input to create partition from file: nan
#> Nodata value recast to float used in partition raster: nan
#> Processes: 1
#> Header read time: 0.002094
#> Data read time: 0.000788
#> Compute time: 0.013749
#> Write time: 0.001609
#> Total time: 0.018240
#> This run may take on the order of 1 minutes to complete.
#> This estimate is very approximate.
#> Run time is highly uncertain as it depends on the complexity of the input data
#> and speed and memory of the computer. This estimate is based on our testing on
#> a dual quad core Dell Xeon E5405 2.0GHz PC with 16GB RAM.
#> ── End of TauDEM output ─────────────────────────────────────────────────────────────
#> ✔ Was able to launch a TauDEM example!
#> ! Double-check above output for serious error messages.
test_dir <- local_tempdir() # temporary directory storing intermediary files created by TauDEM
# Input data for the river Wigger (central Switzerland)
EPSG <- 21781 # EPSG code corresponding to CH1903 (Swiss projected coordinate system)
x_outlet <- 637478 # outlet x coordinate in CH1903 coordinates
y_outlet <- 237413 # outlet y coordinate in CH1903 coordinates
x_ll <- 620000 # x coordinate of the lower-left (SW) limit of the region
x_tr <- 660000 # x coordinate of the top-right (NE) limit of the region
y_ll <- 200000 # y coordinate of the lower-left (SW) limit of the region
y_tr <- 250000 # y coordinate of the top-right (NE) limit of the region
zoom <- 9 # corresponds to cellsize ~= 100 m at the latitude of interest
# use elevatr to download DEM
loc.df <- data.frame(x=c(x_ll, x_tr), y=c(y_ll,y_tr))
r <- rast()
crs(r) <- paste0("epsg:", EPSG)
crs_str <- crs(r)
z <- get_elev_raster(locations = loc.df, prj = crs_str, z = zoom, verbose = F)
#> Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj): PROJ/GDAL PROJ string degradation in workflow
#> repeated warnings suppressed
#> Discarded datum CH1903 in Proj4 definition
#>
#> Accessing raster elevation [-------------------------] 0%
Accessing raster elevation [=====>-------------------] 25%
Accessing raster elevation [===========>-------------] 50%
Accessing raster elevation [==================>------] 75%
Accessing raster elevation [=========================] 100%
#> Mosaicing & Projecting
z <- rast(z)
z <- crop(z, ext(x_ll, x_tr, y_ll, y_tr))
# reclassify: all pixels with elev=NA are set to 0. Then the pit remove algorithm will take care of them
z <- classify(z, matrix(c(NA, NA, 0), 1, 3))
writeRaster(z, filename = file.path(test_dir, "DEM.tif")) # write elevation raster to temporary directory
# apply TauDEM functions
# Remove pits
out_fel <- taudem_pitremove(file.path(test_dir, "DEM.tif"), quiet = TRUE)
# D8 flow directions
out_p <- taudem_d8flowdir(out_fel, quiet = TRUE)
out_p <- out_p$output_d8flowdir_grid # file path of flow direction file
# Contributing area
out_ad8 <- taudem_aread8(out_p, quiet = TRUE)
# Threshold
out_src <- taudem_threshold(out_ad8, quiet = TRUE)
# create shapefile with approximate outlet location
shp.point <- function(x, y, sname = "shape"){ # function to create point shapefile given coordinates
n <- length(x)
dd <- data.frame(Id = 1:n, X = x, Y = y)
ddTable <- data.frame(Id = c(1), Name = paste("Outlet", 1:n, sep = ""))
ddShapefile <- convert.to.shapefile(dd, ddTable, "Id", 1)
write.shapefile(ddShapefile, sname, arcgis = T)
invisible(paste0(sname,".shp"))
}
out_shp <- shp.point(x_outlet, y_outlet, file.path(test_dir, "ApproxOutlet"))
# Move outlet to stream
out_moved.shp <- taudem_moveoutletstostream(out_p, out_src, outlet_file = out_shp,
output_moved_outlets_file = file.path(test_dir, "Outlet.shp"),
quiet = TRUE)
# Contributing area upstream of outlet
out_ssa <- taudem_aread8(out_p, outlet_file = out_moved.shp, quiet = TRUE)
# Derive catchment contour
ssa <- rast(out_ssa)
# reclassify: pixels within the catchment have value 1, others 0
ssa_cont <- classify(ssa, matrix(c(NA, NA, 0, 1, Inf, 1), 2, 3, byrow = T))
# plot DEM of the region and catchment contours
plot(z, col = terrain.colors(1000), range = c(0,1200))
contour(ssa_cont, add = T, labels = "")
title("DEM and catchment contour")
# plot map of drainage area within the catchment
plot(ssa, col = hcl.colors(1000, "viridis"))
title("Drainage area [no. cells]")# total catchment area
ssa_log <- ssa
values(ssa_log) <- log10(values(ssa))
plot(ssa_log, col = hcl.colors(1000, "viridis"))
title("Drainage area [log10(no. cells)]")
cellsize <- sqrt(prod(res(z))) # cellsize is geometric mean of downloaded raster resolution
max(values(ssa),na.rm=T)*cellsize^2 # total catchment area in m^2
#> [1] 363822337The total area of the extracted catchment is 364 km2.