Commit 0a8a78a6 authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

refactor(articles): update texts (airGR v1.7.0)

parent 974cce06
......@@ -24,7 +24,6 @@ airGRteaching <- "[airGR](https://hydrogr.github.io/airGRteaching/)"
```{r, warning=FALSE, include=FALSE, fig.keep='none', results='hide'}
library(airGR)
library(DEoptim)
......@@ -52,11 +51,10 @@ In this article, we use the **GR4J** model to illustrate the different optimizat
In particular, we assume that the R global environment contains input climate data, observed discharge and functions from the [Get Started](page_1_get_started.html) article, as shown below.
Please note that the calibration period is defined in the `CreateRunOptions()` function .
<!--
```{r, warning=FALSE, fig.keep='none', results='hide', fig.height=10, fig.width=10, eval=TRUE, echo=FALSE, message=FALSE}
example("Calibration_Michel", echo = FALSE, ask = FALSE)
```
-->
<!-- ```{r, warning=FALSE, fig.keep='none', results='hide', fig.height=10, fig.width=10, eval=TRUE, echo=FALSE, message=FALSE} -->
<!-- example("Calibration_Michel", echo = FALSE, ask = FALSE) -->
<!-- ``` -->
```{r, echo=TRUE, eval=FALSE, eval=FALSE}
example("Calibration_Michel")
```
......@@ -221,6 +219,7 @@ resGLOB
Multiobjective optimization is used to explore possible trade-offs between different performances criteria.
Here we use the following R implementation of an efficient strategy:
* [caRamel: Automatic Calibration by Evolutionary Multi Objective Algorithm](https://cran.r-project.org/package=caRamel)
Motivated by using the rainfall-runoff model for low flow simulation, we explore the trade-offs between the KGE values obtained without any data transformation and with the inverse transformation.
......
......@@ -29,8 +29,6 @@ library(airGR)
library(coda)
library(FME)
library(ggmcmc)
library(dplyr)
# source("airGR.R")
set.seed(123)
load(system.file("vignettesData/vignetteParamMCMC.rda", package = "airGR"))
```
......@@ -176,7 +174,7 @@ ggmcmc::ggs_traceplot(parDRAM)
The posterior density for each parameter can then be visualised:
```{r, fig.width=6, fig.height=9, warning=FALSE}
burnParDRAM <- dplyr::filter(parDRAM, Iteration > 1000) # to keep only the second half of the series
burnParDRAM <- parDRAM[parDRAM$Iteration > 1000, ] # to keep only the second half of the series
ggmcmc::ggs_density(burnParDRAM)
```
......
......@@ -97,7 +97,7 @@ Ind_Cal <- seq(which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="01/01/
## preparation of the RunOptions object for the calibration period
RunOptions_Cal <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Cal)
InputsModel = InputsModel, IndPeriod_Run = Ind_Cal)
## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit_Cal <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
......@@ -116,8 +116,7 @@ RunOptions_Val <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
## efficiency criterion (Nash-Sutcliffe Efficiency) on the validation period
InputsCrit_Val <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions_Val, Obs = BasinObs$Qmm[Ind_Val])
RunOptions = RunOptions_Val, Obs = BasinObs$Qmm[Ind_Val])
```
# Calibration of the GR4J model with the generalist parameter sets
......@@ -148,7 +147,8 @@ Param_Best <- unlist(Param_Sets_GR4J[which.max(OutputsCrit_Loop), ])
Param_Best
## validation
OutputsModel_Val <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions_Val, Param = Param_Best)
OutputsModel_Val <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions_Val,
Param = Param_Best)
OutputsCrit_Val <- ErrorCrit_NSE(InputsCrit = InputsCrit_Val, OutputsModel = OutputsModel_Val)
```
......@@ -156,7 +156,7 @@ Now we can compute the Nash-Sutcliffe Efficiency criterion on the validation per
# Calibration of the GR4J model with the built-in `Calibration_Michel()` function
# Calibration of the GR4J model with the built-in Calibration_Michel function
As seen above, the Michel's calibration algorithm is based on a local search procedure.
......@@ -177,7 +177,7 @@ CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibrat
## calibration
OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
InputsCrit = InputsCrit_Cal, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J, FUN_CRIT = ErrorCrit_NSE,
FUN_MOD = RunModel_GR4J,
FUN_CALIB = Calibration_Michel)
OutputsModel_Cal <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
Param = OutputsCalib$ParamFinalR, FUN = RunModel_GR4J)
......@@ -214,10 +214,10 @@ CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibrat
## calibration
OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
InputsCrit = InputsCrit_Cal, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J, FUN_CRIT = ErrorCrit_NSE,
FUN_MOD = RunModel_GR4J,
FUN_CALIB = Calibration_Michel)
OutputsModel_Cal <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
Param = OutputsCalib$ParamFinalR, FUN = RunModel_GR4J)
Param = OutputsCalib$ParamFinalR, FUN = RunModel_GR4J)
OutputsCrit_Cal <- ErrorCrit_NSE(InputsCrit = InputsCrit_Cal, OutputsModel = OutputsModel_Cal)
......
......@@ -36,16 +36,18 @@ The `r airGR` package implements semi-distributed model capabilities using a lag
`RunModel_Lag` documentation gives an example of simulating the influence of a reservoir in a lumped model. Try `example(RunModel_Lag)` to get it.
In this vignette, we show how to calibrate 2 sub-catchments in series with a semi-distributed model consisting of 2 GR4J models. For doing this we compare two strategies for calibrating the downstream subcatchment:
In this article, we show how to calibrate 2 sub-catchments in series with a semi-distributed model consisting of 2 GR4J models.
For doing this we compare 3 strategies for calibrating the downstream subcatchment:
- using upstream observed flows
- using upstream simulated flows
- using upstream simulated flows and parameter regularisation [@delavenne_regularization_2019]
We finally compare these calibrations with a theoretical set of parameters.
This comparison is based on the Kling-Gupta Efficiency computed on the root-squared discharges as performance criterion.
## Model description
```{r, warning=FALSE, include=FALSE}
library(airGR)
options(digits = 3)
......@@ -74,6 +76,13 @@ summary(cbind(QObsUp = BasinObs$Qmm, QObsDown))
options(digits = 3)
```
With a delay of 2 days between the 2 gauging stations, the theoretical Velocity parameter should be equal to:
```{r}
Velocity <- 100 * 1e3 / (2 * 86400)
paste("Velocity: ", format(Velocity), "m/s")
```
# Calibration of the upstream subcatchment
The operations are exactly the same as the ones for a GR4J lumped model. So we do exactly the same operations as in the [Get Started](V01_get_started.html) vignette.
......@@ -87,9 +96,11 @@ RunOptionsUp <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModelUp,
IndPeriod_WarmUp = NULL, IndPeriod_Run = Ind_Run,
IniStates = NULL, IniResLevels = NULL)
InputsCritUp <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModelUp,
# Error criterion is KGE computed on the root-squared discharges
InputsCritUp <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModelUp,
RunOptions = RunOptionsUp,
VarObs = "Q", Obs = BasinObs$Qmm[Ind_Run])
VarObs = "Q", Obs = BasinObs$Qmm[Ind_Run],
transfo = "sqrt")
CalibOptionsUp <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel)
OutputsCalibUp <- Calibration_Michel(InputsModel = InputsModelUp, RunOptions = RunOptionsUp,
InputsCrit = InputsCritUp, CalibOptions = CalibOptionsUp,
......@@ -106,7 +117,11 @@ OutputsModelUp <- RunModel_GR4J(InputsModel = InputsModelUp, RunOptions = RunOpt
# Calibration of the downstream subcatchment with upstream flow observations
we need to create the `InputsModel` object completed with upstream information:
# Calibration of the downstream subcatchment
## Creation of the InputsModel objects
We need to create `InputsModel` objects completed with upstream information with upstream observed flow for the calibration of first case and upstream simulated flows for the other cases:
```{r}
InputsModelDown1 <- CreateInputsModel(
......@@ -118,16 +133,38 @@ InputsModelDown1 <- CreateInputsModel(
)
```
And then calibrate the combination of Lag model for upstream flow transfer and GR4J model for the runoff of the downstream subcatchment:
For using upstream simulated flows, we should concatenate a vector with the simulated flows for the entire period of simulation (warm-up + run):
```{r}
Qsim_upstream <- rep(NA, length(BasinObs$DatesR))
# Simulated flow during warm-up period (365 days before run period)
Qsim_upstream[Ind_Run[seq_len(365)] - 365] <- OutputsModelUp$RunOptions$WarmUpQsim
# Simulated flow during run period
Qsim_upstream[Ind_Run] <- OutputsModelUp$Qsim
InputsModelDown2 <- CreateInputsModel(
FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E,
Qupstream = matrix(Qsim_upstream, ncol = 1), # upstream observed flow
LengthHydro = 100, # distance between upstream catchment outlet & the downstream one [km]
BasinAreas = c(180, 180) # upstream and downstream areas [km²]
)
```
## Calibration with upstream flow observations
We calibrate the combination of Lag model for upstream flow transfer and GR4J model for the runoff of the downstream subcatchment:
```{r}
RunOptionsDown <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModelDown1,
IndPeriod_WarmUp = NULL, IndPeriod_Run = Ind_Run,
IniStates = NULL, IniResLevels = NULL)
InputsCritDown <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModelDown1,
InputsCritDown <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModelDown1,
RunOptions = RunOptionsDown,
VarObs = "Q", Obs = QObsDown[Ind_Run])
VarObs = "Q", Obs = QObsDown[Ind_Run],
transfo = "sqrt")
CalibOptionsDown <- CreateCalibOptions(FUN_MOD = RunModel_GR4J,
FUN_CALIB = Calibration_Michel,
IsSD = TRUE) # specify that it's a SD model
......@@ -138,13 +175,6 @@ OutputsCalibDown1 <- Calibration_Michel(InputsModel = InputsModelDown1,
FUN_MOD = RunModel_GR4J)
```
To run the complete model, we should substitute the observed upstream flow by the simulated one:
```{r}
InputsModelDown2 <- InputsModelDown1
InputsModelDown2$Qupstream[Ind_Run] <- OutputsModelUp$Qsim
```
`RunModel` is run in order to automatically combine GR4J and Lag models.
```{r}
......@@ -157,11 +187,11 @@ OutputsModelDown1 <- RunModel(InputsModel = InputsModelDown2,
Performance of the model validation is then:
```{r}
CritDown1 <- ErrorCrit_NSE(InputsCritDown, OutputsModelDown1)
KGE_down1 <- ErrorCrit_KGE(InputsCritDown, OutputsModelDown1)
```
# Calibration of the downstream subcatchment with upstream simulated flow
## Calibration with upstream simulated flow
We calibrate the model with the `InputsModel` object previously created for substituting the observed upstream flow with the simulated one:
......@@ -174,67 +204,105 @@ OutputsCalibDown2 <- Calibration_Michel(InputsModel = InputsModelDown2,
ParamDown2 <- OutputsCalibDown2$ParamFinalR
```
## Calibration with upstream simulated flow and parameter regularisation
# Discussion
The regularisation follow the method proposed by @delavenne_regularization_2019.
## Identification of Velocity parameter
As a priori parameter set, we use the calibrated parameter set of the upstream catchment and the theoretical velocity:
The theoretical Velocity parameter should be equal to:
```{r}
ParamDownTheo <- c(Velocity, OutputsCalibUp$ParamFinalR)
```
The Lavenne criterion is initialised with the a priori parameter set and the value of the KGE of the upstream basin.
```{r}
Velocity <- InputsModelDown1$LengthHydro * 1e3 / (2 * 86400)
paste(format(Velocity), "m/s")
IC_Lavenne <- CreateInputsCrit_Lavenne(InputsModel = InputsModelDown2,
RunOptions = RunOptionsDown,
Obs = QObsDown[Ind_Run],
AprParamR = ParamDownTheo,
AprCrit = OutputsCalibUp$CritFinal)
```
The Lavenne criterion is used instead of the KGE for calibration with regularisation
```{r}
OutputsCalibDown3 <- Calibration_Michel(InputsModel = InputsModelDown2,
RunOptions = RunOptionsDown,
InputsCrit = IC_Lavenne,
CalibOptions = CalibOptionsDown,
FUN_MOD = RunModel_GR4J)
```
The KGE is then calculated for performance comparisons:
```{r}
OutputsModelDown3 <- RunModel(InputsModel = InputsModelDown2,
RunOptions = RunOptionsDown,
Param = OutputsCalibDown3$ParamFinalR,
FUN_MOD = RunModel_GR4J)
KGE_down3 <- ErrorCrit_KGE(InputsCritDown, OutputsModelDown3)
```
# Discussion
## Identification of Velocity parameter
Both calibrations overestimate this parameter:
```{r}
mVelocity <- matrix(c(Velocity,
OutputsCalibDown1$ParamFinalR[1],
OutputsCalibDown2$ParamFinalR[1]),
OutputsCalibDown2$ParamFinalR[1],
OutputsCalibDown3$ParamFinalR[1]),
ncol = 1,
dimnames = list(c("theoretical",
"calibrated with observed upstream flow",
"calibrated with simulated upstream flow"),
"calibrated with simulated upstream flow",
"calibrated with sim upstream flow and regularisation"),
c("Velocity parameter")))
knitr::kable(mVelocity)
```
## Value of the performance criteria with theoretical calibration
Theoretically, the parameters of the downstream GR4J model should be the same as the upstream one and we know the lag time. So this set of parameter should give a better performance criteria:
Theoretically, the parameters of the downstream GR4J model should be the same as the upstream one with the velocity as extra parameter:
```{r}
ParamDownTheo <- c(Velocity, OutputsCalibUp$ParamFinalR)
OutputsModelDownTheo <- RunModel(InputsModel = InputsModelDown2,
RunOptions = RunOptionsDown,
Param = ParamDownTheo,
FUN_MOD = RunModel_GR4J)
CritDownTheo <- ErrorCrit_NSE(InputsCritDown, OutputsModelDownTheo)
KGE_downTheo <- ErrorCrit_KGE(InputsCritDown, OutputsModelDownTheo)
```
## Parameters and performance of each subcatchment for all calibrations
```{r}
comp <- matrix(c(0, OutputsCalibUp$ParamFinalR,
rep(OutputsCalibDown1$ParamFinalR, 2),
OutputsCalibDown2$ParamFinalR,
OutputsCalibDown3$ParamFinalR,
ParamDownTheo),
ncol = 5, byrow = TRUE)
comp <- cbind(comp, c(OutputsCalibUp$CritFinal,
OutputsCalibDown1$CritFinal,
CritDown1$CritValue,
KGE_down1$CritValue,
OutputsCalibDown2$CritFinal,
CritDownTheo$CritValue))
colnames(comp) <- c("Velocity", paste0("X", 1:4), "NSE")
KGE_down3$CritValue,
KGE_downTheo$CritValue))
colnames(comp) <- c("Velocity", paste0("X", 1:4), "KGE(√Q)")
rownames(comp) <- c("Calibration of the upstream subcatchment",
"Calibration 1 with observed upstream flow",
"Validation 1 with simulated upstream flow",
"Calibration 2 with simulated upstream flow",
"Calibration 3 with simulated upstream flow and regularisation",
"Validation theoretical set of parameters")
knitr::kable(comp)
```
Even if calibration with observed upstream flows gives an improved performance criteria, in validation using simulated upstream flows the result is quite similar as the performance obtained with the calibration with upstream simulated flows. The theoretical set of parameters give also an equivalent performance but still underperforming the calibration 2 one.
Even if calibration with observed upstream flows gives an improved performance criteria, in validation using simulated upstream flows the result is quite similar as the performance obtained with the calibration with upstream simulated flows. The theoretical set of parameters give also an equivalent performance but still underperforming the calibration 2 one. Regularisation allows to get similar performance as the one for calibration with simulated flows but with the big advantage of having parameters closer to the theoretical ones (Especially for the velocity parameter).
# References
......@@ -178,4 +178,25 @@
year = {2019},
keywords = {airGRcite},
pages = {70--81}
}
@article{delavenne_regularization_2019,
ids = {lavenneRegularizationApproachImprove2019a},
title = {A {{Regularization Approach}} to {{Improve}} the {{Sequential Calibration}} of a {{Semidistributed Hydrological Model}}},
author = {de Lavenne, Alban and Andréassian, Vazken and Thirel, Guillaume and Ramos, Maria-Helena and Perrin, Charles},
date = {2019},
journaltitle = {Water Resources Research},
volume = {55},
pages = {8821--8839},
issn = {1944-7973},
doi = {10.1029/2018WR024266},
url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018WR024266},
urldate = {2021-06-07},
abstract = {In semidistributed hydrological modeling, sequential calibration usually refers to the calibration of a model by considering not only the flows observed at the outlet of a catchment but also the different gauging points inside the catchment from upstream to downstream. While sequential calibration aims to optimize the performance at these interior gauged points, we show that it generally fails to improve performance at ungauged points. In this paper, we propose a regularization approach for the sequential calibration of semidistributed hydrological models. It consists in adding a priori information on optimal parameter sets for each modeling unit of the semidistributed model. Calibration iterations are then performed by jointly maximizing simulation performance and minimizing drifts from the a priori parameter sets. The combination of these two sources of information is handled by a parameter k to which the method is quite sensitive. The method is applied to 1,305 catchments in France over 30 years. The leave-one-out validation shows that, at locations considered as ungauged, model simulations are significantly improved (over all the catchments, the median KGE criterion is increased from 0.75 to 0.83 and the first quartile from 0.35 to 0.66), while model performance at gauged points is not significantly impacted by the use of the regularization approach. Small catchments benefit most from this calibration strategy. These performances are, however, very similar to the performances obtained with a lumped model based on similar conceptualization.},
annotation = {\_eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2018WR024266},
file = {C\:\\Users\\david.dorchies\\Zotero\\storage\\QNCGSJXB\\Lavenne et al. - 2019 - A Regularization Approach to Improve the Sequentia.pdf;C\:\\Users\\david.dorchies\\Zotero\\storage\\5DMRIV8Y\\2018WR024266.html},
keywords = {regularization,semidistributed model,stepwise calibration,ungauged basins},
langid = {english},
number = {11},
options = {useprefix=true}
}
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment