--- title: "Preparation of input data: geomorphological analysis with whitebox" author: "Alban de Lavenne & Antoine Casquin" date: "`r Sys.Date()`" bibliography: "../inst/REFERENCES.bib" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Preparation of input data: geomorphological analysis with whitebox} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} fig_width: 8 fig_height: 4 --- ```{r setup, include = FALSE} library(transfR) knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` The hydrological modelling of the `transfR` package is based on a geomorphological analysis of the studied catchments. In this vignette, we give some guidance on how to perform this geomorphological analysis. More specifically, we extract the catchment delineation and hydraulic length maps from a digital elevation model (DEM). This analysis is one of the two inputs needed (together with the time series of flow observations) to build a `transfR` object and start using the `transfR` package (see the *Get started with transfR* vignette). Hydraulic length is defined as the distance within the river network along an identified flow path to the outlet. It can be extracted from a DEM in many different ways, such as with the GRASS toolkits [@Jasiewicz2011], Whitebox GAT (see @Lindsay2016 or [WhiteboxTools](https://github.com/jblindsay/whitebox-tools)), TauDEM (D. Tarboton, Utah State University) or online services (@Squividant2015 for catchment delineation only). This vignette presents one possible workflow by making use of two main R packages: * [`elevatr`](https://cran.r-project.org/package=elevatr), which provides elevation data worldwide from various APIs; * [`whitebox`](https://cran.r-project.org/package=whitebox), an R frontend for the 'WhiteboxTools' library, which is an advanced geospatial data analysis platform. Functions in the `whitebox` package generally work with files on disk rather than objects in R memory. A working directory should therefore be defined for these files. ```{r, echo=TRUE, message=FALSE, warning=FALSE, results='hide'} wbt_wd <- tempdir(check = TRUE) ``` ## 1. Retrieve elevation data with `elevatr` The `elevatr` package retrieves grids (rasters) of elevation data worldwide for various zoom levels and data sources. It uses the [Open Topography Global Datasets API](https://opentopography.org/developers) to access Shuttle Radar Topography Mission (SRTM) data. A projected coordinate reference system must be assigned to the catchment delineation and hydraulic length maps used with the `whitebox` package. We use the Lambert-93 projection, the official projection for maps of metropolitan France, whose [EPSG code is 2154](https://epsg.io/2154). ```{r, download_dem, echo=TRUE, message=FALSE, warning=FALSE, eval=FALSE, results='hide'} library(elevatr) library(progress) # Needed by elevatr # Set up a projection (French Lambert-93 projection) EPSG <- 2154 # Define a bbox that will encompass the catchments of the study area blavet_bbox <- st_bbox(c(xmin = -3.3, xmax = -2.7, ymax = 48.11, ymin = 47.77), crs = st_crs(4326)) blavet_loc <- st_as_sfc(blavet_bbox) |> st_sf() # Retrieve elevation data as raster dem_raw <- elevatr::get_elev_raster(blavet_loc, z = 10) # ~76 m resolution # Project and define spatial resolution: dem_100m <- st_warp(st_as_stars(dem_raw), cellsize = 100, crs = st_crs(EPSG)) names(dem_100m) <- "warp" # Set negative values (ocean) to NA dem_100m[dem_100m < 0] <- NA # Write to file write_stars(dem_100m["warp"], file.path(wbt_wd,"dem_100m.tif")) ``` ```{r, echo=FALSE, message=FALSE, warning=TRUE, eval=TRUE, results='hide'} try_chunk <- try({ <> }, silent = TRUE) if(inherits(try_chunk, "try-error")){ warning("\nIssue when downloading elevation data. \nThe vignette will not be fully built.") running <- FALSE }else{ running <- TRUE } ``` ## 2. Retrieve a known river network and burn it into the DEM (optional) The hydrological modelling distinguishes hillslopes from the river network. They have very different transfer dynamics, and the `transfR` package aims to describe only the transfer function of the river network. Defining where this river network begins and which flow path it takes is a key, non-trivial issue for hydrogeomorphologists. The easiest way to draw a drainage network from a DEM is usually to define a minimum drainage area threshold, above which the drainage network is assumed to start. However, defining this threshold can be difficult because it may vary spatially, especially with the geology of the region. It may therefore be better to use a known river network and force the flow paths to follow it. Here we use the [French TOPAGE river network](https://www.data.gouv.fr/datasets/bd-topage-r) as a reference (see the description of the `Blavet` dataset to download it from the "Sandre - Eau France" Web Feature Service (WFS)). We implement a stream-burning technique in the DEM with [`whitebox::wbt_burn_streams_at_roads()`](https://whiteboxr.gishub.org/reference/wbt_burn_streams_at_roads.html), following @Lindsay2016a. ```{r, install_whitebox, echo=FALSE, message=FALSE, warning=TRUE, eval=running, results='hide'} # If the WhiteboxTools executable is not present, install it in the temporary directory library(whitebox) if(!wbt_init()){ wbt_inst <- try(install_whitebox(pkg_dir = wbt_wd), silent = TRUE) if(!inherits(wbt_inst, "try-error")){ exe_path <- file.path(wbt_wd, "WBT", "whitebox_tools") # Unix if(!file.exists(exe_path)) exe_path <- paste0(exe_path,".exe") # Windows if(!file.exists(exe_path)){ warning("WhiteboxTools executable not found") running <- FALSE }else{ wbt_options(exe_path = exe_path, max_procs = 2, wd = wbt_wd, verbose = TRUE) } }else{running <- FALSE} }else{ wbt_options(max_procs = 2, wd = wbt_wd, verbose = TRUE) } ``` ```{r, burn_stream, echo=TRUE, message=FALSE, warning=FALSE, eval=FALSE, results='hide'} library(transfR) library(whitebox) # Get the French Topage river network from the Blavet dataset data(Blavet) CoursEau_Topage2019 <- Blavet$network # Change projection and write files network_topage <- st_transform(CoursEau_Topage2019, EPSG) st_write(network_topage, file.path(wbt_wd, "network_topage.shp"), delete_layer = TRUE, quiet = TRUE) whitebox::wbt_rasterize_streams("network_topage.shp", base = "dem_100m.tif", output = "network_topage.tif", nodata = 0, wd = wbt_wd) # Burn this river network on the DEM # We will neglect the effect of the road embankments at this DEM resolution of 100m # by creating an empty shapefile for roads st_write(st_sfc(st_multilinestring(),crs = EPSG), file.path(wbt_wd,"roads.shp"), delete_layer = TRUE, quiet = TRUE) whitebox::wbt_burn_streams_at_roads(dem = "dem_100m.tif", streams = "network_topage.shp", roads = "roads.shp", output = "dem_100m_burn.tif", wd = wbt_wd) ``` ```{r, echo=FALSE, message=FALSE, warning=TRUE, eval=running, results='hide'} try_chunk <- try({ <> }, silent = TRUE) if(inherits(try_chunk, "try-error")){ warning("\nIssue when burning the river network on the DEM. \nThe vignette will not be fully built.") running <- FALSE } ``` ## 3. Delineate catchments from their outlet coordinates The `whitebox` package provides all the tools required for a standard catchment delineation workflow from a DEM and outlet coordinates. It consists of six steps: * [Fill the depressions](https://whiteboxr.gishub.org/reference/wbt_fill_depressions.html), * [Compute flow direction](https://whiteboxr.gishub.org/reference/wbt_d8_pointer.html) using the D8 algorithm, * [Compute flow accumulation](https://whiteboxr.gishub.org/reference/wbt_d8_flow_accumulation.html) according to flow direction, * [Extract a river network](https://whiteboxr.gishub.org/reference/wbt_extract_streams.html) using a threshold in flow accumulation, * [Snap the outlets](https://whiteboxr.gishub.org/reference/wbt_jenson_snap_pour_points.html) to this river network, * [Delimit the catchments](https://whiteboxr.gishub.org/reference/wbt_watershed.html). ```{r, extract_stream, echo=TRUE, message=FALSE, warning=FALSE, eval=FALSE, results='hide'} # Remove the depressions on the DEM whitebox::wbt_fill_depressions(dem = "dem_100m_burn.tif", output = "dem_fill.tif", wd = wbt_wd) # Flow direction raster whitebox::wbt_d8_pointer(dem = "dem_fill.tif", output = "d8.tif", wd = wbt_wd) # Compute flow accumulation whitebox::wbt_d8_flow_accumulation(input = "d8.tif", pntr = TRUE, output ="facc.tif", wd = wbt_wd) # Extract a stream network (threshold = 1 km^2) consistent with flow direction whitebox::wbt_extract_streams(flow_accum = "facc.tif", threshold = 100, # 100 cells for 1 km^2 output = "network_1km2.tif", zero_background = TRUE, wd = wbt_wd) whitebox::wbt_remove_short_streams(d8_pntr = "d8.tif", streams = "network_1km2.tif", output = "network_d8.tif", min_length= 200, wd = wbt_wd) ``` ```{r, echo=FALSE, message=FALSE, warning=TRUE, eval=running, results='hide'} try_chunk <- try({ <> }, silent = TRUE) if(inherits(try_chunk, "try-error")){ warning("\nIssue when extracting the river network from the DEM. \nThe vignette will not be fully built.") running <- FALSE } ``` Outlet coordinates are retrieved from [hydro.eaufrance.fr](https://www.hydro.eaufrance.fr) and snapped to the river-network pixel that is consistent with the previously defined flow directions. ```{r, delineate_catchments, echo=TRUE, message=FALSE, warning=FALSE, eval=FALSE, results='hide'} # Locate the outlets of the studied catchments (with manual adjustments to facilitate snapping) outlets_coordinates <- data.frame(id = names(Blavet$hl), X = c(254010.612,255940-100,255903,237201,273672,265550), Y = c(6772515.474,6776418-200,6776495,6774304-200,6762681,6783313)) outlets <- st_as_sf(outlets_coordinates, coords = c("X", "Y"), crs=2154) st_write(outlets, dsn = file.path(wbt_wd, "outlets.shp"), delete_layer = TRUE, quiet = TRUE) # Snap the outlets on the stream raster whitebox::wbt_jenson_snap_pour_points(pour_pts = "outlets.shp", streams = "network_d8.tif", output = "outlets_snapped.shp", snap_dist = 200, wd = wbt_wd) outlets_snapped <- st_read(file.path(wbt_wd, "outlets_snapped.shp"), quiet = TRUE) # Delineate catchments catchments <- st_sf(id = character(), geometry = st_sfc(crs = EPSG)) for(id in outlets_snapped$id){ st_write(outlets_snapped[outlets_snapped$id==id,], file.path(wbt_wd, paste0(id, "_outlet.shp")), delete_layer = TRUE, quiet = TRUE) whitebox::wbt_watershed(d8_pntr = "d8.tif", pour_pts = paste0(id, "_outlet.shp"), output = paste0(id, "_catchment.tif"), wd = wbt_wd) # Vectorize catchments drainage <- read_stars(file.path(wbt_wd, paste0(id, "_catchment.tif"))) drainage_polygon <- st_as_sf(drainage, as_points = FALSE, merge = TRUE, na.rm = TRUE) catchment_geometry <- st_union(st_geometry(drainage_polygon)) |> st_cast("POLYGON") # Remove isolated raster fragments and retain the main catchment. catchment_geometry <- catchment_geometry[which.max(st_area(catchment_geometry))] catchments <- rbind(catchments, st_sf(id = id, geometry = catchment_geometry)) } ``` ```{r, echo=FALSE, message=FALSE, warning=TRUE, eval=running, results='hide'} try_chunk <- try({ <> }, silent = TRUE) if(inherits(try_chunk, "try-error")){ warning("\nIssue when delineating catchments. \nThe vignette will not be fully built.") running <- FALSE } ``` The resulting catchment delineations can be plotted and checked. ```{r, echo=TRUE, message=FALSE, warning=FALSE, eval=running} # Compare drainage areas to the dataset provided with transfR compare_areas <- data.frame( name = names(Blavet$hl), expected_area = st_area(st_geometry(Blavet$obs)) |> units::set_units("km^2") |> round(1), computed_area = st_area(catchments) |> units::set_units("km^2") |> round(1) ) print(compare_areas) ``` ```{r, echo=TRUE, message=FALSE, warning=FALSE, eval=running, results='hide', fig.width=7, fig.height=4} # Plot catchment delineation par(oma = c(0, 0, 0, 6)) plot(catchments[,"id"], main = "Catchment delineations", key.pos = 4, reset = FALSE) plot(st_geometry(st_intersection(network_topage,catchments)), col = "white", lwd = 1.5, add = TRUE) plot(outlets_snapped, col = "black", pch = 16, add = TRUE) ``` ## 4. Compute flow path length and hydraulic length for each catchment The `whitebox` package does not provide a function to compute hydraulic length directly. However, it can compute the total [downslope flow path length to the outlet](https://whiteboxr.gishub.org/reference/wbt_downslope_flowpath_length.html) using the `whitebox::wbt_downslope_flowpath_length()` function and the [downslope distance to the stream](https://whiteboxr.gishub.org/reference/wbt_downslope_distance_to_stream.html) using the `whitebox::wbt_downslope_distance_to_stream()` function. Hydraulic length does not include the flow path over hillslopes, so it can simply be calculated as the difference between these two flow distances. Flow path lengths are computed once across the entire study area. For each catchment, this raster is then trimmed (i.e. NA values outside the catchment bounding box are removed), and the flow path lengths are adjusted so that the minimum is equal to half the pixel width or height. Finally, the hydraulic lengths are gathered in a list as expected by the `as_transfr()` function. ```{r, hydraulic_length, echo=TRUE, message=FALSE, warning=FALSE, eval=FALSE, results='hide'} # Compute hydraulic length whitebox::wbt_downslope_flowpath_length(d8_pntr = "d8.tif", output = "fpl.tif", wd = wbt_wd) whitebox::wbt_downslope_distance_to_stream(dem = "dem_fill.tif", streams = "network_topage.tif", output = "d2s.tif", wd = wbt_wd) fpl <- read_stars(file.path(wbt_wd, "fpl.tif")) d2s <- read_stars(file.path(wbt_wd, "d2s.tif")) hl_region <- fpl-d2s names(hl_region) <- "hl" # Crop hydraulic length for each catchment hl <- list() for(id in catchments$id){ crop <- st_crop(hl_region, catchments[catchments$id == id, ]) crop <- crop-min(crop$hl, na.rm = TRUE) crop$hl <- units::set_units(crop$hl, "m") hl[[id]] <- crop } ``` ```{r, echo=FALSE, message=FALSE, warning=TRUE, eval=running, results='hide'} try_chunk <- try({ <> }, silent = TRUE) if(inherits(try_chunk, "try-error")){ warning("\nIssue when computing hydraulic length. \nThe vignette will not be fully built.") running <- FALSE } ``` The resulting hydraulic length maps can be plotted and checked. ```{r, echo=TRUE, message=FALSE, warning=FALSE, eval=running, results='hide', fig.width=7, fig.height=4} i <- 1 network <- st_geometry(st_intersection(network_topage, catchments[i, ])) plot(hl[[i]], main = paste("Hydraulic length of catchment", i,"[m]"), col = hcl.colors(20, palette = "Teal"), key.pos = 1, reset = FALSE) plot(network, col = "white", lwd = 1.5, add = TRUE) ``` ## 5. Creating a transfR object and running a simulation Catchment delineations can be used as the spatial dimension of `stars` objects to georeference the observed flow time series of gauged catchments and locate ungauged catchments (see vignette *Preparation of input data: creation of a stars object* for details). Here we will create a `stars` object using the observed discharge of the `Blavet` dataset and the delineations that we just computed with `whitebox`. ```{r, echo=TRUE, message=FALSE, warning=FALSE, eval=running, results='hide'} obs_st <- st_as_stars(list(Qobs = Blavet$obs$Qobs), dimensions = st_dimensions( time = st_get_dimension_values(Blavet$obs,1), space = st_geometry(catchments))) ``` The hydraulic length maps can then be used to create an object of class `transfR` with the `as_transfr()` function (argument `hl`) and to perform simulations. ```{r, echo=TRUE, message=FALSE, warning=FALSE, eval=running, results='hide'} obs <- as_transfr(st = obs_st, hl = hl) ``` Hydrographs can then be transferred from the gauged catchments to the ungauged catchments using the `quick_transfr()` function. ```{r, echo=TRUE, message=FALSE, warning=FALSE, eval=FALSE} obs <- quick_transfr(obs, velocity = "brittany2013", cv = TRUE, parallel = TRUE, cores = 2) ``` The simulated time series will be available in its stars object as new attributes. ```{r, echo=TRUE, message=FALSE, warning=FALSE, eval=FALSE} obs$st ``` ## References
```{r, echo=FALSE, message=FALSE, warning=FALSE, eval=TRUE, results='hide'} # Cleaning temporary directory unlink(wbt_wd, recursive = TRUE) ```