Skip to content
GitLab

Source

Target

Hide whitespace changes
Inline Side-by-side
Showing
with 1083 additions and 431 deletions
+1083 -431
inst/vignettesData/vignetteParamOptim.rda 0 → 100644
View file @ 4bccfee2
File added
inst/vignettesData/vignetteParamOptimCaramel.rda 0 → 100644
View file @ 4bccfee2
File added
man/BasinInfo.Rd
View file @ 4bccfee2
......@@ -6,7 +6,15 @@
\alias{BasinInfo}
\title{Data sample: characteristics of a fictional catchment (L0123001, L0123002 or L0123003)}
\title{Data sample: characteristics of a different catchments}
\description{
L0123001, L0123002 and L0123003 are fictional catchments. \cr
X0310010 contains actual data from the Durance River at Embrun [La Clapière] (Hautes-Alpes, France). \cr
\cr
R-object containing the code, station's name, area and hypsometric curve of the catchment.
}
\format{List named 'BasinInfo' containing
......@@ -17,19 +25,14 @@
}}
\description{
R-object containing the code, station's name, area and hypsometric curve of the catchment.
}
\seealso{
\code{\link{BasinObs}}.
}
\examples{
library(airGR)
data(L0123001)
str(BasinInfo)
library(airGR)
data(L0123001)
str(BasinInfo)
}
man/BasinObs.Rd
View file @ 4bccfee2
......@@ -7,9 +7,27 @@
\alias{L0123001}
\alias{L0123002}
\alias{L0123003}
\alias{X0310010}
\title{Data sample: time series of observations of a fictional catchment (L0123001, L0123002 or L0123003)}
\title{Data sample: time series of observations of different catchments}
\description{
L0123001, L0123002 or L0123003 are fictional catchments.
\cr
X0310010 contains actual data from the Durance River at Embrun [La Clapière] (Hautes-Alpes, France).
The flows are provided by Electricity of France (EDF) and were retrieved from the Banque Hydro database (http://www.hydro.eaufrance.fr).
The meteorological forcing are derived from the SAFRAN reanalysis from Météo-France (Vidal et al., 2010).
\cr
R-object containing the times series of precipitation, temperature, potential evapotranspiration and discharge.
X0310010 contains in addition MODIS snow cover area (SCA) data retrieved from the National Snow and Ice Data Center (NSIDC) repository (https://nsidc.org/). Five SCA time series are given, corresponding to 5 elevation bands of the CemaNeige model (default configuration). SCA data for days with important cloudiness (> 40 \%) were set to missing values for the sake of data representativeness. .
\cr \cr
Times series for L0123001, L0123002 and X0310010 are at the daily time step for use with daily models such as GR4J, GR5J, GR6J, CemaNeigeGR4J, CemaNeigeGR5J and CemaNeigeGR6J. \cr
Times series for X0310010 are provided in order to test hysteresis version of CemaNeige (see \code{\link{CreateRunOptions}} (Riboust et al., 2019). \cr
Times series for L0123003 are at the hourly time step for use with hourly models such as GR4H or GR5H.
}
\format{Data frame named 'BasinObs' containing
......@@ -19,10 +37,14 @@
}}
\description{
R-object containing the times series of precipitation, temperature, potential evapotranspiration and discharges. \cr
Times series for L0123001 or L0123002 are at the daily time step for use with daily models such as GR4J, GR5J, GR6J, CemaNeigeGR4J, CemaNeigeGR5J and CemaNeigeGR6J.
Times series for L0123003 are at the hourly time step for use with hourly models such as GR4H.
\references{
Riboust, P., Thirel, G., Le Moine, N. and Ribstein P. (2019).
Revisiting a simple degree-day model for integrating satellite data: Implementation of SWE-SCA hystereses.
Journal of Hydrology and Hydromechanics, 67(1), 70–81, \doi{10.2478/johh-2018-0004}.
\cr\cr
Vidal, J.-P., Martin, E., Franchistéguy, L., Baillon, M. and Soubeyroux, J. (2010).
A 50-year high-resolution atmospheric reanalysis over France with the Safran system.
International Journal of Climatology, 30, 1627–1644, \doi{10.1002/joc.2003}.
}
......@@ -32,8 +54,8 @@ Times series for L0123003 are at the hourly time step for use with hourly models
\examples{
library(airGR)
data(L0123001)
str(BasinObs)
library(airGR)
data(L0123001)
str(BasinObs)
}
man/Calibration.Rd
View file @ 4bccfee2
......@@ -5,12 +5,18 @@
\alias{Calibration}
\title{Calibration algorithm which optimises the error criterion selected as objective function using the provided functions}
\title{Calibration algorithm that optimises the error criterion selected as objective function using the provided functions}
\description{
Calibration algorithm that optimises the error criterion selected as objective function using the provided functions.
}
\usage{
Calibration(InputsModel, RunOptions, InputsCrit, CalibOptions, FUN_MOD,
FUN_CRIT, FUN_CALIB = Calibration_Michel, FUN_TRANSFO = NULL,
verbose = TRUE)
Calibration(InputsModel, RunOptions, InputsCrit, CalibOptions,
FUN_MOD, FUN_CRIT, FUN_CALIB = Calibration_Michel,
FUN_TRANSFO = NULL, verbose = TRUE, ...)
}
......@@ -25,13 +31,15 @@ Calibration(InputsModel, RunOptions, InputsCrit, CalibOptions, FUN_MOD,
\item{FUN_MOD}{[function] hydrological model function (e.g. \code{\link{RunModel_GR4J}}, \code{\link{RunModel_CemaNeigeGR4J}})}
\item{FUN_CRIT}{[function] error criterion function (e.g. \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}})}
\item{FUN_CRIT}{(deprecated) [function] error criterion function (e.g. \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}})}
\item{FUN_CALIB}{(optional) [function] calibration algorithm function (e.g. \code{\link{Calibration_Michel}}), \code{default = Calibration_Michel}}
\item{FUN_CALIB}{(deprecated) [function] calibration algorithm function (e.g. \code{\link{Calibration_Michel}}), \code{default = Calibration_Michel}}
\item{FUN_TRANSFO}{(optional) [function] model parameters transformation function, if the \code{FUN_MOD} used is native in the package \code{FUN_TRANSFO} is automatically defined}
\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, default = \code{TRUE}}
\item{...}{Further arguments to pass to \code{\link{RunModel}}}
}
......@@ -40,11 +48,6 @@ Calibration(InputsModel, RunOptions, InputsCrit, CalibOptions, FUN_MOD,
}
\description{
Calibration algorithm which optimises the error criterion selected as objective function using the provided functions.
}
\examples{
library(airGR)
......@@ -52,59 +55,59 @@ library(airGR)
data(L0123001)
## preparation of 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)
## calibration period selection
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="01/01/1990 00:00"),
which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="31/12/1999 00:00"))
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"))
## preparation of RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
## calibration criterion: preparation of the InputsCrit object
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
## preparation of CalibOptions object
CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel)
## calibration
OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions,
InputsCrit = InputsCrit, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J, FUN_CRIT = ErrorCrit_NSE,
InputsCrit = InputsCrit, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J,
FUN_CALIB = Calibration_Michel)
## simulation
Param <- OutputsCalib$ParamFinalR
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions,
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions,
Param = Param, FUN = RunModel_GR4J)
## results preview
plot(OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Kling-Gupta Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
}
\author{
Laurent Coron
Laurent Coron, Olivier Delaigue
}
\seealso{
\code{\link{Calibration_Michel}},
\code{\link{ErrorCrit}}, \code{\link{TransfoParam}},
\code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}},
\code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}.
\code{\link{ErrorCrit}}, \code{\link{TransfoParam}},
\code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}},
\code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}.
}
man/Calibration_Michel.Rd
View file @ 4bccfee2
......@@ -5,12 +5,21 @@
\alias{Calibration_Michel}
\title{Calibration algorithm optimises the error criterion selected as objective function using the Irstea-HBAN procedure described by C. Michel}
\title{Calibration algorithm that optimises the error criterion selected as objective function using the Irstea procedure described by C. Michel}
\description{
Calibration algorithm that optimises the error criterion selected as objective function. \cr
\cr
The algorithm combines a global and a local approach.
First, a screening is performed using either a rough predefined grid or a list of parameter sets.
Then a steepest descent local search algorithm is performed, starting from the result of the screening procedure.
}
\usage{
Calibration_Michel(InputsModel, RunOptions, InputsCrit, CalibOptions,
FUN_MOD, FUN_CRIT, FUN_TRANSFO = NULL, verbose = TRUE)
FUN_MOD, FUN_CRIT, FUN_TRANSFO = NULL, verbose = TRUE, ...)
}
......@@ -25,52 +34,45 @@ Calibration_Michel(InputsModel, RunOptions, InputsCrit, CalibOptions,
\item{FUN_MOD}{[function] hydrological model function (e.g. \code{\link{RunModel_GR4J}}, \code{\link{RunModel_CemaNeigeGR4J}})}
\item{FUN_CRIT}{[function] error criterion function (e.g. \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}})}
\item{FUN_CRIT}{(deprecated) [function] error criterion function (e.g. \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}})}
\item{FUN_TRANSFO}{(optional) [function] model parameters transformation function, if the \code{FUN_MOD} used is native in the package \code{FUN_TRANSFO} is automatically defined}
\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, default = \code{TRUE}}
\item{...}{(optional) arguments to pass to \code{\link{RunModel}}}
}
\value{
[list] list containing the function outputs organised as follows:
\tabular{ll}{
\emph{$ParamFinalR } \tab [numeric] parameter set obtained at the end of the calibration \cr
\emph{$CritFinal } \tab [numeric] error criterion selected as objective function obtained at the end of the calibration \cr
\emph{$NIter } \tab [numeric] number of iterations during the calibration \cr
\emph{$NRuns } \tab [numeric] number of model runs done during the calibration \cr
\emph{$HistParamR } \tab [numeric] table showing the progression steps in the search for optimal set: parameter values \cr
\emph{$HistCrit } \tab [numeric] table showing the progression steps in the search for optimal set: criterion values \cr
\emph{$MatBoolCrit } \tab [boolean] table giving the requested and actual time steps over which the model is calibrated \cr
\emph{$CritName } \tab [character] name of the calibration criterion used as objective function \cr
\emph{$CritBestValue} \tab [numeric] theoretical best criterion value \cr
}
}
\description{
Calibration algorithm optimises the error criterion selected as objective function. \cr
\cr
The algorithm is based on a local search procedure.
First, a screening is performed using either a rough predefined grid or a list of parameter sets
and then a simple steepest descent local search algorithm is performed.
\tabular{ll}{
\emph{$ParamFinalR } \tab [numeric] parameter set obtained at the end of the calibration \cr
\emph{$CritFinal } \tab [numeric] error criterion selected as objective function obtained at the end of the calibration \cr
\emph{$NIter } \tab [numeric] number of iterations during the calibration \cr
\emph{$NRuns } \tab [numeric] number of model runs done during the calibration \cr
\emph{$HistParamR } \tab [numeric] table showing the progression steps in the search for optimal set: parameter values \cr
\emph{$HistCrit } \tab [numeric] table showing the progression steps in the search for optimal set: criterion values \cr
\emph{$MatBoolCrit } \tab [boolean] table giving the requested and actual time steps over which the model is calibrated \cr
\emph{$CritName } \tab [character] name of the calibration criterion used as objective function \cr
\emph{$CritBestValue} \tab [numeric] theoretical best criterion value \cr
}
}
\details{
A screening is first performed either based on a rough predefined grid (considering various initial
A screening is first performed either based on a rough predefined grid (considering various initial
values for each parameter) or from a list of initial parameter sets. \cr
The best set identified in this screening is then used as a starting point for the steepest
The best set identified in this screening is then used as a starting point for the steepest
descent local search algorithm. \cr
For this search, since the ranges of parameter values can be quite different,
simple mathematical transformations are applied to parameters to make them vary
in a similar range and get a similar sensitivity to a predefined search step. This is done using the TransfoParam functions. \cr
During the steepest descent method, at each iteration, we start from a parameter set of NParam values (NParam being the number of
free parameters of the chosen hydrological model) and we determine the 2*NParam-1 new candidates
During the steepest descent method, at each iteration, we start from a parameter set of NParam values (NParam being the number of
free parameters of the chosen hydrological model) and we determine the 2*NParam-1 new candidates
by changing one by one the different parameters (+/- search step). \cr
All these candidates are tested and the best one kept to be the starting point for the next
iteration. At the end of each iteration, the the search step is either increased or decreased to adapt
All these candidates are tested and the best one kept to be the starting point for the next
iteration. At the end of each iteration, the the search step is either increased or decreased to adapt
the progression speed. A composite step can occasionally be done. \cr
The calibration algorithm stops when the search step becomes smaller than a predefined threshold.
}
......@@ -83,28 +85,28 @@ library(airGR)
data(L0123001)
## preparation of 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)
## calibration period selection
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="01/01/1990 00:00"),
which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="31/12/1999 00:00"))
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"))
## preparation of RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
IndPeriod_Run = Ind_Run)
## calibration criterion: preparation of the InputsCrit object
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
## preparation of CalibOptions object
CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel)
## calibration
OutputsCalib <- Calibration_Michel(InputsModel = InputsModel, RunOptions = RunOptions,
InputsCrit = InputsCrit, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J, FUN_CRIT = ErrorCrit_NSE)
OutputsCalib <- Calibration_Michel(InputsModel = InputsModel, RunOptions = RunOptions,
InputsCrit = InputsCrit, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J)
## simulation
Param <- OutputsCalib$ParamFinalR
......@@ -115,32 +117,33 @@ OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
plot(OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Kling-Gupta Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
}
\author{
Laurent Coron, Claude Michel, Olivier Delaigue, Guillaume Thirel
Laurent Coron, Claude Michel, Charles Perrin, Thibault Mathevet, Olivier Delaigue, Guillaume Thirel, David Dorchies
}
\references{
Michel, C. (1991),
Hydrologie appliquée aux petits bassins ruraux, Hydrology handbook (in French), Cemagref, Antony, France.
Hydrologie appliquée aux petits bassins ruraux.
Hydrology handbook (in French), Cemagref, Antony, France.
}
\seealso{
\code{\link{Calibration}},
\code{\link{RunModel_GR4J}}, \code{\link{TransfoParam_GR4J}}, \code{\link{ErrorCrit_RMSE}},
\code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}},
\code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}.
\code{\link{RunModel_GR4J}}, \code{\link{TransfoParam}}, \code{\link{ErrorCrit_RMSE}},
\code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}},
\code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}.
}
man/CreateCalibOptions.Rd
View file @ 4bccfee2
......@@ -5,14 +5,20 @@
\alias{CreateCalibOptions}
\title{Creation of the CalibOptions object required but the Calibration functions}
\title{Creation of the CalibOptions object required but the Calibration* functions}
\description{
Creation of the \emph{CalibOptions} object required by the \code{Calibration*} functions.
}
\usage{
CreateCalibOptions(FUN_MOD, FUN_CALIB = Calibration_Michel,
FUN_TRANSFO = NULL, FixedParam = NULL,
SearchRanges = NULL, StartParamList = NULL,
StartParamDistrib = NULL)
FUN_TRANSFO = NULL, IsHyst = FALSE, IsSD = FALSE,
FixedParam = NULL,
SearchRanges = NULL, StartParamList = NULL,
StartParamDistrib = NULL)
}
......@@ -23,6 +29,10 @@ CreateCalibOptions(FUN_MOD, FUN_CALIB = Calibration_Michel,
\item{FUN_TRANSFO}{(optional) [function] model parameters transformation function, if the \code{FUN_MOD} used is native in the package, \code{FUN_TRANSFO} is automatically defined}
\item{IsHyst}{[boolean] boolean indicating if the hysteresis version of CemaNeige is used. See details}
\item{IsSD}{[boolean] boolean indicating if the semi-distributed lag model is used. See details}
\item{FixedParam}{(optional) [numeric] vector giving the values set for the non-optimised parameter values (NParam columns, 1 line)
\cr Example:
\tabular{llllll}{
......@@ -58,23 +68,27 @@ CreateCalibOptions(FUN_MOD, FUN_CALIB = Calibration_Michel,
}
\value{
[list] object of class \emph{CalibOptions} containing the data required to evaluate the model outputs; it can include the following:
\tabular{ll}{
\emph{$FixedParam } \tab [numeric] vector giving the values to allocate to non-optimised parameter values \cr
\emph{$SearchRanges } \tab [numeric] matrix giving the ranges of raw parameters \cr
\emph{$StartParamList } \tab [numeric] matrix of parameter sets used for grid-screening calibration procedure \cr
\emph{$StartParamDistrib} \tab [numeric] matrix of parameter values used for grid-screening calibration procedure \cr
}
\tabular{ll}{
\emph{$FixedParam } \tab [numeric] vector giving the values to allocate to non-optimised parameter values \cr
\emph{$SearchRanges } \tab [numeric] matrix giving the ranges of raw parameters \cr
\emph{$StartParamList } \tab [numeric] matrix of parameter sets used for grid-screening calibration procedure \cr
\emph{$StartParamDistrib} \tab [numeric] matrix of parameter values used for grid-screening calibration procedure \cr
}
}
\description{
Creation of the \emph{CalibOptions} object required by the \code{Calibration*} functions.
}
\details{
Users wanting to use \code{FUN_MOD}, \code{FUN_CALIB} or \code{FUN_TRANSFO} functions that are not included in
the package must create their own \code{CalibOptions} object accordingly. \cr
## --- CemaNeige version
\details{
Users wanting to use \code{FUN_MOD}, \code{FUN_CALIB} or \code{FUN_TRANSFO} functions that are not included in
the package must create their own \code{CalibOptions} object accordingly.
If \code{IsHyst = FALSE}, the original CemaNeige version from Valéry et al. (2014) is used. \cr
If \code{IsHyst = TRUE}, the CemaNeige version from Riboust et al. (2019) is used. Compared to the original version, this version of CemaNeige needs two more parameters and it includes a representation of the hysteretic relationship between the Snow Cover Area (SCA) and the Snow Water Equivalent (SWE) in the catchment. The hysteresis included in airGR is the Modified Linear hysteresis (LH*); it is represented on panel b) of Fig. 3 in Riboust et al. (2019). Riboust et al. (2019) advise to use the LH* version of CemaNeige with parameters calibrated using an objective function combining 75 \% of KGE calculated on discharge simulated from a rainfall-runoff model compared to observed discharge and 5 \% of KGE calculated on SCA on 5 CemaNeige elevation bands compared to satellite (e.g. MODIS) SCA (see Eq. (18), Table 3 and Fig. 6). Riboust et al. (2019)'s tests were realized with GR4J as the chosen rainfall-runoff model. \cr
## --- Semi-distributed mode
If \code{InputsModel} parameter has been created for using a semi-distributed (SD) model (See \code{\link{CreateInputsModel}}), the parameter \code{IsSD} should be set to \code{TRUE}.
}
......@@ -85,20 +99,20 @@ library(airGR)
data(L0123001)
## preparation of 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)
## calibration period selection
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="01/01/1990 00:00"),
which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="31/12/1999 00:00"))
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"))
## preparation of RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
## calibration criterion: preparation of the InputsCrit object
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
## preparation of CalibOptions object
CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel)
......@@ -106,36 +120,34 @@ CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibrat
## calibration
OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions,
InputsCrit = InputsCrit, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J, FUN_CRIT = ErrorCrit_NSE,
FUN_MOD = RunModel_GR4J,
FUN_CALIB = Calibration_Michel)
## simulation
Param <- OutputsCalib$ParamFinalR
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions,
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions,
Param = Param, FUN = RunModel_GR4J)
## results preview
plot(OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Kling-Gupta Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
}
\author{
Laurent Coron
Laurent Coron, Olivier Delaigue, Guillaume Thirel, David Dorchies
}
\seealso{
\code{\link{Calibration}}, \code{\link{RunModel}}
}
man/CreateErrorCrit_GAPX.Rd 0 → 100644
View file @ 4bccfee2
\encoding{UTF-8}
\name{CreateErrorCrit_GAPX}
\alias{CreateErrorCrit_GAPX}
\title{Creation of the ErrorCrit_GAPX function}
\description{
Function which creates the function \code{ErrorCrit_GAPX}.
The produced function \code{ErrorCrit_GAPX} allows computing an error criterion based on the GAPX formula proposed by Lavenne et al. (2019).
}
\usage{
CreateErrorCrit_GAPX(FUN_TRANSFO)
}
\arguments{
\item{FUN_TRANSFO}{[function] The parameter transformation function used with the model}
}
\value{
[function] function \code{ErrorCrit_GAPX} embedding the parameter transformation function used with the model
}
\details{
In addition to the criterion value, the function outputs include a multiplier (-1 or +1) that allows
the use of the function for model calibration: the product \eqn{CritValue \times Multiplier} is the criterion to be minimised
(\code{Multiplier = -1} for NSE).
}
\examples{
library(airGR)
## loading catchment data
data(L0123001)
## preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
## calibration period selection
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"))
## preparation of RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
IndPeriod_Run = Ind_Run)
## simulation
Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions, Param = Param)
## Creation of the ErrorCrit GAPX function
ErrorCrit_GAPX <- CreateErrorCrit_GAPX(TransfoParam_GR4J)
## The "a priori" parameters for GAPX
AprParamR <- c(X1 = 157, X2 = 0.8, X3 = 100, X4 = 1.5)
AprParamT <- TransfoParam_GR4J(AprParamR, "RT")
## Single efficiency criterion: GAPX with a priori parameters
InputsCrit <- CreateInputsCrit(ErrorCrit_GAPX,
InputsModel,
RunOptions,
Obs = AprParamT,
VarObs = "ParamT")
ErrorCrit <- ErrorCrit_GAPX(InputsCrit, OutputsModel)
str(ErrorCrit)
}
\author{
David Dorchies
}
\references{
de Lavenne, A., Andréassian, V., Thirel, G., Ramos, M.-H. and Perrin, C. (2019).
A Regularization Approach to Improve the Sequential Calibration of a Semidistributed Hydrological Model.
Water Resources Research 55, 8821–8839. \doi{10.1029/2018WR024266}
}
\seealso{
\code{\link{CreateInputsCrit}}, \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}},
\code{\link{ErrorCrit_KGE}}, \code{\link{ErrorCrit_KGE2}}
}
man/CreateIniStates.Rd
View file @ 4bccfee2
......@@ -8,12 +8,20 @@
\title{Creation of the IniStates object possibly required by the CreateRunOptions functions}
\description{
Creation of the \emph{IniStates} object possibly required by the \code{\link{CreateRunOptions}} function.
}
\usage{
CreateIniStates(FUN_MOD, InputsModel,
ProdStore = 350, RoutStore = 90, ExpStore = NULL,
UH1 = NULL, UH2 = NULL,
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = TRUE)
IsHyst = FALSE, IsIntStore = FALSE,
ProdStore = 350, RoutStore = 90,
ExpStore = NULL, IntStore = NULL,
UH1 = NULL, UH2 = NULL,
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
GthrCemaNeigeLayers = NULL, GlocmaxCemaNeigeLayers = NULL,
SD = NULL, verbose = TRUE)
}
......@@ -22,12 +30,18 @@ CreateIniStates(FUN_MOD, InputsModel,
\item{InputsModel}{[object of class \code{InputsModel}] see \code{\link{CreateInputsModel}} for details}
\item{ProdStore}{[numeric] production store level [mm]}
\item{IsHyst}{[boolean] boolean indicating if the hysteresis version of CemaNeige is used. See details}
\item{RoutStore}{[numeric] routing store level [mm]}
\item{IsIntStore}{[boolean] boolean indicating if the interception store is used in GR5H. See details}
\item{ProdStore}{[numeric] production store level [mm] for all GR models except GR1A}
\item{RoutStore}{[numeric] routing store level [mm] for all GR models except GR1A}
\item{ExpStore}{(optional) [numeric] series of exponential store level (negative) [mm] for the GR6J model}
\item{IntStore}{(optional) [numeric] series of rainfall neutralisation or interception store level [mm] for the GR5H model}
\item{UH1}{(optional) [numeric] unit hydrograph 1 levels [mm]}
\item{UH2}{(optional) [numeric] unit hydrograph 2 levels [mm]}
......@@ -36,6 +50,12 @@ CreateIniStates(FUN_MOD, InputsModel,
\item{eTGCemaNeigeLayers}{(optional) [numeric] snow pack thermal state [°C], possibly used to create the CemaNeige model initial state}
\item{GthrCemaNeigeLayers}{(optional) [numeric] melt threshold [mm], possibly used to create the CemaNeige model initial state in case the Linear Hysteresis version is used}
\item{GlocmaxCemaNeigeLayers}{(optional) [numeric] local melt threshold for hysteresis [mm], possibly used to create the CemaNeige model initial state in case the Linear Hysteresis version is used}
\item{SD}{(optional) [list] of [numeric] states of delayed upstream flows for semi-distributed models, the nature of the state and the unit depend on the model and the unit of the upstream flow}
\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, default = \code{TRUE}}
}
......@@ -43,25 +63,21 @@ CreateIniStates(FUN_MOD, InputsModel,
\value{
[list] object of class \code{IniStates} containing the initial model internal states; it always includes the following:
\tabular{ll}{
\emph{$Store } \tab [numeric] list of store levels (\emph{$Prod}, \emph{$Rout} and \emph{$Exp}) \cr
\emph{$UH } \tab [numeric] list of unit hydrographs levels (\emph{$UH1} and \emph{$UH2}) \cr
\emph{$CemaNeigeLayers} \tab [numeric] list of CemaNeige variables (\emph{$G} and \emph{$eTG})
}
\tabular{ll}{
\emph{$Store } \tab [numeric] list of store levels (\emph{$Prod}, \emph{$Rout} and \emph{$Exp}) \cr
\emph{$UH } \tab [numeric] list of unit hydrographs levels (\emph{$UH1} and \emph{$UH2}) \cr
\emph{$CemaNeigeLayers} \tab [numeric] list of CemaNeige variables (\emph{$G}, \emph{$eTG}, \emph{$GthrCemaNeigeLayers} and \emph{$GlocmaxCemaNeigeLayers})
}
}
\description{
Creation of the \emph{IniStates} object possibly required by the \code{\link{CreateRunOptions}} function.
}
\details{
20 numeric values are required for UH1 and 40 numeric values are required for UH2 if GR4J, GR5J or GR6J are used (respectivly 20*24 and 40*24 for the daily model GR4H). \cr
\code{GCemaNeigeLayers} and \code{eTGCemaNeigeLayers} require each numeric values as many as given in \code{\link{CreateInputsModel}} with the \code{NLayersargument}. \code{eTGCemaNeigeLayers} values can be negatives.\cr
The structure of the object of class \code{IniStates} returned is always exactly the same for all models (except for the unit hydrographs levels that contain more values with GR4H), even if some states do nt exist (e.g. \emph{$UH$UH1} for GR2M). \cr
If CemaNeige is not used, \emph{$CemaNeigeLayers$G} and \emph{$CemaNeigeLayers$eTG} are set to \code{NA}. \cr
20 numeric values are required for UH1 and 40 numeric values are required for UH2 if GR4J, GR5J or GR6J are used (respectively 20*24 and 40*24 for the hourly models GR4H and GR5H). Please note that depending on the X4 parameter value that will be provided when running the model, not all the values may be used (only the first int(X4)+1 values are used for UH1 and the first 2*int(X4)+1 for UH2). \cr
\code{GCemaNeigeLayers} and \code{eTGCemaNeigeLayers} require each numeric values as many as given in \code{\link{CreateInputsModel}} with the \code{NLayers} argument. \code{eTGCemaNeigeLayers} values can be negative.\cr
The structure of the object of class \code{IniStates} returned is always exactly the same for all models (except for the unit hydrographs levels that contain more values with GR4H and GR5H), even if some states do not exist (e.g. \emph{$UH$UH1} for GR2M). \cr
If CemaNeige is not used, \emph{$CemaNeigeLayers$G}, \emph{$CemaNeigeLayers$eTG} \emph{$CemaNeigeLayers$GthrCemaNeigeLayers} and \emph{$CemaNeigeLayers$GlocmaxCemaNeigeLayers} are set to \code{NA}. \cr
Nota: the \code{StateEnd} objects from the outputs of \code{RunModel*} functions already respect the format given by the \code{CreateIniStates} function.
Nota: the \code{StateEnd} objects from the outputs of \code{RunModel*} functions already respect the format given by the \code{CreateIniStates} function.
}
......@@ -72,12 +88,12 @@ library(airGR)
data(L0123001)
## 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)
## run period selection
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="01/01/1990 00:00"),
which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="31/12/1991 00:00"))
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"))
### preparation of the IniStates object with low values of ProdStore and RoutStore
......@@ -85,41 +101,45 @@ IniStates <- CreateIniStates(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
ProdStore = 0, RoutStore = 0, ExpStore = NULL,
UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL)
str(IniStates)
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
GthrCemaNeigeLayers = NULL, GlocmaxCemaNeigeLayers = NULL)
str(IniStates)
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
IndPeriod_WarmUp = 0L,
IndPeriod_Run = Ind_Run, IniStates = IniStates)
## simulation
Param <- c(257.238, 1.012, 88.235, 2.208)
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param)
Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions, Param = Param)
## results preview
plot(OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
plot(OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
### preparation of the IniStates object with high values of ProdStore and RoutStore
IniStates <- CreateIniStates(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
ProdStore = 450, RoutStore = 100, ExpStore = NULL,
UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL)
str(IniStates)
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
GthrCemaNeigeLayers = NULL, GlocmaxCemaNeigeLayers = NULL)
str(IniStates)
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
IndPeriod_WarmUp = 0L,
IndPeriod_Run = Ind_Run, IniStates = IniStates)
## simulation
Param <- c(257.238, 1.012, 88.235, 2.208)
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param)
Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions, Param = Param)
## results preview
plot(OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
plot(OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
}
......@@ -131,4 +151,3 @@ Olivier Delaigue
\seealso{
\code{\link{CreateRunOptions}}
}
man/CreateInputsCrit.Rd
View file @ 4bccfee2
......@@ -8,51 +8,96 @@
\title{Creation of the InputsCrit object required to the ErrorCrit functions}
\description{
Creation of the \code{InputsCrit} object required to the \code{ErrorCrit_*} functions. This function is used to define whether the user wants to calculate a single criterion, multiple criteria in the same time, or a composite criterion, which averages several criteria.
}
\usage{
CreateInputsCrit(FUN_CRIT, InputsModel, RunOptions, Qobs, BoolCrit = NULL,
transfo = "", Ind_zeroes = NULL, epsilon = NULL)
CreateInputsCrit(FUN_CRIT, InputsModel, RunOptions,
Obs, VarObs = "Q", BoolCrit = NULL,
transfo = "", Weights = NULL,
epsilon = NULL,
warnings = TRUE)
}
\arguments{
\item{FUN_CRIT}{[function] error criterion function (e.g. ErrorCrit_RMSE, ErrorCrit_NSE)}
\item{FUN_CRIT}{[function (atomic or list)] error criterion function (e.g. \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}})}
\item{InputsModel}{[object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details}
\item{RunOptions}{[object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details}
\item{Qobs}{[numeric] series of observed discharges [mm/time step]}
\item{Obs}{[numeric (atomic or list)] series of observed variable ([mm/time step] for discharge or SWE, [-] for SCA)}
\item{VarObs}{(optional) [character (atomic or list)] names of the observed variable (\code{"Q"} by default, or one of \code{"SCA"}, \code{"SWE"}])}
\item{BoolCrit}{(optional) [boolean] boolean giving the time steps to consider in the computation (all time steps are consider by default)}
\item{BoolCrit}{(optional) [boolean (atomic or list)] boolean (the same length as \code{Obs}) giving the time steps to consider in the computation (all time steps are considered by default. See details)}
\item{transfo}{(optional) [character] name of the transformation (e.g. \code{""}, \code{"sqrt"}, \code{"log"}, \code{"inv"}, \code{"sort"})}
\item{transfo}{(optional) [character (atomic or list)] name of the transformation applied to the variables (e.g. \code{""}, \code{"sqrt"}, \code{"log"}, \code{"inv"}, \code{"sort"}, \code{"boxcox"} or a numeric value for power transformation . See details)}
\item{Ind_zeroes}{(optional) [numeric] indices of the time steps where zeroes are observed}
\item{Weights}{(optional) [numeric (atomic or list)] vector of weights necessary to calculate a composite criterion (the same length as \code{FUN_CRIT}) giving the weights to use for elements of \code{FUN_CRIT} [-]. See details}
\item{epsilon}{(optional) [numeric] epsilon to add to all Qobs and Qsim if \emph{$Ind_zeroes} is not empty}
\item{epsilon}{(optional) [numeric (atomic or list)] small value to add to all observations and simulations when \code{"log"} or \code{"inv"} transformations are used [same unit as \code{Obs}]. See details}
\item{warnings}{(optional) [boolean] boolean indicating if the warning messages are shown, default = \code{TRUE}}
}
\value{
[list] object of class \emph{InputsCrit} containing the data required to evaluate the model outputs; it can include the following:
\tabular{ll}{
\emph{$BoolCrit } \tab [boolean] boolean giving the time steps considered in the computation \cr
\emph{$Qobs } \tab [numeric] series of observed discharges [mm/time step] \cr
\emph{$transfo } \tab [character] name of the transformation (e.g. "", "sqrt", "log", "inv", "sort") \cr
\emph{$Ind_zeroes} \tab [numeric] indices of the time steps where zeroes are observed \cr
\emph{$epsilon } \tab [numeric] epsilon to add to all Qobs and Qsim if \emph{$Ind_zeroes} is not empty \cr
}
\tabular{ll}{
\emph{$FUN_CRIT } \tab [function] error criterion function (e.g. \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}}) \cr
\emph{$Obs } \tab [numeric] series of observed variable(s) ([mm/time step] for discharge or SWE, [-] for SCA) \cr
\emph{$VarObs } \tab [character] names of the observed variable(s) \cr
\emph{$BoolCrit } \tab [boolean] boolean giving the time steps considered in the computation \cr
\emph{$transfo } \tab [character] name of the transformation (e.g. \code{""}, \code{"sqrt"}, \code{"log"}, \code{"inv"}, \code{"sort"}, \code{"boxcox"} or a number for power transformation) \cr
\emph{$epsilon } \tab [numeric] small value to add to all observations and simulations when \code{"log"} or \code{"inv"} transformations are used [same unit as \code{Obs}] \cr
\emph{$Weights } \tab [numeric] vector (same length as \code{VarObs}) giving the weights to use for elements of \code{FUN_CRIT} [-] \cr
}
When \code{Weights = NULL}, \code{CreateInputsCrit} returns an object of class \emph{Single} that is a list such as the one described above. \cr
When \code{Weights} contains at least one \code{NULL} value and \code{Obs} contains a list of observations, \code{CreateInputsCrit} returns an object of class \emph{Multi} that is a list of lists such as the one described above. The \code{\link{ErrorCrit}} function will then compute the different criteria prepared by \code{CreateInputsCrit}. \cr
When \code{Weights} is a list of at least 2 numerical values, \code{CreateInputsCrit} returns an object of class \emph{Compo} that is a list of lists such as the one described above. This object will be useful to compute composite criterion with the \code{\link{ErrorCrit}} function. \cr
To calculate composite or multiple criteria, it is necessary to use the \code{ErrorCrit} function. The other \code{ErrorCrit_*} functions (e.g. \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}}) can only use objects of class \emph{Single} (and not \emph{Multi} or \emph{Compo}). \cr
}
\description{
Creation of the \emph{InputsCrit} object required to the \code{ErrorCrit*} functions.
\details{
Users wanting to use \code{FUN_CRIT} functions that are not included in the package must create their own InputsCrit object accordingly. \cr \cr
## --- Period of calculation
Criteria can be calculated over discontinuous periods (i.e. only over winter periods, or when observed discharge is below a certain threshold). To do so, please indicate in \code{Bool_Crit} which indices must be used for calcullation. Discontinuous periods are allowed in the \code{Bool_Crit} argument.
## --- Transformations
Transformations are simple functions applied to the observed and simulated variables used in order to change their distribution. Transformations are often used in hydrology for modifying the weight put on errors made for high flows or low flows. The following transformations are available: \cr \cr
\itemize{
\item \code{""}: no transformation is used (default case)
\item \code{"sqrt"}: squared root transformation
\item \code{"log"}: logarithmic transformation (see below regarding the specific case of KGE or KGE2)
\item \code{"inv"}: inverse transformation
\item \code{"sort"}: sort transformation (the simulated and observed variables are sorted from lowest to highest)
\item \code{"boxcox"}: Box-Cox transformation (see below for details)
\item numeric: power transformation (see below for details)
}
We do not advise computing KGE or KGE' with log-transformation as it might be wrongly influenced by discharge values close to 0 or 1 and the criterion value is dependent on the discharge unit. See Santos et al. (2018) for more details and alternative solutions (see the references list below). \cr \cr
In order to make sure that KGE and KGE2 remain dimensionless and are not impacted by zero values, the Box-Cox transformation (\code{transfo = "boxcox"}) uses the formulation given in Equation 10 of Santos et al. (2018). Lambda is set to 0.25 accordingly. \cr \cr
The syntax of the power transformation allows a numeric or a string of characters. For example for a squared transformation, the following can be used: \code{transfo = 2}, \code{transfo = "2"} or \code{transfo = "^2"}. Negative values are allowed. Fraction values are not allowed (e.g., \code{"-1/2"} must instead be written \code{"-0.5"}).\cr \cr
## --- The epsilon value
\details{
Users wanting to use \code{FUN_CRIT} functions that are not included in
the package must create their own InputsCrit object accordingly.
The epsilon value is useful when \code{"log"} or \code{"inv"} transformations are used (to avoid calculation of the inverse or of the logarithm of zero). If an epsilon value is provided, then it is added to the observed and simulated variable time series at each time step and before the application of a transformation. The epsilon value has no effect when the \code{"boxcox"} transformation is used. The impact of this value and a recommendation about the epsilon value to use (usually one hundredth of average observation) are discussed in Pushpalatha et al. (2012) for NSE and in Santos et al. (2018) for KGE and KGE'. \cr \cr
## --- Single, multiple or composite criteria calculation
Users can set the following arguments as atomic or list: \code{FUN_CRIT}, \code{Obs}, \code{VarObs}, \code{BoolCrit}, \code{transfo}, \code{Weights}. If the list format is chosen, all the lists must have the same length. \cr
Calculation of a single criterion (e.g. NSE computed on discharge) is prepared by providing to \code{CreateInputsCrit} arguments atomics only. \cr
Calculation of multiple criteria (e.g. NSE computed on discharge and RMSE computed on discharge) is prepared by providing to \code{CreateInputsCrit} arguments lists except for \code{Weights} that must be set as \code{NULL}. \cr
Calculation of a composite criterion (e.g. the average between NSE computed on discharge and NSE computed on log of discharge) is prepared by providing to \code{CreateInputsCrit} arguments lists including \code{Weights}. \cr
\code{\link{ErrorCrit_RMSE}} cannot be used in a composite criterion since it is not a unitless value.
}
......@@ -63,64 +108,70 @@ library(airGR)
data(L0123001)
## 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)
## run period selection
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="01/01/1990 00:00"),
which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="31/12/1999 00:00"))
## calibration period selection
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"))
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
## preparation of RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
IndPeriod_Run = Ind_Run)
## simulation
Param <- c(734.568, -0.840, 109.809, 1.971)
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions,
Param = Param, FUN = RunModel_GR4J)
## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Nash-Sutcliffe Efficiency on log-transformed flows
transfo <- "log"
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run],
transfo = transfo)
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Nash-Sutcliffe Efficiency above a threshold (q75\%)
BoolCrit <- rep(TRUE, length(BasinObs$Qmm[Ind_Run]));
BoolCrit[BasinObs$Qmm[Ind_Run]<quantile(BasinObs$Qmm[Ind_Run], 0.75, na.rm = TRUE)] <- FALSE
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run],
BoolCrit = BoolCrit)
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Kling-Gupta Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Kling-Gupta Efficiency below a threshold (q10\%) on log-transformed flows
transfo <- "log"
BoolCrit <- rep(TRUE, length(BasinObs$Qmm[Ind_Run]));
BoolCrit[BasinObs$Qmm[Ind_Run]>quantile(BasinObs$Qmm[Ind_Run], 0.10, na.rm = TRUE)] <- FALSE
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE,
InputsModel = InputsModel,RunOptions = RunOptions,
Qobs = BasinObs$Qmm[Ind_Run], BoolCrit = BoolCrit, transfo = transfo)
OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param)
## single efficiency criterion: Nash-Sutcliffe Efficiency
InputsCritSingle <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE,
InputsModel = InputsModel, RunOptions = RunOptions,
Obs = list(BasinObs$Qmm[Ind_Run]),
VarObs = "Q", transfo = "",
Weights = NULL)
str(InputsCritSingle)
invisible(ErrorCrit(InputsCrit = InputsCritSingle, OutputsModel = OutputsModel))
## 2 efficiency criteria: RMSE and Nash-Sutcliffe Efficiency
InputsCritMulti <- CreateInputsCrit(FUN_CRIT = list(ErrorCrit_RMSE, ErrorCrit_NSE),
InputsModel = InputsModel, RunOptions = RunOptions,
Obs = list(BasinObs$Qmm[Ind_Run],
BasinObs$Qmm[Ind_Run]),
VarObs = list("Q", "Q"), transfo = list("", "sqrt"),
Weights = NULL)
str(InputsCritMulti)
invisible(ErrorCrit(InputsCrit = InputsCritMulti, OutputsModel = OutputsModel))
## efficiency composite criterion: Nash-Sutcliffe Efficiency mixing
## both raw and log-transformed flows
InputsCritCompo <- CreateInputsCrit(FUN_CRIT = list(ErrorCrit_NSE, ErrorCrit_NSE),
InputsModel = InputsModel, RunOptions = RunOptions,
Obs = list(BasinObs$Qmm[Ind_Run],
BasinObs$Qmm[Ind_Run]),
VarObs = list("Q", "Q"), transfo = list("", "log"),
Weights = list(0.4, 0.6))
str(InputsCritCompo)
invisible(ErrorCrit(InputsCrit = InputsCritCompo, OutputsModel = OutputsModel))
}
\author{
Laurent Coron
Olivier Delaigue, Laurent Coron, Guillaume Thirel
}
\references{
Pushpalatha, R., Perrin, C., Le Moine, N. and Andréassian, V. (2012).
A review of efficiency criteria suitable for evaluating low-flow simulations.
Journal of Hydrology, 420-421, 171-182, doi: 10.1016/j.jhydrol.2011.11.055.
\cr\cr
Santos, L., Thirel, G. and Perrin, C. (2018).
Technical note: Pitfalls in using log-transformed flows within the KGE criterion.
Hydrol. Earth Syst. Sci., 22, 4583-4591, doi: 10.5194/hess-22-4583-2018.
}
\seealso{
\code{\link{RunModel}}, \code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}}, \code{\link{CreateCalibOptions}}
\code{\link{RunModel}}, \code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}}, \code{\link{CreateCalibOptions}}, \code{\link{ErrorCrit}}
}
man/CreateInputsCrit_Lavenne.Rd 0 → 100644
View file @ 4bccfee2
\encoding{UTF-8}
\name{CreateInputsCrit_Lavenne}
\alias{CreateInputsCrit_Lavenne}
\title{Creation of the InputsCrit object for Lavenne Criterion}
\description{
Creation of the \code{InputsCrit} object required to the \code{\link{ErrorCrit}} function. This function defines a composite criterion based on the formula proposed by Lavenne et al. (2019).
}
\usage{
CreateInputsCrit_Lavenne(FUN_CRIT = ErrorCrit_KGE,
InputsModel,
RunOptions,
Obs,
VarObs = "Q",
AprParamR,
AprCrit = 1,
k = 0.15,
BoolCrit = NULL,
transfo = "sqrt",
epsilon = NULL)
}
\arguments{
\item{FUN_CRIT}{[function] error criterion function (e.g. \code{\link{ErrorCrit_KGE}}, \code{\link{ErrorCrit_NSE}}). Default \code{\link{ErrorCrit_KGE}}}
\item{InputsModel}{[object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details}
\item{RunOptions}{[object of class \emph{RunOptions}] see \code{\link{CreateRunOptions}} for details}
\item{Obs}{[numeric (atomic or list)] series of observed variable ([mm/time step] for discharge or SWE, [-] for SCA)}
\item{VarObs}{(optional) [character (atomic or list)] names of the observed variable (\code{"Q"} by default, or one of \code{"SCA"}, \code{"SWE"})}
\item{AprParamR}{[numeric] a priori parameter set from a donor catchment}
\item{AprCrit}{(optional) [numeric] performance criterion of the donor catchment (1 by default)}
\item{k}{(optional) [numeric] coefficient used for the weighted average between \code{FUN_CRIT} and the gap between the optimised parameter set and an a priori parameter set calculated with the function produced by \code{\link{CreateErrorCrit_GAPX}}}
\item{BoolCrit}{(optional) [boolean] boolean (the same length as \code{Obs}) giving the time steps to consider in the computation (all time steps are considered by default. See details)}
\item{transfo}{(optional) [character] name of the transformation applied to the variables (e.g. \code{""}, \code{"sqrt"}, \code{"log"}, \code{"inv"}, \code{"sort"}, \code{"boxcox"} or a numeric value for power transformation for \code{FUN_CRIT}. Default value is \code{"sqrt"}. See details of \code{\link{CreateInputsCrit}}}
\item{epsilon}{(optional) [numeric] small value to add to all observations and simulations for \code{FUN_CRIT} when \code{"log"} or \code{"inv"} transformations are used [same unit as \code{Obs}]. See details of \code{\link{CreateInputsCrit}}}
}
\value{
[list] object of class InputsCrit containing the data required to evaluate the model outputs (see \code{\link{CreateInputsCrit}} for more details).
\code{CreateInputsCrit_Lavenne} returns an object of class \emph{Compo}.
Items \code{Weights} of the criteria are respectively equal to \code{k} and \code{k * max(0, AprCrit)}.
To calculate the Lavenne criterion, it is necessary to use the \code{ErrorCrit} function as for any other composite criterion.
}
\details{
The parameters \code{FUN_CRIT}, \code{Obs}, \code{VarObs}, \code{BoolCrit}, \code{transfo}, and \code{epsilon} must be used as they would be used for \code{\link{CreateInputsCrit}} in the case of a single criterion.
\code{\link{ErrorCrit_RMSE}} cannot be used in a composite criterion since it is not a unitless value.
\code{CreateInputsCrit_Lavenne} creates a composite criterion in respect with Equations 1 and 2 of de Lavenne et al. (2019).
}
\examples{
library(airGR)
## loading catchment data
data(L0123001)
## preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
## calibration period selection
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"))
## preparation of RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
IndPeriod_Run = Ind_Run)
## simulation
Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions, Param = Param)
## The "a priori" parameters for the Lavenne formula
AprParamR <- c(X1 = 157, X2 = 0.8, X3 = 100, X4 = 1.5)
## Single efficiency criterion: GAPX with a priori parameters
IC_DL <- CreateInputsCrit_Lavenne(InputsModel = InputsModel,
RunOptions = RunOptions,
Obs = BasinObs$Qmm[Ind_Run],
AprParamR = AprParamR)
str(ErrorCrit(InputsCrit = IC_DL, OutputsModel = OutputsModel))
}
\author{
David Dorchies
}
\references{
de Lavenne, A., Andréassian, V., Thirel, G., Ramos, M.-H. and Perrin, C. (2019).
A Regularization Approach to Improve the Sequential Calibration of a Semidistributed Hydrological Model.
Water Resources Research 55, 8821–8839. \doi{10.1029/2018WR024266}
}
\seealso{
\code{\link{RunModel}}, \code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}}, \code{\link{CreateCalibOptions}}, \code{\link{ErrorCrit}}
}
man/CreateInputsModel.Rd
View file @ 4bccfee2
......@@ -3,15 +3,25 @@
\name{CreateInputsModel}
\alias{CreateInputsModel}
\alias{[.InputsModel}
\title{Creation of the InputsModel object required to the RunModel functions}
\description{
Creation of the \emph{InputsModel} object required to the \code{RunModel*} functions.
}
\usage{
CreateInputsModel(FUN_MOD, DatesR, Precip, PrecipScale = TRUE, PotEvap = NULL,
TempMean = NULL, TempMin = NULL, TempMax = NULL, ZInputs = NULL, HypsoData = NULL,
NLayers = 5, verbose = TRUE)
TempMean = NULL, TempMin = NULL, TempMax = NULL,
ZInputs = NULL, HypsoData = NULL, NLayers = 5,
Qupstream = NULL, LengthHydro = NULL, BasinAreas = NULL,
QupstrUnit = "mm", verbose = TRUE)
\method{[}{InputsModel}(x, i)
}
......@@ -20,11 +30,11 @@ CreateInputsModel(FUN_MOD, DatesR, Precip, PrecipScale = TRUE, PotEvap = NULL,
\item{DatesR}{[POSIXt] vector of dates required to create the GR model and CemaNeige module inputs}
\item{Precip}{[numeric] time series of total precipitation (catchment average) [mm], required to create the GR model and CemaNeige module inputs}
\item{Precip}{[numeric] time series of total precipitation (catchment average) [mm/time step], required to create the GR model and CemaNeige module inputs}
\item{PrecipScale}{(optional) [boolean] indicating if the mean of the precipitation interpolated on the elevation layers must be kept or not, required to create CemaNeige module inputs, default = \code{TRUE} (the mean of the precipitation is kept to the original value)}
\item{PotEvap}{[numeric] time series of potential evapotranspiration (catchment average) [mm], required to create the GR model inputs}
\item{PotEvap}{[numeric] time series of potential evapotranspiration (catchment average) [mm/time step], required to create the GR model inputs}
\item{TempMean}{(optional) [numeric] time series of mean air temperature [°C], required to create the CemaNeige module inputs}
......@@ -39,33 +49,46 @@ CreateInputsModel(FUN_MOD, DatesR, Precip, PrecipScale = TRUE, PotEvap = NULL,
\item{NLayers}{(optional) [numeric] integer giving the number of elevation layers requested [-], required to create CemaNeige module inputs, default=5}
\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, default = \code{TRUE}}
\item{Qupstream}{(optional) [numerical matrix] time series of upstream flows (catchment average), its unit is defined by the \code{QupstrUnit} parameter, required to create the SD model inputs. See details}
\item{LengthHydro}{(optional) [numeric] real giving the distance between the downstream outlet and each upstream inlet of the sub-catchment [km], required to create the SD model inputs . See details}
\item{BasinAreas}{(optional) [numeric] real giving the area of each upstream sub-catchment [km2] and the area of the downstream sub-catchment in the last item, required to create the SD model inputs . See details}
\item{QupstrUnit}{(optional) [character] unit of the flow in the argument \code{Qupstream}, available units are: "mm" for mm/time-step (default), "m3" for m3/time-step, "m3/s" and "l/s". See details}
\item{x}{[InputsModel] object of class InputsModel}
\item{i}{[integer] of the indices to subset a time series or [character] names of the elements to extract}
}
\value{
[list] object of class \emph{InputsModel} containing the data required to evaluate the model outputs; it can include the following:
\tabular{ll}{
\emph{$DatesR } \tab [POSIXlt] vector of dates \cr
\emph{$Precip } \tab [numeric] time series of total precipitation (catchment average) [mm] \cr
\emph{$PotEvap } \tab [numeric] time series of potential evapotranspiration (catchment average) [mm], \cr\tab defined if FUN_MOD includes GR4H, GR4J, GR5J, GR6J, GR2M or GR1A \cr \cr
\emph{$LayerPrecip } \tab [list] list of time series of precipitation (layer average) [mm], \cr\tab defined if \code{FUN_MOD} includes CemaNeige \cr \cr
\emph{$LayerTempMean } \tab [list] list of time series of mean air temperature (layer average) [°C], \cr\tab defined if \code{FUN_MOD} includes CemaNeige \cr \cr
\emph{$LayerFracSolidPrecip} \tab [list] list of time series of solid precipitation fraction (layer average) [-], \cr\tab defined if \code{FUN_MOD} includes CemaNeige \cr \cr
\tabular{ll}{
\emph{$DatesR } \tab [POSIXlt] vector of dates \cr
\emph{$Precip } \tab [numeric] time series of total precipitation (catchment average) [mm/time step] \cr
\emph{$PotEvap } \tab [numeric] time series of potential evapotranspiration (catchment average) [mm/time step], \cr\tab defined if FUN_MOD includes GR4H, GR5H, GR4J, GR5J, GR6J, GR2M or GR1A \cr \cr
\emph{$LayerPrecip } \tab [list] list of time series of precipitation (layer average) [mm/time step], \cr\tab defined if \code{FUN_MOD} includes CemaNeige \cr \cr
\emph{$LayerTempMean } \tab [list] list of time series of mean air temperature (layer average) [°C], \cr\tab defined if \code{FUN_MOD} includes CemaNeige \cr \cr
\emph{$LayerFracSolidPrecip} \tab [list] list of time series of solid precipitation fraction (layer average) [-], \cr\tab defined if \code{FUN_MOD} includes CemaNeige \cr \cr
}
}
\description{
Creation of the \emph{InputsModel} object required to the \code{RunModel*} functions.
}
\details{
Users wanting to use \code{FUN_MOD} functions that are not included in
the package must create their own InputsModel object accordingly.
Users wanting to use \code{FUN_MOD} functions that are not included in
the package must create their own InputsModel object accordingly. \cr
Please note that if CemaNeige is used, and \code{ZInputs} is different than \code{HypsoData}, then precipitation and temperature are interpolated with the \code{DataAltiExtrapolation_Valery} function.
Users wanting to use a semi-distributed (SD) model should provide valid \code{Qupstream}, \code{LengthHydro}, and \code{BasinAreas} arguments. Each upstream sub-catchment is described by an upstream flow time series (one column in \code{Qupstream} matrix), a distance between the downstream outlet and the upstream inlet (one item in \code{LengthHydro}) and an area (one item in \code{BasinAreas}).
The order of the columns or of the items should be consistent for all these parameters.
\code{BasinAreas} should contain one extra information (stored in the last item) which is the area of the downstream sub-catchment.
Upstream flows that are not related to a sub-catchment such as release or withdraw spots are identified by an area equal to \code{NA}, and if \code{unit="mm"} the upstream flow must be expressed in m3/time step instead of mm/time step which is not possible in absence of a related area.
Please note that the use of SD model requires to use the \code{\link{RunModel}} function instead of \code{\link{RunModel_GR4J}} or the other \code{RunModel_*} functions.
}
\examples{
library(airGR)
......@@ -73,30 +96,30 @@ library(airGR)
data(L0123001)
## 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)
## run period selection
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="01/01/1990 00:00"),
which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="31/12/1999 00:00"))
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"))
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
## simulation
Param <- c(734.568, -0.840, 109.809, 1.971)
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param,
Param <- c(X1 = 734.568, X2 = -0.840, X3 = 109.809, X4 = 1.971)
OutputsModel <- RunModel(InputsModel = InputsModel,
RunOptions = RunOptions, Param = Param,
FUN_MOD = RunModel_GR4J)
## results preview
plot(OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
}
......@@ -106,6 +129,6 @@ Laurent Coron
\seealso{
\code{\link{RunModel}}, \code{\link{CreateRunOptions}}, \code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}, \code{\link{DataAltiExtrapolation_Valery}}
\code{\link{RunModel}}, \code{\link{CreateRunOptions}}, \code{\link{CreateInputsCrit}},
\code{\link{CreateCalibOptions}}, \code{\link{DataAltiExtrapolation_Valery}}
}
man/CreateRunOptions.Rd
View file @ 4bccfee2
......@@ -8,11 +8,18 @@
\title{Creation of the RunOptions object required to the RunModel functions}
\description{
Creation of the RunOptions object required to the \code{RunModel*} functions.
}
\usage{
CreateRunOptions(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run,
IniStates = NULL, IniResLevels = NULL, Outputs_Cal = NULL,
Outputs_Sim = "all", RunSnowModule, MeanAnSolidPrecip = NULL,
verbose = TRUE)
CreateRunOptions(FUN_MOD, InputsModel,
IndPeriod_WarmUp = NULL, IndPeriod_Run,
IniStates = NULL, IniResLevels = NULL, Imax = NULL,
Outputs_Cal = NULL, Outputs_Sim = "all",
MeanAnSolidPrecip = NULL, IsHyst = FALSE,
warnings = TRUE, verbose = TRUE)
}
......@@ -21,52 +28,55 @@ CreateRunOptions(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run,
\item{InputsModel}{[object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details}
\item{IndPeriod_WarmUp}{(optional) [numeric] index of period to be used for the model warm-up [-]}
\item{IndPeriod_WarmUp}{(optional) [numeric] index of period to be used for the model warm-up [-]. See details}
\item{IndPeriod_Run}{[numeric] index of period to be used for the model run [-]}
\item{IndPeriod_Run}{[numeric] index of period to be used for the model run [-]. See details}
\item{IniStates}{(optional) [numeric] object of class \code{IniStates} [mm and °C], see \code{\link{CreateIniStates}} for details}
\item{IniResLevels}{(optional) [numeric] vector of initial fillings for the GR stores (2 or 3 values according to the model) [- and/or mm]; see details}
\item{IniResLevels}{(optional) [numeric] vector of initial fillings for the GR stores (4 values; use NA when not relevant for a given model) [- and/or mm]. See details}
\item{Imax}{(optional) [numeric] an atomic vector of the maximum capacity of the GR5H interception store [mm]; see \code{\link{RunModel_GR5H}}}
\item{Outputs_Cal}{(optional) [character] vector giving the outputs needed for the calibration \cr (e.g. c("Qsim")), the fewer outputs
the faster the calibration}
\item{Outputs_Sim}{(optional) [character] vector giving the requested outputs \cr (e.g. c(\code{"DatesR"}, \code{"Qsim"}, \code{"SnowPack"})), default = \code{"all"}}
\item{RunSnowModule}{(deprecated) [boolean] option indicating whether CemaNeige should be activated. Please adapt \code{FUN_MOD} instead}
\item{MeanAnSolidPrecip}{(optional) [numeric] vector giving the annual mean of average solid precipitation for each layer (computed from InputsModel if not defined) [mm/y]}
\item{IsHyst}{[boolean] boolean indicating if the hysteresis version of CemaNeige is used. See details}
\item{warnings}{(optional) [boolean] boolean indicating if the warning messages are shown, default = \code{TRUE}}
\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, default = \code{TRUE}}
}
\value{
[list] object of class \emph{RunOptions} containing the data required to evaluate the model outputs; it can include the following:
\tabular{ll}{
\emph{IndPeriod_WarmUp } \tab [numeric] index of period to be used for the model warm-up [-] \cr
\emph{IndPeriod_Run } \tab [numeric] index of period to be used for the model run [-] \cr
\emph{IniStates } \tab [numeric] vector of initial model states [mm and °C] \cr
\emph{IniResLevels } \tab [numeric] vector of initial filling rates for production and routing stores [-] and level fothe the exponential store for GR6J [mm]\cr
\emph{Outputs_Cal } \tab [character] character vector giving only the outputs needed for the calibration \cr
\emph{Outputs_Sim } \tab [character] character vector giving the requested outputs \cr
\emph{MeanAnSolidPrecip} \tab [numeric] vector giving the annual mean of average solid precipitation for each layer [mm/y] \cr
}
}
\description{
Creation of the RunOptions object required to the \code{RunModel*} functions.
\tabular{ll}{
\emph{IndPeriod_WarmUp } \tab [numeric] index of period to be used for the model warm-up [-] \cr
\emph{IndPeriod_Run } \tab [numeric] index of period to be used for the model run [-] \cr
\emph{IniStates } \tab [numeric] vector of initial model states [mm and °C] \cr
\emph{IniResLevels } \tab [numeric] vector of initial filling rates for production and routing stores [-] and level for the exponential store for GR6J [mm]\cr
\emph{Outputs_Cal } \tab [character] character vector giving only the outputs needed for the calibration \cr
\emph{Outputs_Sim } \tab [character] character vector giving the requested outputs \cr
\emph{Imax } \tab [numeric] vector giving the maximal capacity of the GR5H interception store \cr
\emph{MeanAnSolidPrecip} \tab [numeric] vector giving the annual mean of average solid precipitation for each layer [mm/y] \cr
}
}
\details{
Users wanting to use \code{FUN_MOD} functions that are not included in
Users wanting to use \code{FUN_MOD} functions that are not included in
the package must create their own \code{RunOptions} object accordingly.
##### Initialisation options #####
## --- IndPeriod_WarmUp and IndPeriod_Run
Since the hydrological models included in airGR are continuous models, meaning that internal states of the models are propagated to the next time step, \code{IndPeriod_WarmUp} and \code{IndPeriod_Run} must be continuous periods, represented by continuous indices values; no gaps are allowed. To calculate criteria or to calibrate a model over discontinuous periods, please see the \code{Bool_Crit} argument of the \code{\link{CreateInputsCrit}} function.
## --- Initialisation options
The model initialisation options can either be set to a default configuration or be defined by the user.
......@@ -76,27 +86,35 @@ A default configuration is used for initialisation if these vectors are not defi
(1) Default initialisation options:
\itemize{
\item \code{IndPeriod_WarmUp} default setting ensures a one-year warm-up using the time steps preceding the \code{IndPeriod_Run}.
\item \code{IndPeriod_WarmUp} default setting ensures a one-year warm-up using the time steps preceding the \code{IndPeriod_Run}.
The actual length of this warm-up might be shorter depending on data availability (no missing value of climate inputs being allowed in model input series).
\item \code{IniStates} and \code{IniResLevels} are automatically set to initialise all the model states at 0, except for the production and routing stores which are respectively initialised at 30 \% and 50 \% of their capacity. In case GR6J is used, the exponential store is initialised by default with 0 mm. This initialisation is made at the very beginning of the model call (i.e. at the beginning of \code{IndPeriod_WarmUp} or at the beginning of \code{IndPeriod_Run} if the warm-up period is disabled).
\item \code{IniStates} and \code{IniResLevels} are automatically set to initialise all the model states at 0, except for the production and routing stores levels which are respectively initialised at 30 \% and 50 \% of their capacity. In case GR5H is used with an interception store, the intercetion store level is initialised by default with 0 mm. In case GR6J is used, the exponential store level is initialised by default with 0 mm. This initialisation is made at the very beginning of the model call (i.e. at the beginning of \code{IndPeriod_WarmUp} or at the beginning of \code{IndPeriod_Run} if the warm-up period is disabled).
}
(2) Customisation of initialisation options:
\itemize{
\item \code{IndPeriod_WarmUp} can be used to specify the indices of the warm-up period (within the time series prepared in InputsModel). \cr
- remark 1: for most common cases, indices corresponding to one or several years preceding \code{IndPeriod_Run} are used (e.g. \code{IndPeriod_WarmUp = 1000:1365} and \code{IndPeriod_Run = 1366:5000)}. \cr
\item \code{IndPeriod_WarmUp} can be used to specify the indices of the warm-up period (within the time series prepared in InputsModel).
\itemize{
\item remark 1: for most common cases, indices corresponding to one or several years preceding \code{IndPeriod_Run} are used (e.g. \code{IndPeriod_WarmUp = 1000:1365} and \code{IndPeriod_Run = 1366:5000)}. \cr
However, it is also possible to perform a long-term initialisation if other indices than the warm-up ones are set in \code{IndPeriod_WarmUp} (e.g. \code{IndPeriod_WarmUp = c(1:5000, 1:5000, 1:5000, 1000:1365)}). \cr
- remark 2: it is also possible to completely disable the warm-up period when using \code{IndPeriod_WarmUp = 0L}. This is necessary if you want \code{IniStates} and / or \code{IniResLevels} to be the actual initial values of the model variables from your simulation (e.g. to perform a forecast form a given initial state).
\item \code{IniStates} and \code{IniResLevels} can be used to specify the initial model states. \cr
- remark 1: \code{IniStates} and \code{IniResLevels} can not be used with GR1A. \cr
- remark 2: if \code{IniStates} is used, two possibilities are offered:
- \code{IniStates} can be set to the \code{$StateEnd} output of a previous \code{RunModel} call, as \code{$StateEnd} already respects the correct format; \cr
- \code{IniStates} can be created with the \code{\link{CreateIniStates}} function.
- remark 3: in addition to \code{IniStates}, \code{IniResLevels} allows to set the filling rate of the production and routing stores for the GR models. For instance for GR4J and , GR5J: \code{IniResLevels = c(0.3, 0.5)} should be used to obtain initial fillings of 30 \% and 50 \% for the production and routing stores, respectively. For GR6J, \code{IniResLevels = c(0.3, 0.5, 0)} shold be use to obtain initial fillings of 30 \% and 50 \% for the production, routing stores and 0 mm for the exponential store, respectively. \code{IniResLevels} is optional and can only be used if \code{IniStates} is also defined (the state values corresponding to these two other stores in \code{IniStates} are not used in such case). \cr \cr
\item remark 2: it is also possible to completely disable the warm-up period when using \code{IndPeriod_WarmUp = 0L}. This is necessary if you want \code{IniStates} and/or \code{IniResLevels} to be the actual initial values of the model variables from your simulation (e.g. to perform a forecast form a given initial state).
}
\item \code{IniStates} and \code{IniResLevels} can be used to specify the initial model states.
\itemize{
\item remark 1: \code{IniStates} and \code{IniResLevels} can not be used with GR1A. \cr
\item remark 2: if \code{IniStates} is used, two possibilities are offered:\cr
- \code{IniStates} can be set to the \emph{$StateEnd} output of a previous \code{RunModel} call, as \emph{$StateEnd} already respects the correct format; \cr
- \code{IniStates} can be created with the \code{\link{CreateIniStates}} function.
\item remark 3: in addition to \code{IniStates}, \code{IniResLevels} allows to set the filling rate of the production and routing stores for the GR models. For instance for GR4J and GR5J: \code{IniResLevels = c(0.3, 0.5, NA, NA)} should be used to obtain initial fillings of 30 \% and 50 \% for the production and routing stores, respectively. For GR6J, \code{IniResLevels = c(0.3, 0.5, 0, NA)} should be used to obtain initial fillings of 30 \% and 50 \% for the production and routing stores levels and 0 mm for the exponential store level, respectively. For GR5H with an interception store, \code{IniResLevels = c(0.3, 0.5, NA, 0.4)} should be used to obtain initial fillings of 30 \%, 50 \% and 40 \% for the production, routing and interception stores levels, respectively. \code{IniResLevels} is optional and can only be used if \code{IniStates} is also defined (the state values corresponding to these two other stores in \code{IniStates} are not used in such case).
}
}
## --- CemaNeige version
If \code{IsHyst = FALSE}, the original CemaNeige version from Valéry et al. (2014) is used. \cr
If \code{IsHyst = TRUE}, the CemaNeige version from Riboust et al. (2019) is used. Compared to the original version, this version of CemaNeige needs two more parameters and it includes a representation of the hysteretic relationship between the Snow Cover Area (SCA) and the Snow Water Equivalent (SWE) in the catchment. The hysteresis included in airGR is the Modified Linear hysteresis (LH*); it is represented on panel b) of Fig. 3 in Riboust et al. (2019). Riboust et al. (2019) advise to use the LH* version of CemaNeige with parameters calibrated using an objective function combining 75 \% of KGE calculated on discharge simulated from a rainfall-runoff model compared to observed discharge and 5 \% of KGE calculated on SCA on 5 CemaNeige elevation bands compared to satellite (e.g. MODIS) SCA (see Eq. (18), Table 3 and Fig. 6). Riboust et al. (2019)'s tests were realized with GR4J as the chosen rainfall-runoff model. \cr
}
......@@ -107,28 +125,30 @@ library(airGR)
data(L0123001)
## 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)
## run period selection
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="01/01/1990 00:00"),
which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="31/12/1999 00:00"))
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"))
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
## simulation
Param <- c(734.568, -0.840, 109.809, 1.971)
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param,
Param <- c(X1 = 734.568, X2 = -0.840, X3 = 109.809, X4 = 1.971)
OutputsModel <- RunModel(InputsModel = InputsModel,
RunOptions = RunOptions, Param = Param,
FUN_MOD = RunModel_GR4J)
## results preview
plot(OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions,
Obs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
}
......@@ -139,6 +159,6 @@ Laurent Coron, Olivier Delaigue, Guillaume Thirel
\seealso{
\code{\link{RunModel}}, \code{\link{CreateInputsModel}}, \code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}, \code{\link{CreateIniStates}}
\code{\link{RunModel}}, \code{\link{CreateInputsModel}}, \code{\link{CreateInputsCrit}},
\code{\link{CreateCalibOptions}}, \code{\link{CreateIniStates}}, \code{\link{Imax}}
}
man/DataAltiExtrapolation_Valery.Rd
View file @ 4bccfee2
......@@ -8,17 +8,22 @@
\title{Altitudinal extrapolation of precipitation and temperature series described by A. Valery}
\description{
Function which extrapolates the precipitation and air temperature series for different elevation layers (method from Valéry, 2010).
}
\usage{
DataAltiExtrapolation_Valery(DatesR, Precip, PrecipScale = TRUE,
TempMean, TempMin = NULL, TempMax = NULL,
ZInputs, HypsoData, NLayers, verbose = TRUE)
TempMean, TempMin = NULL, TempMax = NULL,
ZInputs, HypsoData, NLayers, verbose = TRUE)
}
\arguments{
\item{DatesR}{[POSIXt] vector of dates}
\item{Precip}{[numeric] time series of daily total precipitation (catchment average) [mm]}
\item{Precip}{[numeric] time series of daily total precipitation (catchment average) [mm/time step]}
\item{PrecipScale}{(optional) [boolean] indicating if the mean of the precipitation interpolated on the elevation layers must be kept or not, required to create CemaNeige module inputs, default = \code{TRUE} (the mean of the precipitation is kept to the original value)}
......@@ -38,43 +43,43 @@ DataAltiExtrapolation_Valery(DatesR, Precip, PrecipScale = TRUE,
}
\value{
list containing the extrapolated series of precip. and air temp. on each elevation layer
\tabular{ll}{
\emph{$LayerPrecip } \tab [list] list of time series of daily precipitation (layer average) [mm] \cr
\emph{$LayerTempMean } \tab [list] list of time series of daily mean air temperature (layer average) [°C] \cr
\emph{$LayerTempMin } \tab [list] list of time series of daily min air temperature (layer average) [°C] \cr
\emph{$LayerTempMax } \tab [list] list of time series of daily max air temperature (layer average) [°C] \cr
\emph{$LayerFracSolidPrecip} \tab [list] list of time series of daily solid precip. fract. (layer average) [-] \cr
\emph{$ZLayers } \tab [numeric] vector of median elevation for each layer \cr
\tabular{ll}{
\emph{$LayerPrecip } \tab [list] list of time series of daily precipitation (layer average) [mm/time step] \cr
\emph{$LayerTempMean } \tab [list] list of time series of daily mean air temperature (layer average) [°C] \cr
\emph{$LayerTempMin } \tab [list] list of time series of daily min air temperature (layer average) [°C] \cr
\emph{$LayerTempMax } \tab [list] list of time series of daily max air temperature (layer average) [°C] \cr
\emph{$LayerFracSolidPrecip} \tab [list] list of time series of daily solid precip. fract. (layer average) [-] \cr
\emph{$ZLayers } \tab [numeric] vector of median elevation for each layer \cr
}
}
\description{
Function which extrapolates the precipitation and air temperature series for different elevation layers (method from Valéry, 2010).
}
\details{
Elevation layers of equal surface are created the 101 elevation quantiles (\code{HypsoData})
Elevation layers of equal surface are created the 101 elevation quantiles (\code{HypsoData})
and the number requested elevation layers (\code{NLayers}). \cr
Forcing data (precipitation and air temperature) are extrapolated using gradients from Valery (2010).
(e.g. gradP = 0.0004 [m-1] for France and gradT = 0.434 [°C/100m] for January, 1st). \cr
This function is used by the \code{CreateInputsModel} function.
This function is used by the \code{\link{CreateInputsModel}} function.
}
\author{
Laurent Coron, Audrey Valéry, Pierre Brigode, Olivier Delaigue, Guillaume Thirel
Laurent Coron, Audrey Valéry, Olivier Delaigue, Pierre Brigode, Guillaume Thirel
}
\references{
Turcotte, R., L.-G. Fortin, V. Fortin, J.-P. Fortin and J.-P. Villeneuve (2007),
Operational analysis of the spatial distribution and the temporal evolution of the snowpack water equivalent
in southern Quebec, Canada, Nordic Hydrology, 38(3), 211, doi:10.2166/nh.2007.009. \cr
Valéry, A. (2010), Modélisation précipitations-débit sous influence nivale ? : Elaboration d'un module neige
et évaluation sur 380 bassins versants, PhD thesis (in french), AgroParisTech, Paris, France. \cr
USACE (1956), Snow Hydrology, pp. 437, U.S. Army Coprs of Engineers (USACE) North Pacific Division, Portland, Oregon, USA.
Turcotte, R., Fortin, L.-G., Fortin, V., Fortin, J.-P. and Villeneuve, J.-P. (2007).
Operational analysis of the spatial distribution and the temporal evolution of the snowpack water equivalent in southern Quebec, Canada.
Nordic Hydrology, 38(3), 211, \doi{10.2166/nh.2007.009}.
\cr\cr
Valéry, A. (2010),
Modélisation précipitations-débit sous influence nivale ? : Elaboration d'un module neige et évaluation sur 380 bassins versants.
PhD thesis (in French), AgroParisTech - Cemagref Antony, Paris, France.
\cr\cr
USACE (1956),
Snow Hydrology, pp. 437.
U.S. Army Coprs of Engineers (USACE) North Pacific Division, Portland, Oregon, USA.
}
......
man/ErrorCrit.Rd
View file @ 4bccfee2
......@@ -8,8 +8,13 @@
\title{Error criterion using the provided function}
\description{
Function which computes an error criterion with the provided function.
}
\usage{
ErrorCrit(InputsCrit, OutputsModel, FUN_CRIT, warnings = TRUE, verbose = TRUE)
ErrorCrit(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE)
}
......@@ -18,19 +23,35 @@ ErrorCrit(InputsCrit, OutputsModel, FUN_CRIT, warnings = TRUE, verbose = TRUE)
\item{OutputsModel}{[object of class \emph{OutputsModel}] see \code{\link{RunModel_GR4J}} or \code{\link{RunModel_CemaNeigeGR4J}} for details}
\item{FUN_CRIT}{[function] error criterion function (e.g. \code{ErrorCrit_RMSE}, \code{ErrorCrit_NSE})}
\item{warnings}{(optional) [boolean] boolean indicating if the warning messages are shown, default = \code{TRUE}}
\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, default = \code{TRUE}}
}
\value{
[list] list containing the function outputs, see \code{\link{ErrorCrit_RMSE}} or \code{\link{ErrorCrit_NSE}} for details
}
\description{
Function which computes an error criterion with the provided function.
\value{
If \code{InputsCrit} is of class \emph{Single}:
\tabular{ll}{
[list] containing the \code{ErrorCrit_*} functions outputs, see \code{\link{ErrorCrit_RMSE}} or \code{\link{ErrorCrit_NSE}} for details \cr
}
If \code{InputsCrit} is of class \emph{Multi}:
\tabular{ll}{
[list] of list containing the \code{ErrorCrit_*} functions outputs, see \code{\link{ErrorCrit_RMSE}} or \code{\link{ErrorCrit_NSE}} for details \cr
}
If \code{InputsCrit} is of class \emph{Compo}:
\tabular{ll}{
\emph{$CritValue } \tab [numeric] value of the composite criterion \cr
\emph{$CritName } \tab [character] name of the composite criterion \cr
\emph{$CritBestValue} \tab [numeric] theoretical best criterion value \cr
\emph{$Multiplier } \tab [numeric] integer indicating whether the criterion is indeed an error (+1) or an efficiency (-1) \cr
\emph{$CritCompo$MultiCritValues} \tab [numeric] values of the sub-criteria \cr
\emph{$CritCompo$MultiCritNames} \tab [numeric] names of the sub-criteria \cr
\emph{$CritCompo$MultiCritWeights} \tab [character] weighted values of the sub-criteria \cr
\emph{$MultiCrit } \tab [list] of list containing the \code{ErrorCrit_*} functions outputs, see \code{\link{ErrorCrit_NSE}} or \code{\link{ErrorCrit_KGE}} for details \cr
}
}
......@@ -41,64 +62,58 @@ library(airGR)
data(L0123001)
## 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)
## run period selection
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="01/01/1990 00:00"),
which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="31/12/1999 00:00"))
## calibration period selection
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"))
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
## preparation of RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
IndPeriod_Run = Ind_Run)
## simulation
Param <- c(734.568, -0.840, 109.809, 1.971)
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions,
Param = Param, FUN = RunModel_GR4J)
## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Nash-Sutcliffe Efficiency on log-transformed flows
transfo <- "log"
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run],
transfo = transfo)
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Nash-Sutcliffe Efficiency above a threshold (q75\%)
BoolCrit <- rep(TRUE, length(BasinObs$Qmm[Ind_Run]))
BoolCrit[BasinObs$Qmm[Ind_Run]<quantile(BasinObs$Qmm[Ind_Run], 0.75, na.rm = TRUE)] <- FALSE
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run],
BoolCrit = BoolCrit)
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Kling-Gupta Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Kling-Gupta Efficiency below a threshold (q10\%) on log-transformed flows
transfo <- "log"
BoolCrit <- rep(TRUE, length(BasinObs$Qmm[Ind_Run]))
BoolCrit[BasinObs$Qmm[Ind_Run]>quantile(BasinObs$Qmm[Ind_Run], 0.10, na.rm = TRUE)] <- FALSE
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run],
BoolCrit = BoolCrit, transfo = transfo)
OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
Param <- c(X1 = 257.238, X2 = 1.012, X3 = 88.235, X4 = 2.208)
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions, Param = Param)
## single efficiency criterion: Nash-Sutcliffe Efficiency
InputsCritSingle <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE,
InputsModel = InputsModel, RunOptions = RunOptions,
Obs = list(BasinObs$Qmm[Ind_Run]),
VarObs = "Q", transfo = "",
Weights = NULL)
str(ErrorCrit(InputsCrit = InputsCritSingle, OutputsModel = OutputsModel))
## 2 efficiency critera: RMSE and the Nash-Sutcliffe Efficiency
InputsCritMulti <- CreateInputsCrit(FUN_CRIT = list(ErrorCrit_RMSE, ErrorCrit_NSE),
InputsModel = InputsModel, RunOptions = RunOptions,
Obs = list(BasinObs$Qmm[Ind_Run],
BasinObs$Qmm[Ind_Run]),
VarObs = list("Q", "Q"), transfo = list("", "sqrt"),
Weights = NULL)
str(ErrorCrit(InputsCrit = InputsCritMulti, OutputsModel = OutputsModel))
## efficiency composite criterion: Nash-Sutcliffe Efficiency mixing
## both raw and log-transformed flows
InputsCritCompo <- CreateInputsCrit(FUN_CRIT = list(ErrorCrit_NSE, ErrorCrit_NSE),
InputsModel = InputsModel, RunOptions = RunOptions,
Obs = list(BasinObs$Qmm[Ind_Run],
BasinObs$Qmm[Ind_Run]),
VarObs = list("Q", "Q"), transfo = list("", "log"),
Weights = list(0.4, 0.6))
str(ErrorCrit(InputsCrit = InputsCritCompo, OutputsModel = OutputsModel))
}
\author{
Laurent Coron
Olivier Delaigue
}
\seealso{
\code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}}, \code{\link{ErrorCrit_KGE}}
\code{\link{CreateInputsCrit}}, \code{\link{ErrorCrit_RMSE}}, \code{\link{ErrorCrit_NSE}},
\code{\link{ErrorCrit_KGE}}, \code{\link{ErrorCrit_KGE2}}
}
man/ErrorCrit_KGE.Rd
View file @ 4bccfee2
......@@ -8,6 +8,11 @@
\title{Error criterion based on the KGE formula}
\description{
Function which computes an error criterion based on the KGE formula proposed by Gupta et al. (2009).
}
\usage{
ErrorCrit_KGE(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE)
}
......@@ -38,14 +43,9 @@ ErrorCrit_KGE(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE)
}
\description{
Function which computes an error criterion based on the KGE formula proposed by Gupta et al. (2009).
}
\details{
In addition to the criterion value, the function outputs include a multiplier (-1 or +1) which allows
the use of the function for model calibration: the product CritValue * Multiplier is the criterion to be minimised (Multiplier = -1 for KGE).\cr\cr
In addition to the criterion value, the function outputs include a multiplier (-1 or +1) which allows
the use of the function for model calibration: the product \eqn{CritValue \times Multiplier} is the criterion to be minimised (Multiplier = -1 for KGE).\cr\cr
The KGE formula is
\deqn{KGE = 1 - \sqrt{(r - 1)^2 + (\alpha - 1)^2 + (\beta - 1)^2}}{KGE = 1 - sqrt((r - 1)² + (\alpha - 1)² + (\beta - 1)²)}
with the following sub-criteria:
......@@ -56,19 +56,58 @@ with the following sub-criteria:
\examples{
## see example of the ErrorCrit function
library(airGR)
## loading catchment data
data(L0123001)
## preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
## run period selection
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"))
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
## simulation
Param <- c(X1 = 734.568, X2 = -0.840, X3 = 109.809, X4 = 1.971)
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions,
Param = Param, FUN = RunModel_GR4J)
## efficiency criterion: Kling-Gupta Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Kling-Gupta Efficiency on square-root-transformed flows
transfo <- "sqrt"
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run],
transfo = transfo)
OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Kling-Gupta Efficiency above a threshold (quant. 75 \%)
BoolCrit <- BasinObs$Qmm[Ind_Run] >= quantile(BasinObs$Qmm[Ind_Run], 0.75, na.rm = TRUE)
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run],
BoolCrit = BoolCrit)
OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
}
\author{
Laurent Coron
Laurent Coron, Olivier Delaigue
}
\references{
Gupta, H. V., Kling, H., Yilmaz, K. K. and Martinez, G. F. (2009),
Decomposition of the mean squared error and NSE performance criteria: Implications
for improving hydrological modelling, Journal of Hydrology, 377(1-2), 80-91, doi:10.1016/j.jhydrol.2009.08.003. \cr
Gupta, H. V., Kling, H., Yilmaz, K. K. and Martinez, G. F. (2009).
Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling.
Journal of Hydrology, 377(1-2), 80-91, \doi{10.1016/j.jhydrol.2009.08.003}.
}
......
man/ErrorCrit_KGE2.Rd
View file @ 4bccfee2
......@@ -8,6 +8,11 @@
\title{Error criterion based on the KGE' formula}
\description{
Function which computes an error criterion based on the KGE' formula proposed by Kling et al. (2012).
}
\usage{
ErrorCrit_KGE2(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE)
}
......@@ -38,14 +43,9 @@ ErrorCrit_KGE2(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE)
}
\description{
Function which computes an error criterion based on the KGE' formula proposed by Kling et al. (2012).
}
\details{
In addition to the criterion value, the function outputs include a multiplier (-1 or +1) which allows
the use of the function for model calibration: the product CritValue * Multiplier is the criterion to be minimised (Multiplier = -1 for KGE2).\cr\cr
In addition to the criterion value, the function outputs include a multiplier (-1 or +1) which allows
the use of the function for model calibration: the product \eqn{CritValue \times Multiplier} is the criterion to be minimised (Multiplier = -1 for KGE2).\cr\cr
The KGE' formula is
\deqn{KGE' = 1 - \sqrt{(r - 1)^2 + (\gamma - 1)^2 + (\beta - 1)^2}}{KGE' = 1 - sqrt((r - 1)² + (\gamma - 1)² + (\beta - 1)²)}
with the following sub-criteria:
......@@ -56,22 +56,62 @@ with the following sub-criteria:
\examples{
## see example of the ErrorCrit function
library(airGR)
## loading catchment data
data(L0123001)
## preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
## run period selection
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"))
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
## simulation
Param <- c(X1 = 734.568, X2 = -0.840, X3 = 109.809, X4 = 1.971)
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions,
Param = Param, FUN = RunModel_GR4J)
## efficiency criterion: Kling-Gupta Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE2, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_KGE2(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Kling-Gupta Efficiency on square-root-transformed flows
transfo <- "sqrt"
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE2, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run],
transfo = transfo)
OutputsCrit <- ErrorCrit_KGE2(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Kling-Gupta Efficiency above a threshold (quant. 75 \%)
BoolCrit <- BasinObs$Qmm[Ind_Run] >= quantile(BasinObs$Qmm[Ind_Run], 0.75, na.rm = TRUE)
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE2, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run],
BoolCrit = BoolCrit)
OutputsCrit <- ErrorCrit_KGE2(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
}
\author{
Laurent Coron
Laurent Coron, Olivier Delaigue
}
\references{
Gupta, H. V., Kling, H., Yilmaz, K. K. and Martinez, G. F. (2009),
Decomposition of the mean squared error and NSE performance criteria: Implications
for improving hydrological modelling, Journal of Hydrology, 377(1-2), 80-91, doi:10.1016/j.jhydrol.2009.08.003. \cr
Kling, H., Fuchs, M. and Paulin, M. (2012),
Runoff conditions in the upper Danube basin under an ensemble of climate change scenarios,
Journal of Hydrology, 424-425, 264-277, doi:10.1016/j.jhydrol.2012.01.011.
Gupta, H. V., Kling, H., Yilmaz, K. K. and Martinez, G. F. (2009).
Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling.
Journal of Hydrology, 377(1-2), 80-91, \doi{10.1016/j.jhydrol.2009.08.003}.
\cr\cr
Kling, H., Fuchs, M. and Paulin, M. (2012).
Runoff conditions in the upper Danube basin under an ensemble of climate change scenarios.
Journal of Hydrology, 424-425, 264-277, \doi{10.1016/j.jhydrol.2012.01.011}.
}
......
man/ErrorCrit_NSE.Rd
View file @ 4bccfee2
......@@ -13,6 +13,11 @@ ErrorCrit_NSE(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE)
}
\description{
Function which computes an error criterion based on the NSE formula proposed by Nash & Sutcliffe (1970).
}
\arguments{
\item{InputsCrit}{[object of class \emph{InputsCrit}] see \code{\link{CreateInputsCrit}} for details}
......@@ -36,30 +41,65 @@ ErrorCrit_NSE(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE)
}
\description{
Function which computes an error criterion based on the NSE formula proposed by Nash & Sutcliffe (1970).
}
\details{
In addition to the criterion value, the function outputs include a multiplier (-1 or +1) which allows
the use of the function for model calibration: the product CritValue * Multiplier is the criterion to be minimised
In addition to the criterion value, the function outputs include a multiplier (-1 or +1) which allows
the use of the function for model calibration: the product \eqn{CritValue \times Multiplier} is the criterion to be minimised
(Multiplier = -1 for NSE).
}
\examples{
## see example of the ErrorCrit function
library(airGR)
## loading catchment data
data(L0123001)
## preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
## run period selection
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"))
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
## simulation
Param <- c(X1 = 734.568, X2 = -0.840, X3 = 109.809, X4 = 1.971)
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions,
Param = Param, FUN = RunModel_GR4J)
## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Nash-Sutcliffe Efficiency on log-transformed flows
transfo <- "log"
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run],
transfo = transfo)
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Kling-Gupta Efficiency above a threshold (quant. 75 \%)
BoolCrit <- BasinObs$Qmm[Ind_Run] >= quantile(BasinObs$Qmm[Ind_Run], 0.75, na.rm = TRUE)
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run],
BoolCrit = BoolCrit)
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
}
\author{
Laurent Coron
Laurent Coron, Olivier Delaigue
}
\references{
Nash, J.E. and Sutcliffe, J.V. (1970),
River flow forecasting through conceptual models part 1.
A discussion of principles, Journal of Hydrology, 10(3), 282-290, doi:10.1016/0022-1694(70)90255-6. \cr
Nash, J. E. and Sutcliffe, J. V. (1970).
River flow forecasting through conceptual models.
Part 1 - A discussion of principles.
Journal of Hydrology, 10(3), 282-290, \doi{10.1016/0022-1694(70)90255-6}.
}
......
man/ErrorCrit_RMSE.Rd
View file @ 4bccfee2
......@@ -3,6 +3,8 @@
\name{ErrorCrit_RMSE}
\alias{ErrorCrit_RMSE}
\title{Error criterion based on the RMSE}
......@@ -11,6 +13,11 @@ ErrorCrit_RMSE(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE)
}
\description{
Function which computes an error criterion based on the root-mean-square error (RMSE).
}
\arguments{
\item{InputsCrit}{[object of class \emph{InputsCrit}] see \code{\link{CreateInputsCrit}} for details}
......@@ -25,34 +32,68 @@ ErrorCrit_RMSE(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE)
\value{
[list] list containing the function outputs organised as follows:
\tabular{ll}{
\emph{$CritValue } \tab [numeric] value of the criterion \cr
\emph{$CritName } \tab [character] name of the criterion \cr
\emph{$CritBestValue } \tab [numeric] theoretical best criterion value \cr
\emph{$Multiplier } \tab [numeric] integer indicating whether the criterion is indeed an error (+1) or an efficiency (-1) \cr
\emph{$Ind_notcomputed} \tab [numeric] indices of the time steps where \emph{InputsCrit$BoolCrit} = \code{FALSE} or no data is available \cr
\emph{$CritValue } \tab [numeric] value of the criterion \cr
\emph{$CritName } \tab [character] name of the criterion \cr
\emph{$CritBestValue } \tab [numeric] theoretical best criterion value \cr
\emph{$Multiplier } \tab [numeric] integer indicating whether the criterion is indeed an error (+1) or an efficiency (-1) \cr
\emph{$Ind_notcomputed} \tab [numeric] indices of the time steps where \emph{InputsCrit$BoolCrit} = \code{FALSE} or no data is available \cr
}
}
\description{
Function which computes an error criterion based on the root mean square error (RMSE).
}
\details{
In addition to the criterion value, the function outputs include a multiplier (-1 or +1) which allows
the use of the function for model calibration: the product CritValue * Multiplier is the criterion to be minimised
the use of the function for model calibration: the product \eqn{CritValue \times Multiplier} is the criterion to be minimised
(Multiplier = +1 for RMSE).
}
\examples{
## see example of the ErrorCrit function
library(airGR)
## loading catchment data
data(L0123001)
## preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
## run period selection
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"))
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
## simulation
Param <- c(X1 = 734.568, X2 = -0.840, X3 = 109.809, X4 = 1.971)
OutputsModel <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions,
Param = Param, FUN = RunModel_GR4J)
## efficiency criterion: root-mean-square error
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_RMSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_RMSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: root-mean-square error on log-transformed flows
transfo <- "log"
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_RMSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run],
transfo = transfo)
OutputsCrit <- ErrorCrit_RMSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
## efficiency criterion: Kling-Gupta Efficiency above a threshold (quant. 75 \%)
BoolCrit <- BasinObs$Qmm[Ind_Run] >= quantile(BasinObs$Qmm[Ind_Run], 0.75, na.rm = TRUE)
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_RMSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run],
BoolCrit = BoolCrit)
OutputsCrit <- ErrorCrit_RMSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
}
\author{
Laurent Coron, Ludovic Oudin
Laurent Coron, Ludovic Oudin, Olivier Delaigue
}
......
man/Imax.Rd 0 → 100644
View file @ 4bccfee2
\encoding{UTF-8}
\name{Imax}
\alias{Imax}
\title{Computation of the maximum capacity of the GR5H interception store}
\description{
Function which calculates the maximal capacity of the GR5H interception store. This function compares the interception evapotranspiration from the GR5H interception store for different maximal capacity values with the interception evapotranspiration classically used in the daily GR models (e.g. GR4J). Among all the \code{TestedValues}, the value that gives the closest interception evapotranspiration flux over the whole period is kept.
}
\usage{
Imax(InputsModel, IndPeriod_Run,
TestedValues = seq(from = 0.1, to = 3, by = 0.1))
}
\arguments{
\item{InputsModel}{[object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details}
\item{IndPeriod_Run}{[numeric] index of period to be used for the model run [-]}
\item{TestedValues}{[numeric] vector of tested Imax values [mm]}
}
\value{
Optimal Imax value [mm].
}
\examples{
library(airGR)
## loading catchment data
data(L0123003)
## preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR5H, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
## run period selection
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d \%H")=="2006-01-01 00"),
which(format(BasinObs$DatesR, format = "\%Y-\%m-\%d \%H")=="2006-12-31 23"))
## Imax computation
Imax <- Imax(InputsModel = InputsModel, IndPeriod_Run = Ind_Run,
TestedValues = seq(from = 0, to = 3, by = 0.2))
## preparation of the RunOptions object
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR5H, Imax = Imax,
InputsModel = InputsModel, IndPeriod_Run = Ind_Run)
## simulation
Param <- c(X1 = 706.912, X2 = -0.163, X3 = 188.880, X4 = 2.575, X5 = 0.104)
OutputsModel <- RunModel_GR5H(InputsModel = InputsModel,
RunOptions = RunOptions, Param = Param)
## results preview
plot(OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
}
\author{
Guillaume Thirel, Olivier Delaigue
}
\references{
Ficchi, A. (2017).
An adaptive hydrological model for multiple time-steps:
Diagnostics and improvements based on fluxes consistency.
PhD thesis, UPMC - Irstea Antony, Paris, France.
\cr\cr
Ficchi, A., Perrin, C. and Andréassian, V. (2019).
Hydrological modelling at multiple sub-daily time steps: model improvement via flux-matching.
Journal of Hydrology, 575, 1308-1327, \doi{10.1016/j.jhydrol.2019.05.084}.
}
\seealso{
\code{\link{RunModel_GR5H}},
\code{\link{CreateInputsModel}}, \code{\link{CreateRunOptions}}.
}