---
title: "V03 - Usages data preprocessing"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{V03 - Usages data preprocessing}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, 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 = "70%",
  fig.align = "center"
  # <<< Graphic settings
)
```

```{r setup}
library(tidyverse)
library(airGRccia)
```

## Exploring available data on hubeau

Documentation:

-  https://inrae.github.io/hubeau/reference/get_prelevements.html for
R functions retrieving data on hubeau API.
- https://hubeau.eaufrance.fr/page/api-prelevements-eau for available
filters and variables available in the database

```{r}
hubeau::list_endpoints("prelevements")
hubeau::list_params("prelevements", "chroniques")
```

The "bbbox" parameter is convenient for retrieving data from GPS coordinates (WGS84).

From hubeau's documentation:

> Rectangle d'emprise de l'objet demandé, emprise au format : min longitude,
> min latitude, max longitude, max latitude avec les coordonnées en WGS84
> (EPSG:4326), le point doit être utilisé comme séparateur décimal, exemple :
> 1.6194,47.7965,2.1910,47.9988

## Search for withdrawals data in the Herault catchment area

Get the bounding box of the catchment area

```{r}
sf_sb <- readRDS(getCachePath("V01", "sb", "RDS"))
sf_sb
```

Coordinate system used here is RGF93. We need to convert bbox to WGS84.

```{r}
sf_sb <- sf::st_transform(sf_sb, crs = "WGS84")
bbox <- sf::st_bbox(sf_sb)
bbox
```

Retrieving data

```{r}
dfWD <- hubeau::get_prelevements_chroniques(bbox = paste(round(bbox, 4), collapse = ","))
dfWD <- dfWD |> filter(volume > 0, annee >= 2012)
str(dfWD)
```

```{r, out.width="100%"}
sf_WD <- sf::st_as_sf(dfWD, crs = "WGS84", coords = c("geometry.coordinates1", "geometry.coordinates2"))
library(tmap)
tmap_mode("view")
plot_WD_map <- function(selected_year = 2021) {
  tm_basemap("OpenStreetMap.France") +
    tm_shape(sf_sb) + tm_borders("black", lwd = 2) +
    tm_shape(sf_WD |> filter(annee == selected_year)) +
    tm_symbols(col = "libelle_usage", size = "volume", scale = 5, popup.vars = c(
      "libelle_usage", "code_ouvrage", "nom_ouvrage", "volume"
    ))
}
plot_WD_map()
```

## Cleaning and validation of yearly data

Lets remove Water flowing through turbines of run of river hydro facilities and
only keep withdrawals located in the catchment area.

```{r, out.width="100%"}
sf_WD <- sf_WD |> filter(code_usage != "BAR")
sf_WD <- sf::st_intersection(sf_WD, sf_sb |> select(id))
plot_WD_map()
```



```{r}
lapply(setNames(nm = unique(sf_WD$libelle_usage)), function(x) {
  unique(sf_WD |> filter(libelle_usage == x) |> pull(nom_ouvrage))
})
```

Volumes registered in "CANAUX" correspond to the withdrawn at the irrigation canal
intake.
Volumes registered in "IRRIGATION" correspond to the amount of water used from
the canal for irrigating the fields, which means that these volumes are counted
twice.
The difference of volume between "CANAUX" and "IRRIGATION" correspond to the amount
of water that has been released in the river all along the irrigation system.

```{r}
knitr::include_graphics(
  rprojroot::find_package_root_file("vignettes/scheme_water_fluxes_irrigation_system.png")
)
```

As first approximation, we consider that the impact on the river correspond to
only water consumption, so we are not considering "CANAUX" withdrawn.

```{r}
sf_WD <- sf_WD |> filter(code_usage != "CAN")
```

```{r}
plot_withdraw_year <- function(df) {
  ggplot(
    df |> group_by(annee, code_usage) |> summarise(total_volume = sum(volume)),
    aes(x = annee, y = total_volume, fill = code_usage)
  ) +
    geom_col() +
    theme(legend.position = "bottom")
}
plot_withdraw_year(sf_WD)
```

The calibration period begins in August 2013 and ends at the 31th December 2024.
Years after `r max(sf_WD$annee)` are missing. By default, we complete the gap with the data
of 2021.

```{r}
# Reference year to replicate is the last year with data
ref_year <- max(sf_WD$annee)

# Target years to fill
target_years <- (ref_year + 1):2024

# Generate rows for missing years by replicating 2021
sf_WD <- bind_rows(
  sf_WD,
  purrr::map_dfr(target_years, ~ sf_WD %>%
            filter(annee == ref_year) %>%
            mutate(annee = .x))
)
plot_withdraw_year(sf_WD)
```


And then we apply a a water consumption ratio to all raw withdrawals in order to
get the real amount of water removed from the river.

Aater consumption ratio are extracted from:

Eric Sauquet, Jean-Pierre Vergnes, Guillaume Thirel, Laurent Strohmenger.
Identification de stations hydrométriques et piézomètres pour l'évaluation des
modélisations hydrologiques et hydrogéologiques. INRAE; BRGM. 2022. ([hal-03940233](https://hal.inrae.fr/hal-03940233))
https://professionnels.ofb.fr/sites/default/files/pdf/projets/Explore2/Rapport%20Explore2%20reseauReference%20Evaluation%20-%20VF.pdf

```{r}
cons_ratio <- c(AEP = 0.2, IRR = 0.7, IND = 0.07)
sf_WD$cons_ratio <- cons_ratio[sf_WD$code_usage]
sf_WD$raw_volume <- sf_WD$volume
sf_WD$volume <- sf_WD$raw_volume * sf_WD$cons_ratio
plot_withdraw_year(sf_WD)
```

We aggregate the data at sub-basin scale:

```{r}
wd_yearly <-  as.data.frame(sf_WD) |>
  group_by(id, annee, code_usage) |>
  summarise(volume = sum(volume), .groups = 'drop')
head(wd_yearly)
```

## Usage data deaggregation at daily time step

Sauquet et al. (2022) consider a constant withdrawal throughout the calendar
year for Drinking water and industry and a concentration of the irrigation between
June and September.

```{r}
code_usages <- c("AEP", "IND", "IRR")
monthly_ratio <- data.frame(
  code_usage = rep(code_usages, each = 12),
  month = rep(1:12, 3),
  ratio = c(rep(1/12, 2 * 12), rep(0, 5), rep(1/4, 4), rep(0, 3))
)
monthly_ratio
```

More precise monthly distribution can be taken into account given other sources,
notably on https://fleuve-herault.fr website.

```{r}
get_daily_wd <- function(yearly) {
  # Create vector of daily dates for the current year
  bound_dates <- as.POSIXlt(c("2000-01-01", "2000-12-31"))
  lubridate::year(bound_dates) <- yearly$annee
  DatesR <- seq(bound_dates[1], bound_dates[2], "1 day")
  # Preparation of the table
  daily <- data.frame(DatesR = DatesR,
                      month = lubridate::month(DatesR),
                      id = yearly$id,
                      code_usage = yearly$code_usage)
  # Computation of the ratio to apply by day
  day_ratio <- 1 / 365.25 * 12 # Average ratio for a monthly volume
  month_ratio <- monthly_ratio |>
    filter(code_usage == yearly$code_usage) |>
    arrange(month) |>
    pull(ratio)
  daily$day_ratio <- day_ratio * month_ratio[daily$month]
  # Computation of the net withdrawal
  daily$volume <- yearly$volume * daily$day_ratio
  return(daily)
}

l_daily <- lapply(seq(nrow(wd_yearly)), function(i) {
  get_daily_wd(yearly = wd_yearly[i, , drop = FALSE])
})
wd_daily <- do.call(rbind, l_daily)
summary(wd_daily)
```

We aggregate the data and format it for airGRiwrm:

```{r}
wd_daily <- wd_daily |> group_by(DatesR, id) |> summarise(volume = sum(volume), .groups = "drop")
Qinf <- tidyr::pivot_wider(wd_daily, names_from = id, values_from = volume, values_fill = 0)
summary(Qinf)
```

```{r, out.width="100%"}
TSstudio::ts_plot(Qinf, slider = TRUE)
```

```{r}
file_Qinf <- getCachePath("V03", "Qinf_historical", ".xlsx")
if(!file.exists(file_Qinf)) {
  openxlsx2::write_xlsx(Qinf, file_Qinf)
}
```