---
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.
