---
title: "V02 - Hydroclimatic data preprocessing"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{V02 - Hydroclimatic 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)
library(stars)
```

The final goal of this vignette is to create a *BasinsObs* object containing
both hydrometric and climatic observations.

## Retrieving and format climatic data

### Climatic data: French SAFRAN reanalysis

First, we retrieve climatic data from the SAFRAN reanalysis produced by Meteo France (todo: add some details about this database),
on the study area.
We proceed by querying the RADIS API with the function
[RADIS package](https://umr-g-eau.pages-forge.inrae.fr/radis/) developed by
the G-EAU Mixte Research Unit.

```{r}
# Load sub-basins shapefile for getting the extent of the data retrieval
sf_sb <- readRDS(getCachePath("V01", "sb", "RDS"))

# Get the SAFRAN stars object with RADIS
s_safran <- get_climate_safran(sf_sb,
                               date_start = "2012-08-01",
                               date_end = "2024-12-31",
                               overlay_mode = "extent")
```

The `s_safran` object is a `stars` object containing the climatic variables
in grid format:

```{r, fig.asp=1}
plot_climate_bvi <- function(s, sf_sb, variable = "TempMean", i_time, limits = c(-20, 40)) {
  ggplot() +
    geom_stars(data = s[variable, , , i_time]) +
    geom_sf(
      data = sf_sb,
      fill = NA,
      color = "grey20",
      linewidth = 1
    ) +
    #coord_equal() +
    coord_sf(datum = sf::st_crs(s)) +
    geom_sf_text(data = sf_sb, aes(label = id), size = 1.5) +
    scale_fill_gradient2(
      limits = limits,
      low = scales::muted("blue"),
      high = scales::muted("red"),
      mid  = "white"
    ) +
    ggtitle(sprintf(
      "Variable '%s' at date %s",
      variable,
      format(stars::st_get_dimension_values(s_safran, "time")[i_time])
    ))
}

# Every layer should be in same CRS in order to be plotted together
plot_climate_bvi(s_safran, sf_sb, i_time = 1)
```

Semi-distrbuted hydrologixcal models such as airGRiwrm require mean climate data
at sub-basin scale.
Such data are available using the agrument `overlay_mode = "aggregate"` when
retrieving the data:

```{r, out.width="100%"}
s_sb <- get_climate_safran(
  sf_sb,
  date_start = "2012-08-01",
  date_end = "2024-12-31",
  overlay_mode = "aggregate"
)
plot(s_sb["TempMean", , 1:6])
```


## Creation of the dataset for hydrometric data

We retrieve observed flow time series from site codes of 8 character-lengths used
in the sub-basin delineation.

```{r}
Qobs <- get_hubeau_flows(
  x = sf_sb$id,  # Use sub-basin ids
  date_start = "2012-08-01",
  date_end = "2024-12-31"
)
summary(Qobs)
```

We get the observed flows by station code (10 character-lengths) and there
are potentially more than one station by site.

Nota Bene: no observed flows are available for the gauging site Y2230010
(Le Salagou à Clermont-l'Hérault [Mas Audran]).

We can plot the data availability for each station:

```{r}
# Plot data availability for each station
ggplot(
  Qobs %>% tidyr::pivot_longer(-DatesR, names_to = "code_station", values_to = "Q") %>% filter(!is.na(Q)),
  aes(x = as.POSIXct(DatesR), y = code_station)
) +
  geom_point() +
  ggtitle("Gauging station data availabilities")
```

As station `Y214001001` has definitely been replaced by station `Y214001002` in 2020
we remove the latter for this site.

```{r}
# Remove
Qobs <- Qobs |> select(-Y214001001)
```

Then We can convert station code into site codes to be compliant with the code
used with sub-basins.

```{r}
names(Qobs) <- substr(names(Qobs), 1, 8)  # Keep only the first 8 characters
```

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


## Create the BasinsObs object

```{r}
BasinsObs <- CreateBasinsObs(s_sb, ids = sf_sb$id, Qobs = Qobs)
str(BasinsObs)
```

Save the BasinsObs object on the cloud

```{r}
file_BasinsObs <- getCachePath("V02", "BasinsObs_safran_historical", "RDS")
if (!file.exists(file_BasinsObs)) saveRDS(BasinsObs, file_BasinsObs)
```