Skip to content
GitLab

Source

Target

Hide whitespace changes
Inline Side-by-side
Showing
with 956 additions and 117 deletions
+956 -117
tests/testthat/test-vignettes.R 0 → 100644
View file @ 4bccfee2
context("Test vignette chunks")
test_that("V01_get_started works", {
skip_on_cran()
rm(list = ls())
expect_true(RunVignetteChunks("V01_get_started"))
TestQmmQlsConversion(BasinObs, BasinInfo$BasinArea)
})
test_that("V02.1_param_optim works", {
skip_on_cran()
rm(list = ls())
load(system.file("vignettesData/vignetteParamOptim.rda", package = "airGR"))
load(system.file("vignettesData/vignetteParamOptimCaramel.rda", package = "airGR"))
rda_resGLOB <- resGLOB
rda_resPORT <- resPORT
rda_optMO <- optMO
expect_true(RunVignetteChunks("V02.1_param_optim"))
expect_equal(summary(resGLOB), summary(rda_resGLOB), tolerance = 1e-7)
expect_equal(resGLOB[, -1], rda_resGLOB[, -1], tolerance = 1e-2) # High tolerance due to randomisation in optimisations
expect_equal(summary(optMO$parameters), summary(rda_optMO$parameters), tolerance = 1e-7)
})
test_that("V02.2_param_mcmc works", {
skip_on_cran()
rm(list = ls())
load(system.file("vignettesData/vignetteParamMCMC.rda", package = "airGR"))
rda_gelRub <- gelRub
rda_multDRAM <- multDRAM
expect_true(RunVignetteChunks("V02.2_param_mcmc"))
expect_equal(gelRub, rda_gelRub, tolerance = 1e-7)
expect_equal(multDRAM, rda_multDRAM, tolerance = 1e-7)
})
test_that("V03_param_sets_GR4J works", {
skip_on_cran()
rm(list = ls())
expect_true(RunVignetteChunks("V03_param_sets_GR4J"))
})
test_that("V04_cemaneige_hysteresis works", {
skip_on_cran()
rm(list = ls())
load(system.file("vignettesData/vignetteCNHysteresis.rda", package = "airGR"))
rda_OutputsCrit_Cal <- OutputsCrit_Cal
rda_OutputsCrit_Cal_NoHyst <- OutputsCrit_Cal_NoHyst
rda_OutputsCrit_Val <- OutputsCrit_Val
rda_OutputsCrit_Val_NoHyst <- OutputsCrit_Val_NoHyst
expect_true(RunVignetteChunks("V04_cemaneige_hysteresis"))
TestQmmQlsConversion(BasinObs, BasinInfo$BasinArea)
expect_equal(OutputsCrit_Cal, rda_OutputsCrit_Cal, tolerance = 1e-7)
expect_equal(OutputsCrit_Cal_NoHyst, rda_OutputsCrit_Cal_NoHyst, tolerance = 1e-7)
expect_equal(OutputsCrit_Val, rda_OutputsCrit_Val, tolerance = 1e-7)
expect_equal(OutputsCrit_Val_NoHyst, rda_OutputsCrit_Val_NoHyst, tolerance = 1e-7)
})
vignettes/V00_airgr_ref.bib
View file @ 4bccfee2
@phdthesis{ficchi_adaptive_2017,
title = {An adaptive hydrological model for multiple time-steps: diagnostics and improvements based on fluxes consistency},
shorttitle = {An adaptive hydrological model for multiple time-steps},
url = {https://www.theses.fr/2017PA066097},
abstract = {Cette thèse vise à explorer la question du changement d'échelle temporelle en modélisation hydrologique conceptuelle. Les principaux objectifs sont : (i) étudier les effets du changement du pas de temps sur les performances, les paramètres et la structure des modèles hydrologiques ; (ii) mettre au point un modèle pluie-débit applicable à différents pas de temps. Notre point de départ est le modèle global journalier GR4J, développé à Irstea. Ce modèle a été choisi comme le modèle de référence à adapter à d'autres résolutions plus fines, jusqu'à des pas de temps infra-horaires, en suivant une approche descendante. Pour nos tests, nous avons construit une base de données de 240 bassins versants non influencés en France, à différents pas de temps allant de 6 minutes à 1 jour, en utilisant: (i) les données pluviométriques à 6 minutes et la réanalyse des lames d'eau journalières à plus haute résolution spatiale ; (ii) les données de température journalière pour le calcul de l'évapotranspiration potentielle ; (iii) les données hydrométriques à pas de temps variable. Nous avons étudié l'impact de la distribution temporelle des entrées sur les performances du modèle en se focalisant sur la simulation de crue, sur la base de 2400 événements. Ensuite, notre évaluation du modèle a porté sur l'analyse de la cohérence des flux internes du modèle à différents pas de temps, afin d'assurer une performance satisfaisante à travers un fonctionnement du modèle cohérent. Notre diagnostic du modèle nous a permis d'identifier une amélioration de la structure du modèle à différents pas de temps infra-journaliers basée sur la complexification de la composante d'interception du modèle.},
urldate = {2017-11-24},
school = {Université Pierre et Marie Curie, Paris 6},
author = {Ficchi, Andrea},
month = feb,
year = {2017},
keywords = {Rainfall-runoff modelling, airGRcite, Cohérence structurelle, Événements de crue, Interception, Modèles GR, Modélisation pluie-Débit, Pas de temps fin, Short time step}
}
@article{ficchi_hydrological_2019,
title = {Hydrological modelling at multiple sub-daily time steps: model improvement via flux-matching},
issn = {00221694},
shorttitle = {Hydrological modelling at multiple sub-daily time steps},
url = {https://linkinghub.elsevier.com/retrieve/pii/S0022169419305281},
doi = {10.1016/j.jhydrol.2019.05.084},
language = {en},
urldate = {2019-06-11},
journal = {Journal of Hydrology},
author = {Ficchì, Andrea and Perrin, Charles and Andréassian, Vazken},
month = jun,
year = {2019},
keywords = {airGR, airGRcite}
}
@phdthesis{mathevet_quels_2005, @phdthesis{mathevet_quels_2005,
address = {Paris}, address = {Paris},
title = {Quels modèles pluie-débit globaux au pas de temps horaire ? {Développements} empiriques et comparaison de modèles sur un large échantillon de bassins versants}, title = {Quels modèles pluie-débit globaux au pas de temps horaire ? {Développements} empiriques et comparaison de modèles sur un large échantillon de bassins versants},
...@@ -133,4 +161,42 @@ ...@@ -133,4 +161,42 @@
month = oct, month = oct,
year = {2014}, year = {2014},
pages = {8356--8366} pages = {8356--8366}
} }
\ No newline at end of file
@article{riboust_revisiting_2019,
title = {Revisiting a {Simple} {Degree}-{Day} {Model} for {Integrating} {Satellite} {Data}: {Implementation} of {Swe}-{Sca} {Hystereses}},
volume = {67},
issn = {0042-790X},
shorttitle = {Revisiting a {Simple} {Degree}-{Day} {Model} for {Integrating} {Satellite} {Data}},
url = {http://content.sciendo.com/view/journals/johh/67/1/article-p70.xml},
doi = {10.2478/johh-2018-0004},
number = {1},
urldate = {2019-02-18},
journal = {Journal of Hydrology and Hydromechanics},
author = {Riboust, Philippe and Thirel, Guillaume and Moine, Nicolas Le and Ribstein, Pierre},
month = mar,
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}
}
vignettes/V01_get_started.Rmd
View file @ 4bccfee2
--- ---
title: "Get Started with airGR" title: "Get Started with airGR"
author: "Guillaume Thirel, Olivier Delaigue, Laurent Coron"
bibliography: V00_airgr_ref.bib bibliography: V00_airgr_ref.bib
output: rmarkdown::html_vignette output: rmarkdown::html_vignette
vignette: > vignette: >
...@@ -10,29 +11,32 @@ vignette: > ...@@ -10,29 +11,32 @@ vignette: >
# Introduction # Introduction
**airGR** is a package that brings into the [**R software**](https://cran.r-project.org/) the hydrological modelling tools used and developed at the [Catchment Hydrology Research Group](https://webgr.irstea.fr/en/) at [Irstea (France)](http://www.irstea.fr/en/), including the [**GR rainfall-runoff models**](https://webgr.irstea.fr/en/modeles/) and a snowmelt and accumulation model, [**CemaNeige**](https://webgr.irstea.fr/en/modeles/modele-de-neige/). Each model core is coded in **Fortran** to ensure low computational time. The other package functions (i.e. mainly the calibration algorithm and the efficiency criteria calculation) are coded in **R**. **airGR** is a package that brings into the [**R software**](https://cran.r-project.org/) the hydrological modelling tools used and developed at the [Catchment Hydrology Research Group](https://webgr.inrae.fr/home/) at [INRAE (France)](https://www.inrae.fr/en), including the [**GR rainfall-runoff models**](https://webgr.inrae.fr/models/) that can be applied either on a **lumped** or **semi-distributed** way. A snow accumulation and melt model ([**CemaNeige**](https://webgr.inrae.fr/models/snow-model/)) and the associated functions for the calibration and evaluation of models are also included. Each model core is coded in **Fortran** to ensure low computational time. The other package functions (i.e. mainly the calibration algorithm and the efficiency criteria calculation) are coded in **R**.
The **airGR** package has been designed to fulfill two major requirements: to facilitate the use by non-expert users and to allow flexibility regarding the addition of external criteria, models or calibration algorithms. The names of the functions and their arguments were chosen to this end. **airGR** also contains basics plotting facilities. The **airGR** package has been designed to fulfill two major requirements: to facilitate the use by non-expert users and to allow flexibility regarding the addition of external criteria, models or calibration algorithms. The names of the functions and their arguments were chosen to this end. **airGR** also contains basics plotting facilities.
Six hydrological models and one snowmelt and accumulation model are implemented in **airGR**. The snow model can be used alone or together with the daily hydrological models. Seven hydrological models and one snow melt and accumulation model are implemented in **airGR**. The hydrological models can be applied either on a lumped way or on a semi-distributed way (on sub-catchments). The snow model can either be used alone or with the daily or hourly hydrological models. Naturally each hydrological model can also be used alone.
The models can be called within **airGR** using the following functions: The models can be called within **airGR** using the following functions:
* `RunModel_GR4H()`: four-parameter hourly lumped hydrological model [@mathevet_quels_2005] * `RunModel_GR4H()`: four-parameter hourly lumped hydrological model [@mathevet_quels_2005]
* `RunModel_GR5H()`: five-parameter hourly lumped hydrological model [@ficchi_adaptive_2017; @ficchi_hydrological_2019]
* `RunModel_GR4J()`: four-parameter daily lumped hydrological model [@perrin_improvement_2003] * `RunModel_GR4J()`: four-parameter daily lumped hydrological model [@perrin_improvement_2003]
* `RunModel_GR5J()`: five-parameter daily lumped hydrological model [@le_moine_bassin_2008] * `RunModel_GR5J()`: five-parameter daily lumped hydrological model [@le_moine_bassin_2008]
* `RunModel_GR6J()`: six-parameter daily lumped hydrological model [@pushpalatha_downward_2011] * `RunModel_GR6J()`: six-parameter daily lumped hydrological model [@pushpalatha_downward_2011]
* `RunModel_GR2M()`: two-parameter monthly lumped hydrological model [@mouelhi_vers_2003; @mouelhi_stepwise_2006] * `RunModel_GR2M()`: two-parameter monthly lumped hydrological model [@mouelhi_vers_2003; @mouelhi_stepwise_2006]
* `RunModel_GR1A()`: one-parameter yearly lumped hydrological model [@mouelhi_vers_2003; @mouelhi_linking_2006] * `RunModel_GR1A()`: one-parameter yearly lumped hydrological model [@mouelhi_vers_2003; @mouelhi_linking_2006]
* `RunModel_CemaNeige()`: two-parameter degree-day snowmelt and accumulation model [@valery_as_2014] * `RunModel_CemaNeige()`: two-parameter degree-day snowmelt and accumulation model [@valery_as_2014; @riboust_revisiting_2019]
* `RunModel_CemaNeigeGR4H()`: combined use of **GR4H** and **CemaNeige**
* `RunModel_CemaNeigeGR5H()`: combined use of **GR5H** and **CemaNeige**
* `RunModel_CemaNeigeGR4J()`: combined use of **GR4J** and **CemaNeige** * `RunModel_CemaNeigeGR4J()`: combined use of **GR4J** and **CemaNeige**
* `RunModel_CemaNeigeGR5J()`: combined use of **GR5J** and **CemaNeige** * `RunModel_CemaNeigeGR5J()`: combined use of **GR5J** and **CemaNeige**
* `RunModel_CemaNeigeGR6J()`: combined use of **GR6J** and **CemaNeige** * `RunModel_CemaNeigeGR6J()`: combined use of **GR6J** and **CemaNeige**
The [**GRP**](https://webgr.irstea.fr/en/modeles/modele-de-prevision-grp/) forecasting model and the [**Otamin**](https://webgr.irstea.fr/en/modeles/otamin/) predictive uncertainty tool are not available in **airGR**. The [**GRP**](https://webgr.inrae.fr/models/hydrological-forecasting-model-grp/) forecasting model and the [**Otamin**](https://webgr.inrae.fr/software/otamin/) predictive uncertainty tool are not available in **airGR**.
In this vignette, we show how to prepare and run a calibration and a simulation with airGR hydrological models. In this vignette, we show how to prepare and run a calibration and a simulation with airGR hydrological models.
...@@ -44,7 +48,7 @@ In the following example, we use a data sample contained in the package. For rea ...@@ -44,7 +48,7 @@ In the following example, we use a data sample contained in the package. For rea
First, it is necessary to load the **airGR** package: First, it is necessary to load the **airGR** package:
```{r} ```{r, warning=FALSE}
library(airGR) library(airGR)
``` ```
...@@ -59,7 +63,7 @@ Below is presented an example of a `data.frame` of daily hydrometeorological obs ...@@ -59,7 +63,7 @@ Below is presented an example of a `data.frame` of daily hydrometeorological obs
```{r} ```{r}
data(L0123001) data(L0123001)
summary(BasinObs) summary(BasinObs, digits = 2)
``` ```
The usual functions (e.g. `read.table()`) can be used to load real-case data sets. The usual functions (e.g. `read.table()`) can be used to load real-case data sets.
...@@ -82,7 +86,7 @@ To facilitate the use of the package, there are several functions dedicated to t ...@@ -82,7 +86,7 @@ To facilitate the use of the package, there are several functions dedicated to t
To run a GR hydrological model or CemaNeige, the user has to prepare the input data with the `CreateInputsModel()` function. To run a GR hydrological model or CemaNeige, the user has to prepare the input data with the `CreateInputsModel()` function.
As arguments, this function needs the function name corresponding to the model the user wants to run, a vector of dates, a vector of precipitation and a vector of potential evapotranspiration. As arguments, this function needs the function name corresponding to the model the user wants to run, a vector of dates, a vector of precipitation and a vector of potential evapotranspiration.
In the example below, we already have the potential evapotranspiration. If the user does not have these data, it is possible to compute it with the [Oudin's formula](http://dx.doi.org/10.1016/j.jhydrol.2004.08.026) with the `PEdaily_Oudin()` function (this function only needs Julian days, daily average air temperature and latitude). In the example below, we already have the potential evapotranspiration. If the user does not have these data, it is possible to compute it with the [Oudin's formula](http://dx.doi.org/10.1016/j.jhydrol.2004.08.026) with the `PE_Oudin()` function (this function only needs Julian days, daily average air temperature and latitude).
Missing values (`NA`) of precipitation (or potential evapotranspiration) are **not allowed**. Missing values (`NA`) of precipitation (or potential evapotranspiration) are **not allowed**.
...@@ -94,7 +98,6 @@ str(InputsModel) ...@@ -94,7 +98,6 @@ str(InputsModel)
``` ```
## RunOptions object ## RunOptions object
The `CreateRunOptions()` function allows to prepare the options required to the `RunModel*()` functions, which are the actual models functions. The `CreateRunOptions()` function allows to prepare the options required to the `RunModel*()` functions, which are the actual models functions.
...@@ -102,14 +105,14 @@ The `CreateRunOptions()` function allows to prepare the options required to the ...@@ -102,14 +105,14 @@ The `CreateRunOptions()` function allows to prepare the options required to the
The user must at least define the following arguments: The user must at least define the following arguments:
* `FUN_MOD`: the name of the model function to run * `FUN_MOD`: the name of the model function to run
* `InputsModel`: the associated inputs data * `InputsModel`: the associated input data
* `IndPeriod_Run`: the period on which the model is run * `IndPeriod_Run`: the period on which the model is run
To select a period for which the user wants to run the model, select the corresponding indexes for different time periods (not the POSIXt dates), as follows: To select a period for which the user wants to run the model, select the corresponding indexes for different time periods (not the POSIXt dates), as follows:
```{r} ```{r}
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%d/%m/%Y")=="01/01/1990"), Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"),
which(format(BasinObs$DatesR, format = "%d/%m/%Y")=="31/12/1999")) which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1999-12-31"))
str(Ind_Run) str(Ind_Run)
``` ```
...@@ -136,22 +139,25 @@ The `CreateRunOptions()` function returns warnings if the default initialization ...@@ -136,22 +139,25 @@ The `CreateRunOptions()` function returns warnings if the default initialization
## InputsCrit object ## InputsCrit object
The `CreateInputsCrit()` function allows to prepare the input in order to calculate a criterion. It is possible to define the following arguments: The `CreateInputsCrit()` function allows to prepare the input in order to calculate a criterion. It is possible to define the following arguments:
* `FUN_CRIT`: the name of the error criterion function (the available functions are introduced later on) * `FUN_CRIT`: the name of the error criterion function (the available functions are introduced later on)
* `InputsModel`: the inputs of the hydrological model previously prepared by the `CreateInputsModel()` function * `InputsModel`: the inputs of the hydrological model previously prepared by the `CreateInputsModel()` function
* `RunOptions`: the options of the hydrological model previously prepared by the `CreateRunOptions()` function * `RunOptions`: the options of the hydrological model previously prepared by the `CreateRunOptions()` function
* `Qobs`: the observed discharge expressed in *mm/time step* * `VarObs`: the name of the considered variable (by default `"Q"` for the discharge)
* `Obs`: the observed variable time serie (e.g. the discharge expressed in *mm/time step*)
Missing values (`NA`) are **allowed** for observed discharge. Missing values (`NA`) are **allowed** for observed discharge.
It is possible to compute a composite criterion (e.g. the average between NSE computed on discharge and NSE computed on log of discharge). In this case, users have to provide lists to the following arguments (some of the are optional): `FUN_CRIT`, `Obs`, `VarObs`, `BoolCrit`, `transfo`, `Weights.`
```{r} ```{r}
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run]) RunOptions = RunOptions, VarObs = "Q", Obs = BasinObs$Qmm[Ind_Run])
str(InputsCrit) str(InputsCrit)
``` ```
## CalibOptions object ## CalibOptions object
Before using the automatic calibration tool, the user needs to prepare the calibration options with the `CreateCalibOptions()` function. For that, it is necessary to define the following arguments: Before using the automatic calibration tool, the user needs to prepare the calibration options with the `CreateCalibOptions()` function. For that, it is necessary to define the following arguments:
...@@ -195,7 +201,7 @@ The calibration algorithm optimizes the error criterion selected as objective-fu ...@@ -195,7 +201,7 @@ The calibration algorithm optimizes the error criterion selected as objective-fu
```{r} ```{r}
OutputsCalib <- Calibration_Michel(InputsModel = InputsModel, RunOptions = RunOptions, OutputsCalib <- Calibration_Michel(InputsModel = InputsModel, RunOptions = RunOptions,
InputsCrit = InputsCrit, CalibOptions = CalibOptions, InputsCrit = InputsCrit, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J, FUN_CRIT = ErrorCrit_NSE) FUN_MOD = RunModel_GR4J)
Param <- OutputsCalib$ParamFinalR Param <- OutputsCalib$ParamFinalR
Param Param
``` ```
...@@ -217,7 +223,6 @@ Performing a model control with **airGR** is similar to running a simulation (se ...@@ -217,7 +223,6 @@ Performing a model control with **airGR** is similar to running a simulation (se
# Simulation # Simulation
## Simulation run ## Simulation run
To run a model, the user has to use the `RunModel*()` functions (`InputsModel`, `RunOptions` and parameters). To run a model, the user has to use the `RunModel*()` functions (`InputsModel`, `RunOptions` and parameters).
...@@ -229,10 +234,9 @@ str(OutputsModel) ...@@ -229,10 +234,9 @@ str(OutputsModel)
``` ```
## Results preview ## Results preview
Although it is possible for the user to design its own graphics from the outputs of the `RunModel*()` functions, the **airGR** package offers the possibility to make use of the `plot.OutputsModel()` function (or `plot()` with an `OutputsModel` object). This function returns a dashboard of results including various graphs (depending on the model used): Although it is possible for the user to design its own graphics from the outputs of the `RunModel*()` functions, the **airGR** package offers the possibility to make use of the `plot()` function. This function returns a dashboard of results including various graphs (depending on the model used):
* time series of total precipitation and simulated discharge (and observed discharge if provided) * time series of total precipitation and simulated discharge (and observed discharge if provided)
* interannual average daily simulated discharge (and daily observed discharge if provided) and interannual average monthly precipitation * interannual average daily simulated discharge (and daily observed discharge if provided) and interannual average monthly precipitation
...@@ -252,11 +256,19 @@ To evaluate the efficiency of the model, it is possible to use the same criterio ...@@ -252,11 +256,19 @@ To evaluate the efficiency of the model, it is possible to use the same criterio
```{r} ```{r}
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel) OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
str(OutputsCrit)
``` ```
```{r} ```{r}
OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel) OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
str(OutputsCrit)
``` ```
## Features diagram
The figure below presents the implementation flow of the **airGR** functions.
<br><br>
<img src="fig/airGR_features_diagram.svg" alt="airGR features diagram" width=75% />
# References # References
vignettes/V02.1_param_optim.Rmd
View file @ 4bccfee2
--- ---
title: "Plugging in new calibration algorithms in airGR" title: "Plugging in new calibration algorithms in airGR"
author: "François Bourgin" author: "François Bourgin, Guillaume Thirel"
output: rmarkdown::html_vignette output: rmarkdown::html_vignette
vignette: > vignette: >
%\VignetteEngine{knitr::rmarkdown} %\VignetteEngine{knitr::rmarkdown}
...@@ -10,13 +10,18 @@ vignette: > ...@@ -10,13 +10,18 @@ vignette: >
```{r, warning=FALSE, include=FALSE, fig.keep='none', results='hide'} ```{r setup, warning=FALSE, include=FALSE, fig.keep='none', results='hide'}
library(airGR) library(airGR)
library(DEoptim) library(DEoptim)
library(hydroPSO) # library(hydroPSO) # Needs R version >= 3.6 or latticeExtra <= 0.6-28 on R 3.5. Archived on 2023-10-16 as requires archived packages 'hydroTSM' and 'hydroGOF'.
library(Rmalschains) library(Rmalschains)
library(caRamel)
library(ggplot2)
library(GGally)
# source("airGR.R") # source("airGR.R")
set.seed(321) set.seed(321)
load(system.file("vignettesData/vignetteParamOptim.rda", package = "airGR"))
load(system.file("vignettesData/vignetteParamOptimCaramel.rda", package = "airGR"))
``` ```
...@@ -32,18 +37,27 @@ In this vignette, we use the **GR4J** model to illustrate the different optimiza ...@@ -32,18 +37,27 @@ In this vignette, we use the **GR4J** model to illustrate the different optimiza
In particular, we assume that the R global environment contains input climate data, observed discharge and functions from the [Get Started](V01_get_started.html) vignette, as shown below. In particular, we assume that the R global environment contains input climate data, observed discharge and functions from the [Get Started](V01_get_started.html) vignette, as shown below.
Please note that the calibration period is defined in the `CreateRunOptions()` function . 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} <!-- ```{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) <!-- example("Calibration_Michel", echo = FALSE, ask = FALSE) -->
``` <!-- ``` -->
```{r, echo=TRUE, eval=FALSE}
```{r Calibration_Michel, echo=TRUE, eval=FALSE}
example("Calibration_Michel") example("Calibration_Michel")
``` ```
In order for the `RunModel_*()` functions to run faster during the parameter estimation process, it is recommended that the outputs contain only the simulated flows (see the `Outputs_Sim` argument in the `CreateRunOptions()` help page).
```{r RunOptions, results='hide', eval=FALSE}
RunOptions <- airGR::CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
IndPeriod_Run = Ind_Run,
Outputs_Sim = "Qsim")
```
Regarding the different optimization strategies presented here, we refer to each package for in-depth information about the description of the methods used. Regarding the different optimization strategies presented here, we refer to each package for in-depth information about the description of the methods used.
Please note that this vignette is only for illustration purposes and does not provide any guidance about which optimization strategies is recommended for the family of the **GR** models. Please note that this vignette is only for illustration purposes and does not provide any guidance about which optimization strategies is recommended for the family of the **GR** models.
## Definition of the necessary function ## Definition of the necessary function
Parameter estimation can be performed by defining a function that takes a parameter set as input and returns the value of the performance criterion. Parameter estimation can be performed by defining a function that takes a parameter set as input and returns the value of the performance criterion.
...@@ -52,16 +66,15 @@ Here we choose to minimize the root mean square error. ...@@ -52,16 +66,15 @@ Here we choose to minimize the root mean square error.
The change of the repository from the "real" parameter space to a "transformed" space ensures homogeneity of displacement in the different dimensions of the parameter space during the step-by-step procedure of the calibration algorithm of the model. The change of the repository from the "real" parameter space to a "transformed" space ensures homogeneity of displacement in the different dimensions of the parameter space during the step-by-step procedure of the calibration algorithm of the model.
```{r OptimGR4J, warning=FALSE, results='hide', eval=FALSE}
```{r, warning=FALSE, results='hide'} OptimGR4J <- function(ParamOptim) {
OptimGR4J <- function(Param_Optim) {
## Transformation of the parameter set to real space ## Transformation of the parameter set to real space
Param_Optim_Vre <- airGR::TransfoParam_GR4J(ParamIn = Param_Optim, RawParamOptim <- airGR::TransfoParam_GR4J(ParamIn = ParamOptim,
Direction = "TR") Direction = "TR")
## Simulation given a parameter set ## Simulation given a parameter set
OutputsModel <- airGR::RunModel_GR4J(InputsModel = InputsModel, OutputsModel <- airGR::RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions, RunOptions = RunOptions,
Param = Param_Optim_Vre) Param = RawParamOptim)
## Computation of the value of the performance criteria ## Computation of the value of the performance criteria
OutputsCrit <- airGR::ErrorCrit_RMSE(InputsCrit = InputsCrit, OutputsCrit <- airGR::ErrorCrit_RMSE(InputsCrit = InputsCrit,
OutputsModel = OutputsModel, OutputsModel = OutputsModel,
...@@ -72,52 +85,63 @@ OptimGR4J <- function(Param_Optim) { ...@@ -72,52 +85,63 @@ OptimGR4J <- function(Param_Optim) {
In addition, we need to define the lower and upper bounds of the four **GR4J** parameters in the transformed parameter space: In addition, we need to define the lower and upper bounds of the four **GR4J** parameters in the transformed parameter space:
```{r, warning=FALSE, results='hide'}
```{r boundsGR4J, warning=FALSE, results='hide', eval=FALSE}
lowerGR4J <- rep(-9.99, times = 4) lowerGR4J <- rep(-9.99, times = 4)
upperGR4J <- rep(+9.99, times = 4) upperGR4J <- rep(+9.99, times = 4)
``` ```
# Local optimization # Local optimization
We start with a local optimization strategy by using the PORT routines (using the `nlminb()` of the `stats` package) and by setting a starting point in the transformed parameter space: We start with a local optimization strategy by using the PORT routines (using the `nlminb()` of the `stats` package) and by setting a starting point in the transformed parameter space:
```{r, warning=FALSE, results='hide'}
```{r local1, warning=FALSE, results='hide', eval=FALSE}
startGR4J <- c(4.1, 3.9, -0.9, -8.7) startGR4J <- c(4.1, 3.9, -0.9, -8.7)
optPORT <- stats::nlminb(start = startGR4J, optPORT <- stats::nlminb(start = startGR4J,
objective = OptimGR4J, objective = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J, lower = lowerGR4J, upper = upperGR4J,
control = list(trace = 1)) control = list(trace = 1))
``` ```
The RMSE value reaches a local minimum value after 35 iterations. The RMSE value reaches a local minimum value after 35 iterations.
We can also try a multi-start approach to test the consistency of the local optimization. We can also try a multi-start approach to test the consistency of the local optimization.
Here we use the same grid used for the filtering step of the Michel's calibration strategy (`Calibration_Michel()` function). Here we use the same grid used for the filtering step of the Michel's calibration strategy (`Calibration_Michel()` function).
For each starting point, a local optimization is performed. For each starting point, a local optimization is performed.
```{r, warning=FALSE, results='hide'}
startGR4J <- expand.grid(data.frame(CalibOptions$StartParamDistrib)) ```{r local2, warning=FALSE, results='hide', eval=FALSE}
startGR4JDistrib <- TransfoParam_GR4J(ParamIn = CalibOptions$StartParamDistrib,
Direction = "RT")
startGR4J <- expand.grid(data.frame(startGR4JDistrib))
optPORT_ <- function(x) { optPORT_ <- function(x) {
opt <- stats::nlminb(start = x, opt <- stats::nlminb(start = x,
objective = OptimGR4J, objective = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J, lower = lowerGR4J, upper = upperGR4J,
control = list(trace = 1)) control = list(trace = 1))
} }
list_opt <- apply(startGR4J, 1, optPORT_) listOptPORT <- apply(startGR4J, MARGIN = 1, FUN = optPORT_)
``` ```
We can then extract the best parameter sets and the value of the performance criteria: We can then extract the best parameter sets and the value of the performance criteria:
```{r, warning=FALSE, results='hide'}
list_par <- t(sapply(list_opt, function(x) x$par)) ```{r local3, warning=FALSE, results='hide', eval=FALSE}
list_obj <- sapply(list_opt, function(x) x$objective) parPORT <- t(sapply(listOptPORT, function(x) x$par))
df_port <- data.frame(list_par, RMSE = list_obj) objPORT <- sapply(listOptPORT, function(x) x$objective)
resPORT <- data.frame(parPORT, RMSE = objPORT)
``` ```
As can be seen below, the optimum performance criterion values (column *objective*) can differ from the global optimum value in many cases, resulting in various parameter sets. As can be seen below, the optimum performance criterion values (column *objective*) can differ from the global optimum value in many cases, resulting in various parameter sets.
```{r, warning=FALSE}
summary(df_port) ```{r local4, warning=FALSE}
summary(resPORT)
``` ```
The existence of several local minima illustrates the importance of defining an appropriate starting point or of using a multi-start strategy or a global optimization strategy. The existence of several local minima illustrates the importance of defining an appropriate starting point or of using a multi-start strategy or a global optimization strategy.
# Global optimization # Global optimization
Global optimization is most often used when facing a complex response surface, with multiple local mimina. Global optimization is most often used when facing a complex response surface, with multiple local mimina.
...@@ -127,8 +151,10 @@ Here we use the following R implementation of some popular strategies: ...@@ -127,8 +151,10 @@ Here we use the following R implementation of some popular strategies:
* [hydroPSO: particle swarm](https://cran.r-project.org/package=hydroPSO) * [hydroPSO: particle swarm](https://cran.r-project.org/package=hydroPSO)
* [Rmalschains: memetic algorithms](https://cran.r-project.org/package=Rmalschains) * [Rmalschains: memetic algorithms](https://cran.r-project.org/package=Rmalschains)
## Differential Evolution ## Differential Evolution
```{r, warning=FALSE, results='hide'}
```{r optDE, warning=FALSE, results='hide', eval=FALSE}
optDE <- DEoptim::DEoptim(fn = OptimGR4J, optDE <- DEoptim::DEoptim(fn = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J, lower = lowerGR4J, upper = upperGR4J,
control = DEoptim::DEoptim.control(NP = 40, trace = 10)) control = DEoptim::DEoptim.control(NP = 40, trace = 10))
...@@ -136,36 +162,146 @@ optDE <- DEoptim::DEoptim(fn = OptimGR4J, ...@@ -136,36 +162,146 @@ optDE <- DEoptim::DEoptim(fn = OptimGR4J,
## Particle Swarm ## Particle Swarm
```{r, warning=FALSE, results='hide', message=FALSE}
```{r hydroPSO1, warning=FALSE, results='hide', message=FALSE, eval=FALSE}
# to install the package temporary removed from CRAN
# Rtools needed (windows : https://cran.r-project.org/bin/windows/Rtools/)
# install.packages("https://cran.r-project.org/src/contrib/Archive/hydroPSO/hydroPSO_0.5-1.tar.gz",
# repos = NULL, type = "source", dependencies = TRUE)
```
```{r hydroPSO2, warning=FALSE, results='hide', message=FALSE, eval=FALSE}
optPSO <- hydroPSO::hydroPSO(fn = OptimGR4J, optPSO <- hydroPSO::hydroPSO(fn = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J, lower = lowerGR4J, upper = upperGR4J,
control = list(write2disk = FALSE, verbose = FALSE)) control = list(write2disk = FALSE, verbose = FALSE))
``` ```
## MA-LS-Chains ## MA-LS-Chains
```{r, warning=FALSE, results='hide'}
```{r optMALS, warning=FALSE, results='hide', eval=FALSE}
optMALS <- Rmalschains::malschains(fn = OptimGR4J, optMALS <- Rmalschains::malschains(fn = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J, lower = lowerGR4J, upper = upperGR4J,
maxEvals = 2000) maxEvals = 2000)
``` ```
# Results # Results
As it can be seen in the table below, the four additional optimization strategies tested lead to very close optima. As it can be seen in the table below, the four additional optimization strategies tested lead to very close optima.
```{r, warning=FALSE, echo=FALSE} ```{r resGLOB1, warning=FALSE, echo=FALSE, eval=FALSE}
data.frame(Algo = c("airGR", "PORT", "DE", "PSO", "MA-LS"), resGLOB <- data.frame(Algo = c("airGR", "PORT", "DE", "PSO", "MA-LS"),
round(rbind( round(rbind(
OutputsCalib$ParamFinalR , OutputsCalib$ParamFinalR,
airGR::TransfoParam_GR4J(ParamIn = optPORT$par , Direction = "TR"), airGR::TransfoParam_GR4J(ParamIn = optPORT$par , Direction = "TR"),
airGR::TransfoParam_GR4J(ParamIn = as.numeric(optDE$optim$bestmem), Direction = "TR"), airGR::TransfoParam_GR4J(ParamIn = as.numeric(optDE$optim$bestmem), Direction = "TR"),
airGR::TransfoParam_GR4J(ParamIn = as.numeric(optPSO$par) , Direction = "TR"), airGR::TransfoParam_GR4J(ParamIn = as.numeric(optPSO$par) , Direction = "TR"),
airGR::TransfoParam_GR4J(ParamIn = optMALS$sol , Direction = "TR")), airGR::TransfoParam_GR4J(ParamIn = optMALS$sol , Direction = "TR")),
digits = 3)) digits = 3))
```
```{r resGLOB2, warning=FALSE, echo=FALSE}
resGLOB
``` ```
<!-- This is an expected result because the response surface for quadratic performance criteria of the **GR4J** model is generally sufficiently well defined in the transformed parameter space to allow using a local optimization strategy instead of a more time consuming global optimization strategy. --> <!-- This is an expected result because the response surface for quadratic performance criteria of the **GR4J** model is generally sufficiently well defined in the transformed parameter space to allow using a local optimization strategy instead of a more time consuming global optimization strategy. -->
# Multiobjective optimization
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.
First, the `OptimGR4J()` function previously used is modified to return two values.
```{r, warning=FALSE, results='hide', eval=FALSE}
InputsCrit_inv <- InputsCrit
InputsCrit_inv$transfo <- "inv"
MOptimGR4J <- function(i) {
if (algo == "caRamel") {
ParamOptim <- x[i, ]
}
## Transformation of the parameter set to real space
RawParamOptim <- airGR::TransfoParam_GR4J(ParamIn = ParamOptim,
Direction = "TR")
## Simulation given a parameter set
OutputsModel <- airGR::RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions,
Param = RawParamOptim)
## Computation of the value of the performance criteria
OutputsCrit1 <- airGR::ErrorCrit_KGE(InputsCrit = InputsCrit,
OutputsModel = OutputsModel,
verbose = FALSE)
## Computation of the value of the performance criteria
OutputsCrit2 <- airGR::ErrorCrit_KGE(InputsCrit = InputsCrit_inv,
OutputsModel = OutputsModel,
verbose = FALSE)
return(c(OutputsCrit1$CritValue, OutputsCrit2$CritValue))
}
```
## caRamel
caRamel is a multiobjective evolutionary algorithm combining the MEAS algorithm and the NGSA-II algorithm.
```{r, warning=FALSE, results='hide', eval=FALSE}
algo <- "caRamel"
optMO <- caRamel::caRamel(nobj = 2,
nvar = 4,
minmax = rep(TRUE, 2),
bounds = matrix(c(lowerGR4J, upperGR4J), ncol = 2),
func = MOptimGR4J,
popsize = 100,
archsize = 100,
maxrun = 15000,
prec = rep(1.e-3, 2),
carallel = FALSE,
graph = FALSE)
```
The algorithm returns parameter sets that describe the pareto front, illustrating the trade-off between overall good performance and good performance for low flow.
```{r, fig.width=6, fig.height=6, warning=FALSE}
ggplot() +
geom_point(aes(optMO$objectives[, 1], optMO$objectives[, 2])) +
coord_equal(xlim = c(0.4, 0.9), ylim = c(0.4, 0.9)) +
xlab("KGE") + ylab("KGE [1/Q]") +
theme_bw()
```
The parameter sets can be viewed in the parameter space, illustrating different populations.
```{r fig.height=6, fig.width=6, message=FALSE, warning=FALSE}
param_optMO <- apply(optMO$parameters, MARGIN = 1, FUN = function(x) {
airGR::TransfoParam_GR4J(x, Direction = "TR")
})
GGally::ggpairs(data.frame(t(param_optMO)), diag = NULL) + theme_bw()
```
```{r fig.height=6, fig.width=12, message=FALSE, warning=FALSE}
RunOptions$Outputs_Sim <- "Qsim"
run_optMO <- apply(optMO$parameters, MARGIN = 1, FUN = function(x) {
airGR::RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions,
Param = x)
}$Qsim)
run_optMO <- data.frame(run_optMO)
ggplot() +
geom_line(aes(x = as.POSIXct(InputsModel$DatesR[Ind_Run]),
y = run_optMO$X1)) +
geom_line(aes(x = as.POSIXct(InputsModel$DatesR[Ind_Run]),
y = run_optMO$X54),
colour = "darkred") +
scale_x_datetime(limits = c(as.POSIXct("1998-01-01"), NA)) +
ylab("Discharge [mm/d]") + xlab("Date") +
scale_y_sqrt() +
theme_bw()
```
vignettes/V02.2_param_mcmc.Rmd
View file @ 4bccfee2
...@@ -15,9 +15,8 @@ library(airGR) ...@@ -15,9 +15,8 @@ library(airGR)
library(coda) library(coda)
library(FME) library(FME)
library(ggmcmc) library(ggmcmc)
library(dplyr)
# source("airGR.R")
set.seed(123) set.seed(123)
load(system.file("vignettesData/vignetteParamMCMC.rda", package = "airGR"))
``` ```
...@@ -30,13 +29,21 @@ In this vignette, we give an example of parameter estimation within a Bayesian M ...@@ -30,13 +29,21 @@ In this vignette, we give an example of parameter estimation within a Bayesian M
We use the **GR4J** model and we assume that the R global environment contains data and functions from the **airGR** [Get Started](V01_get_started.html) vignette. We use the **GR4J** model and we assume that the R global environment contains data and functions from the **airGR** [Get Started](V01_get_started.html) vignette.
```{r, warning=FALSE, fig.keep='none', results='hide', fig.height=10, fig.width=10, eval=TRUE, echo=FALSE, message=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) <!-- example("Calibration_Michel", echo = FALSE, ask = FALSE) -->
``` <!-- ``` -->
```{r, echo=TRUE, eval=FALSE} ```{r, echo=TRUE, eval=FALSE, eval=FALSE}
example("Calibration_Michel") example("Calibration_Michel")
``` ```
In order for the `RunModel_*()` functions to run faster during the parameter estimation process, it is recommended that the outputs contain only the simulated flows (see the `Outputs_Sim` argument in the `CreateRunOptions()` help page).
```{r, results='hide', eval=FALSE}
RunOptions <- airGR::CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
IndPeriod_Run = Ind_Run,
Outputs_Sim = "Qsim")
```
Please refer to the [2.1 Plugging in new calibration](V02.1_param_optim.html) vignette for explanations on how to plug in a parameter estimation algorithm to **airGR**. Please refer to the [2.1 Plugging in new calibration](V02.1_param_optim.html) vignette for explanations on how to plug in a parameter estimation algorithm to **airGR**.
...@@ -47,30 +54,30 @@ Please note that this vignette is only for illustration purposes and does not pr ...@@ -47,30 +54,30 @@ Please note that this vignette is only for illustration purposes and does not pr
We show how to use the DRAM algorithm for SLS Bayesian inference, with the `modMCMC()` function of the [FME](https://cran.r-project.org/package=FME) package. We show how to use the DRAM algorithm for SLS Bayesian inference, with the `modMCMC()` function of the [FME](https://cran.r-project.org/package=FME) package.
First, we need to define a function that returns twice the opposite of the log-likelihood for a given parameter set. First, we need to define a function that returns twice the opposite of the log-likelihood for a given parameter set.
Nota: in the `RunAirGR4J()` function, the computation of the log-likelihood is simplified in order to ensure a good computing performance. It corresponds to a translation of the two following lines. Nota: in the `LogLikeGR4J()` function, the computation of the log-likelihood is simplified in order to ensure a good computing performance. It corresponds to a translation of the two following lines.
```{r, echo=TRUE, eval=FALSE}
Likelihood <- sum((ObsY - ModY)^2, na.rm = TRUE)^(-sum(!is.na(ObsY)) / 2)
LL <- -2 * log(Likelihood)
```{r, echo=TRUE, eval=FALSE, purl=FALSE}
Likelihood <- sum((ObsY - ModY)^2, na.rm = TRUE)^(-sum(!is.na(ObsY)) / 2)
LogLike <- -2 * log(Likelihood)
``` ```
In our simplified setting of Gaussian likelihood with measurement error integrated out, the log of the sum of squared error is related to the log-likelihood. In our simplified setting of Gaussian likelihood with measurement error integrated out, the log of the sum of squared error is related to the log-likelihood.
Note that we do not use here any discharge transformation, although applying Box-Cox transformation is quite popular in hydrological modelling. Note that we do not use here any discharge transformation, although applying Box-Cox transformation is quite popular in hydrological modelling.
```{r, results='hide'} ```{r, results='hide', eval=FALSE}
RunAirGR4J <- function(Param_Optim) { LogLikeGR4J <- function(ParamOptim) {
## Transformation to real space ## Transformation to real space
Param_Optim_Vre <- airGR::TransfoParam_GR4J(ParamIn = Param_Optim, RawParamOptim <- airGR::TransfoParam_GR4J(ParamIn = ParamOptim,
Direction = "TR") Direction = "TR")
## Simulation given a parameter set ## Simulation given a parameter set
OutputsModel <- airGR::RunModel_GR4J(InputsModel = InputsModel, OutputsModel <- airGR::RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions, RunOptions = RunOptions,
Param = Param_Optim_Vre) Param = RawParamOptim)
## Computation of the log-likelihood: N * log(SS) ## Computation of the log-likelihood: N * log(SS)
ObsY <- InputsCrit$Qobs ObsY <- InputsCrit$Obs
ModY <- OutputsModel$Qsim ModY <- OutputsModel$Qsim
LL <- sum(!is.na(ObsY)) * log(sum((ObsY - ModY)^2, na.rm = TRUE)) LogLike <- sum(!is.na(ObsY)) * log(sum((ObsY - ModY)^2, na.rm = TRUE))
} }
``` ```
...@@ -78,12 +85,12 @@ RunAirGR4J <- function(Param_Optim) { ...@@ -78,12 +85,12 @@ RunAirGR4J <- function(Param_Optim) {
## Estimation of the best-fit parameters as a starting point ## Estimation of the best-fit parameters as a starting point
We start by using the PORT optimization routine to estimate the best-fit parameters. We start by using the PORT optimization routine to estimate the best-fit parameters.
```{r, results='hide'} ```{r, results='hide', eval=FALSE}
optPORT <- stats::nlminb(start = c(4.1, 3.9, -0.9, -8.7), optPORT <- stats::nlminb(start = c(4.1, 3.9, -0.9, -8.7),
objective = RunAirGR4J, objective = LogLikeGR4J,
lower = rep(-9.9, times = 4), upper = rep(9.9, times = 4), lower = rep(-9.9, times = 4), upper = rep(9.9, times = 4),
control = list(trace = 1)) control = list(trace = 1))
IniParam <- optPORT$par iniParPORT <- optPORT$par
``` ```
## Running 3 chains for convergence assessment ## Running 3 chains for convergence assessment
...@@ -95,24 +102,26 @@ Nota: in this example, there are relatively few iterations (2000), in order to l ...@@ -95,24 +102,26 @@ Nota: in this example, there are relatively few iterations (2000), in order to l
With the DRAM algorithm, the covariance of the proposal is updated every 100 runs and delayed rejection is applied. With the DRAM algorithm, the covariance of the proposal is updated every 100 runs and delayed rejection is applied.
```{r, results='hide'} ```{r, results='hide', eval=FALSE}
ListIniParam <- data.frame(Chain1 = IniParam, Chain2 = IniParam, Chain3 = IniParam, iniParPORT <- data.frame(Chain1 = iniParPORT,
row.names = paste0("X", 1:4)) Chain2 = iniParPORT,
ListIniParam <- sweep(ListIniParam, MARGIN = 2, STATS = c(1, 0.9, 1.1), FUN = "*") Chain3 = iniParPORT,
ListIniParam[ListIniParam < -9.9] <- -9.9 row.names = paste0("X", 1:4))
ListIniParam[ListIniParam > +9.9] <- +9.9 iniParPORT <- sweep(iniParPORT, MARGIN = 2, STATS = c(1, 0.9, 1.1), FUN = "*")
iniParPORT[iniParPORT < -9.9] <- -9.9
mcmcDRAM <- apply(ListIniParam, 2, function(iListIniParam) { iniParPORT[iniParPORT > +9.9] <- +9.9
FME::modMCMC(f = RunAirGR4J,
p = iListIniParam, mcmcDRAM <- apply(iniParPORT, MARGIN = 2, FUN = function(iIniParPORT) {
FME::modMCMC(f = LogLikeGR4J,
p = iIniParPORT,
lower = rep(-9.9, times = 4), ## lower bounds for GR4J lower = rep(-9.9, times = 4), ## lower bounds for GR4J
upper = rep(+9.9, times = 4), ## upper bounds for GR4J upper = rep(+9.9, times = 4), ## upper bounds for GR4J
niter = 2000, niter = 2000,
jump = 0.01, jump = 0.01,
outputlength = 2000, outputlength = 2000,
burninlength = 0, burninlength = 0,
updatecov = 100, ## Adaptative Metropolis updatecov = 100, ## adaptative Metropolis (AM)
ntrydr = 2) ## Delayed Rejection ntrydr = 2) ## delayed rejection (RD)
}) })
``` ```
...@@ -126,12 +135,14 @@ The result will be better with more iterations than 2000. As we kept the iterati ...@@ -126,12 +135,14 @@ The result will be better with more iterations than 2000. As we kept the iterati
Note that we rescale model parameter sets of the GR4J model from the transformed space to the real space. Note that we rescale model parameter sets of the GR4J model from the transformed space to the real space.
```{r, results='hide'} ```{r, results='hide', eval=FALSE}
MultDRAM <- coda::as.mcmc(lapply(mcmcDRAM, function(x) { multDRAM <- coda::as.mcmc.list(lapply(mcmcDRAM, FUN = function(x) {
coda::as.mcmc(airGR::TransfoParam_GR4J(as.matrix(x$pars), Direction = "TR")) coda::as.mcmc(airGR::TransfoParam_GR4J(as.matrix(x$pars), Direction = "TR"))
})) }))
GelRub <- coda::gelman.diag(MultDRAM, autoburnin = TRUE)$mpsrf gelRub <- coda::gelman.diag(multDRAM, autoburnin = TRUE)$mpsrf
GelRub ```
```{r}
gelRub
``` ```
In addition, graphical tools can be used, with for example the [ggmcmc](https://cran.r-project.org/package=ggmcmc) package. In addition, graphical tools can be used, with for example the [ggmcmc](https://cran.r-project.org/package=ggmcmc) package.
...@@ -141,21 +152,20 @@ In addition, graphical tools can be used, with for example the [ggmcmc](https:// ...@@ -141,21 +152,20 @@ In addition, graphical tools can be used, with for example the [ggmcmc](https://
First, the evolution of the Markov chains can be seen with a traceplot: First, the evolution of the Markov chains can be seen with a traceplot:
```{r, fig.width=6, fig.height=9, warning=FALSE} ```{r, fig.width=6, fig.height=9, warning=FALSE}
ParamDRAM <- ggmcmc::ggs(MultDRAM) ## to convert objet for using by all ggs_* graphical functions parDRAM <- ggmcmc::ggs(multDRAM) ## to convert object for using by all ggs_* graphical functions
ggmcmc::ggs_traceplot(ParamDRAM) ggmcmc::ggs_traceplot(parDRAM)
``` ```
The posterior density for each parameter can then be visualised: The posterior density for each parameter can then be visualised:
```{r, fig.width=6, fig.height=9, warning=FALSE} ```{r, fig.width=6, fig.height=9, warning=FALSE}
ParamDRAM_burn <- dplyr::filter(ParamDRAM, 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(ParamDRAM_burn) ggmcmc::ggs_density(burnParDRAM)
``` ```
Finally, a paired plot can be useful to look at the correlation between parameters. Finally, a paired plot can be useful to look at the correlation between parameters.
Highly correlated parameters can be seen as an indication for a structural deficit of the model used. Highly correlated parameters can be seen as an indication for a structural deficit of the model used.
```{r, fig.width=6, fig.height=6} ```{r, fig.width=6, fig.height=6, results='hide'}
ggmcmc::ggs_pairs(ParamDRAM_burn, lower = list(continuous = "density")) ggmcmc::ggs_pairs(burnParDRAM, lower = list(continuous = "density"))
``` ```
## Exploring further possibilities ## Exploring further possibilities
......
vignettes/V03_param_sets_GR4J.Rmd
View file @ 4bccfee2
--- ---
title: "Generalist parameter sets for the GR4J model" title: "Generalist parameter sets for the GR4J model"
author: "Olivier Delaigue, Guillaume Thirel"
bibliography: V00_airgr_ref.bib bibliography: V00_airgr_ref.bib
output: rmarkdown::html_vignette output: rmarkdown::html_vignette
vignette: > vignette: >
...@@ -36,7 +37,6 @@ In this vignette, we show a method using a known parameter set that can be used ...@@ -36,7 +37,6 @@ In this vignette, we show a method using a known parameter set that can be used
We load an example data set from the package and the GR4J parameter sets. We load an example data set from the package and the GR4J parameter sets.
```{r, warning=FALSE} ```{r, warning=FALSE}
## loading catchment data ## loading catchment data
data(L0123001) data(L0123001)
...@@ -70,7 +70,6 @@ The calibration period has been defined from **1990-01-01** to **1990-02-28**, a ...@@ -70,7 +70,6 @@ The calibration period has been defined from **1990-01-01** to **1990-02-28**, a
As a consequence, in this example the calibration period is very short, less than 6 months. As a consequence, in this example the calibration period is very short, less than 6 months.
```{r, warning=FALSE, include=TRUE} ```{r, warning=FALSE, include=TRUE}
## preparation of the InputsModel object ## preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR, InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E) Precip = BasinObs$P, PotEvap = BasinObs$E)
...@@ -83,11 +82,11 @@ Ind_Cal <- seq(which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="01/01/ ...@@ -83,11 +82,11 @@ Ind_Cal <- seq(which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="01/01/
## preparation of the RunOptions object for the calibration period ## preparation of the RunOptions object for the calibration period
RunOptions_Cal <- CreateRunOptions(FUN_MOD = RunModel_GR4J, RunOptions_Cal <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Cal) InputsModel = InputsModel, IndPeriod_Run = Ind_Cal)
## efficiency criterion: Nash-Sutcliffe Efficiency ## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit_Cal <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, InputsCrit_Cal <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions_Cal, Qobs = BasinObs$Qmm[Ind_Cal]) RunOptions = RunOptions_Cal, Obs = BasinObs$Qmm[Ind_Cal])
## ---- validation step ## ---- validation step
...@@ -102,8 +101,7 @@ RunOptions_Val <- CreateRunOptions(FUN_MOD = RunModel_GR4J, ...@@ -102,8 +101,7 @@ RunOptions_Val <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
## efficiency criterion (Nash-Sutcliffe Efficiency) on the validation period ## efficiency criterion (Nash-Sutcliffe Efficiency) on the validation period
InputsCrit_Val <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, InputsCrit_Val <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions_Val, Qobs = BasinObs$Qmm[Ind_Val]) RunOptions = RunOptions_Val, Obs = BasinObs$Qmm[Ind_Val])
``` ```
# Calibration of the GR4J model with the generalist parameter sets # Calibration of the GR4J model with the generalist parameter sets
...@@ -134,7 +132,8 @@ Param_Best <- unlist(Param_Sets_GR4J[which.max(OutputsCrit_Loop), ]) ...@@ -134,7 +132,8 @@ Param_Best <- unlist(Param_Sets_GR4J[which.max(OutputsCrit_Loop), ])
Param_Best Param_Best
## validation ## 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) OutputsCrit_Val <- ErrorCrit_NSE(InputsCrit = InputsCrit_Val, OutputsModel = OutputsModel_Val)
``` ```
...@@ -142,7 +141,7 @@ Now we can compute the Nash-Sutcliffe Efficiency criterion on the validation per ...@@ -142,7 +141,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. As seen above, the Michel's calibration algorithm is based on a local search procedure.
...@@ -156,22 +155,22 @@ By default, the predefined grid used by the `Calibration_Michel()` function cont ...@@ -156,22 +155,22 @@ By default, the predefined grid used by the `Calibration_Michel()` function cont
```{r, warning=FALSE, message=FALSE} ```{r, warning=FALSE, message=FALSE}
CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel) CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel)
``` ```
```{r, warning=FALSE, message=FALSE, include=FALSE} ```{r, warning=FALSE, message=FALSE, include=FALSE}
## calibration ## calibration
OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions_Cal, OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
InputsCrit = InputsCrit_Cal, CalibOptions = CalibOptions, InputsCrit = InputsCrit_Cal, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J, FUN_CRIT = ErrorCrit_NSE, FUN_MOD = RunModel_GR4J,
FUN_CALIB = Calibration_Michel) FUN_CALIB = Calibration_Michel)
OutputsModel_Cal <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions_Cal, OutputsModel_Cal <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
Param = OutputsCalib$ParamFinalR, FUN = RunModel_GR4J) Param = OutputsCalib$ParamFinalR, FUN_MOD = RunModel_GR4J)
OutputsCrit_Cal <- ErrorCrit_NSE(InputsCrit = InputsCrit_Cal, OutputsModel = OutputsModel_Cal) OutputsCrit_Cal <- ErrorCrit_NSE(InputsCrit = InputsCrit_Cal, OutputsModel = OutputsModel_Cal)
## validation ## validation
OutputsModel_Val <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions_Val, Param = OutputsCalib$ParamFinalR) OutputsModel_Val <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions_Val,
Param = OutputsCalib$ParamFinalR)
OutputsCrit_Val <- ErrorCrit_NSE(InputsCrit = InputsCrit_Val, OutputsModel = OutputsModel_Val) OutputsCrit_Val <- ErrorCrit_NSE(InputsCrit = InputsCrit_Val, OutputsModel = OutputsModel_Val)
``` ```
...@@ -199,10 +198,10 @@ CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibrat ...@@ -199,10 +198,10 @@ CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibrat
## calibration ## calibration
OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions_Cal, OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
InputsCrit = InputsCrit_Cal, CalibOptions = CalibOptions, InputsCrit = InputsCrit_Cal, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J, FUN_CRIT = ErrorCrit_NSE, FUN_MOD = RunModel_GR4J,
FUN_CALIB = Calibration_Michel) FUN_CALIB = Calibration_Michel)
OutputsModel_Cal <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions_Cal, OutputsModel_Cal <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
Param = OutputsCalib$ParamFinalR, FUN = RunModel_GR4J) Param = OutputsCalib$ParamFinalR, FUN_MOD = RunModel_GR4J)
OutputsCrit_Cal <- ErrorCrit_NSE(InputsCrit = InputsCrit_Cal, OutputsModel = OutputsModel_Cal) OutputsCrit_Cal <- ErrorCrit_NSE(InputsCrit = InputsCrit_Cal, OutputsModel = OutputsModel_Cal)
......
vignettes/V04_cemaneige_hysteresis.Rmd 0 → 100644
View file @ 4bccfee2
---
title: "Using satellite snow cover area data for calibrating and improving CemaNeige"
author: "Guillaume Thirel, Olivier Delaigue"
bibliography: V00_airgr_ref.bib
output: rmarkdown::html_vignette
vignette: >
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{Using CemaNeige with Hysteresis}
%\VignetteEncoding{UTF-8}
---
```{r, warning=FALSE, include=FALSE}
library(airGR)
load(system.file("vignettesData/vignetteCNHysteresis.rda", package = "airGR"))
```
# Introduction
## Scope
Rainfall-runoff models that include a snow accumulation and melt module are still often calibrated using only discharge observations. This can result in poor snow modeling as the swnow module parameters can distorted in order to allow skilful discharge simulation.
After the work of @riboust_revisiting_2019, we propose now in **airGR** an improved version of the degree-day CemaNeige snow and accumulation module. This new version is based on a more accurate representation of the relationship that exists at the basin scale between the Snow Water Equivalent (SWE) and the Snow Cover Area (SCA). To do so, a linear SWE-SCA hysteresis, which represents the fact that snow accumulation is rather homogeneous on the basin and snow melt is more heterogeneous, was implemented.
This new CemaNeige version has two more parameters to calibrate. It also presents the advantage of allowing using satellite snow data to constrain the calibration in addition to discharge.
@riboust_revisiting_2019 show that while the simulated discharge is not significantly improved, the snow simulation is much improved. In addition, they show that the model is more robust (i.e. transferable in time) in terms of discharge, which has many implications for climate change impact studies.
The configuration that was identified as optimal by @riboust_revisiting_2019 includes a CemaNeige module run on 5 elevation bands and an objective function determine by a composite function of KGE' calculated on discharge (75 % weight) and KGE' calculated on each elevation band (5 % for each).
In this page, we show how to use and calibrate this new CemaNeige version.
## Data preparation
We load an example data set from the package. Please note that this data set includes MODIS data that was pre-calculated for 5 elevation bands and for which days with few data (more than 40 % cloud coverage) were assigned as missing values.
## loading catchment data
```{r, warning=FALSE, eval=FALSE}
data(X0310010)
summary(BasinObs)
```
## Object model preparation
We assume that the R global environment contains data and functions from the [Get Started](V01_get_started.html) vignette.
The calibration period has been defined from 2000-09-01 to 2005-08-31, and the validation period from 2005-09-01 to 2010-07-31. CemaNeige will be used in coupling with GR4J in this vignette.
```{r, warning=FALSE, eval=FALSE}
## preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_CemaNeigeGR4J,
DatesR = BasinObs$DatesR, Precip = BasinObs$P,
PotEvap = BasinObs$E, TempMean = BasinObs$T,
ZInputs = median(BasinInfo$HypsoData),
HypsoData = BasinInfo$HypsoData, NLayers = 5)
## ---- calibration step
## calibration period selection
Ind_Cal <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "2000-09-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "2005-08-31"))
## ---- validation step
## validation period selection
Ind_Val <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "2005-09-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "2010-07-31"))
```
# Calibration and evaluation of the new CemaNeige module
In order to use the Linear Hysteresis, a new argument (`IsHyst`) is added in the `CreateRunOptions()` and `CreateCalibOptions()` functions and has to be set to `TRUE`.
```{r, warning=FALSE, eval=FALSE}
## preparation of the RunOptions object for the calibration period
RunOptions_Cal <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Cal,
IsHyst = TRUE)
## preparation of the RunOptions object for the validation period
RunOptions_Val <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Val,
IsHyst = TRUE)
## preparation of the CalibOptions object
CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
FUN_CALIB = Calibration_Michel,
IsHyst = TRUE)
```
In order to calibrate and assess the model performance, we will follow the recommendations of @riboust_revisiting_2019. This is now possible in **airGR** with the added functionality that permits calculating composite criteria by combining different metrics.
```{r, warning=FALSE, eval=FALSE}
## efficiency criterion: 75 % KGE'(Q) + 5 % KGE'(SCA) on each of the 5 layers
InputsCrit_Cal <- CreateInputsCrit(FUN_CRIT = rep("ErrorCrit_KGE2", 6),
InputsModel = InputsModel, RunOptions = RunOptions_Cal,
Obs = list(BasinObs$Qmm[Ind_Cal],
BasinObs$SCA1[Ind_Cal],
BasinObs$SCA2[Ind_Cal],
BasinObs$SCA3[Ind_Cal],
BasinObs$SCA4[Ind_Cal],
BasinObs$SCA5[Ind_Cal]),
VarObs = list("Q", "SCA", "SCA", "SCA", "SCA", "SCA"),
Weights = list(0.75, 0.05, 0.05, 0.05, 0.05, 0.05))
InputsCrit_Val <- CreateInputsCrit(FUN_CRIT = rep("ErrorCrit_KGE2", 6),
InputsModel = InputsModel, RunOptions = RunOptions_Val,
Obs = list(BasinObs$Qmm[Ind_Val],
BasinObs$SCA1[Ind_Val],
BasinObs$SCA2[Ind_Val],
BasinObs$SCA3[Ind_Val],
BasinObs$SCA4[Ind_Val],
BasinObs$SCA5[Ind_Val]),
VarObs = list("Q", "SCA", "SCA", "SCA", "SCA", "SCA"),
Weights = list(0.75, 0.05, 0.05, 0.05, 0.05, 0.05))
```
We can now calibrate the model.
```{r, warning=FALSE, eval=FALSE}
## calibration
OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
InputsCrit = InputsCrit_Cal, CalibOptions = CalibOptions,
FUN_MOD = RunModel_CemaNeigeGR4J,
FUN_CALIB = Calibration_Michel)
```
Now we can run it on the calibration period and assess it.
```{r, warning=FALSE, message=FALSE, eval=FALSE}
## run on the calibration period
OutputsModel_Cal <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
RunOptions = RunOptions_Cal,
Param = OutputsCalib$ParamFinalR)
## evaluation
OutputsCrit_Cal <- ErrorCrit(InputsCrit = InputsCrit_Cal, OutputsModel = OutputsModel_Cal)
```
Find below the performance of the model over the calibration period.
```{r, warning=FALSE}
str(OutputsCrit_Cal, max.level = 2)
```
Now we can run the model on the validation period and assess it.
```{r, warning=FALSE, message=FALSE, eval=FALSE}
## run on the validation period
OutputsModel_Val <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
RunOptions = RunOptions_Val,
Param = OutputsCalib$ParamFinalR)
## evaluation
OutputsCrit_Val <- ErrorCrit(InputsCrit = InputsCrit_Val, OutputsModel = OutputsModel_Val)
```
Find below the performance of the model over the validation period.
```{r, warning=FALSE}
str(OutputsCrit_Val, max.level = 2)
```
# Comparison with the performance of the initial CemaNeige version
Here we use the same InputsModel object and calibration and validation periods. However, we have to redefine the way we run the model (`RunOptions` argument), calibrate and assess it (`InputsCrit` argument). The objective function is only based on KGE'(Q). Note how we set the `IsHyst` argument to `FALSE` in the `CreateRunOptions()` and the `CreateCalibOptions()` functions.
```{r, warning=FALSE, eval=FALSE}
## preparation of RunOptions object
RunOptions_Cal_NoHyst <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
InputsModel = InputsModel,
IndPeriod_Run = Ind_Cal,
IsHyst = FALSE)
RunOptions_Val_NoHyst <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
InputsModel = InputsModel,
IndPeriod_Run = Ind_Val,
IsHyst = FALSE)
InputsCrit_Cal_NoHyst <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE2,
InputsModel = InputsModel,
RunOptions = RunOptions_Cal_NoHyst,
Obs = BasinObs$Qmm[Ind_Cal], VarObs = "Q")
## preparation of CalibOptions object
CalibOptions_NoHyst <- CreateCalibOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
FUN_CALIB = Calibration_Michel,
IsHyst = FALSE)
```
We can now calibrate the model.
```{r, warning=FALSE, eval=FALSE}
## calibration
OutputsCalib_NoHyst <- Calibration(InputsModel = InputsModel,
InputsCrit = InputsCrit_Cal_NoHyst,
RunOptions = RunOptions_Cal_NoHyst,
CalibOptions = CalibOptions_NoHyst,
FUN_MOD = RunModel_CemaNeigeGR4J,
FUN_CALIB = Calibration_Michel)
```
And run it over the calibration and validation periods.
```{r, warning=FALSE, eval=FALSE}
OutputsModel_Cal_NoHyst <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
RunOptions = RunOptions_Cal_NoHyst,
Param = OutputsCalib_NoHyst$ParamFinalR)
OutputsModel_Val_NoHyst <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
RunOptions = RunOptions_Val_NoHyst,
Param = OutputsCalib_NoHyst$ParamFinalR)
```
In order to assess the model performance over the two periods, we will use the InputsCrit objects prepared before, which allow assessing also the performance in terms of snow simulation.
```{r, warning=FALSE, message=FALSE, eval=FALSE}
OutputsCrit_Cal_NoHyst <- ErrorCrit(InputsCrit = InputsCrit_Cal,
OutputsModel = OutputsModel_Cal_NoHyst)
OutputsCrit_Val_NoHyst <- ErrorCrit(InputsCrit = InputsCrit_Val,
OutputsModel = OutputsModel_Val_NoHyst)
```
We can check the performance over the calibration and the validation period.
```{r, warning=FALSE}
str(OutputsCrit_Cal_NoHyst, max.level = 2)
str(OutputsCrit_Val_NoHyst, max.level = 2)
```
We can see above that the performance of the initial model is slightly better than the new one over the calibration period in terms of discharge, but also that the new version calibrated using SCA provides better performance in terms of snow.
However, over the validation period, we see that the discharge simulated by the new version brings better performance (in addition to improved SCA also). This shows the interests of the combined use of a linear hysteresis and of SCA data for calibration in CemaNeige.
# References
vignettes/V05_sd_model.Rmd 0 → 100644
View file @ 4bccfee2
---
title: "Simulated vs observed upstream flows in calibration of semi-distributed GR4J model"
author: "David Dorchies"
bibliography: V00_airgr_ref.bib
output: rmarkdown::html_vignette
vignette: >
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{Simulated vs observed upstream flows in calibration of semi-distributed GR4J model}
%\VignetteEncoding{UTF-8}
---
```{r, include=FALSE, fig.keep='none', results='hide'}
library(airGR)
options(digits = 3)
```
# Introduction
## Scope
The **airGR** package implements semi-distributed model capabilities using a lag model between subcatchments. It allows to chain together several lumped models as well as integrating anthropogenic influence such as reservoirs or withdrawals.
Here we explain how to implement the semi-distribution with **airGR**. For everyday use, however, it is easier to use the [**airGRiwrm**](https://cran.r-project.org/package=airGRiwrm) package.
`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 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)
```
We use an example data set from the package that unfortunately contains data for only one catchment.
```{r, warning=FALSE}
## loading catchment data
data(L0123001)
```
Let's imagine that this catchment of 360 km² is divided into 2 subcatchments:
- An upstream subcatchment of 180 km²
- 100 km downstream another subcatchment of 180 km²
We consider that meteorological data are homogeneous on the whole catchment, so we use the same pluviometry `BasinObs$P` and the same evapotranspiration `BasinObs$E` for the 2 subcatchments.
For the observed flow at the downstream outlet, we generate it with the assumption that the upstream flow arrives at downstream with a constant delay of 2 days.
```{r}
QObsDown <- (BasinObs$Qmm + c(0, 0, BasinObs$Qmm[1:(length(BasinObs$Qmm)-2)])) / 2
options(digits = 5)
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.
```{r}
InputsModelUp <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1999-12-31"))
RunOptionsUp <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModelUp,
IndPeriod_WarmUp = NULL, IndPeriod_Run = Ind_Run,
IniStates = NULL, IniResLevels = NULL)
# 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],
transfo = "sqrt")
CalibOptionsUp <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel)
OutputsCalibUp <- Calibration_Michel(InputsModel = InputsModelUp, RunOptions = RunOptionsUp,
InputsCrit = InputsCritUp, CalibOptions = CalibOptionsUp,
FUN_MOD = RunModel_GR4J)
```
And see the result of the simulation:
```{r}
OutputsModelUp <- RunModel_GR4J(InputsModel = InputsModelUp, RunOptions = RunOptionsUp,
Param = OutputsCalibUp$ParamFinalR)
```
# 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(
FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E,
Qupstream = matrix(BasinObs$Qmm, 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²]
)
```
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_KGE, InputsModel = InputsModelDown1,
RunOptions = RunOptionsDown,
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
OutputsCalibDown1 <- Calibration_Michel(InputsModel = InputsModelDown1,
RunOptions = RunOptionsDown,
InputsCrit = InputsCritDown,
CalibOptions = CalibOptionsDown,
FUN_MOD = RunModel_GR4J)
```
`RunModel` is run in order to automatically combine GR4J and Lag models.
```{r}
OutputsModelDown1 <- RunModel(InputsModel = InputsModelDown2,
RunOptions = RunOptionsDown,
Param = OutputsCalibDown1$ParamFinalR,
FUN_MOD = RunModel_GR4J)
```
Performance of the model validation is then:
```{r}
KGE_down1 <- ErrorCrit_KGE(InputsCritDown, OutputsModelDown1)
```
## 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:
```{r}
OutputsCalibDown2 <- Calibration_Michel(InputsModel = InputsModelDown2,
RunOptions = RunOptionsDown,
InputsCrit = InputsCritDown,
CalibOptions = CalibOptionsDown,
FUN_MOD = RunModel_GR4J)
ParamDown2 <- OutputsCalibDown2$ParamFinalR
```
## Calibration with upstream simulated flow and parameter regularisation
The regularisation follow the method proposed by @delavenne_regularization_2019.
As a priori parameter set, we use the calibrated parameter set of the upstream catchment and the theoretical velocity:
```{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}
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],
OutputsCalibDown3$ParamFinalR[1]),
ncol = 1,
dimnames = list(c("theoretical",
"calibrated with observed 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 with the velocity as extra parameter:
```{r}
OutputsModelDownTheo <- RunModel(InputsModel = InputsModelDown2,
RunOptions = RunOptionsDown,
Param = ParamDownTheo,
FUN_MOD = RunModel_GR4J)
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,
KGE_down1$CritValue,
OutputsCalibDown2$CritFinal,
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. 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
vignettes/fig/airGR_features_diagram.svg 0 → 100644
View file @ 4bccfee2
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!-- Created with Inkscape (http://www.inkscape.org/) -->
<svg xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:cc="http://creativecommons.org/ns#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:svg="http://www.w3.org/2000/svg" xmlns="http://www.w3.org/2000/svg" xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd" xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape" id="svg2" version="1.1" inkscape:version="0.91 r13725" xml:space="preserve" width="900" height="675" viewBox="0 0 900 675" sodipodi:docname="airGR_features.svg"><metadata id="metadata8"><rdf:RDF><cc:Work rdf:about=""><dc:format>image/svg+xml</dc:format><dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/><dc:title/></cc:Work></rdf:RDF></metadata><defs id="defs6"><clipPath clipPathUnits="userSpaceOnUse" id="clipPath16"><path d="m 0,1.2207e-4 720,0 0,539.99999793 -720,0 L 0,1.2207e-4 Z" id="path18" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath><clipPath clipPathUnits="userSpaceOnUse" id="clipPath26"><path d="M 3.08e-7,1.2207e-4 720,1.2207e-4 720,540.00012 3.08e-7,540.00012 l 0,-539.99999793 z" id="path28" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath><clipPath clipPathUnits="userSpaceOnUse" id="clipPath102"><path d="M -2.645e-6,540 720,540 720,1.2207e-4 -6.582e-6,1.2207e-4" id="path104" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath><clipPath clipPathUnits="userSpaceOnUse" id="clipPath116"><path d="M -1.161e-6,540 720,540 720,1.2207e-4 l -720.000001161,0" id="path118" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath><clipPath clipPathUnits="userSpaceOnUse" id="clipPath132"><path d="M -2.728e-6,540 720,540 720,1.2207e-4 l -720.000006665,0" id="path134" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath><clipPath clipPathUnits="userSpaceOnUse" id="clipPath178"><path d="M -1.52e-6,540 720,540 720,6.104e-6 -1.52e-6,6.104e-6" id="path180" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath></defs><sodipodi:namedview pagecolor="#ffffff" bordercolor="#666666" borderopacity="1" objecttolerance="10" gridtolerance="10" guidetolerance="10" inkscape:pageopacity="0" inkscape:pageshadow="2" inkscape:window-width="1920" inkscape:window-height="1132" id="namedview4" showgrid="false" inkscape:zoom="0.98890193" inkscape:cx="333.64663" inkscape:cy="425.84952" inkscape:window-x="-8" inkscape:window-y="-8" inkscape:window-maximized="1" inkscape:current-layer="g10"/><g id="g10" inkscape:groupmode="layer" inkscape:label="airGR_features" transform="matrix(1.25,0,0,-1.25,0,675)"><g id="g12"><g id="g14" clip-path="url(#clipPath16)"><path d="m 0,6.104e-6 720,0 L 720,540.00001 l -720,0 L 0,6.104e-6 Z" style="fill:#ffffff;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path20" inkscape:connector-curvature="0"/></g></g><g id="g22"><g id="g24" clip-path="url(#clipPath26)"><path d="M 391.5,540 390.85,-2.1562" style="fill:none;stroke:#7f7f7f;stroke-width:2;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:8, 6;stroke-dashoffset:0;stroke-opacity:1" id="path30" inkscape:connector-curvature="0"/></g></g><path d="m 308.97,304.02 c 0,3.65 2.96,6.61 6.62,6.61 l 151.19,0 c 3.66,0 6.62,-2.96 6.62,-6.61 l 0,-26.46 c 0,-3.65 -2.96,-6.62 -6.62,-6.62 l -151.19,0 c -3.66,0 -6.62,2.97 -6.62,6.62 l 0,26.46 z" style="fill:#00b0f0;fill-opacity:0.2;fill-rule:evenodd;stroke:none" id="path32" inkscape:connector-curvature="0"/><g id="g34"><path d="m 308.97,304.02 c 0,3.65 2.96,6.61 6.62,6.61 l 151.19,0 c 3.66,0 6.62,-2.96 6.62,-6.61 l 0,-26.46 c 0,-3.65 -2.96,-6.62 -6.62,-6.62 l -151.19,0 c -3.66,0 -6.62,2.97 -6.62,6.62 l 0,26.46 z" style="fill:none;stroke:#00b0f0;stroke-width:3;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:9, 3;stroke-dashoffset:0;stroke-opacity:1" id="path36" inkscape:connector-curvature="0"/></g><text transform="matrix(1,0,0,-1,320.57,284.74)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#000000;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text38"><tspan x="0 9.5939999 15.606 24.57 33.119999 38.880001 47.844002 52.416 61.866001 71.351997 80.802002 86.832001 93.870003 109.26 118.746 128.196 137.214" y="0" sodipodi:role="line" id="tspan40">CreateInputsModel</tspan></text>
<path d="m 481.47,473.17 c 0,3.66 2.96,6.62 6.61,6.62 l 151.2,0 c 3.65,0 6.61,-2.96 6.61,-6.62 l 0,-26.46 c 0,-3.65 -2.96,-6.61 -6.61,-6.61 l -151.2,0 c -3.65,0 -6.61,2.96 -6.61,6.61 l 0,26.46 z" style="fill:#cc00cc;fill-opacity:0.2;fill-rule:evenodd;stroke:none" id="path42" inkscape:connector-curvature="0"/><g id="g44"><path d="m 481.47,473.17 c 0,3.66 2.96,6.62 6.61,6.62 l 151.2,0 c 3.65,0 6.61,-2.96 6.61,-6.62 l 0,-26.46 c 0,-3.65 -2.96,-6.61 -6.61,-6.61 l -151.2,0 c -3.65,0 -6.61,2.96 -6.61,6.61 l 0,26.46 z" style="fill:none;stroke:#cc00cc;stroke-width:3;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:9, 3;stroke-dashoffset:0;stroke-opacity:1" id="path46" inkscape:connector-curvature="0"/></g><text transform="matrix(1,0,0,-1,493.2,453.94)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#000000;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text48"><tspan x="0 9.5939999 15.606 24.57 33.119999 38.880001 47.844002 57.438 66.096001 70.181999 74.267998 83.718002 95.634003 104.994 111.024 115.074 124.56 134.00999" y="0" sodipodi:role="line" id="tspan50">CreateCalibOptions</tspan></text>
<path d="m 481.5,400.38 c 0,3.66 2.96,6.62 6.61,6.62 l 151.2,0 c 3.66,0 6.62,-2.96 6.62,-6.62 l 0,-26.45 c 0,-3.66 -2.96,-6.62 -6.62,-6.62 l -151.2,0 c -3.65,0 -6.61,2.96 -6.61,6.62 l 0,26.45 z" style="fill:#cc00cc;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path52" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,494.81,381.12)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text54"><tspan x="0 9.5939999 18.216 22.356001 26.406 35.855999 41.759998 50.273998 56.304001 60.354 69.839996 79.290001 88.290001 103.68 107.73 115.29 124.74 133.758" y="0" sodipodi:role="line" id="tspan56">Calibration_Michel</tspan></text>
<path d="m 308.97,131.92 c 0,3.66 2.96,6.62 6.61,6.62 l 151.2,0 c 3.66,0 6.62,-2.96 6.62,-6.62 l 0,-26.45 c 0,-3.66 -2.96,-6.618 -6.62,-6.618 l -151.2,0 c -3.65,0 -6.61,2.958 -6.61,6.618 l 0,26.45 z" style="fill:#92d050;fill-opacity:0.30196001;fill-rule:evenodd;stroke:none" id="path58" inkscape:connector-curvature="0"/><g id="g60"><path d="m 308.97,131.92 c 0,3.66 2.96,6.62 6.61,6.62 l 151.2,0 c 3.66,0 6.62,-2.96 6.62,-6.62 l 0,-26.45 c 0,-3.66 -2.96,-6.618 -6.62,-6.618 l -151.2,0 c -3.65,0 -6.61,2.958 -6.61,6.618 l 0,26.45 z" style="fill:none;stroke:#92d050;stroke-width:3;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:9, 3;stroke-dashoffset:0;stroke-opacity:1" id="path62" inkscape:connector-curvature="0"/></g><text transform="matrix(1,0,0,-1,331.25,112.61)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#000000;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text64"><tspan x="0 9.5939999 15.606 24.57 33.119999 38.880001 47.844002 52.416 61.866001 71.351997 80.802002 86.832001 93.870003 103.464 109.746 113.832" y="0" sodipodi:role="line" id="tspan66">CreateInputsCrit</tspan></text>
<path d="m 31.145,132.03 c 0,3.66 2.961,6.62 6.615,6.62 l 49.139,0 c 3.653,0 6.615,-2.96 6.615,-6.62 l 0,-26.46 c 0,-3.65 -2.962,-6.61 -6.615,-6.61 l -49.139,0 c -3.654,0 -6.615,2.96 -6.615,6.61 l 0,26.46 z" style="fill:#92d050;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path68" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,47.784,112.73)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text70"><tspan x="0 9.4499998 13.59 23.076" y="0" sodipodi:role="line" id="tspan72">plot</tspan></text>
<path d="m 481.5,385.53 -23.9,0 c -0.89,0 -1.62,0.73 -1.62,1.62 l 0,0.03 1.62,-1.63 -15.77,0 0,3.25 15.77,0 c 0.9,0 1.63,-0.72 1.63,-1.62 l 0,-0.03 -1.63,1.63 23.9,0 0,-3.25 z m -38.04,-3.23 -9.75,4.88 9.75,4.87 0,-9.75 z" style="fill:#cc00cc;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path74" inkscape:connector-curvature="0"/><path d="m 308.97,120.32 -14.17,0 c -0.9,0 -1.63,-0.73 -1.63,-1.63 l 0,-0.09 1.63,1.62 -6.05,0 0,-3.25 6.05,0 c 0.89,0 1.62,0.73 1.62,1.63 l 0,0.09 -1.62,-1.62 14.17,0 0,3.25 z m -18.6,3.15 -9.75,-4.87 9.75,-4.88 0,9.75 z" style="fill:#92d050;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path76" inkscape:connector-curvature="0"/><path d="m 120.22,388.8 -57.891,0 c -0.897,0 -1.625,-0.72 -1.625,-1.62 l 0,-240.41 3.25,0 0,240.41 -1.625,-1.63 57.891,0 0,3.25 z m -62.766,-240.4 4.875,-9.75 4.875,9.75 -9.75,0 z" style="fill:#7030a0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path78" inkscape:connector-curvature="0"/><path d="m 308.97,289.16 -77.96,0 c -0.9,0 -1.62,0.73 -1.62,1.63 l 0,68.39 3.25,0 0,-68.39 -1.63,1.62 77.96,0 0,-3.25 z m -82.83,68.4 4.9,9.74 4.85,-9.76 -9.75,0.02 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path80" inkscape:connector-curvature="0"/><path d="m 308.97,220.52 -78.19,0 c -0.89,0 -1.62,0.72 -1.62,1.62 l 0,137.07 3.25,0 0,-137.07 -1.63,1.63 78.19,0 0,-3.25 z m -83.06,137.06 4.87,9.75 4.88,-9.75 -9.75,0 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path82" inkscape:connector-curvature="0"/><path d="m 473.4,220.52 90.31,0 c 0.9,0 1.63,0.72 1.63,1.62 l 0,137.05 -3.25,0 0,-137.05 1.62,1.63 -90.31,0 0,-3.25 z m 95.19,137.04 -4.88,9.75 -4.87,-9.75 9.75,0 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path84" inkscape:connector-curvature="0"/><path d="m 365.67,51.027 -303.341,0 c -0.897,0 -1.625,0.728 -1.625,1.625 l 0,38.183 3.25,0 0,-38.183 -1.625,1.625 303.341,0 0,-3.25 z M 57.454,89.21 l 4.875,9.75 4.875,-9.75 -9.75,0 z" style="fill:#e46c0a;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path86" inkscape:connector-curvature="0"/><path d="m 473.4,289.16 90.31,0 c 0.9,0 1.63,0.73 1.63,1.63 l 0,68.4 -3.25,0 0,-68.4 1.62,1.62 -90.31,0 0,-3.25 z m 95.19,68.4 -4.88,9.75 -4.87,-9.75 9.75,0 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path88" inkscape:connector-curvature="0"/><path d="m 342.99,387.18 9.92,19.84 70.88,0 9.92,-19.84 -9.92,-19.85 -70.88,0 -9.92,19.85 z" style="fill:#00508c;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path90" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,365.11,381.14)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text92"><tspan x="0 8.8739996 17.496 23.4 32.021999" y="0" sodipodi:role="line" id="tspan94">Param</tspan></text>
<path d="m 564.12,173.36 -172.94,0 c -0.89,0 -1.62,-0.73 -1.62,-1.63 l 0,-25.07 3.25,0 0,25.07 -1.63,-1.62 172.94,0 0,3.25 z m -177.81,-25.07 4.87,-9.75 4.88,9.75 -9.75,0 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path96" inkscape:connector-curvature="0"/><g id="g98"><g id="g100" clip-path="url(#clipPath102)"><path d="m 562.06,440.1 0,-16.55 c 0,-0.9 0.72,-1.63 1.62,-1.63 l 0.03,0 -1.62,1.63 0,-8.43 3.25,0 0,8.43 c 0,0.9 -0.73,1.62 -1.63,1.62 l -0.03,0 1.63,-1.62 0,16.55 -3.25,0 z m -3.22,-23.35 4.87,-9.75 4.88,9.75 -9.75,0 z" style="fill:#cc00cc;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path106" inkscape:connector-curvature="0"/></g></g><path d="m 473.4,117.07 124.74,0 c 0.89,0 1.62,0.73 1.62,1.62 l 0,240.49 -3.25,0 0,-240.49 1.63,1.63 -124.74,0 0,-3.25 z m 129.61,240.48 -4.87,9.75 -4.88,-9.75 9.75,0 z" style="fill:#92d050;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path108" inkscape:connector-curvature="0"/><path d="m 354.33,388.8 -34.84,0 c -0.9,0 -1.63,-0.72 -1.63,-1.62 l 0,0 1.63,1.62 -26.72,0 0,-3.25 26.72,0 c 0.9,0 1.62,0.73 1.62,1.63 l 0,0 -1.62,-1.63 34.84,0 0,3.25 z m -59.93,3.25 -9.75,-4.87 9.75,-4.88 0,9.75 z" style="fill:#00508c;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path110" inkscape:connector-curvature="0"/><g id="g112"><g id="g114" clip-path="url(#clipPath116)"><path d="m 200.03,367.33 0,-114.44 c 0,-0.9 -0.73,-1.63 -1.62,-1.63 l 0,0 1.62,1.63 0,-106.32 -3.25,0 0,106.32 c 0,0.89 0.73,1.62 1.63,1.62 l 0,0 -1.63,-1.62 0,114.44 3.25,0 z m 3.25,-219.14 -4.87,-9.75 -4.88,9.75 9.75,0 z" style="fill:#7030a0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path120" inkscape:connector-curvature="0"/></g></g><path d="m 120.22,400.41 c 0,3.65 2.96,6.61 6.62,6.61 l 151.19,0 c 3.66,0 6.62,-2.96 6.62,-6.61 l 0,-26.46 c 0,-3.66 -2.96,-6.62 -6.62,-6.62 l -151.19,0 c -3.66,0 -6.62,2.96 -6.62,6.62 l 0,26.46 z" style="fill:#7030a0;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path122" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,164.42,381.14)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text124"><tspan x="0 9.7740002 19.224001 28.674 44.063999 53.549999 63 71.963997" y="0" sodipodi:role="line" id="tspan126">RunModel</tspan></text>
<g id="g128"><g id="g130" clip-path="url(#clipPath132)"><path d="m 389.56,270.94 0,-14.47 c 0,-0.9 0.73,-1.63 1.62,-1.63 l 0,0 -1.62,1.63 0,-6.36 3.25,0 0,6.36 c 0,0.89 -0.73,1.62 -1.63,1.62 l 0,0 1.63,-1.62 0,14.47 -3.25,0 z m -3.25,-19.2 4.87,-9.75 4.88,9.75 -9.75,0 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path136" inkscape:connector-curvature="0"/></g></g><text transform="matrix(1,0,0,-1,509.9,509.71)" style="font-variant:normal;font-weight:bold;font-size:24px;font-family:Calibri;-inkscape-font-specification:'Calibri Bold';writing-mode:lr-tb;fill:#5f5f5f;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text138"><tspan x="0 12.696 24.552 30.455999 36.360001 49.248001 57.240002 68.879997 77.208 83.112 96.024002" y="0" sodipodi:role="line" id="tspan140">Calibration</tspan></text>
<text transform="matrix(1,0,0,-1,142.75,509.38)" style="font-variant:normal;font-weight:bold;font-size:24px;font-family:Calibri;-inkscape-font-specification:'Calibri Bold';writing-mode:lr-tb;fill:#5f5f5f;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text142"><tspan x="0 11.352 17.256001 36.84 49.728001 55.632 67.199997 75.528 81.431999 94.344002" y="0" sodipodi:role="line" id="tspan144">Simulation</tspan></text>
<path d="m 308.97,235.37 c 0,3.65 2.96,6.62 6.62,6.62 l 151.19,0 c 3.66,0 6.62,-2.97 6.62,-6.62 l 0,-26.46 c 0,-3.65 -2.96,-6.61 -6.62,-6.61 l -151.19,0 c -3.66,0 -6.62,2.96 -6.62,6.61 l 0,26.46 z" style="fill:#00b0f0;fill-opacity:0.2;fill-rule:evenodd;stroke:none" id="path146" inkscape:connector-curvature="0"/><g id="g148"><path d="m 308.97,235.37 c 0,3.65 2.96,6.62 6.62,6.62 l 151.19,0 c 3.66,0 6.62,-2.97 6.62,-6.62 l 0,-26.46 c 0,-3.65 -2.96,-6.61 -6.62,-6.61 l -151.19,0 c -3.66,0 -6.62,2.96 -6.62,6.61 l 0,26.46 z" style="fill:none;stroke:#00b0f0;stroke-width:3;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:9, 3;stroke-dashoffset:0;stroke-opacity:1" id="path150" inkscape:connector-curvature="0"/></g><text transform="matrix(1,0,0,-1,324.29,216.07)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#000000;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text152"><tspan x="0 9.5939999 15.606 24.57 33.119999 38.880001 47.844002 57.618 67.068001 76.517998 88.433998 97.793999 103.824 107.874 117.36 126.81" y="0" sodipodi:role="line" id="tspan154">CreateRunOptions</tspan></text>
<path d="m 116.19,400.41 c 0,3.65 2.96,6.61 6.62,6.61 l 151.2,0 c 3.65,0 6.61,-2.96 6.61,-6.61 l 0,-26.46 c 0,-3.66 -2.96,-6.62 -6.61,-6.62 l -151.2,0 c -3.66,0 -6.62,2.96 -6.62,6.62 l 0,26.46 z" style="fill:#7030a0;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path156" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,160.39,381.14)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text158"><tspan x="0 9.7740002 19.224001 28.674 44.063999 53.549999 63 71.963997" y="0" sodipodi:role="line" id="tspan160">RunModel</tspan></text>
<path d="m 116.19,131.83 c 0,3.65 2.96,6.61 6.62,6.61 l 151.2,0 c 3.65,0 6.61,-2.96 6.61,-6.61 l 0,-26.46 c 0,-3.66 -2.96,-6.618 -6.61,-6.618 l -151.2,0 c -3.66,0 -6.62,2.958 -6.62,6.618 l 0,26.46 z" style="fill:#92d050;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path162" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,167.11,112.51)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text164"><tspan x="0 8.7840004 14.994 20.988001 30.474001 36.756001 46.349998 52.542 56.627998" y="0" sodipodi:role="line" id="tspan166">ErrorCrit</tspan></text>
<path d="m 351.85,52.652 9.92,19.845 59.53,0 9.93,-19.845 -9.93,-19.844 -59.53,0 -9.92,19.844 z" style="fill:#e46c0a;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path168" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,372.5,46.56)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text170"><tspan x="0 12.114 21.6 30.959999" y="0" sodipodi:role="line" id="tspan172">Qobs</tspan></text>
<g id="g174"><g id="g176" clip-path="url(#clipPath178)"><path d="m 392.81,71.77 0,13.541 c 0,0.897 -0.73,1.625 -1.63,1.625 l 0,0 1.63,-1.625 0,5.416 -3.25,0 0,-5.416 c 0,-0.897 0.73,-1.625 1.62,-1.625 l 0,0 -1.62,1.625 0,-13.541 3.25,0 z m 3.25,17.332 -4.88,9.75 -4.87,-9.75 9.75,0 z" style="fill:#e46c0a;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path182" inkscape:connector-curvature="0"/></g></g><path d="m 563.71,222.75 0.62,-50.29" style="fill:none;stroke:#00b0f0;stroke-width:3.25;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:none;stroke-opacity:1" id="path184" inkscape:connector-curvature="0"/></g></svg>
\ No newline at end of file