Nothing
## ----setup, include = FALSE---------------------------------------------------
library(transfR)
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
## ----echo=TRUE, message=FALSE, warning=FALSE, results='hide'------------------
wbt_wd <- tempdir(check = TRUE)
## ----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) # ~76m 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"))
## ----echo=FALSE, message=FALSE, warning=TRUE, eval=TRUE, results='hide'-------
try_chunk <- try({
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) # ~76m 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"))
}, 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
}
## ----install_whitebox, echo=FALSE, message=FALSE, warning=TRUE, eval=running, results='hide'----
# If WhiteboxTools executable are 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)
}
## ----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)
## ----echo=FALSE, message=FALSE, warning=TRUE, eval=running, results='hide'----
try_chunk <- try({
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)
}, 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
}
## ----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 km2) consistent with flow direction
# whitebox::wbt_extract_streams(flow_accum = "facc.tif",
# threshold = 100, # 100 cells for 1 km2
# 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)
## ----echo=FALSE, message=FALSE, warning=TRUE, eval=running, results='hide'----
try_chunk <- try({
# 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 km2) consistent with flow direction
whitebox::wbt_extract_streams(flow_accum = "facc.tif",
threshold = 100, # 100 cells for 1 km2
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)
}, 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
}
## ----delineate_catchments, echo=TRUE, message=FALSE, warning=FALSE, eval=FALSE, results='hide'----
# # Localize the outlets of the studied catchments (with manual adjustments to help 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_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")))
# contours <- st_contour(drainage, breaks = 1) |> st_geometry() |> st_cast("POLYGON")
# contours <- contours[which.max(st_area(contours))]
# catchments <- rbind(catchments, st_sf(data.frame(id, geom = contours)))
# }
## ----echo=FALSE, message=FALSE, warning=TRUE, eval=running, results='hide'----
try_chunk <- try({
# Localize the outlets of the studied catchments (with manual adjustments to help 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_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")))
contours <- st_contour(drainage, breaks = 1) |> st_geometry() |> st_cast("POLYGON")
contours <- contours[which.max(st_area(contours))]
catchments <- rbind(catchments, st_sf(data.frame(id, geom = contours)))
}
}, silent = TRUE)
if(inherits(try_chunk, "try-error")){
warning("\nIssue when delineating catchments. \nThe vignette will not be fully built.")
running <- FALSE
}
## ----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)
## ----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 = "Catchments delineation", 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)
## ----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
# }
#
## ----echo=FALSE, message=FALSE, warning=TRUE, eval=running, results='hide'----
try_chunk <- try({
# 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
}
}, silent = TRUE)
if(inherits(try_chunk, "try-error")){
warning("\nIssue when computing hydraulic length. \nThe vignette will not be fully built.")
running <- FALSE
}
## ----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)
## ----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)))
## ----echo=TRUE, message=FALSE, warning=FALSE, eval=running, results='hide'----
obs <- as_transfr(st = obs_st, hl = hl)
## ----echo=TRUE, message=FALSE, warning=FALSE, eval=FALSE----------------------
# obs <- quick_transfr(obs, velocity = "brittany2013", parallel = TRUE, cores = 2, cv = TRUE)
## ----echo=TRUE, message=FALSE, warning=FALSE, eval=FALSE----------------------
# obs$st
## ----echo=FALSE, message=FALSE, warning=FALSE, eval=TRUE, results='hide'------
# Cleaning temporary directory
unlink(wbt_wd, recursive = TRUE)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.