An error occurred while loading the file. Please try again.
-
Delaigue Olivier authored
Signed-off-by:Olivier Delaigue <olivier.delaigue@irstea.fr>
ec2371ca
#*************************************************************************************************#' Creation of the RunOptions object required to the RunModel functions.#'#' Users wanting to use FUN_MOD functions that are not included in #' the package must create their own RunOptions object accordingly.#' #' ##### Initialisation options ######' #' The model initialisation options can either be set to a default configuration or be defined by the user.#' #' This is done via three vectors: \cr \emph{IndPeriod_WarmUp}, \emph{IniStates}, \emph{IniResLevels}. \cr#' A default configuration is used for initialisation if these vectors are not defined.#' #' (1) Default initialisation options:#' #' \itemize{#' \item \emph{IndPeriod_WarmUp} default setting ensures a one-year warm-up using the time-steps preceding the \emph{IndPeriod_Run}. #' The actual length of this warm-up might be shorter depending on data availability (no missing value being allowed on model input series).#' #' \item \emph{IniStates} and \emph{IniResLevels} are automatically set to initialise all the model states at 0, except for the production and routing stores which are initialised at 50\% of their capacity. This initialisation is made at the very beginning of the model call (i.e. at the beginning of \emph{IndPeriod_WarmUp} or at the beginning of IndPeriod_Run if the warm-up period is disabled).#' }#' #' (2) Customisation of initialisation options:#' #' \itemize{#' \item \emph{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 \emph{IndPeriod_Run} are used (e.g. \emph{IndPeriod_WarmUp <- 1000:1365} and \emph{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 \emph{IndPeriod_WarmUp} (e.g. \emph{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 \emph{IndPeriod_WarmUp <- 0}.#' #' \item \emph{IniStates} and \emph{IniResLevels} can be used to specify the initial model states. \cr#' - remark 1: if \emph{IniStates} is used, all model states must be provided (e.g. 60 floats [mm] are required for GR4J, GR5J and GR6J; 60+2*NLayers floats [mm] are required for CemaNeigeGR4J, CemaNeigeGR5J and CemaNeigeGR6J; see fortran source code for details). \cr#' - remark 2: in addition to \emph{IniStates}, \emph{IniResLevels} allows to set the filling rate of the production and routing stores for the GR models. For instance for GR4J, GR5J and GR6J: \emph{IniResLevels <- c(0.3,0.5)} should be used to obtain initial fillings of 30\% and 50\% for the production and routing stores, respectively. \emph{IniResLevels} is optional and can only be used if \emph{IniStates} is also defined (the state values corresponding to these two stores in \emph{IniStates} are not used in such case). \cr \cr#' }#*************************************************************************************************#' @title Creation of the RunOptions object required to the RunModel functions#' @author Laurent Coron (June 2014)#' @seealso \code{\link{RunModel}}, \code{\link{CreateInputsModel}}, \code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}#' @example tests/example_RunModel.R#' @encoding UTF-8#' @export#_FunctionInputs__________________________________________________________________________________#' @param FUN_MOD [function] hydrological model function (e.g. RunModel_GR4J, RunModel_CemaNeigeGR4J)#' @param InputsModel [object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details#' @param IndPeriod_WarmUp (optional) [numeric] index of period to be used for the model warm-up [-]#' @param IndPeriod_Run [numeric] index of period to be used for the model run [-]#' @param IniStates (optional) [numeric] vector of initial model states [mm]#' @param IniResLevels (optional) [numeric] vector of initial filling rates for production and routing stores (2 values between 0 and 1) [-]#' @param Outputs_Cal (optional) [character] vector giving the outputs needed for the calibration \cr (e.g. c("Qsim")), the least outputs the fastest the calibration#' @param Outputs_Sim (optional) [character] vector giving the requested outputs \cr (e.g. c("DatesR","Qsim","SnowPack")), default="all"#' @param RunSnowModule (optional) [boolean] option indicating whether CemaNeige should be activated, default=TRUE#' @param MeanAnSolidPrecip (optional) [numeric] vector giving the annual mean of average solid precipitation for each layer (computed from InputsModel if not defined) [mm/y]#' @param quiet (optional) [boolean] boolean indicating if the function is run in quiet mode or not, default=FALSE#_FunctionOutputs_________________________________________________________________________________#' @return [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] \cr#' \emph{IniResLevels } \tab [numeric] vector of initial filling rates for production and routing stores [-] \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{RunSnowModule } \tab [boolean] option indicating whether CemaNeige should be activated \cr#' \emph{MeanAnSolidPrecip} \tab [numeric] vector giving the annual mean of average solid precipitation for each layer [mm/y] \cr#' }#**************************************************************************************************'CreateRunOptions <- function(FUN_MOD,InputsModel,IndPeriod_WarmUp=NULL,IndPeriod_Run,IniStates=NULL,IniResLevels=NULL, Outputs_Cal=NULL,Outputs_Sim="all",RunSnowModule=TRUE,MeanAnSolidPrecip=NULL,quiet=FALSE){7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
ObjectClass <- NULL;
##check_FUN_MOD
BOOL <- FALSE;
if(identical(FUN_MOD,RunModel_GR4J) | identical(FUN_MOD,RunModel_GR5J) | identical(FUN_MOD,RunModel_GR6J)){
ObjectClass <- c(ObjectClass,"GR","daily");
BOOL <- TRUE;
}
if(identical(FUN_MOD,RunModel_CemaNeige)){
ObjectClass <- c(ObjectClass,"CemaNeige","daily");
BOOL <- TRUE;
}
if(identical(FUN_MOD,RunModel_CemaNeigeGR4J) | identical(FUN_MOD,RunModel_CemaNeigeGR5J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)){
ObjectClass <- c(ObjectClass,"GR","CemaNeige","daily");
BOOL <- TRUE;
}
if(!BOOL){ stop("incorrect FUN_MOD for use in CreateRunOptions \n"); return(NULL); }
##check_InputsModel
if(!inherits(InputsModel,"InputsModel")){
stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }
if("daily" %in% ObjectClass & !inherits(InputsModel,"daily")){
stop("InputsModel must be of class 'daily' \n"); return(NULL); }
##check_IndPeriod_Run
if(!is.vector( IndPeriod_Run)){ stop("IndPeriod_Run must be a vector of numeric values \n"); return(NULL); }
if(!is.numeric(IndPeriod_Run)){ stop("IndPeriod_Run must be a vector of numeric values \n"); return(NULL); }
if(identical(as.integer(IndPeriod_Run),as.integer(seq(from=IndPeriod_Run[1],to=tail(IndPeriod_Run,1),by=1)))==FALSE){
stop("IndPeriod_Run must be a continuous sequence of integers \n"); return(NULL); }
if(storage.mode(IndPeriod_Run)!="integer"){ stop("IndPeriod_Run should be of type integer \n"); return(NULL); }
##check_IndPeriod_WarmUp
WTxt <- NULL;
if(is.null(IndPeriod_WarmUp)){
WTxt <- paste(WTxt,"\t Model warm-up period not defined -> default configuration used \n",sep="");
##If_the_run_period_starts_at_the_very_beginning_of_the_time_series
if(IndPeriod_Run[1]==as.integer(1)){
IndPeriod_WarmUp <- as.integer(0);
WTxt <- paste(WTxt,"\t No data were found for model warm-up! \n",sep="");
##We_look_for_the_longest_period_preceeding_the_run_period_with_a_maximum_of_one_year
} else {
TmpDateR <- InputsModel$DatesR[IndPeriod_Run[1]] - 365*24*60*60; ### minimal date to start the warmup
IndPeriod_WarmUp <- which(InputsModel$DatesR==max(InputsModel$DatesR[1],TmpDateR)) : (IndPeriod_Run[1]-1);
if("daily" %in% ObjectClass){ TimeStep <- as.integer( 24*60*60); }
if(length(IndPeriod_WarmUp)*TimeStep/(365*24*60*60)>=1){
WTxt <- paste(WTxt,"\t The year preceding the run period is used \n",sep="");
} else {
WTxt <- paste(WTxt,"\t Less than a year (without missing values) was found for model warm-up: \n",sep="");
WTxt <- paste(WTxt,"\t Only ",length(IndPeriod_WarmUp)," time-steps are used! \n",sep="");
}
}
}
if(!is.null(IndPeriod_WarmUp)){
if(!is.vector( IndPeriod_WarmUp)){ stop("IndPeriod_Run must be a vector of numeric values \n"); return(NULL); }
if(!is.numeric(IndPeriod_WarmUp)){ stop("IndPeriod_Run must be a vector of numeric values \n"); return(NULL); }
if(storage.mode(IndPeriod_WarmUp)!="integer"){ stop("IndPeriod_Run should be of type integer \n"); return(NULL); }
if(identical(IndPeriod_WarmUp,as.integer(0))){
WTxt <- paste(WTxt,"\t No warm-up period is used! \n",sep=""); }
if((IndPeriod_Run[1]-1)!=tail(IndPeriod_WarmUp,1)){
WTxt <- paste(WTxt,"\t Model warm-up period is not directly before the model run period \n",sep=""); }
}
if(!is.null(WTxt) & !quiet){ warning(WTxt); }
##check_IniStates_and_IniResLevels
if(is.null(IniStates) & is.null(IniResLevels) & !quiet){
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
warning("\t Model states initialisation not defined -> default configuration used \n"); }
if("GR" %in% ObjectClass){
if("daily" %in% ObjectClass){ NH <- 20; }
} else {
NH <- 0;
}
if("CemaNeige" %in% ObjectClass){ NLayers <- length(InputsModel$LayerPrecip); } else { NLayers <- 0; }
NState <- 3*NH + 2*NLayers;
if(!is.null(IniStates)){
if(!is.vector( IniStates) ){ stop("IniStates must be a vector of numeric values \n"); return(NULL); }
if(!is.numeric(IniStates) ){ stop("IniStates must be a vector of numeric values \n"); return(NULL); }
if(length(IniStates)!=NState){ stop(paste("the length of IniStates must be ",NState," for the chosen FUN_MOD \n",sep="")); return(NULL); }
} else {
IniStates <- as.double(rep(0.0,NState));
}
if(!is.null(IniResLevels)){
if(!is.vector(IniResLevels) ){ stop("IniResLevels must be a vector of numeric values \n"); return(NULL); }
if(!is.numeric(IniResLevels)){ stop("IniResLevels must be a vector of numeric values \n"); return(NULL); }
if(length(IniResLevels)!=2 ) { stop("the length of IniStates must be 2 for the chosen FUN_MOD \n"); return(NULL); }
} else {
if("GR" %in% ObjectClass){ IniResLevels <- as.double(c(0.3,0.5)); }
}
##check_Outputs_Cal_and_Sim
##Outputs_all
Outputs_all <- NULL;
if(identical(FUN_MOD,RunModel_GR4J) | identical(FUN_MOD,RunModel_CemaNeigeGR4J)){
Outputs_all <- c(Outputs_all,"PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim"); }
if(identical(FUN_MOD,RunModel_GR5J) | identical(FUN_MOD,RunModel_CemaNeigeGR5J)){
Outputs_all <- c(Outputs_all,"PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim"); }
if(identical(FUN_MOD,RunModel_GR6J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)){
Outputs_all <- c(Outputs_all,"PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QR1","Exp","QD","Qsim"); }
if("CemaNeige" %in% ObjectClass){
Outputs_all <- c(Outputs_all,"Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt"); }
##check_Outputs_Sim
if(!is.vector( Outputs_Sim)){ stop("Outputs_Sim must be a vector of characters \n"); return(NULL); }
if(!is.character(Outputs_Sim)){ stop("Outputs_Sim must be a vector of characters \n"); return(NULL); }
if(sum(is.na(Outputs_Sim))!=0){ stop("Outputs_Sim must not contain NA \n"); return(NULL); }
if("all" %in% Outputs_Sim){ Outputs_Sim <- c("DatesR",Outputs_all,"StateEnd"); }
Test <- which(Outputs_Sim %in% c("DatesR",Outputs_all,"StateEnd") == FALSE); if(length(Test)!=0){
stop(paste("Outputs_Sim is incorrectly defined: ",paste(Outputs_Sim[Test],collapse=", ")," not found \n",sep="")); return(NULL); }
Outputs_Sim <- Outputs_Sim[!duplicated(Outputs_Sim)];
##check_Outputs_Cal
if(is.null(Outputs_Cal)){
if("GR" %in% ObjectClass ){ Outputs_Cal <- c("Qsim"); }
if("CemaNeige" %in% ObjectClass ){ Outputs_Cal <- c("all"); }
if("GR" %in% ObjectClass & "CemaNeige" %in% ObjectClass){ Outputs_Cal <- c("PliqAndMelt","Qsim"); }
} else {
if(!is.vector( Outputs_Cal)){ stop("Outputs_Cal must be a vector of characters \n"); return(NULL); }
if(!is.character(Outputs_Cal)){ stop("Outputs_Cal must be a vector of characters \n"); return(NULL); }
if(sum(is.na(Outputs_Cal))!=0){ stop("Outputs_Cal must not contain NA \n"); return(NULL); }
}
if("all" %in% Outputs_Cal){ Outputs_Cal <- c("DatesR",Outputs_all,"StateEnd"); }
Test <- which(Outputs_Cal %in% c("DatesR",Outputs_all,"StateEnd") == FALSE); if(length(Test)!=0){
stop(paste("Outputs_Cal is incorrectly defined: ",paste(Outputs_Cal[Test],collapse=", ")," not found \n",sep="")); return(NULL); }
Outputs_Cal <- Outputs_Cal[!duplicated(Outputs_Cal)];
##check_RunSnowModule
if("CemaNeige" %in% ObjectClass){
if(!is.vector( RunSnowModule)){ stop("RunSnowModule must be a single boolean \n"); return(NULL); }
if(!is.logical(RunSnowModule)){ stop("RunSnowModule must be either TRUE or FALSE \n"); return(NULL); }
if(length(RunSnowModule)!=1 ){ stop("RunSnowModule must be either TRUE or FALSE \n"); return(NULL); }
}
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261
##check_MeanAnSolidPrecip
if("CemaNeige" %in% ObjectClass & is.null(MeanAnSolidPrecip)){
NLayers <- length(InputsModel$LayerPrecip);
SolidPrecip <- NULL; for(iLayer in 1:NLayers){
if(iLayer==1){ SolidPrecip <- InputsModel$LayerFracSolidPrecip[[1]]*InputsModel$LayerPrecip[[iLayer]]/NLayers;
} else { SolidPrecip <- SolidPrecip + InputsModel$LayerFracSolidPrecip[[iLayer]]*InputsModel$LayerPrecip[[iLayer]]/NLayers; } }
Factor <- NULL;
if(inherits(InputsModel,"hourly" )){ Factor <- 365.25*24; }
if(inherits(InputsModel,"daily" )){ Factor <- 365.25; }
if(inherits(InputsModel,"monthly")){ Factor <- 12; }
if(inherits(InputsModel,"yearly" )){ Factor <- 1; }
if(is.null(Factor)){ stop("InputsModel must be of class 'hourly', 'daily', 'monthly' or 'yearly' \n"); return(NULL); }
MeanAnSolidPrecip <- rep(mean(SolidPrecip)*Factor,NLayers); ### default value: same Gseuil for all layers
if(!quiet){ warning("\t MeanAnSolidPrecip not defined -> it was automatically set to c(",paste(round(MeanAnSolidPrecip),collapse=","),") \n"); }
}
if("CemaNeige" %in% ObjectClass & !is.null(MeanAnSolidPrecip)){
if(!is.vector( MeanAnSolidPrecip) ){ stop(paste("MeanAnSolidPrecip must be a vector of numeric values \n",sep="")); return(NULL); }
if(!is.numeric(MeanAnSolidPrecip) ){ stop(paste("MeanAnSolidPrecip must be a vector of numeric values \n",sep="")); return(NULL); }
if(length(MeanAnSolidPrecip)!=NLayers){ stop(paste("MeanAnSolidPrecip must be a numeric vector of length ",NLayers," \n",sep="")); return(NULL); }
}
##check_PliqAndMelt
if(RunSnowModule & "GR" %in% ObjectClass & "CemaNeige" %in% ObjectClass){
if("PliqAndMelt" %in% Outputs_Cal == FALSE & "all" %in% Outputs_Cal == FALSE){
WTxt <- NULL;
WTxt <- paste(WTxt,"\t PliqAndMelt was not defined in Outputs_Cal but is needed to feed the hydrological model with the snow module outputs \n",sep="");
WTxt <- paste(WTxt,"\t -> it was automatically added \n",sep="");
if(!is.null(WTxt) & !quiet){ warning(WTxt); }
Outputs_Cal <- c(Outputs_Cal,"PliqAndMelt"); }
if("PliqAndMelt" %in% Outputs_Sim == FALSE & "all" %in% Outputs_Sim == FALSE){
WTxt <- NULL;
WTxt <- paste(WTxt,"\t PliqAndMelt was not defined in Outputs_Sim but is needed to feed the hydrological model with the snow module outputs \n",sep="");
WTxt <- paste(WTxt,"\t -> it was automatically added \n",sep="");
if(!is.null(WTxt) & !quiet){ warning(WTxt); }
Outputs_Sim <- c(Outputs_Sim,"PliqAndMelt"); }
}
##Create_RunOptions
RunOptions <- list(IndPeriod_WarmUp=IndPeriod_WarmUp,IndPeriod_Run=IndPeriod_Run,IniStates=IniStates,IniResLevels=IniResLevels,
Outputs_Cal=Outputs_Cal,Outputs_Sim=Outputs_Sim);
if("CemaNeige" %in% ObjectClass){
RunOptions <- c(RunOptions,list(RunSnowModule=RunSnowModule,MeanAnSolidPrecip=MeanAnSolidPrecip)); }
class(RunOptions) <- c("RunOptions",ObjectClass);
return(RunOptions);
}