--- title: "V01 - Watershed Delineation" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{V01 - Watershed Delineation} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r knitr_setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", # >>> Graphic settings # from https://r4ds.had.co.nz/graphics-for-communication.html?q=fig.asp#figure-sizing fig.width = 6, fig.asp = 0.618, out.width = "100%", fig.align = "center" # <<< Graphic settings ) ``` # 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): ```{r setup} library(tidyverse) library(sf) library(terra) library(tmap) library(airGRccia) ``` ## 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*. ```{r bounding_box, out.width="100%", fig.asp=1} # 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") 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. ```{r dem, fig.asp=1.2, out.width = "60%"} dem <- get_ign_dem(bbox) plot(dem) ``` ### Obtaining the Stream Network We download the stream network from the *IGN API*. ```{r stream_network, fig.asp = 1} # Define cache path for the stream network streams <- get_ign_streams(bbox) # Plot (not in "view" mode because it is too large for a vignette) original_mode <- tmap_mode("plot") tm_basemap("OpenStreetMap.France") + tm_shape(streams) + tm_lines(col = "nature") tmap_mode(original_mode) ``` ### 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: ```{r stations} stations <- hubeau::get_hydrometrie_sites( code_site = c( "Y2140010", "Y2140020", "Y2230010", "Y2210010", "Y2300020", "Y2100020" ) ) stations ``` ## Sub-Basins Delineation The `airGRciaa` package contains functions that automate the watershed delineation and flow path calculation using the library *WhiteboxTools*. ```{r sub_basins, fig.width=9, out.width="100%"} # 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 ) sb plot(terra::vect(sb), c("id", "down")) ``` ## Flow Path Calculation The flow path calculation requires the sub-basin object produced by `extract_sub_basins`. ```{r reaches} reaches <- extract_reaches(sb = sb) ``` One can then map the details of the sub-basin delineations with flow path between outlets: ```{r flow_path_map, fig.asp=1.1} 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 ```{r hypsometric} # 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) ``` ## Creating Candidate Data Frame for GRiwrm Network ```{r network_nodes} 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 ``` ## Saving the Object for Next Vignettes ```{r save_network} 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.