V01 - Watershed Delineation

Watershed Delineation for Semi-Distributed Hydrological Modeling

This vignette demonstrates how to delineate a watershed area for a semi-distributed hydrological model. We will use the airGRccia package along with other spatial data processing tools to achieve this.

Prerequisites

Before starting, ensure you have the following packages installed (See README):

library(tidyverse)
#> ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
#> ✔ dplyr     1.2.1     ✔ readr     2.2.0
#> ✔ forcats   1.0.1     ✔ stringr   1.6.0
#> ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
#> ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
#> ✔ purrr     1.2.2     
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag()    masks stats::lag()
#> ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)
#> Linking to GEOS 3.14.1, GDAL 3.12.2, PROJ 9.7.1; sf_use_s2() is TRUE
library(terra)
#> terra 1.9.34
#> 
#> Attaching package: 'terra'
#> 
#> The following object is masked from 'package:tidyr':
#> 
#>     extract
library(tmap)
library(airGRccia)
#> Loading required package: airGRiwrm
#> Loading required package: airGR
#> 
#> Attaching package: 'airGRiwrm'
#> 
#> The following objects are masked from 'package:airGR':
#> 
#>     Calibration, CreateCalibOptions, CreateInputsCrit,
#>     CreateInputsModel, CreateRunOptions, RunModel

Data Setup

Obtaining the Digital Elevation Model (DEM)

First, we need to define the bounding box of our watershed area. We can obtain the coordinates from Google Earth or Geoportail.

# Define the bounding box coordinates in RGF93
bbox <- sf::st_bbox(
  c(
    xmin = 716000,
    ymin = 6238000,
    xmax = 774000,
    ymax = 6355000
  ),
  crs = 2154
)

# Create an sf object from the bounding box
bbox_sf <- sf::st_as_sfc(bbox) %>% st_as_sf()

tmap_mode("view")
#> ℹ tmap modes "plot" - "view"
#> ℹ toggle with `tmap::ttm()`
tm_basemap("OpenTopoMap") +
  tm_shape(bbox_sf) +
  tm_polygons(lwd = 3, col = "red", fill = NA, fill_alpha = 0)

Next, we download the DEM for the defined bounding box.

dem <- get_ign_dem(bbox)
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
#> Warp executed successfully.
plot(dem)

Obtaining the Stream Network

We download the stream network from the IGN API.

# Define cache path for the stream network
streams <- get_ign_streams(bbox)
#> iterating ■■■                                7% | ETA: 17s
#> iterating ■■■■■                             15% | ETA: 19s
#> iterating ■■■■■■■■                          24% | ETA: 19s
#> iterating ■■■■■■■■■■■■                      37% | ETA: 16s
#> iterating ■■■■■■■■■■■■■■■                   46% | ETA: 14s
#> iterating ■■■■■■■■■■■■■■■■■■■               59% | ETA: 11s
#> iterating ■■■■■■■■■■■■■■■■■■■■■■■■          78% | ETA:  5s
#> iterating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s

# Plot (not in "view" mode because it is too large for a vignette)
original_mode <- tmap_mode("plot")
#> ℹ tmap modes "plot" - "view"
tm_basemap("OpenStreetMap.France") +
  tm_shape(streams) +
  tm_lines(col = "nature")

tmap_mode(original_mode)
#> ℹ tmap modes "plot" - "view"

Points of Interest Locations

The final goal is to create the GRiwrm object which describes the network for the airGRiwrm hydrological semi-distributed model. The semi-distributed model cuts the watershed into sub-basins delineated by Points of Interest (POIs).

If we refer to createGRiwrm function, we need to get: - POIs list sorted in upstream-downstream order - Area of each basin/sub-basin - Hydraulic distance between POIs

Gauging stations are the first important POIs to consider. They can be selected from: - Search engine on Hydroportail: https://www.hydro.eaufrance.fr/rechercher/entites-hydrometriques - Hubeau’s hydrometry viewer: https://hubeau.eaufrance.fr/sites/default/files/api/demo/hydro_tr.htm

After searching on the River Herault and its affluents, the following gauging sites can be selected:

stations <- hubeau::get_hydrometrie_sites(
  code_site = c(
    "Y2140010",
    "Y2140020",
    "Y2230010",
    "Y2210010",
    "Y2300020",
    "Y2100020"
  )
)
stations
#> # A tibble: 6 × 39
#>   code_site libelle_site           type_site coordonnee_x_site coordonnee_y_site
#>   <chr>     <chr>                  <chr>                 <dbl>             <dbl>
#> 1 Y2100020  L'Hérault à Laroque    STANDARD            759105           6313209 
#> 2 Y2140010  L'Hérault à Gignac     STANDARD            743142           6283755 
#> 3 Y2140020  L'Hérault à Saint-Gui… STANDARD            748274.          6294640.
#> 4 Y2210010  La Lergue à Lodève     STANDARD            726110           6292216 
#> 5 Y2230010  Le Salagou à Clermont… STANDARD            733408           6285655 
#> 6 Y2300020  L'Hérault à Aspiran    STANDARD            737989           6274704 
#> # ℹ 34 more variables: code_projection <int>, longitude_site <dbl>,
#> #   latitude_site <dbl>, altitude_site <dbl>, code_systeme_alti_site <int>,
#> #   surface_bv <dbl>, statut_site <int>, premier_mois_etiage_site <int>,
#> #   premier_mois_annee_hydro_site <int>, influence_generale_site <int>,
#> #   code_entite_hydro_site <chr>, code_troncon_hydro_site <chr>,
#> #   code_commune_site <chr>, code_zone_hydro_site <chr>, libelle_commune <chr>,
#> #   code_departement <chr>, code_region <chr>, libelle_region <chr>, …

Sub-Basins Delineation

The airGRciaa package contains functions that automate the watershed delineation and flow path calculation using the library WhiteboxTools.

# Define cache path for Sub-basin polygons
file_sb <- getCachePath("V01", "sb", "RDS")

# Format the POIs location for `extract_sub_basins`
# Option 1: Using a data.frame (original approach)
outlets <- stations[, c(
  "code_site",
  "coordonnee_x_site",
  "coordonnee_y_site",
  "surface_bv"
)]
colnames(outlets) <- c("id", "x", "y", "area")
outlets$area <- outlets$area * 1E6 # Convert from km² to m²

# Option 2: Using an sf object (coordinates extracted from geometry)
# outlets_sf <- sf::st_as_sf(
#   stations[, c("code_site", "coordonnee_x_site", "coordonnee_y_site", "surface_bv")],
#   coords = c("coordonnee_x_site", "coordonnee_y_site"),
#   crs = 2154
# )
# colnames(outlets_sf)[1] <- "id"
# outlets_sf$area <- outlets_sf$surface_bv * 1E6  # Convert from km² to m²

# Areas from hydroportail are wrong in this area
# E.g.: https://hydrogr.github.io/BDD-HydroClim/images/fr/Y210002001_INRAE_BDD-HydroClim_fact_sheet_FR.png
# or worse: https://hydrogr.github.io/BDD-HydroClim/images/fr/Y221001001_INRAE_BDD-HydroClim_fact_sheet_FR.png
# So we set area_tolerance = 0.2 for avoiding changing outlet position for a wrong reason
sb <- extract_sub_basins(
  outlets = outlets,
  dem = dem,
  streams = streams,
  filename = file_sb,
  area_tolerance = 0.2
)
#> Performing one-time download of WhiteboxTools binary from
#>   https://www.whiteboxgeo.com/WBT_Linux/WhiteboxTools_linux_amd64.zip 
#> (This could take a few minutes, please be patient...)
#> WhiteboxTools binary is located here:  /github/home/.local/share/R/whitebox/WBT/whitebox_tools 
#> You can now start using whitebox
#>     library(whitebox)
#>     wbt_version()
#> Writing layer `streams_16bfea658b4c260b3707df96d3034ab1' to data source 
#>   `/github/home/.cache/airGRccia/streams_16bfea658b4c260b3707df96d3034ab1.shp' using driver `ESRI Shapefile'
#> Writing 20133 features with 1 fields and geometry type Line String.
sb
#> Simple feature collection with 6 features and 10 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 718450 ymin: 6273325 xmax: 770000 ymax: 6336175
#> Projected CRS: RGF93 v1 / Lambert-93
#> # A tibble: 6 × 11
#>   id            area mvd_outlet_x mvd_outlet_y apriori_area   tot_area
#> * <chr>        [m^2]        <dbl>        <dbl>        <dbl>      [m^2]
#> 1 Y2230010  78768750      733388.     6285662.           NA   78768750
#> 2 Y2210010 182417500      725862.     6291762.    228000000  182417500
#> 3 Y2100020 910304375      759112.     6313212.    760000000  910304375
#> 4 Y2140020 316295625      748262.     6294638.   1200000000 1226600000
#> 5 Y2140010 180587500      743138.     6283762.   1312000000 1407188750
#> 6 Y2300020 286473750      738012.     6274712.           NA 1954851875
#> # ℹ 5 more variables: geometry <POLYGON [m]>, outlet_x <dbl>, outlet_y <dbl>,
#> #   mvd_dist <dbl>, down <chr>

plot(terra::vect(sb), c("id", "down"))

Flow Path Calculation

The flow path calculation requires the sub-basin object produced by extract_sub_basins.

reaches <- extract_reaches(sb = sb)
#> Reading layer `streams_16bfea658b4c260b3707df96d3034ab1' from data source 
#>   `/github/home/.cache/airGRccia/streams_16bfea658b4c260b3707df96d3034ab1.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 20133 features and 1 field
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 714437 ymin: 6236544 xmax: 778416.6 ymax: 6356142
#> Projected CRS: RGF93 v1 / Lambert-93
#> Rasterization of streams...
#> Buffering around streams step 1/2...
#> Buffering around streams step 2/2...
#> Running extract_reach on 3 cores... This can take a while.
#> Search From Y2230010 to Y2300020
#> Reach Upstream (2774, 696) - Downstream (3212, 881). Fly distance = 475 cells
#> Found downstream outlet with dist = 0 | min dist = 0 | reach length = 791 (unit: cells)
#> Search From Y2210010 to Y2300020
#> Reach Upstream (2530, 395) - Downstream (3212, 881). Fly distance = 837 cells
#> Found downstream outlet with dist = 0 | min dist = 0 | reach length = 1273 (unit: cells)
#> Search From Y2100020 to Y2140020
#> Reach Upstream (1672, 1725) - Downstream (2415, 1291). Fly distance = 860 cells
#> Found downstream outlet with dist = 0 | min dist = 0 | reach length = 1131 (unit: cells)
#> Search From Y2140020 to Y2140010
#> Reach Upstream (2415, 1291) - Downstream (2850, 1086). Fly distance = 481 cells
#> Found downstream outlet with dist = 0 | min dist = 0 | reach length = 684 (unit: cells)
#> Search From Y2140010 to Y2300020
#> Reach Upstream (2850, 1086) - Downstream (3212, 881). Fly distance = 416 cells
#> Found downstream outlet with dist = 0 | min dist = 0 | reach length = 515 (unit: cells)
#> Writing layer `reaches_ddbe1cd2ea7a8763f3f29772e584b83e' to data source 
#>   `/github/home/.cache/airGRccia/reaches_ddbe1cd2ea7a8763f3f29772e584b83e.gpkg' using driver `GPKG'
#> Writing 5 features with 5 fields and geometry type Line String.

One can then map the details of the sub-basin delineations with flow path between outlets:

moved_outlets <- sf::st_as_sf(
  sf::st_drop_geometry(sb)[, c("id", "mvd_outlet_x", "mvd_outlet_y")],
  coords = c("mvd_outlet_x", "mvd_outlet_y"),
  crs = sf::st_crs(sb)
)

# Temporary fix due to https://github.com/r-tmap/tmap/issues/1259
sb_plot <- sb
sb_plot$apriori_area <- NULL
tm_basemap("OpenStreetMap.France") +
  tm_shape(sb_plot) +
  tm_polygons(fill = "id", fill_alpha = 0.3) +
  tm_shape(reaches) +
  tm_lines(col = "blue", lwd = 3) +
  tm_shape(moved_outlets) +
  tm_symbols(fill = "id")

Hypsometric analysis

# Compute hypsometric curve for each sub-basin from DEM
file_hypso <- getCachePath("V01", "hypso", "xlsx")
if (!file.exists(file_hypso)) {
  HypsoData <- calc_hypso(sb, dem)
  openxlsx2::write_xlsx(HypsoData, file_hypso)
} else {
  HypsoData <- openxlsx2::read_xlsx(file_hypso)
}
str(HypsoData)
#>  num [1:101, 1:6] 130 170 194 212 227 ...
#>  - attr(*, "dimnames")=List of 2
#>   ..$ : chr [1:101] "0%" "1%" "2%" "3%" ...
#>   ..$ : chr [1:6] "Y2100020" "Y2140010" "Y2140020" "Y2210010" ...

Creating Candidate Data Frame for GRiwrm Network

nodes <- sf::st_drop_geometry(sb)

# Add distance between node and downstream node
dist <- c(units::drop_units(sf::st_length(reaches)) / 1000, NA_real_)
names(dist) <- c(reaches$id, nodes$id[is.na(nodes$down)])
nodes$length <- dist[nodes$id]

libelle <- setNames(object = stations$libelle_site, nm = stations$code_site)
nodes$libelle_site <- libelle[nodes$id]
nodes$tot_area <- units::drop_units(units::set_units(nodes$tot_area, km^2))
nodes$area <- units::drop_units(units::set_units(nodes$area, km^2))
nodes
#> # A tibble: 6 × 12
#>   id      area mvd_outlet_x mvd_outlet_y apriori_area tot_area outlet_x outlet_y
#>   <chr>  <dbl>        <dbl>        <dbl>        <dbl>    <dbl>    <dbl>    <dbl>
#> 1 Y2230…  78.8      733388.     6285662.           NA     78.8  733408  6285655 
#> 2 Y2210… 182.       725862.     6291762.    228000000    182.   726110  6292216 
#> 3 Y2100… 910.       759112.     6313212.    760000000    910.   759105  6313209 
#> 4 Y2140… 316.       748262.     6294638.   1200000000   1227.   748274. 6294640.
#> 5 Y2140… 181.       743138.     6283762.   1312000000   1407.   743142  6283755 
#> 6 Y2300… 286.       738012.     6274712.           NA   1955.   737989  6274704 
#> # ℹ 4 more variables: mvd_dist <dbl>, down <chr>, length <dbl>,
#> #   libelle_site <chr>

Saving the Object for Next Vignettes

file_network <- getCachePath("V01", "network", "xlsx")

if (!file.exists(file_network)) {
  openxlsx2::write_xlsx(nodes, file_network)
}

Conclusion

In this vignette, we have demonstrated how to delineate a watershed area for a semi-distributed hydrological model. We encourage you to experiment with different parameters and data sources to see how they affect the results.