---
title: "V05 - Hydrologic model calibration"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{V05 - Hydrologic model calibration}
  %\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)
```

## Network definition: GRiwrm object

Load network from V01:

```{r}
nodes <- openxlsx2::read_xlsx(getCachePath("V01", "network", "xlsx"))
# Remove unecessary columns
nodes <- nodes |>
  select(id, down, length, tot_area, libelle_site, outlet_x, outlet_y)
```

Adding withdrawals as Direct Injection nodes:

```{r}
nodes <- rbind(
  nodes,
  data.frame(
    id = paste0(nodes$id, "_WD"),
    down = nodes$id,
    length = 0,
    tot_area = NA_real_,
    libelle_site = paste("Withdrawals on", nodes$libelle_site),
    outlet_x = nodes$outlet_x,
    outlet_y = nodes$outlet_y
  )
)
nodes
```


Load the observations:

```{r}
BasinsObs <- readRDS(getCachePath("V02", "BasinsObs_safran_historical", "RDS"))
str(BasinsObs)
```

No observed flow is available at station Y2230010 which corresponds to location
of the Salagou dam.

So we need to:

- Add a reservoir in the model corresponding to the Salagou dam
- model the reservoir inflows using an ungauged model with parameters
transferred from a nearby basin.

```{r}
# Add the dam
nodes <- rbind(
  nodes,
  data.frame(
    id = "DAM",
    down = "Y2300020",
    length = 48.50,
    tot_area = as.numeric(NA),
    libelle_site = "Barrage du Salagou",
    outlet_x = nodes$outlet_x[nodes$id == "Y2230010"],
    outlet_y = nodes$outlet_y[nodes$id == "Y2230010"]
  )
)
# Adjust Dam inflow basin
nodes$down[nodes$id == "Y2230010"] <- "DAM"
nodes$length[nodes$id == "Y2230010"] <- 0
```

Choice of the models used:

```{r}
nodes$model <- "RunModel_GR4J"
nodes$model[nodes$id == "Y2230010"] <- "Ungauged"
nodes$model[nodes$id == "DAM"] <- "RunModel_Reservoir"
nodes$model[grep("_WD", nodes$id)] <- NA_character_
```

Definition of the basin given its parameters for the ungauged basin Y2230010.
The basin "La Lergue à Lodève (Y2210010)" is the closest one.

```{r}
nodes$donor <- as.character(NA)
nodes$donor[nodes$id == "Y2230010"] <- "Y2210010"
```

Creation of the GRiwrm object:

```{r}
griwrm <- CreateGRiwrm(nodes, cols = list(area = "tot_area"))
plot(griwrm)
```



## Inputs of the model: `createInputsModel`

Format all the inputs:

```{r}
Qinf <- openxlsx2::read_xlsx(getCachePath(
  "V03",
  "Qinf_historical",
  ".xlsx"
)) %>%
  mutate(DatesR = as.Date(DatesR))
Qinf <- Qinf %>%
  filter(DatesR >= min(BasinsObs$DatesR) & DatesR <= max(BasinsObs$DatesR))
names(Qinf)[-1] <- paste0(names(Qinf)[-1], "_WD")
str(Qinf)
```

```{r}
Qrelease <- openxlsx2::read_xlsx(getCachePath(
  "V04",
  "Qrelease_historical",
  "xlsx"
))
str(Qrelease)
```

Remove extra columns and convert flows into m<sup>3</sup>/day:

```{r}
Qinf$DatesR <- NULL
Qrelease$DatesR <- NULL
Qinf <- -Qinf # withdrawals are negative flows
Qrelease <- Qrelease * 86400 # m3/s -> m3/day
```

Call CreateInputsModel:

```{r}
IM <- CreateInputsModel(
  griwrm,
  DatesR = BasinsObs$DatesR,
  Precip = BasinsObs$Precip,
  PotEvap = BasinsObs$PotEvap,
  Qinf = Qinf,
  Qrelease = Qrelease
)
```

## `CreateRunoptions`

Definition of the calibration period

```{r}
Ind_Run <- seq(
  which(BasinsObs$DatesR == as.POSIXlt("2013-08-01", tz = "UTC")),
  which(BasinsObs$DatesR == as.POSIXlt("2023-07-31", tz = "UTC"))
)
```

```{r}
RO <- CreateRunOptions(IM, IndPeriod_Run = Ind_Run)
```

It seems that in 2012, the dam was barely full.

```{r}
RO[["DAM"]]$IniStates <- c("Reservoir.V" = 100E6)
```



## `CreateInputsCrit`

Observed flows need to be converted from m3/s to mm/day.

```{r}
Qobs <- BasinsObs$Qobs[Ind_Run, ]
areas <- sapply(names(Qobs), function(x) griwrm$area[griwrm$id == x])
k <- 86400 * 1E3 / (areas * 1E6) # seconds in a day * m3->L / area in m2 = L/m2/day = mm/day
Qobs <- t(replicate(nrow(Qobs), k)) * Qobs
summary(Qobs)
```

```{r}
IC <- CreateInputsCrit(
  InputsModel = IM,
  FUN_CRIT = ErrorCrit_KGE2,
  RunOptions = RO,
  Obs = Qobs
)
```


## `CreateCalibOptions`

```{r}
CO <- CreateCalibOptions(
  IM,
  FixedParam = list(DAM = c(Vmax = 102.2E6, celerity = 1))
)
```

## `Calibration`

```{r}
OC <- Calibration(
  InputsModel = IM,
  RunOptions = RO,
  InputsCrit = IC,
  CalibOptions = CO
)
```

```{r}
params <- extractParam(OC)
str(params)
```

# Check calibration result

```{r}
OM <- RunModel(
  IM,
  RunOptions = RO,
  Param = params
)
```

```{r, out.width="100%", fig.asp=1}
plot(OM, Qobs = Qobs[Ind_Run, ])
```