inst/doc/V03_inputs_preparation_whitebox.R

In transfR: Transfer of Hydrograph from Gauged to Ungauged Catchments

## ----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)