---
title: "Severn_02: Calibration of a GR4J semi-distributed model network"
author: "David Dorchies"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Severn_02: Calibration of a GR4J semi-distributed model network}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  comment = "#>",
  fig.width = 6,
  fig.asp = 0.68,
  out.width = "70%",
  fig.align = "center"
)
```

**airGRiwrm** automates the execution of **airGR** semi-distributed models. The steps for running or calibrating the model are the same as the ones of 'airGR'.

## Load library

```{r}
library(airGRiwrm)
```

## Preparation of function inputs

To run a model, as for **airGR**, the functions of the  **airGRiwrm** package (e.g. the models, calibration and criteria calculation functions) require data and options with specific formats.

To facilitate the use of the package, there are several functions dedicated to the creation of these objects:

  * `CreateInputsModel()`: prepares the inputs for the different hydrological models (times series of dates, precipitation, observed discharge, etc.)
  * `CreateRunOptions()`: prepares the options for the hydrological model run (warm up period, calibration period, etc.)
  * `CreateInputsCrit()`: prepares the options in order to compute the efficiency criterion (choice of the criterion, choice of the transformation on discharge: "log", "sqrt", etc.)
  * `CreateCalibOptions()`: prepares the options for the hydrological model calibration algorithm (choice of parameters to optimize, predefined values for uncalibrated parameters, etc.)

### GRiwrmInputsModel object

The method used for producing the `GRiwrmInputsModel` object is detailed in the vignette "V01_Structure_SD_model" of the package. The following code chunk resumes all the steps of this vignette:

```{r}
data(Severn)
nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
nodes$model <- "RunModel_GR4J"
griwrm <- CreateGRiwrm(nodes, list(id = "gauge_id", down = "downstream_id", length = "distance_downstream"))
BasinsObs <- Severn$BasinsObs
DatesR <- BasinsObs[[1]]$DatesR
PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
Qobs <- cbind(sapply(BasinsObs, function(x) {x$discharge_spec}))
Precip <- ConvertMeteoSD(griwrm, PrecipTot)
PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap)
str(InputsModel)
```

## GRiwrmRunOptions object

The `CreateRunOptions()` function allows to prepare the options required for the `RunModel()` function.

The user must at least define the following arguments:

* `InputsModel`: the associated input data
* `IndPeriod_Run`: the period on which the model is run

Below, we define a one-year warm up period and we start the run period just after the warm up period.

```{r}
IndPeriod_Run <- seq(
  which(InputsModel[[1]]$DatesR == (InputsModel[[1]]$DatesR[1] + 365*24*60*60)), # Set aside warm-up period
  length(InputsModel[[1]]$DatesR) # Until the end of the time series
)
IndPeriod_WarmUp <- seq(1, IndPeriod_Run[1] - 1)
```

Arguments of the `CreateRunOptions` function for **airGRiwrm** are the same as for the function in **airGR** and are copied for each node running a rainfall-runoff model.

```{r}
RunOptions <- CreateRunOptions(
  InputsModel,
  IndPeriod_WarmUp = IndPeriod_WarmUp,
  IndPeriod_Run = IndPeriod_Run
)
```

## GRiwrmInputsCrit object

The `CreateInputsCrit()` function allows to prepare the input in order to calculate a criterion.
We use composed criterion with a parameter regularization based on @delavenneRegularizationApproachImprove2019.

It needs the following arguments:

* `InputsModel`: the inputs of the `GRiwrm` network previously prepared by the `CreateInputsModel()` function
* `FUN_CRIT`: the name of the error criterion function (see the available functions description in the **airGR** package)
* `RunOptions`: the options of the `GRiwrm` network previously prepared by the `CreateRunOptions()` function
* `Qobs`: the observed variable time series (e.g. the discharge expressed in *mm/time step*)
* `AprioriIds`: the list of the sub-catchments IDs where to apply a parameter regularization based on the parameters of an upstream sub-catchment (e.g. here below the parameters of the sub-catchment "54057" is regulated by the parameters of the sub-catchment "54032")
* `transfo`: a transformation function applied on the flow before calculation of the criterion (square-root transformation is recommended for the De Lavenne regularization)
* `k`: coefficient used for the weighted average between the performance criterion and the gap between the optimized parameter set and an a priori parameter set (a value equal to 0.15 is recommended for the De Lavenne regularization)

```{r InputsCrit}
InputsCrit <- CreateInputsCrit(
  InputsModel = InputsModel,
  FUN_CRIT = ErrorCrit_KGE2,
  RunOptions = RunOptions,
  Obs = Qobs[IndPeriod_Run, ],
  AprioriIds = c(
      "54057" = "54032",
      "54032" = "54001",
      "54001" = "54095"
  ),
  transfo = "sqrt",
  k = 0.15
)
str(InputsCrit)
```


### GRiwrmCalibOptions object

Before using the automatic calibration tool, the user needs to prepare the calibration options with the `CreateCalibOptions()` function. The `GRiwrmInputsModel` argument contains all the necessary information:

```{r CalibOption}
CalibOptions <- CreateCalibOptions(InputsModel)
```

## Calibration

The **airGR** calibration process is applied on each node of the `GRiwrm` network from upstream nodes to downstream nodes.

```{r Calibration}
OutputsCalib <- suppressWarnings(
  Calibration(InputsModel, RunOptions, InputsCrit, CalibOptions))
```

## Run the model with the optimized model parameters

```{r RunModel}
OutputsModels <- RunModel(
  InputsModel,
  RunOptions = RunOptions,
  Param = extractParam(OutputsCalib)
)
```

## Plot the results for each basin

```{r, fig.height = 5, fig.width = 8}
plot(OutputsModels, Qobs = Qobs[IndPeriod_Run,])
```

The resulting flows of each node in m<sup>3</sup>/s are directly available and can be plotted with these commands:

```{r}
Qm3s <- attr(OutputsModels, "Qm3s")
plot(Qm3s[1:150,])
```