Skip to content
GitLab

Source

Target

Hide whitespace changes
Inline Side-by-side
Showing
with 1795 additions and 571 deletions
+1795 -571
R/RunModel_GR5J.R
View file @ 4bccfee2
RunModel_GR5J <- function(InputsModel,RunOptions,Param){
RunModel_GR5J <- function(InputsModel, RunOptions, Param) {
NParam <- 5;
FortranOutputs <- .FortranOutputs(GR = "GR5J")$GR
.ArgumentsCheckGR(InputsModel, RunOptions, Param)
##Arguments_check
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("'InputsModel' must be of class 'InputsModel'") }
if(inherits(InputsModel,"daily" )==FALSE){ stop("'InputsModel' must be of class 'daily' ") }
if(inherits(InputsModel,"GR" )==FALSE){ stop("'InputsModel' must be of class 'GR' ") }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("'RunOptions' must be of class 'RunOptions' ") }
if(inherits(RunOptions,"GR" )==FALSE){ stop("'RunOptions' must be of class 'GR' ") }
if(!is.vector(Param) | !is.numeric(Param)){ stop("'Param' must be a numeric vector") }
if(sum(!is.na(Param))!=NParam){ stop(paste("'Param' must be a vector of length ",NParam," and contain no NA",sep="")) }
Param <- as.double(Param);
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
##Input_data_preparation
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
LInputSeries <- as.integer(length(IndPeriod1))
if("all" %in% RunOptions$Outputs_Sim){ IndOutputs <- as.integer(1:length(FortranOutputs));
} else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim); }
Param <- as.double(Param)
##Output_data_preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim;
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
##Use_of_IniResLevels
if(!is.null(RunOptions$IniResLevels)){
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1]; ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3]; ### routing store level (mm)
}
Param_X1X3_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[1L] <- Param_X1X3_threshold
}
if (Param[3L] < Param_X1X3_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3_threshold, Param_X1X3_threshold))
Param[3L] <- Param_X1X3_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
##Call_fortan
RESULTS <- .Fortran("frun_GR5J",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
InputsPrecip=InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/d]
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam=as.integer(length(Param)), ### number of model parameter
Param=Param, ### parameter set
NStates=as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart=RunOptions$IniStates, ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputs)), ### number of output series
IndOutputs=IndOutputs, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputs)), ### output series [mm]
StateEnd=rep(as.double(-999.999),length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
if (ExportStateEnd) {
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5J, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = NULL, UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
##Output_data_preparation
##OutputsModel_only
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(OutputsModel) <- FortranOutputs[IndOutputs]; }
##DatesR_and_OutputsModel_only
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs]); }
##OutputsModel_and_SateEnd_only
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputs[IndOutputs],"StateEnd"); }
##DatesR_and_OutputsModel_and_SateEnd
if((ExportDatesR==TRUE & ExportStateEnd==TRUE) | "all" %in% RunOptions$Outputs_Sim){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs],"StateEnd"); }
## Output data preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
##End
rm(RESULTS);
class(OutputsModel) <- c("OutputsModel","daily","GR");
return(OutputsModel);
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
}
## Call GR model Fortan
RESULTS <- .Fortran("frun_gr5j", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/d]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam = as.integer(length(Param)), ### number of model parameter
Param = Param, ### parameter set
NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart = RunOptions$IniStates, ### state variables used when the model run starts
NOutputs = as.integer(length(IndOutputs)), ### number of output series
IndOutputs = IndOutputs, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm or mm/d]
StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5J, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
UH1 = NULL,
UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param)
}
R/RunModel_GR6J.R
View file @ 4bccfee2
RunModel_GR6J <- function(InputsModel,RunOptions,Param){
RunModel_GR6J <- function(InputsModel, RunOptions, Param) {
NParam <- 6;
FortranOutputs <- .FortranOutputs(GR = "GR6J")$GR
.ArgumentsCheckGR(InputsModel, RunOptions, Param)
##Arguments_check
if(inherits(InputsModel,"InputsModel")==FALSE){ stop("'InputsModel' must be of class 'InputsModel'") }
if(inherits(InputsModel,"daily" )==FALSE){ stop("'InputsModel' must be of class 'daily' ") }
if(inherits(InputsModel,"GR" )==FALSE){ stop("'InputsModel' must be of class 'GR' ") }
if(inherits(RunOptions,"RunOptions" )==FALSE){ stop("'RunOptions' must be of class 'RunOptions' ") }
if(inherits(RunOptions,"GR" )==FALSE){ stop("'RunOptions' must be of class 'GR' ") }
if(!is.vector(Param) | !is.numeric(Param)){ stop("'Param' must be a numeric vector") }
if(sum(!is.na(Param))!=NParam){ stop(paste("'Param' must be a vector of length ",NParam," and contain no NA",sep="")) }
Param <- as.double(Param);
Param_X1X3X6_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[1L] <- Param_X1X3X6_threshold
}
if (Param[3L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[3L] <- Param_X1X3X6_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
if (Param[6L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[6] (X6: coefficient for emptying exponential store [mm]) < %.2f\n X6 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[6L] <- Param_X1X3X6_threshold
}
Param <- as.double(Param)
##Input_data_preparation
if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; }
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
LInputSeries <- as.integer(length(IndPeriod1))
if("all" %in% RunOptions$Outputs_Sim){ IndOutputs <- as.integer(1:length(FortranOutputs));
} else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim); }
Param_X1X3X6_threshold <- 1e-2
Param_X4_threshold <- 0.5
if (Param[1L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[1] (X1: production store capacity [mm]) < %.2f\n X1 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[1L] <- Param_X1X3X6_threshold
}
if (Param[3L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[3] (X3: routing store capacity [mm]) < %.2f\n X3 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[3L] <- Param_X1X3X6_threshold
}
if (Param[4L] < Param_X4_threshold) {
warning(sprintf("Param[4] (X4: unit hydrograph time constant [d]) < %.2f\n X4 set to %.2f", Param_X4_threshold, Param_X4_threshold))
Param[4L] <- Param_X4_threshold
}
if (Param[6L] < Param_X1X3X6_threshold) {
warning(sprintf("Param[6] (X6: coefficient for emptying exponential store [mm]) < %.2f\n X6 set to %.2f", Param_X1X3X6_threshold, Param_X1X3X6_threshold))
Param[6L] <- Param_X1X3X6_threshold
}
##Output_data_preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim;
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim;
##Use_of_IniResLevels
if(!is.null(RunOptions$IniResLevels)){
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1]; ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3]; ### routing store level (mm)
RunOptions$IniStates[3] <- RunOptions$IniResLevels[3] ### exponential store level (mm)
}
## Input data preparation
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
if ("all" %in% RunOptions$Outputs_Sim) {
IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
} else {
IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
}
##Call_fortan
RESULTS <- .Fortran("frun_GR6J",PACKAGE="airGR",
##inputs
LInputs=LInputSeries, ### length of input and output series
InputsPrecip=InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/d]
InputsPE=InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam=as.integer(length(Param)), ### number of model parameter
Param=Param, ### parameter set
NStates=as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart=RunOptions$IniStates, ### state variables used when the model run starts
NOutputs=as.integer(length(IndOutputs)), ### number of output series
IndOutputs=IndOutputs, ### indices of output series
##outputs
Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputs)), ### output series [mm]
StateEnd=rep(as.double(-999.999),length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
if (ExportStateEnd) {
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR6J, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = RESULTS$StateEnd[3L],
UH1 = RESULTS$StateEnd[(1:20)+7], UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## Output data preparation
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
ExportDatesR <- "DatesR" %in% RunOptions$Outputs_Sim
ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
##Output_data_preparation
##OutputsModel_only
if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
names(OutputsModel) <- FortranOutputs[IndOutputs]; }
##DatesR_and_OutputsModel_only
if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs]); }
##OutputsModel_and_SateEnd_only
if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c(FortranOutputs[IndOutputs],"StateEnd"); }
##DatesR_and_OutputsModel_and_SateEnd
if((ExportDatesR==TRUE & ExportStateEnd==TRUE) | "all" %in% RunOptions$Outputs_Sim){
OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
list(RESULTS$StateEnd) );
names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs],"StateEnd"); }
## Use of IniResLevels
if (!is.null(RunOptions$IniResLevels)) {
RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * Param[3] ### routing store level (mm)
RunOptions$IniStates[3] <- RunOptions$IniResLevels[3] ### exponential store level (mm)
}
##End
rm(RESULTS);
class(OutputsModel) <- c("OutputsModel","daily","GR");
return(OutputsModel);
## Call GR model Fortan
RESULTS <- .Fortran("frun_gr6j", PACKAGE = "airGR",
## inputs
LInputs = LInputSeries, ### length of input and output series
InputsPrecip = InputsModel$Precip[IndPeriod1], ### input series of total precipitation [mm/d]
InputsPE = InputsModel$PotEvap[IndPeriod1], ### input series potential evapotranspiration [mm/d]
NParam = as.integer(length(Param)), ### number of model parameter
Param = Param, ### parameter set
NStates = as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
StateStart = RunOptions$IniStates, ### state variables used when the model run starts
NOutputs = as.integer(length(IndOutputs)), ### number of output series
IndOutputs = IndOutputs, ### indices of output series
## outputs
Outputs = matrix(as.double(-99e9), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm or mm/d]
StateEnd = rep(as.double(-99e9), length(RunOptions$IniStates)) ### state variables at the end of the model run
)
RESULTS$Outputs[RESULTS$Outputs <= -99e8] <- NA
RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
if (ExportStateEnd) {
RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR6J, InputsModel = InputsModel,
ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = RESULTS$StateEnd[3L],
UH1 = RESULTS$StateEnd[(1:20) + 7],
UH2 = RESULTS$StateEnd[(1:40) + (7+20)],
GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
verbose = FALSE)
}
## OutputsModel generation
.GetOutputsModelGR(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param)
}
\ No newline at end of file
R/RunModel_Lag.R 0 → 100644
View file @ 4bccfee2
RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
NParam <- 1
## argument check
if (!inherits(InputsModel, "InputsModel")) {
stop("'InputsModel' must be of class 'InputsModel'")
}
if (!inherits(InputsModel, "SD")) {
stop("'InputsModel' must be of class 'SD'")
}
if (!inherits(RunOptions, "RunOptions")) {
stop("'RunOptions' must be of class 'RunOptions'")
}
if (!is.vector(Param) | !is.numeric(Param)) {
stop("'Param' must be a numeric vector")
}
if (sum(!is.na(Param)) != NParam) {
stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
}
if (inherits(QcontribDown, "OutputsModel")) {
if (is.null(QcontribDown$Qsim)) {
stop("'QcontribDown' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment")
}
if (length(QcontribDown$Qsim) != length(RunOptions$IndPeriod_Run)) {
stop("Time series Qsim in 'QcontribDown' should have the same length as 'RunOptions$IndPeriod_Run'")
}
if (!identical(RunOptions$IndPeriod_WarmUp, 0L) && !identical(RunOptions$Outputs_Sim, RunOptions$Outputs_Cal)) {
# This test is not necessary during calibration but usefull in other cases because
# WarmUpQsim is then used for downstream sub-basins because of the delay in Qupstream
if (is.null(QcontribDown$RunOptions$WarmUpQsim) ||
length(QcontribDown$RunOptions$WarmUpQsim) != length(RunOptions$IndPeriod_WarmUp)) {
stop("Time series WarmUpQsim in 'QcontribDown' should have the same length as 'RunOptions$IndPeriod_WarmUp'")
}
}
} else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
if (length(QcontribDown) != length(RunOptions$IndPeriod_Run)) {
stop("'QcontribDown' should have the same length as 'RunOptions$IndPeriod_Run'")
}
} else {
stop("'QcontribDown' must be a numeric vector or a 'OutputsModel' object")
}
# data set up
NbUpBasins <- length(InputsModel$LengthHydro)
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
if (inherits(QcontribDown, "OutputsModel")) {
OutputsModel <- QcontribDown
if (is.null(OutputsModel$RunOptions$WarmUpQsim)) {
OutputsModel$RunOptions$WarmUpQsim <- rep(NA, length(RunOptions$IndPeriod_WarmUp))
}
QsimDown <- c(OutputsModel$RunOptions$WarmUpQsim, OutputsModel$Qsim)
} else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
OutputsModel <- list()
class(OutputsModel) <- c("OutputsModel", class(RunOptions)[-1])
QsimDown <- c(rep(NA, length(RunOptions$IndPeriod_WarmUp)),
QcontribDown)
}
## propagation time from upstream meshes to outlet
PT <- InputsModel$LengthHydro * 1e3 / Param[1L] / RunOptions$FeatFUN_MOD$TimeStep
HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))
## set up initial states
IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
if (length(IniSD) > 0) {
if (sum(floor(PT)) + NbUpBasins != length(IniSD)) {
stop(
sprintf(
"SD initial states has a length of %i and a length of %i is required",
length(IniSD),
sum(floor(PT)) + NbUpBasins
)
)
}
IniStates <- lapply(seq_len(NbUpBasins), function(x) {
iStart <- 1
if (x > 1) {
iStart <- iStart + sum(floor(PT[1:x - 1]) + 1)
}
as.vector(IniSD[iStart:(iStart + PT[x])])
})
} else {
IniStates <- lapply(
seq_len(NbUpBasins),
function(iUpBasins) {
iWarmUp <- seq.int(
from = max(1, IndPeriod1[1] - floor(PT[iUpBasins]) - 1),
to = max(1, IndPeriod1[1] - 1)
)
ini <- InputsModel$Qupstream[iWarmUp, iUpBasins]
if (length(ini) != floor(PT[iUpBasins] + 1)) {
# If warm-up period is not enough long complete beginning with first value
ini <- c(rep(ini[1], floor(PT[iUpBasins] + 1) - length(ini)), ini)
}
return(as.vector(ini))
}
)
}
# message("IniStates: ", paste(IniStates, collapse = ", "))
## Lag model computation
Qsim_m3 <- QsimDown *
InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
for (upstream_basin in seq_len(NbUpBasins)) {
Qupstream <- c(IniStates[[upstream_basin]],
InputsModel$Qupstream[IndPeriod1, upstream_basin])
# message("Qupstream[", upstream_basin, "]: ", paste(Qupstream, collapse = ", "))
Qsim_m3 <- Qsim_m3 +
Qupstream[2:(1 + LInputSeries)] * HUTRANS[1, upstream_basin] +
Qupstream[1:LInputSeries] * HUTRANS[2, upstream_basin]
}
## OutputsModel
if ("Qsim_m3" %in% RunOptions$Outputs_Sim) {
OutputsModel$Qsim_m3 <- Qsim_m3[IndPeriod2]
}
if ("Qsim" %in% RunOptions$Outputs_Sim) {
# Convert back Qsim to mm
OutputsModel$Qsim <- Qsim_m3[IndPeriod2] / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
# message("Qsim: ", paste(OutputsModel$Qsim, collapse = ", "))
}
if ("QsimDown" %in% RunOptions$Outputs_Sim) {
# Convert back Qsim to mm
OutputsModel$QsimDown <- QsimDown[IndPeriod2]
}
# Warning for negative flows or NAs only in extended outputs
if (length(RunOptions$Outputs_Sim) > 2) {
if (any(OutputsModel$Qsim[!is.na(OutputsModel$Qsim)] < 0)) {
warning(length(which(OutputsModel$Qsim < 0)), " time steps with negative flow, set to zero.")
OutputsModel$Qsim[OutputsModel$Qsim < 0] <- 0
}
# Warning for NAs
if (any(is.na(OutputsModel$Qsim))) {
warning(length(which(is.na(OutputsModel$Qsim))), " time steps with NA values")
}
}
if ("StateEnd" %in% RunOptions$Outputs_Sim) {
SD <- lapply(seq(NbUpBasins), function(x) {
lastTS <- RunOptions$IndPeriod_Run[length(RunOptions$IndPeriod_Run)]
InputsModel$Qupstream[(lastTS - floor(PT[x])):lastTS, x]
})
if (is.null(OutputsModel$StateEnd)) {
OutputsModel$StateEnd <- CreateIniStates(RunModel_Lag, InputsModel, SD = SD)
} else {
OutputsModel$StateEnd$SD <- SD
}
# message("StateEnd: ", paste(OutputsModel$StateEnd$SD, collapse = ", "))
}
if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) {
OutputsModel$RunOptions$WarmUpQsim <- Qsim_m3[seq_len(length(RunOptions$IndPeriod_WarmUp))] / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
}
if ("Param" %in% RunOptions$Outputs_Sim) {
OutputsModel$RunOptions$Param <- c(Param, OutputsModel$RunOptions$Param)
}
class(OutputsModel) <- c(class(OutputsModel), "SD")
return(OutputsModel)
}
R/SeriesAggreg.InputsModel.R 0 → 100644
View file @ 4bccfee2
SeriesAggreg.InputsModel <- function(x, Format, ...) {
SeriesAggreg.list(x,
Format,
ConvertFun = NA,
except = c("ZLayers", "LengthHydro", "BasinAreas"),
...)
}
R/SeriesAggreg.OutputsModel.R 0 → 100644
View file @ 4bccfee2
SeriesAggreg.OutputsModel <- function(x, Format, ...) {
SeriesAggreg.list(x,
Format,
ConvertFun = NA,
except = c("RunOptions", "StateEnd"),
...)
}
R/SeriesAggreg.R
View file @ 4bccfee2
SeriesAggreg <- function(TabSeries,
TimeFormat, NewTimeFormat,
ConvertFun,
YearFirstMonth = 1, TimeLag = 0,
verbose = TRUE) {
##_Arguments_check
##check_TabSeries
if (is.null(TabSeries) ) {
stop("'TabSeries' must be a data.frame containing the dates and data to be converted")
}
if (!is.data.frame(TabSeries)) {
stop("'TabSeries' must be a data.frame containing the dates and data to be converted")
}
if (ncol(TabSeries) < 2) {
stop("'TabSeries' must contain at least two columns (including the coulmn of dates")
}
##check_TimeFormat
if (!any(class(TabSeries[, 1]) %in% "POSIXt")) {
stop("'TabSeries' first column must be a vector of class 'POSIXlt' or 'POSIXct'")
}
if (any(class(TabSeries[, 1]) %in% "POSIXlt")) {
TabSeries[, 1] <- as.POSIXct(TabSeries[, 1])
}
for (iCol in 2:ncol(TabSeries)) {
if (!is.numeric(TabSeries[,iCol])) {
stop("'TabSeries' columns (other than the first one) be of numeric class")
}
}
if (is.null(TimeFormat)) {
stop("'TimeFormat' must be one of 'hourly', 'daily', 'monthly' or 'yearly'")
}
if (!is.vector(TimeFormat)) {
stop("'TimeFormat' must be one of 'hourly', 'daily', 'monthly' or 'yearly'")
}
if (!is.character(TimeFormat)) {
stop("'TimeFormat' must be one of 'hourly', 'daily', 'monthly' or 'yearly'")
}
if (length(TimeFormat) != 1) {
stop("'TimeFormat' must be one of 'hourly', 'daily', 'monthly' or 'yearly'")
}
if (! TimeFormat %in% c("hourly", "daily", "monthly", "yearly")) {
stop("'TimeFormat' must be one of 'hourly', 'daily', 'monthly' or 'yearly'")
}
##check_NewTimeFormat
if (is.null(NewTimeFormat)) {
stop("'NewTimeFormat' must be one of 'hourly', 'daily', 'monthly' or 'yearly'")
}
if (!is.vector(NewTimeFormat)) {
stop("'NewTimeFormat' must be one of 'hourly', 'daily', 'monthly' or 'yearly'")
}
if (!is.character(NewTimeFormat)) {
stop("'NewTimeFormat' must be one of 'hourly', 'daily', 'monthly' or 'yearly'")
}
if (length(NewTimeFormat) != 1) {
stop("'NewTimeFormat' must be one of 'hourly', 'daily', 'monthly' or 'yearly'")
}
if (! NewTimeFormat %in% c("hourly", "daily", "monthly", "yearly")) {
stop("'NewTimeFormat' must be one of 'hourly', 'daily', 'monthly' or 'yearly'")
}
##check_ConvertFun
if (is.null(ConvertFun)) {
stop("'ConvertFun' must be a vector of character")
}
if (!is.vector(ConvertFun)) {
stop("'ConvertFun' must be a vector of character")
}
if (!is.character(ConvertFun)) {
stop("'ConvertFun' must be a vector of character")
}
if (length(ConvertFun) != (ncol(TabSeries) - 1)) {
stop(
paste("'ConvertFun' must be of length", ncol(TabSeries) - 1, "(length=ncol(TabSeries)-1)")
)
}
if (sum(ConvertFun %in% c("sum", "mean") == FALSE) != 0) {
stop("'ConvertFun' elements must be one of 'sum' or 'mean'")
}
##check_YearFirstMonth
if (is.null(YearFirstMonth)) {
stop("'YearFirstMonth' must be an integer between 1 and 12")
}
if (!is.vector(YearFirstMonth)) {
stop("'YearFirstMonth' must be an integer between 1 and 12")
}
if (!is.numeric(YearFirstMonth)) {
stop("'YearFirstMonth' must be an integer between 1 and 12")
}
YearFirstMonth <- as.integer(YearFirstMonth)
if (length(YearFirstMonth) != 1) {
stop("'YearFirstMonth' must be only one integer between 1 and 12")
}
if (YearFirstMonth %in% (1:12) == FALSE) {
stop("'YearFirstMonth' must be only one integer between 1 and 12")
}
##check_DatesR_integrity
if (TimeFormat == "hourly") {
by <- "hours"
}
if (TimeFormat == "daily") {
by <- "days"
}
if (TimeFormat == "monthly") {
by <- "months"
}
if (TimeFormat == "yearly") {
by <- "years"
}
TmpDatesR <- seq(from = TabSeries[1, 1], to = tail(TabSeries[, 1], 1), by = by)
if (!identical(TabSeries[, 1], TmpDatesR)) {
stop("some dates might not be ordered or are missing in 'TabSeries'")
}
##check_conversion_direction
if ((TimeFormat == "daily" & NewTimeFormat %in% c("hourly") ) |
(TimeFormat == "monthly" & NewTimeFormat %in% c("hourly","daily") ) |
(TimeFormat == "yearly" & NewTimeFormat %in% c("hourly","daily","monthly"))) {
stop("only time aggregation can be performed")
}
##check_if_conversion_not_needed
if ((TimeFormat == "hourly" & NewTimeFormat == "hourly" ) |
(TimeFormat == "daily" & NewTimeFormat == "daily" ) |
(TimeFormat == "monthly" & NewTimeFormat == "monthly") |
(TimeFormat == "yearly" & NewTimeFormat == "yearly" )) {
if (verbose) {
warning("the old and new format are identical \n\t -> no time-step conversion was performed")
return(TabSeries)
}
}
##_Time_step_conversion
##_Handle_conventional_difference_between_hourly_series_and_others
TmpDatesR <- TabSeries[, 1]
#if (TimeFormat=="hourly") { TmpDatesR <- TmpDatesR - 60*60; }
TmpDatesR <- TmpDatesR + TimeLag
Hmax <- "00"
if (TimeFormat == "hourly") {
Hmax <- "23"
}
##_Identify_the_part_of_the_series_to_be_aggregated
NDaysInMonth <- list("31", c("28", "29"), "31", "30", "31", "30", "31", "31", "30", "31", "30", "31")
YearLastMonth <- YearFirstMonth + 11
if (YearLastMonth > 12) {
YearLastMonth <- YearLastMonth - 12
}
YearFirstMonthTxt <- formatC(YearFirstMonth, format = "d", width = 2, flag = "0")
YearLastMonthTxt <- formatC(YearLastMonth , format = "d", width = 2, flag = "0")
if (NewTimeFormat == "daily") {
Ind1 <- which(format(TmpDatesR, "%H") == "00")
Ind2 <- which(format(TmpDatesR, "%H") == Hmax)
if (Ind2[1] < Ind1[1]) {
Ind2 <- Ind2[2:length(Ind2)]
}
if (tail(Ind1, 1) > tail(Ind2, 1)) {
Ind1 <- Ind1[1:(length(Ind1) - 1)]
}
### Aggr <- NULL; iii <- 0; for(kkk in 1:length(Ind1)) {
### iii <- iii+1; Aggr <- c(Aggr,rep(iii,length(Ind1[kkk]:Ind2[kkk]))); }
Aggr <- as.numeric(format(TmpDatesR[min(Ind1):max(Ind2)], "%Y%m%d"))
### more efficient
NewDatesR <- data.frame(seq(from = TmpDatesR[min(Ind1)], to = TmpDatesR[max(Ind2)], by = "days"))
}
if (NewTimeFormat=="monthly") {
Ind1 <- which(format(TmpDatesR, "%d%H") == "0100")
Ind2 <- which(format(TmpDatesR,"%m%d%H") %in% paste0(c("0131", "0228", "0229", "0331", "0430", "0531", "0630", "0731", "0831", "0930", "1031", "1130", "1231"), Hmax))
Ind2[1:(length(Ind2) - 1)][diff(Ind2) == 1] <- NA
Ind2 <- Ind2[!is.na(Ind2)] ### to keep only feb 29 if both feb 28 and feb 29 exists
if (Ind2[1] < Ind1[1]) {
Ind2 <- Ind2[2:length(Ind2)]
}
if (tail(Ind1, 1) > tail(Ind2, 1)) {
Ind1 <- Ind1[1:(length(Ind1) - 1)]
}
### Aggr <- NULL; iii <- 0; for(kkk in 1:length(Ind1)) {
### iii <- iii+1; Aggr <- c(Aggr,rep(iii,length(Ind1[kkk]:Ind2[kkk]))); }
Aggr <- as.numeric(format(TmpDatesR[min(Ind1):max(Ind2)],"%Y%m")); ### more efficient
NewDatesR <- data.frame(seq(from=TmpDatesR[min(Ind1)],to=TmpDatesR[max(Ind2)],by="months"))
}
if (NewTimeFormat == "yearly") {
Ind1 <- which(format(TmpDatesR, "%m%d%H") %in% paste0(YearFirstMonthTxt, "0100"))
Ind2 <- which(format(TmpDatesR, "%m%d%H") %in% paste0(YearLastMonthTxt, NDaysInMonth[[YearLastMonth]], Hmax))
Ind2[1:(length(Ind2) - 1)][diff(Ind2) == 1] <- NA
Ind2 <- Ind2[!is.na(Ind2)]
### to keep only feb 29 if both feb 28 and feb 29 exists
if (Ind2[1] < Ind1[1]) {
Ind2 <- Ind2[2:length(Ind2)]
}
if (tail(Ind1, 1) > tail(Ind2, 1)) {
Ind1 <- Ind1[1:(length(Ind1) - 1)]
}
Aggr <- NULL
iii <- 0
for (kkk in 1:length(Ind1)) {
iii <- iii + 1
Aggr <- c(Aggr, rep(iii, length(Ind1[kkk]:Ind2[kkk])))
}
### Aggr <- as.numeric(format(TmpDatesR[min(Ind1):max(Ind2)],"%Y")); ### not working if YearFirstMonth != 01
NewDatesR <- data.frame(seq(from = TmpDatesR[min(Ind1)], to = TmpDatesR[max(Ind2)], by = "years"))
}
##_Aggreation_and_export
NewTabSeries <- data.frame(NewDatesR)
for (iCol in 2:ncol(TabSeries)) {
AggregData <- aggregate(TabSeries[min(Ind1):max(Ind2), iCol],
by = list(Aggr),
FUN = ConvertFun[iCol - 1], na.rm = FALSE)[, 2]
NewTabSeries <- data.frame(NewTabSeries, AggregData)
}
names(NewTabSeries) <- names(TabSeries)
return(NewTabSeries)
}
\ No newline at end of file
SeriesAggreg <- function(x, Format, ...) {
UseMethod("SeriesAggreg")
}
R/SeriesAggreg.data.frame.R 0 → 100644
View file @ 4bccfee2
SeriesAggreg.data.frame <- function(x,
Format,
ConvertFun,
TimeFormat = NULL, # deprecated
NewTimeFormat = NULL, # deprecated
YearFirstMonth = 1,
TimeLag = 0,
...) {
## Arguments checks
if (!is.null(TimeFormat)) {
warning("deprecated 'TimeFormat' argument", call. = FALSE)
}
if (missing(Format)) {
Format <- .GetSeriesAggregFormat(NewTimeFormat)
} else if (!is.null(NewTimeFormat)) {
warning("deprecated 'NewTimeFormat' argument: 'Format' argument is used instead",
call. = FALSE)
}
if (is.null(Format)) {
stop("argument 'Format' is missing")
}
## check x
if (!is.data.frame(x)) {
stop("'x' must be a data.frame containing the dates and data to be aggregated")
}
if (ncol(x) < 2) {
stop("'x' must contain at least two columns (including the column of dates)")
}
## check x date column
if (!inherits(x[[1L]], "POSIXt")) {
stop("'x' first column must be a vector of class 'POSIXlt' or 'POSIXct'")
}
if (inherits(x[[1L]], "POSIXlt")) {
x[[1L]] <- as.POSIXct(x[[1L]])
}
## check x other columns (boolean converted to numeric)
apply(x[, -1L, drop = FALSE],
MARGIN = 2,
FUN = function(iCol) {
if (!is.numeric(iCol)) {
stop("'x' columns (other than the first one) must be of numeric class")
}
})
## check Format
listFormat <- c("%Y%m%d", "%Y%m", "%Y", "%m", "%d")
Format <- gsub(pattern = "[[:punct:]]+", replacement = "%", Format)
Format <- match.arg(Format, choices = listFormat)
## check ConvertFun
if (length(ConvertFun) != (ncol(x) - 1)) {
stop(sprintf("'ConvertFun' must be of length %i (ncol(x)-1)", ncol(x) - 1))
}
listConvertFun <- lapply(unique(ConvertFun), function(y) {
if (!grepl("^q\\d+$", y, ignore.case = TRUE)) {
return(match.fun(y))
}
})
names(listConvertFun) <- unique(ConvertFun)
lapply(ConvertFun, function(y) {
if (!grepl("^q\\d+$", y, ignore.case = TRUE)) {
TestOutput <- listConvertFun[[y]](1:10)
if (!is.numeric(TestOutput)) {
stop(sprintf("Returned value of '%s' function should be numeric", y))
}
if (length(TestOutput) != 1) {
stop(sprintf("Returned value of '%s' function should be of length 1", y))
}
}
})
## check YearFirstMonth
msgYearFirstMonth <- "'YearFirstMonth' should be a single vector of numeric value between 1 and 12"
YearFirstMonth <- match(YearFirstMonth, 1:12)
if (anyNA(YearFirstMonth)) {
stop(msgYearFirstMonth)
}
if (length(YearFirstMonth) != 1) {
stop(msgYearFirstMonth)
}
if (YearFirstMonth != 1 & Format != "%Y") {
warning("'YearFirstMonth' is ignored because Format != '%Y'")
}
## check TimeLag
msgTimeLag <- "'TimeLag' should be a single vector of a positive numeric value"
if (!is.vector(TimeLag)) {
stop(msgTimeLag)
}
if (!is.numeric(TimeLag)) {
stop(msgTimeLag)
}
if (length(TimeLag) != 1 | !any(TimeLag >= 0)) {
stop(msgTimeLag)
}
TabSeries0 <- x
colnames(TabSeries0)[1L] <- "DatesR"
TabSeries0$DatesR <- TabSeries0$DatesR + TimeLag
TabSeries2 <- TabSeries0
if (!Format %in% c("%d", "%m")) {
start <- sprintf("%i-01-01 00:00:00",
as.numeric(format(TabSeries2$DatesR[1L], format = "%Y")) - 1)
stop <- sprintf("%i-12-31 00:00:00",
as.numeric(format(TabSeries2$DatesR[nrow(TabSeries2)], format = "%Y")) + 1)
unitTs <- format(diff(x[1:2, 1]))
if (gsub("[0-9]+ ", "", unitTs) == "hours") {
byTs <- "hours"
} else {
if (gsub(" days$", "", unitTs) == "1") {
byTs <- "days"
} else {
byTs <- "months"
}
}
fakeTs <- data.frame(DatesR = seq(from = as.POSIXct(start, tz = "UTC"),
to = as.POSIXct(stop , tz = "UTC"),
by = byTs) + TimeLag)
TabSeries2 <- merge(fakeTs, TabSeries2, by = "DatesR", all.x = TRUE)
}
TabSeries2$DatesRini <- TabSeries2$DatesR %in% TabSeries0$DatesR
TabSeries2$Selec2 <- format(TabSeries2$DatesR, Format)
if (nchar(Format) > 2 | Format == "%Y") {
# Compute aggregation
TabSeries2$Selec <- !duplicated(TabSeries2$Selec2)
if (all(TabSeries2$Selec)) {
warning("the requested time 'Format' is the same as the one in 'x'. No time-step conversion was performed")
return(x)
}
if (Format == "%Y") {
yfm <- sprintf("%02.f", YearFirstMonth)
spF1 <- "%m"
spF2 <- "%Y-%m"
TabSeries2$Selec1 <- format(TabSeries2$DatesR, spF1)
TabSeries2$Selec2 <- format(TabSeries2$DatesR, spF2)
TabSeries2$Selec <- !duplicated(TabSeries2$Selec2) & TabSeries2$Selec1 == yfm
}
TabSeries2$Fac2 <- cumsum(TabSeries2$Selec)
} else {
# Compute regime
if (Format == "%d") {
spF2 <- "%m-%d"
TabSeries2$Selec2 <- format(TabSeries2$DatesR, spF2)
}
TabSeries2$Fac2 <- TabSeries2$Selec2
TabSeries2$Selec <- !duplicated(TabSeries2$Selec2)
}
listTsAggreg <- lapply(names(listConvertFun), function(y) {
if (any(ConvertFun == y)) {
colTsAggreg <- c("Fac2", colnames(x)[-1L][ConvertFun == y])
if (grepl("^q\\d+$", y, ignore.case = TRUE)) {
probs <- as.numeric(gsub("^q", "", y, ignore.case = TRUE)) / 100
if (probs < 0 || probs > 1) {
stop("'Q...' format of argument 'ConvertFun' must be an integer between 0 and 100")
}
aggregate(. ~ Fac2,
data = TabSeries2[, colTsAggreg],
FUN = quantile,
na.action = na.pass,
probs = probs,
type = 8,
na.rm = TRUE)
} else {
aggregate(. ~ Fac2,
data = TabSeries2[, colTsAggreg],
FUN = listConvertFun[[y]],
na.action = na.pass)
}
} else {
NULL
}
})
listTsAggreg <- listTsAggreg[!sapply(listTsAggreg, is.null)]
tsAggreg <- do.call(cbind, listTsAggreg)
tsAggreg <- tsAggreg[, !duplicated(colnames(tsAggreg))]
tsAggreg <- merge(tsAggreg,
TabSeries2[, c("Fac2", "DatesR", "DatesRini", "Selec")],
by = "Fac2",
all.x = TRUE,
all.y = FALSE)
tsAggreg <- tsAggreg[tsAggreg$Selec & tsAggreg$DatesRini, ]
tsAggreg <- tsAggreg[, colnames(TabSeries0)]
tsAggreg <- tsAggreg[order(tsAggreg$DatesR), ] # reorder by date especially for regime time series
colnames(tsAggreg)[1L] <- colnames(x)[1L] # keep original column names
return(tsAggreg)
}
R/SeriesAggreg.list.R 0 → 100644
View file @ 4bccfee2
SeriesAggreg.list <- function(x,
Format,
ConvertFun,
NewTimeFormat = NULL,
simplify = FALSE,
except = NULL,
recursive = TRUE,
...) {
classIni <- class(x)
class(x) <- "list" # in order to avoid the use of '['.InputsModel' or '['.OutputsModel' when x[i] is used
if (missing(Format)) {
Format <- .GetSeriesAggregFormat(NewTimeFormat)
} else if (!is.null(NewTimeFormat)) {
warning("deprecated 'NewTimeFormat' argument: 'Format' argument is used instead",
call. = FALSE)
}
# Check ConvertFun
if (any(classIni %in% c("InputsModel", "OutputsModel"))) {
if (!all(is.na(ConvertFun))) {
warning("Argument 'ConvertFun' is ignored on 'InputsModel' and 'OutputsModel' objects")
}
} else if (length(ConvertFun) != 1) {
stop("Argument 'ConvertFun' must be of length 1 with 'list' object")
} else if (!is.character(ConvertFun)) {
stop("Argument 'ConvertFun' must be a character")
}
# Determination of DatesR
if (!is.null(x$DatesR)) {
if (!inherits(x$DatesR, "POSIXt")) {
stop("'x$DatesR' should be of class 'POSIXt'")
}
DatesR <- x$DatesR
} else {
# Auto-detection of POSIXt item in Tabseries
itemPOSIXt <- which(sapply(x, function(x) {
inherits(x, "POSIXt")
}, simplify = TRUE))[1]
if (is.na(itemPOSIXt)) {
stop("At least one item of argument 'x' should be of class 'POSIXt'")
}
warning("Item 'DatesR' not found in 'x' argument: the item ",
names(x)[itemPOSIXt],
" has been automatically chosen")
DatesR <- x[[names(x)[itemPOSIXt]]]
}
# Selection of numeric items for aggregation
numericCols <- names(which(sapply(x, inherits, "numeric")))
arrayCols <- names(which(sapply(x, inherits, "array")))
numericCols <- setdiff(numericCols, arrayCols)
if (!is.null(except)) {
if (!inherits(except, "character")) {
stop("Argument 'except' should be a 'character'")
}
numericCols <- setdiff(numericCols, except)
}
cols <- x[numericCols]
lengthCols <- sapply(cols, length, simplify = TRUE)
if (any(lengthCols != length(DatesR))) {
sErr <- paste0(names(lengthCols)[lengthCols != length(DatesR)],
" (", lengthCols[lengthCols != length(DatesR)], ")",
collapse = ", ")
warning("The length of the following `numeric` items in 'x' ",
"is different than the length of 'DatesR (",
length(DatesR),
")': they will be ignored in the aggregation: ",
sErr)
cols <- cols[lengthCols == length(DatesR)]
}
dfOut <- NULL
if (length(cols)) {
# Treating aggregation at root level
if (is.na(ConvertFun)) {
ConvertFun2 <- .GetAggregConvertFun(names(cols), Format)
} else {
ConvertFun2 <- rep(ConvertFun, length(cols))
}
dfOut <- SeriesAggreg(cbind(DatesR, as.data.frame(cols)),
Format,
...,
ConvertFun = ConvertFun2)
}
if (simplify) {
# Returns data.frame of numeric found in the first level of the list
return(dfOut)
} else {
res <- list()
# Convert aggregated data.frame into list
if (!is.null(dfOut)) {
res <- as.list(dfOut)
## To be consistent with InputsModel class and because plot.OutputsModel use the POSIXlt class
res$DatesR <- as.POSIXlt(res$DatesR)
}
# Exploration of embedded lists and data.frames
if (is.null(recursive) || recursive) {
listCols <- x[!names(x) %in% except]
listCols <- listCols[sapply(listCols, inherits, "list")]
dfCols <- x[sapply(x, inherits, "data.frame")]
dfCols <- c(dfCols, x[sapply(x, inherits, "matrix")])
listCols <- listCols[setdiff(names(listCols), names(dfCols))]
if (length(listCols) > 0) {
if (is.na(ConvertFun)) {
# Check for predefined ConvertFun for all sub-elements
listConvertFun <- .GetAggregConvertFun(names(listCols), Format)
}
# Run SeriesAggreg for each embedded list
listRes <- lapply(names(listCols), function(y) {
listCols[[y]]$DatesR <- DatesR
if (is.na(ConvertFun)) {
SeriesAggreg.list(listCols[[y]],
Format = Format,
recursive = NULL,
...,
ConvertFun = listConvertFun[y])
} else {
SeriesAggreg.list(listCols[[y]],
Format = Format,
recursive = NULL,
...)
}
})
names(listRes) <- names(listCols)
if (is.null(res$DatesR)) {
# Copy DatesR in top level list
res$DatesR <- listRes[[1]]$DatesR
}
# Remove DatesR in embedded lists
lapply(names(listRes), function(x) {
listRes[[x]]$DatesR <<- NULL
})
res <- c(res, listRes)
}
if (length(dfCols) > 0) {
# Processing matrix and dataframes
for (i in length(dfCols)) {
key <- names(dfCols)[i]
m <- dfCols[[i]]
if (nrow(m) != length(DatesR)) {
warning(
"The number of rows of the 'matrix' item ",
key, " (", nrow(m),
") is different than the length of 'DatesR ('", length(DatesR),
"), it will be ignored in the aggregation"
)
} else {
if (is.na(ConvertFun)) {
ConvertFun2 <- rep(.GetAggregConvertFun(key, Format), ncol(m))
} else {
ConvertFun2 <- rep(ConvertFun, ncol(m))
}
res[[key]] <- SeriesAggreg.data.frame(data.frame(DatesR, m),
Format = Format,
ConvertFun = ConvertFun2)
}
}
}
}
# Store elements that are not aggregated
res <- c(res, x[setdiff(names(x), names(res))])
class(res) <- gsub("hourly|daily|monthly|yearly",
.GetSeriesAggregClass(Format),
classIni)
return(res)
}
}
R/TransfoParam_CemaNeige.R
View file @ 4bccfee2
......@@ -28,7 +28,7 @@ TransfoParam_CemaNeige <- function(ParamIn, Direction) {
}
if (Direction == "RT") {
ParamOut <- ParamIn
ParamOut[, 1] <- ParamIn[, 1] * 19.98 - 9.99 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
ParamOut[, 1] <- ParamIn[, 1] * 19.98 - 9.99 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
ParamOut[, 2] <- log(ParamIn[, 2] * 200) ### CemaNeige X2 (degree-day melt coefficient)
}
......
R/TransfoParam_CemaNeigeHyst.R
View file @ 4bccfee2
......@@ -22,15 +22,15 @@ TransfoParam_CemaNeigeHyst <- function(ParamIn, Direction) {
## transformation
if (Direction == "TR") {
ParamOut <- ParamIn
ParamOut <- ParamIn
ParamOut[, 1] <- (ParamIn[, 1] + 9.99) / 19.98 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
ParamOut[, 2] <- exp(ParamIn[, 2]) / 200 ### CemaNeige X2 (degree-day melt coefficient)
ParamOut[, 3] <- (ParamIn[, 3] * 5) + 50 ### Hyst Gaccum
ParamOut[, 4] <- (ParamIn[, 4] / 19.98) + 0.5 ### Hyst CV
}
if (Direction == "RT") {
ParamOut <- ParamIn
ParamOut[, 1] <- ParamIn[, 1] * 19.98 - 9.99 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
ParamOut <- ParamIn
ParamOut[, 1] <- ParamIn[, 1] * 19.98 - 9.99 ### CemaNeige X1 (weighting coefficient for snow pack thermal state)
ParamOut[, 2] <- log(ParamIn[, 2] * 200) ### CemaNeige X2 (degree-day melt coefficient)
ParamOut[, 3] <- (ParamIn[, 3] - 50) / 5 ### Hyst Gaccum
ParamOut[, 4] <- (ParamIn[, 4] - 0.5) * 19.98 ### Hyst CV
......
R/TransfoParam_GR1A.R
View file @ 4bccfee2
......@@ -2,7 +2,6 @@ TransfoParam_GR1A <- function(ParamIn, Direction) {
## number of model parameters
NParam <- 1L
NParam <- 2L
## check arguments
......@@ -26,7 +25,7 @@ TransfoParam_GR1A <- function(ParamIn, Direction) {
ParamOut <- (ParamIn + 10.0) / 8
}
if (Direction == "RT") {
ParamOut <- ParamIn * 8 - 10.0
ParamOut <- ParamIn * 8 - 10.0
}
......
R/TransfoParam_GR5H.R 0 → 100644
View file @ 4bccfee2
TransfoParam_GR5H <- function(ParamIn, Direction) {
## number of model parameters
NParam <- 5L
## check arguments
isVecParamIn <- is.vector(ParamIn)
if (isVecParamIn) {
ParamIn <- matrix(ParamIn, nrow = 1)
}
if (!inherits(ParamIn, "matrix")) {
stop("'ParamIn' must be of class 'matrix'")
}
if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
}
if (ncol(ParamIn) != NParam) {
stop(sprintf("the GR5H model requires %i parameters", NParam))
}
## transformation
if (Direction == "TR") {
ParamOut <- ParamIn
ParamOut[, 1] <- exp(ParamIn[, 1]) ### GR5H X1 (production store capacity)
ParamOut[, 2] <- sinh(ParamIn[, 2]) ### GR5H X2 (groundwater exchange coefficient)
ParamOut[, 3] <- exp(ParamIn[, 3]) ### GR5H X3 (routing store capacity)
ParamOut[, 4] <- 480 + (480 - 0.01) * (ParamIn[, 4] - 10) / 20 ### GR5H X4 (unit hydrograph time constant)
ParamOut[, 5] <- (ParamIn[, 5] + 10) / 20 ### GR5H X5 (groundwater exchange coefficient 2)
}
if (Direction == "RT") {
ParamOut <- ParamIn
ParamOut[, 1] <- log(ParamIn[, 1]) ### GR5H X1 (production store capacity)
ParamOut[, 2] <- asinh(ParamIn[, 2]) ### GR5H X2 (groundwater exchange coefficient)
ParamOut[, 3] <- log(ParamIn[, 3]) ### GR5H X3 (routing store capacity)
ParamOut[, 4] <- (ParamIn[, 4] - 480) * 20 / (480 - 0.01) + 10 ### GR5H X4 (unit hydrograph time constant)
ParamOut[, 5] <- ParamIn[, 5] * 20 - 10 ### GR5H X5 (groundwater exchange coefficient 2)
}
## check output
if (isVecParamIn) {
ParamOut <- as.vector(ParamOut)
}
return(ParamOut)
}
R/TransfoParam_GR5J.R
View file @ 4bccfee2
TransfoParam_GR5J <- function(ParamIn, Direction) {
## number of model parameters
NParam <- 5L
## check arguments
isVecParamIn <- is.vector(ParamIn)
if (isVecParamIn) {
ParamIn <- matrix(ParamIn, nrow = 1)
}
}
if (!inherits(ParamIn, "matrix")) {
stop("'ParamIn' must be of class 'matrix'")
}
......@@ -17,10 +17,10 @@ TransfoParam_GR5J <- function(ParamIn, Direction) {
stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
}
if (ncol(ParamIn) != NParam) {
stop(sprintf("the GR4J model requires %i parameters", NParam))
stop(sprintf("the GR5J model requires %i parameters", NParam))
}
## transformation
if (Direction == "TR") {
ParamOut <- ParamIn
......@@ -29,7 +29,7 @@ TransfoParam_GR5J <- function(ParamIn, Direction) {
ParamOut[, 3] <- exp(ParamIn[, 3]) ### GR5J X3 (routing store capacity)
ParamOut[, 4] <- 20 + 19.5 * (ParamIn[, 4] - 9.99) / 19.98 ### GR5J X4 (unit hydrograph time constant)
ParamOut[, 5] <- (ParamIn[, 5] + 9.99) / 19.98 ### GR5J X5 (groundwater exchange coefficient 2)
}
}
if (Direction == "RT") {
ParamOut <- ParamIn
ParamOut[, 1] <- log(ParamIn[, 1]) ### GR5J X1 (production store capacity)
......@@ -37,14 +37,14 @@ TransfoParam_GR5J <- function(ParamIn, Direction) {
ParamOut[, 3] <- log(ParamIn[, 3]) ### GR5J X3 (routing store capacity)
ParamOut[, 4] <- 9.99 + 19.98 * (ParamIn[, 4] - 20) / 19.5 ### GR5J X4 (unit hydrograph time constant)
ParamOut[, 5] <- ParamIn[, 5] * 19.98 - 9.99 ### GR5J X5 (groundwater exchange coefficient 2)
}
}
## check output
if (isVecParamIn) {
ParamOut <- as.vector(ParamOut)
}
return(ParamOut)
}
......
R/TransfoParam_Lag.R 0 → 100644
View file @ 4bccfee2
TransfoParam_Lag <- function(ParamIn, Direction) {
## number of model parameters
NParam <- 1L
## check arguments
isVecParamIn <- is.vector(ParamIn)
if (isVecParamIn) {
ParamIn <- matrix(ParamIn, nrow = 1)
}
if (!inherits(ParamIn, "matrix")) {
stop("'ParamIn' must be of class 'matrix'")
}
if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
}
if (ncol(ParamIn) != NParam) {
stop(sprintf("the Lag model requires %i parameters", NParam))
}
## transformation
if (Direction == "TR") {
ParamOut <- 20 * (ParamIn + 10) / 20.0
}
if (Direction == "RT") {
ParamOut <- ParamIn * 20.0 / 20 - 10
}
## check output
if (isVecParamIn) {
ParamOut <- as.vector(ParamOut)
}
return(ParamOut)
}
R/Utils.R
View file @ 4bccfee2
......@@ -3,7 +3,7 @@
## function to check
## =================================================================================
# .onLoad <- function(libname, pkgname){
# .onLoad <- function(libname, pkgname) {
# if (requireNamespace("airGRteaching", quietly = TRUE)) {
# if (packageVersion("airGRteaching") %in% package_version(c("0.2.0.9", "0.2.2.2", "0.2.3.2"))) {
# packageStartupMessage("In order to be compatible with the present version of 'airGR', please update your version of the 'airGRteaching' package.")
......@@ -12,15 +12,82 @@
# }
## =================================================================================
## function to extract model features
## =================================================================================
## table of feature models
.FeatModels <- function() {
path <- system.file("modelsFeatures/FeatModelsGR.csv", package = "airGR")
read.table(path, header = TRUE, sep = ";", stringsAsFactors = FALSE)
}
## function to extract model features
.GetFeatModel <- function(FUN_MOD, DatesR = NULL) {
FeatMod <- .FeatModels()
NameFunMod <- ifelse(test = FeatMod$Pkg %in% "airGR",
yes = paste("RunModel", FeatMod$NameMod, sep = "_"),
no = FeatMod$NameMod)
FunMod <- lapply(NameFunMod, FUN = match.fun)
IdMod <- which(sapply(FunMod, FUN = function(x) identical(FUN_MOD, x)))
if (length(IdMod) < 1) {
stop("'FUN_MOD' must be one of ", paste(NameFunMod, collapse = ", "))
} else {
res <- as.list(FeatMod[IdMod, ])
res$NameFunMod <- NameFunMod[IdMod]
if (!is.null(DatesR)) {
DiffTimeStep <- as.numeric(difftime(DatesR[length(DatesR)],
DatesR[length(DatesR)-1],
units = "secs"))
if (is.na(res$TimeUnit)) {
if (any(DiffTimeStep %in% 3600:3601)) { # 3601: leap second
res$TimeUnit <- "hourly"
} else {
res$TimeUnit <- "daily"
}
}
}
res$TimeStep <- switch(res$TimeUnit,
hourly = 1,
daily = 1 * 24,
monthly = 28:31 * 24,
yearly = 365:366 * 24)
res$TimeStepMean <- switch(res$TimeUnit,
hourly = 1,
daily = 1 * 24,
monthly = 365.25 / 12 * 24,
yearly = 365.25 * 24)
res$TimeStep <- res$TimeStep * 3600
res$TimeStepMean <- as.integer(res$TimeStepMean * 3600)
res$Class <- c(res$TimeUnit, res$Class)
res$CodeModHydro <- res$CodeMod
if (grepl("CemaNeige", res$NameFunMod)) {
res$Class <- c(res$Class, "CemaNeige")
res$CodeModHydro <- gsub("CemaNeige", "", res$CodeMod)
}
res$Class <- res$Class[!is.na(res$Class)]
if (!is.null(DatesR)) {
if (all(DiffTimeStep != res$TimeStep)) {
stop("the time step of the model inputs must be ", res$TimeUnit)
}
}
return(res)
}
}
## =================================================================================
## function to manage Fortran outputs
## =================================================================================
.FortranOutputs <- function(GR = NULL, isCN = FALSE) {
outGR <- NULL
outCN <- NULL
if (is.null(GR)) {
GR <- ""
}
......@@ -28,26 +95,28 @@
outGR <- c("PotEvap", "Precip",
"Qsim")
} else if (GR == "GR2M") {
outGR <- c("PotEvap", "Precip", "Prod", "Pn",
outGR <- c("PotEvap", "Precip", "Prod", "Pn", "Ps",
"AE",
"Perc", "PR",
"Rout", "Exch",
"Rout",
"AExch",
"Qsim")
} else if (GR == "GR4H") {
outGR <- c("PotEvap", "Precip", "Prod",
"AE",
} else if (GR == "GR5H") {
outGR <- c("PotEvap", "Precip", "Interc", "Prod", "Pn", "Ps",
"AE", "EI", "ES",
"Perc", "PR",
"Q9", "Q1",
"Rout", "Exch",
"AExch1", "AExch2",
"AExch", "QR",
"QD",
"Qsim")
} else if (GR %in% c("GR4J", "GR5J")) {
} else if (GR %in% c("GR4J", "GR5J", "GR4H")) {
outGR <- c("PotEvap", "Precip", "Prod", "Pn", "Ps",
"AE",
"Perc", "PR",
"Q9", "Q1",
"Rout", "Exch",
"Rout", "Exch",
"AExch1", "AExch2",
"AExch", "QR",
"QD",
......@@ -65,163 +134,130 @@
"Qsim")
}
if (isCN) {
outCN <- c("Pliq", "Psol",
"SnowPack", "ThermalState", "Gratio",
"PotMelt", "Melt", "PliqAndMelt", "Temp",
outCN <- c("Pliq", "Psol",
"SnowPack", "ThermalState", "Gratio",
"PotMelt", "Melt", "PliqAndMelt", "Temp",
"Gthreshold", "Glocalmax")
}
res <- list(GR = outGR, CN = outCN)
}
## =================================================================================
## function to manage inputs of specific ErrorCrit_*() functions
## functions to extract parts of InputsModel or OutputsModel objects
## =================================================================================
.ErrorCrit <- function(InputsCrit, Crit, OutputsModel, warnings) {
## Arguments check
if (!inherits(InputsCrit, "InputsCrit")) {
stop("'InputsCrit' must be of class 'InputsCrit'", call. = FALSE)
}
if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) {
if (Crit == "RMSE") {
stop("'InputsCrit' must be of class 'Single'. Use the 'ErrorCrit' function on objects of class 'Multi' with RMSE", call. = FALSE)
} else {
stop(paste0("'InputsCrit' must be of class 'Single'. Use the 'ErrorCrit' function on objects of class 'Multi' or 'Compo' with ", Crit), call. = FALSE)
## InputsModel
.ExtractInputsModel <- function(x, i) {
res <- lapply(x, function(x) {
if (is.matrix(x)) {
res0 <- x[i, , drop = FALSE]
}
}
## Initialisation
CritName <- NA
CritVar <- InputsCrit$VarObs
if (InputsCrit$transfo == "") {
CritName <- paste0(Crit, "[CritVar]")
}
if (InputsCrit$transfo %in% c("sqrt", "log", "sort", "boxcox")) {
CritName <- paste0(Crit, "[", InputsCrit$transfo, "(CritVar)]")
}
if (InputsCrit$transfo == "inv") {
CritName <- paste0(Crit, "[1/CritVar]")
}
if (grepl("\\^", InputsCrit$transfo)) {
transfoPow <- suppressWarnings(as.numeric(gsub("\\^", "", InputsCrit$transfo)))
CritName <- paste0(Crit, "[CritVar^", transfoPow, "]")
}
CritName <- gsub(pattern = "CritVar", replacement = CritVar, x = CritName)
CritValue <- NA
if (Crit %in% c("RMSE")) {
CritBestValue <- +1
Multiplier <- +1
}
if (Crit %in% c("NSE", "KGE", "KGE2")) {
CritBestValue <- +1
Multiplier <- -1
}
## Data preparation
VarObs <- InputsCrit$Obs
VarObs[!InputsCrit$BoolCrit] <- NA
if (InputsCrit$VarObs == "Q") {
VarSim <- OutputsModel$Qsim
}
if (InputsCrit$VarObs == "SCA") {
VarSim <- rowMeans(sapply(OutputsModel$CemaNeigeLayers[InputsCrit$idLayer], FUN = "[[", "Gratio"))
}
if (InputsCrit$VarObs == "SWE") {
VarSim <- rowMeans(sapply(OutputsModel$CemaNeigeLayers[InputsCrit$idLayer], FUN = "[[", "SnowPack"))
}
VarSim[!InputsCrit$BoolCrit] <- NA
## Data transformation
if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & warnings) {
if (any(VarObs %in% 0)) {
warning("zeroes detected in 'Qobs': the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL", call. = FALSE)
if (is.vector(x) | inherits(x, "POSIXt")) {
res0 <- x[i]
}
if (any(VarSim %in% 0)) {
warning("zeroes detected in 'Qsim': the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL", call. = FALSE)
}
}
if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon)) {
VarObs <- VarObs + InputsCrit$epsilon
VarSim <- VarSim + InputsCrit$epsilon
}
if (InputsCrit$transfo == "sqrt") {
VarObs <- sqrt(VarObs)
VarSim <- sqrt(VarSim)
}
if (InputsCrit$transfo == "log") {
VarObs <- log(VarObs)
VarSim <- log(VarSim)
VarSim[VarSim < -1e100] <- NA
}
if (InputsCrit$transfo == "inv") {
VarObs <- 1 / VarObs
VarSim <- 1 / VarSim
VarSim[abs(VarSim) > 1e+100] <- NA
}
if (InputsCrit$transfo == "sort") {
VarSim[is.na(VarObs)] <- NA
VarSim <- sort(VarSim, na.last = TRUE)
VarObs <- sort(VarObs, na.last = TRUE)
InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE)
if (is.list(x) & !inherits(x, "POSIXt")) {
if (inherits(x, "OutputsModel")) {
res0 <- .ExtractOutputsModel(x = x, i = i)
} else {
res0 <- .ExtractInputsModel(x = x, i = i)
}
}
return(res0)
})
if (!is.null(x$ZLayers)) {
res$ZLayers <- x$ZLayers
}
if (InputsCrit$transfo == "boxcox") {
VarSim <- (VarSim^0.25 - 0.01 * mean(VarSim, na.rm = TRUE)) / 0.25
VarObs <- (VarObs^0.25 - 0.01 * mean(VarObs, na.rm = TRUE)) / 0.25
if (inherits(x, "SD")) {
res$LengthHydro <- x$LengthHydro
res$BasinAreas <- x$BasinAreas
}
if (grepl("\\^", InputsCrit$transfo)) {
VarObs <- VarObs^transfoPow
VarSim <- VarSim^transfoPow
class(res) <- class(x)
res
}
'[.InputsModel' <- function(x, i) {
if (!inherits(x, "InputsModel")) {
stop("'x' must be of class 'InputsModel'")
}
## TS_ignore
TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit
Ind_TS_ignore <- which(TS_ignore)
if (length(Ind_TS_ignore) == 0) {
Ind_TS_ignore <- NULL
if (is.factor(i)) {
i <- as.character(i)
}
if (sum(!TS_ignore) == 0 | (sum(!TS_ignore) == 1 & Crit %in% c("KGE", "KGE2"))) {
CritCompute <- FALSE
if (is.numeric(i)) {
.ExtractInputsModel(x, i)
} else {
CritCompute <- TRUE
NextMethod()
}
if (inherits(OutputsModel, "hourly")) {
WarningTS <- 365
}
## OutputsModel
.ExtractOutputsModel <- function(x, i) {
res <- lapply(x, function(x) {
if (is.matrix(x) && length(dim(x)) == 2L) {
res0 <- x[i, ]
}
if (is.array(x) && length(dim(x)) == 3L) {
res0 <- x[i, , ]
}
if (is.vector(x) | inherits(x, "POSIXt")) {
res0 <- x[i]
}
if (is.list(x) & !inherits(x, "POSIXt")) {
res0 <- .ExtractOutputsModel(x = x, i = i)
}
return(res0)
})
if (!is.null(x$RunOptions)) {
res$RunOptions <- x$RunOptions
}
if (inherits(OutputsModel, "daily")) {
WarningTS <- 365
if (!is.null(x$StateEnd)) {
res$StateEnd <- x$StateEnd
}
if (inherits(OutputsModel, "monthly")) {
WarningTS <- 12
class(res) <- class(x)
res
}
.IndexOutputsModel <- function(x, i) {
# '[.OutputsModel' <- function(x, i) {
if (!inherits(x, "OutputsModel")) {
stop("'x' must be of class 'OutputsModel'")
}
if (inherits(OutputsModel, "yearly")) {
WarningTS <- 3
if (is.factor(i)) {
i <- as.character(i)
}
if (sum(!TS_ignore) < WarningTS & warnings) {
warning("\t criterion computed on less than ", WarningTS, " time-steps", call. = FALSE)
if (is.numeric(i)) {
.ExtractOutputsModel(x, i)
} else {
NextMethod()
}
## Outputs
OutputsCritCheck <- list(WarningTS = WarningTS,
VarObs = VarObs,
VarSim = VarSim,
CritBestValue = CritBestValue,
Multiplier = Multiplier,
CritName = CritName,
CritVar = CritVar,
CritCompute = CritCompute,
TS_ignore = TS_ignore,
Ind_TS_ignore = Ind_TS_ignore)
}
## =================================================================================
## function to try to set local time in English
## =================================================================================
.TrySetLcTimeEN <- function() {
locale <- list("English_United Kingdom",
"en_US",
"en_US.UTF-8",
"en_US.utf8",
"en")
dateTest <- as.POSIXct("2000-02-15", tz = "UTC", format = "%Y-%m-%d")
monthTestTarget <- "February"
monthTest <- function() {
format(dateTest, format = "%B")
}
lapply(locale, function(x) {
if (monthTest() != monthTestTarget) {
Sys.setlocale(category = "LC_TIME", locale = x)
}
})
}
R/UtilsCalibOptions.R 0 → 100644
View file @ 4bccfee2
.FunTransfo <- function(FeatFUN_MOD) {
IsHyst <- FeatFUN_MOD$IsHyst
IsSD <- FeatFUN_MOD$IsSD
## set FUN_GR
if (FeatFUN_MOD$NameFunMod == "Cemaneige") {
if (IsHyst) {
FUN_GR <- TransfoParam_CemaNeigeHyst
} else {
FUN_GR <- TransfoParam_CemaNeige
}
} else {
# fatal error if the TransfoParam function does not exist
FUN_GR <- match.fun(sprintf("TransfoParam_%s", FeatFUN_MOD$CodeModHydro))
}
## set FUN_SNOW
if ("CemaNeige" %in% FeatFUN_MOD$Class) {
if (IsHyst) {
FUN_SNOW <- TransfoParam_CemaNeigeHyst
} else {
FUN_SNOW <- TransfoParam_CemaNeige
}
}
## set FUN_LAG
if (IsSD) {
FUN_LAG <- TransfoParam_Lag
}
## set FUN_TRANSFO
if (! "CemaNeige" %in% FeatFUN_MOD$Class) {
if (!IsSD) {
FUN_TRANSFO <- FUN_GR
} else {
FUN_TRANSFO <- function(ParamIn, Direction) {
Bool <- is.matrix(ParamIn)
if (!Bool) {
ParamIn <- rbind(ParamIn)
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
ParamOut[, 2:NParam] <- FUN_GR(ParamIn[, 2:NParam], Direction)
ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
return(ParamOut)
}
}
} else {
if (IsHyst & !IsSD) {
FUN_TRANSFO <- function(ParamIn, Direction) {
Bool <- is.matrix(ParamIn)
if (!Bool) {
ParamIn <- rbind(ParamIn)
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
ParamOut[, 1:(NParam - 4) ] <- FUN_GR(ParamIn[, 1:(NParam - 4)], Direction)
ParamOut[, (NParam - 3):NParam] <- FUN_SNOW(ParamIn[, (NParam - 3):NParam], Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
return(ParamOut)
}
}
if (!IsHyst & !IsSD) {
FUN_TRANSFO <- function(ParamIn, Direction) {
Bool <- is.matrix(ParamIn)
if (!Bool) {
ParamIn <- rbind(ParamIn)
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
if (NParam <= 3) {
ParamOut[, 1:(NParam - 2)] <- FUN_GR(cbind(ParamIn[, 1:(NParam - 2)]), Direction)
} else {
ParamOut[, 1:(NParam - 2)] <- FUN_GR(ParamIn[, 1:(NParam - 2)], Direction)
}
ParamOut[, (NParam - 1):NParam] <- FUN_SNOW(ParamIn[, (NParam - 1):NParam], Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
return(ParamOut)
}
}
if (IsHyst & IsSD) {
FUN_TRANSFO <- function(ParamIn, Direction) {
Bool <- is.matrix(ParamIn)
if (!Bool) {
ParamIn <- rbind(ParamIn)
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
ParamOut[, 2:(NParam - 4) ] <- FUN_GR(ParamIn[, 2:(NParam - 4)], Direction)
ParamOut[, (NParam - 3):NParam] <- FUN_SNOW(ParamIn[, (NParam - 3):NParam], Direction)
ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
return(ParamOut)
}
}
if (!IsHyst & IsSD) {
FUN_TRANSFO <- function(ParamIn, Direction) {
Bool <- is.matrix(ParamIn)
if (!Bool) {
ParamIn <- rbind(ParamIn)
}
ParamOut <- NA * ParamIn
NParam <- ncol(ParamIn)
if (NParam <= 3) {
ParamOut[, 2:(NParam - 2)] <- FUN_GR(cbind(ParamIn[, 2:(NParam - 2)]), Direction)
} else {
ParamOut[, 2:(NParam - 2)] <- FUN_GR(ParamIn[, 2:(NParam - 2)], Direction)
}
ParamOut[, (NParam - 1):NParam] <- FUN_SNOW(ParamIn[, (NParam - 1):NParam], Direction)
ParamOut[, 1 ] <- FUN_LAG(as.matrix(ParamIn[, 1]), Direction)
if (!Bool) {
ParamOut <- ParamOut[1, ]
}
return(ParamOut)
}
}
}
if (is.null(FUN_TRANSFO)) {
stop("'FUN_TRANSFO' was not found")
}
return(FUN_TRANSFO)
}
R/UtilsCemaNeige.R 0 → 100644
View file @ 4bccfee2
## daily gradients for mean, min and max air temperature
.GradT_Valery2010 <- as.data.frame(matrix(c(
01, 01, 0.434, 0.366, 0.498,
02, 01, 0.434, 0.366, 0.500,
03, 01, 0.435, 0.367, 0.501,
04, 01, 0.436, 0.367, 0.503,
05, 01, 0.437, 0.367, 0.504,
06, 01, 0.439, 0.367, 0.506,
07, 01, 0.440, 0.367, 0.508,
08, 01, 0.441, 0.368, 0.510,
09, 01, 0.442, 0.368, 0.512,
10, 01, 0.444, 0.368, 0.514,
11, 01, 0.445, 0.368, 0.517,
12, 01, 0.446, 0.368, 0.519,
13, 01, 0.448, 0.369, 0.522,
14, 01, 0.450, 0.369, 0.525,
15, 01, 0.451, 0.369, 0.527,
16, 01, 0.453, 0.370, 0.530,
17, 01, 0.455, 0.370, 0.533,
18, 01, 0.456, 0.370, 0.537,
19, 01, 0.458, 0.371, 0.540,
20, 01, 0.460, 0.371, 0.543,
21, 01, 0.462, 0.371, 0.547,
22, 01, 0.464, 0.372, 0.550,
23, 01, 0.466, 0.372, 0.554,
24, 01, 0.468, 0.373, 0.558,
25, 01, 0.470, 0.373, 0.561,
26, 01, 0.472, 0.374, 0.565,
27, 01, 0.474, 0.374, 0.569,
28, 01, 0.476, 0.375, 0.573,
29, 01, 0.478, 0.375, 0.577,
30, 01, 0.480, 0.376, 0.582,
31, 01, 0.483, 0.376, 0.586,
01, 02, 0.485, 0.377, 0.590,
02, 02, 0.487, 0.377, 0.594,
03, 02, 0.489, 0.378, 0.599,
04, 02, 0.492, 0.379, 0.603,
05, 02, 0.494, 0.379, 0.607,
06, 02, 0.496, 0.380, 0.612,
07, 02, 0.498, 0.381, 0.616,
08, 02, 0.501, 0.381, 0.621,
09, 02, 0.503, 0.382, 0.625,
10, 02, 0.505, 0.383, 0.630,
11, 02, 0.508, 0.384, 0.634,
12, 02, 0.510, 0.384, 0.639,
13, 02, 0.512, 0.385, 0.643,
14, 02, 0.515, 0.386, 0.648,
15, 02, 0.517, 0.387, 0.652,
16, 02, 0.519, 0.387, 0.657,
17, 02, 0.522, 0.388, 0.661,
18, 02, 0.524, 0.389, 0.666,
19, 02, 0.526, 0.390, 0.670,
20, 02, 0.528, 0.391, 0.674,
21, 02, 0.530, 0.392, 0.679,
22, 02, 0.533, 0.393, 0.683,
23, 02, 0.535, 0.393, 0.687,
24, 02, 0.537, 0.394, 0.691,
25, 02, 0.539, 0.395, 0.695,
26, 02, 0.541, 0.396, 0.699,
27, 02, 0.543, 0.397, 0.703,
28, 02, 0.545, 0.398, 0.707,
29, 02, 0.546, 0.399, 0.709,
01, 03, 0.547, 0.399, 0.711,
02, 03, 0.549, 0.400, 0.715,
03, 03, 0.551, 0.401, 0.718,
04, 03, 0.553, 0.402, 0.722,
05, 03, 0.555, 0.403, 0.726,
06, 03, 0.557, 0.404, 0.729,
07, 03, 0.559, 0.405, 0.732,
08, 03, 0.560, 0.406, 0.736,
09, 03, 0.562, 0.406, 0.739,
10, 03, 0.564, 0.407, 0.742,
11, 03, 0.566, 0.408, 0.745,
12, 03, 0.567, 0.409, 0.748,
13, 03, 0.569, 0.410, 0.750,
14, 03, 0.570, 0.411, 0.753,
15, 03, 0.572, 0.412, 0.756,
16, 03, 0.573, 0.413, 0.758,
17, 03, 0.575, 0.414, 0.761,
18, 03, 0.576, 0.415, 0.763,
19, 03, 0.577, 0.416, 0.765,
20, 03, 0.579, 0.417, 0.767,
21, 03, 0.580, 0.417, 0.769,
22, 03, 0.581, 0.418, 0.771,
23, 03, 0.582, 0.419, 0.773,
24, 03, 0.583, 0.420, 0.774,
25, 03, 0.584, 0.421, 0.776,
26, 03, 0.585, 0.422, 0.777,
27, 03, 0.586, 0.422, 0.779,
28, 03, 0.587, 0.423, 0.780,
29, 03, 0.588, 0.424, 0.781,
30, 03, 0.589, 0.425, 0.782,
31, 03, 0.590, 0.425, 0.783,
01, 04, 0.591, 0.426, 0.784,
02, 04, 0.591, 0.427, 0.785,
03, 04, 0.592, 0.427, 0.785,
04, 04, 0.593, 0.428, 0.786,
05, 04, 0.593, 0.429, 0.787,
06, 04, 0.594, 0.429, 0.787,
07, 04, 0.595, 0.430, 0.787,
08, 04, 0.595, 0.431, 0.788,
09, 04, 0.596, 0.431, 0.788,
10, 04, 0.596, 0.432, 0.788,
11, 04, 0.597, 0.432, 0.788,
12, 04, 0.597, 0.433, 0.788,
13, 04, 0.597, 0.433, 0.788,
14, 04, 0.598, 0.434, 0.788,
15, 04, 0.598, 0.434, 0.788,
16, 04, 0.598, 0.435, 0.787,
17, 04, 0.599, 0.435, 0.787,
18, 04, 0.599, 0.436, 0.787,
19, 04, 0.599, 0.436, 0.786,
20, 04, 0.599, 0.436, 0.786,
21, 04, 0.600, 0.437, 0.785,
22, 04, 0.600, 0.437, 0.785,
23, 04, 0.600, 0.437, 0.784,
24, 04, 0.600, 0.438, 0.784,
25, 04, 0.600, 0.438, 0.783,
26, 04, 0.601, 0.438, 0.783,
27, 04, 0.601, 0.438, 0.782,
28, 04, 0.601, 0.439, 0.781,
29, 04, 0.601, 0.439, 0.781,
30, 04, 0.601, 0.439, 0.780,
01, 05, 0.601, 0.439, 0.779,
02, 05, 0.601, 0.439, 0.778,
03, 05, 0.601, 0.439, 0.778,
04, 05, 0.601, 0.440, 0.777,
05, 05, 0.601, 0.440, 0.776,
06, 05, 0.601, 0.440, 0.775,
07, 05, 0.601, 0.440, 0.775,
08, 05, 0.601, 0.440, 0.774,
09, 05, 0.601, 0.440, 0.773,
10, 05, 0.602, 0.440, 0.772,
11, 05, 0.602, 0.440, 0.772,
12, 05, 0.602, 0.440, 0.771,
13, 05, 0.602, 0.440, 0.770,
14, 05, 0.602, 0.440, 0.770,
15, 05, 0.602, 0.440, 0.769,
16, 05, 0.602, 0.440, 0.768,
17, 05, 0.602, 0.440, 0.768,
18, 05, 0.602, 0.440, 0.767,
19, 05, 0.602, 0.440, 0.767,
20, 05, 0.602, 0.440, 0.766,
21, 05, 0.602, 0.440, 0.766,
22, 05, 0.602, 0.440, 0.765,
23, 05, 0.602, 0.440, 0.765,
24, 05, 0.602, 0.440, 0.764,
25, 05, 0.602, 0.440, 0.764,
26, 05, 0.602, 0.440, 0.764,
27, 05, 0.602, 0.439, 0.763,
28, 05, 0.602, 0.439, 0.763,
29, 05, 0.602, 0.439, 0.763,
30, 05, 0.602, 0.439, 0.762,
31, 05, 0.602, 0.439, 0.762,
01, 06, 0.602, 0.439, 0.762,
02, 06, 0.602, 0.439, 0.762,
03, 06, 0.602, 0.439, 0.762,
04, 06, 0.602, 0.439, 0.762,
05, 06, 0.602, 0.439, 0.762,
06, 06, 0.602, 0.438, 0.761,
07, 06, 0.602, 0.438, 0.761,
08, 06, 0.602, 0.438, 0.761,
09, 06, 0.602, 0.438, 0.761,
10, 06, 0.602, 0.438, 0.761,
11, 06, 0.602, 0.438, 0.762,
12, 06, 0.602, 0.438, 0.762,
13, 06, 0.602, 0.438, 0.762,
14, 06, 0.602, 0.438, 0.762,
15, 06, 0.602, 0.437, 0.762,
16, 06, 0.602, 0.437, 0.762,
17, 06, 0.602, 0.437, 0.762,
18, 06, 0.602, 0.437, 0.762,
19, 06, 0.602, 0.437, 0.763,
20, 06, 0.602, 0.437, 0.763,
21, 06, 0.602, 0.437, 0.763,
22, 06, 0.602, 0.436, 0.763,
23, 06, 0.602, 0.436, 0.763,
24, 06, 0.602, 0.436, 0.764,
25, 06, 0.602, 0.436, 0.764,
26, 06, 0.601, 0.436, 0.764,
27, 06, 0.601, 0.436, 0.764,
28, 06, 0.601, 0.436, 0.764,
29, 06, 0.601, 0.435, 0.765,
30, 06, 0.601, 0.435, 0.765,
01, 07, 0.601, 0.435, 0.765,
02, 07, 0.600, 0.435, 0.765,
03, 07, 0.600, 0.435, 0.765,
04, 07, 0.600, 0.434, 0.766,
05, 07, 0.600, 0.434, 0.766,
06, 07, 0.599, 0.434, 0.766,
07, 07, 0.599, 0.434, 0.766,
08, 07, 0.599, 0.434, 0.766,
09, 07, 0.598, 0.433, 0.766,
10, 07, 0.598, 0.433, 0.766,
11, 07, 0.598, 0.433, 0.766,
12, 07, 0.597, 0.433, 0.766,
13, 07, 0.597, 0.432, 0.767,
14, 07, 0.597, 0.432, 0.767,
15, 07, 0.596, 0.432, 0.767,
16, 07, 0.596, 0.432, 0.766,
17, 07, 0.595, 0.431, 0.766,
18, 07, 0.595, 0.431, 0.766,
19, 07, 0.594, 0.431, 0.766,
20, 07, 0.594, 0.430, 0.766,
21, 07, 0.593, 0.430, 0.766,
22, 07, 0.593, 0.430, 0.766,
23, 07, 0.592, 0.429, 0.765,
24, 07, 0.592, 0.429, 0.765,
25, 07, 0.591, 0.428, 0.765,
26, 07, 0.590, 0.428, 0.765,
27, 07, 0.590, 0.428, 0.764,
28, 07, 0.589, 0.427, 0.764,
29, 07, 0.588, 0.427, 0.764,
30, 07, 0.588, 0.426, 0.763,
31, 07, 0.587, 0.426, 0.763,
01, 08, 0.586, 0.425, 0.762,
02, 08, 0.586, 0.425, 0.762,
03, 08, 0.585, 0.424, 0.761,
04, 08, 0.584, 0.424, 0.761,
05, 08, 0.583, 0.423, 0.760,
06, 08, 0.583, 0.423, 0.760,
07, 08, 0.582, 0.422, 0.759,
08, 08, 0.581, 0.421, 0.758,
09, 08, 0.580, 0.421, 0.758,
10, 08, 0.579, 0.420, 0.757,
11, 08, 0.578, 0.420, 0.756,
12, 08, 0.578, 0.419, 0.755,
13, 08, 0.577, 0.418, 0.754,
14, 08, 0.576, 0.418, 0.754,
15, 08, 0.575, 0.417, 0.753,
16, 08, 0.574, 0.416, 0.752,
17, 08, 0.573, 0.415, 0.751,
18, 08, 0.572, 0.415, 0.750,
19, 08, 0.571, 0.414, 0.749,
20, 08, 0.570, 0.413, 0.748,
21, 08, 0.569, 0.413, 0.747,
22, 08, 0.569, 0.412, 0.746,
23, 08, 0.568, 0.411, 0.745,
24, 08, 0.567, 0.410, 0.744,
25, 08, 0.566, 0.409, 0.743,
26, 08, 0.565, 0.409, 0.742,
27, 08, 0.564, 0.408, 0.741,
28, 08, 0.563, 0.407, 0.740,
29, 08, 0.562, 0.406, 0.738,
30, 08, 0.561, 0.405, 0.737,
31, 08, 0.560, 0.405, 0.736,
01, 09, 0.558, 0.404, 0.735,
02, 09, 0.557, 0.403, 0.734,
03, 09, 0.556, 0.402, 0.732,
04, 09, 0.555, 0.401, 0.731,
05, 09, 0.554, 0.401, 0.730,
06, 09, 0.553, 0.400, 0.728,
07, 09, 0.552, 0.399, 0.727,
08, 09, 0.551, 0.398, 0.725,
09, 09, 0.550, 0.397, 0.724,
10, 09, 0.549, 0.396, 0.723,
11, 09, 0.548, 0.396, 0.721,
12, 09, 0.546, 0.395, 0.720,
13, 09, 0.545, 0.394, 0.718,
14, 09, 0.544, 0.393, 0.717,
15, 09, 0.543, 0.392, 0.715,
16, 09, 0.542, 0.391, 0.713,
17, 09, 0.541, 0.391, 0.712,
18, 09, 0.540, 0.390, 0.710,
19, 09, 0.538, 0.389, 0.709,
20, 09, 0.537, 0.388, 0.707,
21, 09, 0.536, 0.388, 0.705,
22, 09, 0.535, 0.387, 0.703,
23, 09, 0.533, 0.386, 0.702,
24, 09, 0.532, 0.385, 0.700,
25, 09, 0.531, 0.385, 0.698,
26, 09, 0.530, 0.384, 0.696,
27, 09, 0.528, 0.383, 0.694,
28, 09, 0.527, 0.383, 0.692,
29, 09, 0.526, 0.382, 0.690,
30, 09, 0.525, 0.381, 0.688,
01, 10, 0.523, 0.381, 0.686,
02, 10, 0.522, 0.380, 0.684,
03, 10, 0.521, 0.379, 0.682,
04, 10, 0.519, 0.379, 0.680,
05, 10, 0.518, 0.378, 0.678,
06, 10, 0.517, 0.377, 0.676,
07, 10, 0.515, 0.377, 0.674,
08, 10, 0.514, 0.376, 0.671,
09, 10, 0.512, 0.376, 0.669,
10, 10, 0.511, 0.375, 0.667,
11, 10, 0.510, 0.375, 0.664,
12, 10, 0.508, 0.374, 0.662,
13, 10, 0.507, 0.374, 0.659,
14, 10, 0.505, 0.373, 0.657,
15, 10, 0.504, 0.373, 0.654,
16, 10, 0.502, 0.372, 0.652,
17, 10, 0.501, 0.372, 0.649,
18, 10, 0.499, 0.372, 0.647,
19, 10, 0.498, 0.371, 0.644,
20, 10, 0.496, 0.371, 0.641,
21, 10, 0.495, 0.371, 0.639,
22, 10, 0.493, 0.370, 0.636,
23, 10, 0.492, 0.370, 0.633,
24, 10, 0.490, 0.370, 0.630,
25, 10, 0.489, 0.369, 0.628,
26, 10, 0.487, 0.369, 0.625,
27, 10, 0.485, 0.369, 0.622,
28, 10, 0.484, 0.368, 0.619,
29, 10, 0.482, 0.368, 0.616,
30, 10, 0.481, 0.368, 0.613,
31, 10, 0.479, 0.368, 0.610,
01, 11, 0.478, 0.368, 0.607,
02, 11, 0.476, 0.367, 0.604,
03, 11, 0.475, 0.367, 0.601,
04, 11, 0.473, 0.367, 0.598,
05, 11, 0.471, 0.367, 0.595,
06, 11, 0.470, 0.367, 0.592,
07, 11, 0.468, 0.367, 0.589,
08, 11, 0.467, 0.366, 0.586,
09, 11, 0.465, 0.366, 0.583,
10, 11, 0.464, 0.366, 0.580,
11, 11, 0.462, 0.366, 0.577,
12, 11, 0.461, 0.366, 0.574,
13, 11, 0.459, 0.366, 0.571,
14, 11, 0.458, 0.366, 0.568,
15, 11, 0.456, 0.366, 0.565,
16, 11, 0.455, 0.366, 0.562,
17, 11, 0.454, 0.366, 0.559,
18, 11, 0.452, 0.365, 0.556,
19, 11, 0.451, 0.365, 0.553,
20, 11, 0.450, 0.365, 0.550,
21, 11, 0.448, 0.365, 0.547,
22, 11, 0.447, 0.365, 0.544,
23, 11, 0.446, 0.365, 0.542,
24, 11, 0.445, 0.365, 0.539,
25, 11, 0.443, 0.365, 0.536,
26, 11, 0.442, 0.365, 0.533,
27, 11, 0.441, 0.365, 0.531,
28, 11, 0.440, 0.365, 0.528,
29, 11, 0.439, 0.365, 0.526,
30, 11, 0.438, 0.365, 0.523,
01, 12, 0.437, 0.365, 0.521,
02, 12, 0.436, 0.365, 0.519,
03, 12, 0.435, 0.365, 0.517,
04, 12, 0.434, 0.365, 0.515,
05, 12, 0.434, 0.365, 0.513,
06, 12, 0.433, 0.365, 0.511,
07, 12, 0.432, 0.365, 0.509,
08, 12, 0.431, 0.365, 0.507,
09, 12, 0.431, 0.365, 0.505,
10, 12, 0.430, 0.365, 0.504,
11, 12, 0.430, 0.365, 0.502,
12, 12, 0.429, 0.365, 0.501,
13, 12, 0.429, 0.365, 0.500,
14, 12, 0.429, 0.365, 0.498,
15, 12, 0.428, 0.365, 0.497,
16, 12, 0.428, 0.365, 0.496,
17, 12, 0.428, 0.365, 0.496,
18, 12, 0.428, 0.365, 0.495,
19, 12, 0.428, 0.365, 0.494,
20, 12, 0.428, 0.365, 0.494,
21, 12, 0.428, 0.365, 0.494,
22, 12, 0.428, 0.365, 0.493,
23, 12, 0.429, 0.365, 0.493,
24, 12, 0.429, 0.366, 0.493,
25, 12, 0.429, 0.366, 0.493,
26, 12, 0.430, 0.366, 0.494,
27, 12, 0.430, 0.366, 0.494,
28, 12, 0.431, 0.366, 0.495,
29, 12, 0.431, 0.366, 0.495,
30, 12, 0.432, 0.366, 0.496,
31, 12, 0.433, 0.366, 0.497),
ncol = 5, byrow = TRUE,
dimnames = list(NULL, c("day", "month", "grad_Tmean", "grad_Tmin", "grad_Tmax"))))
R/UtilsErrorCrit.R 0 → 100644
View file @ 4bccfee2
## =================================================================================
## function to manage inputs of specific ErrorCrit_*() functions
## =================================================================================
.ErrorCrit <- function(InputsCrit, Crit, OutputsModel, warnings) {
## Arguments check
if (!inherits(InputsCrit, "InputsCrit")) {
stop("'InputsCrit' must be of class 'InputsCrit'", call. = FALSE)
}
if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) {
stop(paste0("'InputsCrit' must be of class 'Single'. Use the 'ErrorCrit' function on objects of class 'Multi' or 'Compo' with ", Crit), call. = FALSE)
}
## Initialisation
CritName <- NA
CritVar <- InputsCrit$VarObs
if (InputsCrit$transfo == "") {
CritName <- paste0(Crit, "[CritVar]")
}
if (InputsCrit$transfo %in% c("sqrt", "log", "sort", "boxcox")) {
CritName <- paste0(Crit, "[", InputsCrit$transfo, "(CritVar)]")
}
if (InputsCrit$transfo == "inv") {
CritName <- paste0(Crit, "[1/CritVar]")
}
if (grepl("\\^", InputsCrit$transfo)) {
transfoPow <- suppressWarnings(as.numeric(gsub("\\^", "", InputsCrit$transfo)))
CritName <- paste0(Crit, "[CritVar^", transfoPow, "]")
}
CritName <- gsub(pattern = "CritVar", replacement = CritVar, x = CritName)
CritValue <- NA
if (Crit %in% c("RMSE")) {
CritBestValue <- +1
Multiplier <- +1
}
if (Crit %in% c("NSE", "KGE", "KGE2", "GAPX")) {
CritBestValue <- +1
Multiplier <- -1
}
## Data preparation
VarObs <- InputsCrit$Obs
VarObs[!InputsCrit$BoolCrit] <- NA
VarSim <- switch(
InputsCrit$VarObs,
Q = OutputsModel$Qsim,
SCA = rowMeans(sapply(OutputsModel$CemaNeigeLayers[InputsCrit$idLayer], FUN = "[[", "Gratio")),
SWE = rowMeans(sapply(OutputsModel$CemaNeigeLayers[InputsCrit$idLayer], FUN = "[[", "SnowPack")),
ParamT = OutputsModel$RunOptions$ParamT
)
VarSim[!InputsCrit$BoolCrit] <- NA
## Data transformation
if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & warnings) {
if (any(VarObs %in% 0)) {
warning("zeroes detected in 'Qobs': the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL", call. = FALSE)
}
if (any(VarSim %in% 0)) {
warning("zeroes detected in 'Qsim': the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL", call. = FALSE)
}
}
if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon) & !(InputsCrit$transfo == "boxcox")) {
VarObs <- VarObs + InputsCrit$epsilon
VarSim <- VarSim + InputsCrit$epsilon
}
if (InputsCrit$transfo == "sqrt") {
VarObs <- sqrt(VarObs)
VarSim <- sqrt(VarSim)
}
if (InputsCrit$transfo == "log") {
VarObs <- log(VarObs)
VarSim <- log(VarSim)
VarSim[VarSim < -1e100] <- NA
}
if (InputsCrit$transfo == "inv") {
VarObs <- 1 / VarObs
VarSim <- 1 / VarSim
VarSim[abs(VarSim) > 1e+100] <- NA
}
if (InputsCrit$transfo == "sort") {
VarSim[is.na(VarObs)] <- NA
VarSim <- sort(VarSim, na.last = TRUE)
VarObs <- sort(VarObs, na.last = TRUE)
InputsCrit$BoolCrit <- sort(InputsCrit$BoolCrit, decreasing = TRUE)
}
if (InputsCrit$transfo == "boxcox") {
muTransfoVarObs <- (0.01 * mean(VarObs, na.rm = TRUE))^0.25
VarSim <- (VarSim^0.25 - muTransfoVarObs) / 0.25
VarObs <- (VarObs^0.25 - muTransfoVarObs) / 0.25
}
if (grepl("\\^", InputsCrit$transfo)) {
VarObs <- VarObs^transfoPow
VarSim <- VarSim^transfoPow
}
## TS_ignore
TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit
Ind_TS_ignore <- which(TS_ignore)
if (length(Ind_TS_ignore) == 0) {
Ind_TS_ignore <- NULL
}
if (sum(!TS_ignore) == 0 | (sum(!TS_ignore) == 1 & Crit %in% c("KGE", "KGE2"))) {
CritCompute <- FALSE
} else {
CritCompute <- TRUE
}
if (Crit != "GAPX") {
WarningTS <- 10
if (sum(!TS_ignore) < WarningTS & warnings) {
warning("\t criterion computed on less than ", WarningTS, " time-steps", call. = FALSE)
}
} else {
WarningTS <- 4 # For at least daily time step models (GR4J)
if (sum(!TS_ignore) < WarningTS & warnings) {
warning("\t criterion GAPX computed on less than ", WarningTS, " parameters", call. = FALSE)
}
}
## Outputs
OutputsCritCheck <- list(WarningTS = WarningTS,
VarObs = VarObs,
VarSim = VarSim,
CritBestValue = CritBestValue,
Multiplier = Multiplier,
CritName = CritName,
CritVar = CritVar,
CritCompute = CritCompute,
TS_ignore = TS_ignore,
Ind_TS_ignore = Ind_TS_ignore)
}
R/UtilsRunModel.R 0 → 100644
View file @ 4bccfee2
#' Create `OutputsModel` for GR non-Cemaneige models
#'
#' @param InputsModel output of [CreateInputsModel]
#' @param RunOptions output of [CreateRunOptions]
#' @param RESULTS outputs of [.Fortran]
#' @param LInputSeries number of time steps of warm-up + run periods
#' @param Param [numeric] vector of model parameters
#' @param CemaNeigeLayers outputs of Cemaneige pre-process
#'
#' @return OutputsModel object
#' @noRd
#'
.GetOutputsModelGR <- function(InputsModel,
RunOptions,
RESULTS,
LInputSeries,
Param,
CemaNeigeLayers = NULL) {
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
FortranOutputs <- RunOptions$FortranOutputs$GR
IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim)
OutputsModel <- list()
if ("DatesR" %in% RunOptions$Outputs_Sim) {
OutputsModel$DatesR <- InputsModel$DatesR[RunOptions$IndPeriod_Run]
}
seqOutputs <- seq_len(RESULTS$NOutputs)
names(seqOutputs) <- FortranOutputs[IndOutputs]
OutputsModel <- c(OutputsModel,
lapply(seqOutputs, function(i) RESULTS$Outputs[IndPeriod2, i]))
if (!is.null(CemaNeigeLayers)) {
OutputsModel$CemaNeigeLayers <- CemaNeigeLayers
}
if ("WarmUpQsim" %in% RunOptions$Outputs_Sim && !identical(RunOptions$IndPeriod_WarmUp, 0L)) {
OutputsModel$RunOptions$WarmUpQsim <- RESULTS$Outputs[seq_len(length(RunOptions$IndPeriod_WarmUp)),
which(FortranOutputs == "Qsim")]
# class(OutputsModel$RunOptions$WarmUpQsim) <- c("WarmUpOutputsModelItem", class(OutputsModel$RunOptions$WarmUpQsim))
}
if ("Param" %in% RunOptions$Outputs_Sim) {
OutputsModel$RunOptions$Param <- Param
}
if ("StateEnd" %in% RunOptions$Outputs_Sim) {
OutputsModel$StateEnd <- RESULTS$StateEnd
}
class(OutputsModel) <- c("OutputsModel", class(RunOptions)[-1])
return(OutputsModel)
}
#' Check arguments of `RunModel_*GR*` functions
#'
#' @param InputsModel see [CreateInputsModel]
#' @param RunOptions see [CreateRunOptions]
#' @param Param [numeric] [vector] model calibration parameters
#'
#' @return [NULL]
#' @noRd
#'
.ArgumentsCheckGR <- function(InputsModel, RunOptions, Param) {
if (!inherits(InputsModel, "InputsModel")) {
stop("'InputsModel' must be of class 'InputsModel'")
}
if (!inherits(InputsModel, RunOptions$FeatFUN_MOD$TimeUnit)) {
stop("'InputsModel' must be of class '", RunOptions$FeatFUN_MOD$TimeUnit, "'")
}
if (!inherits(InputsModel, "GR")) {
stop("'InputsModel' must be of class 'GR'")
}
if (class(RunOptions)[1] != "RunOptions") {
if (!inherits(RunOptions, "RunOptions")) {
stop("'RunOptions' must be of class 'RunOptions'")
} else {
stop("'RunOptions' class of 'RunOptions' must be in first position")
}
}
if (!inherits(RunOptions, "GR")) {
stop("'RunOptions' must be of class 'GR'")
}
if ("CemaNeige" %in% RunOptions$FeatFUN_MOD$Class) {
if (!inherits(InputsModel, "CemaNeige")) {
stop("'InputsModel' must be of class 'CemaNeige'")
}
if (!inherits(RunOptions, "CemaNeige")) {
stop("'RunOptions' must be of class 'CemaNeige'")
}
}
if (!is.vector(Param) | !is.numeric(Param)) {
stop("'Param' must be a numeric vector")
}
if (sum(!is.na(Param)) != RunOptions$FeatFUN_MOD$NbParam) {
stop(paste("'Param' must be a vector of length", RunOptions$FeatFUN_MOD$NbParam, "and contain no NA"))
}
}
R/UtilsSeriesAggreg.R 0 → 100644
View file @ 4bccfee2
.GetSeriesAggregFormat <- function(NewTimeFormat) {
errNewTimeFormat <- FALSE
if (missing(NewTimeFormat)) {
errNewTimeFormat <- TRUE
} else if (is.null(NewTimeFormat)) {
errNewTimeFormat <- TRUE
}
if (errNewTimeFormat) {
stop("Argument `Format` is missing")
}
if (!is.null(NewTimeFormat)) {
TimeStep <- c("hourly", "daily", "monthly", "yearly")
NewTimeFormat <- match.arg(NewTimeFormat, choices = TimeStep)
Format <- switch(NewTimeFormat,
hourly = "%Y%m%d%h",
daily = "%Y%m%d",
monthly = "%Y%m",
yearly = "%Y")
msgNewTimeFormat <- sprintf("'Format' automatically set to %s", sQuote(Format))
warning("deprecated 'NewTimeFormat' argument: please use 'Format' instead.",
msgNewTimeFormat,
call. = FALSE)
return(Format)
}
return(NULL)
}
.GetSeriesAggregClass <- function(Format) {
Format <- substr(Format,
start = nchar(Format),
stop = nchar(Format))
switch(Format,
h = "hourly",
d = "daily",
m = "monthly",
Y = "yearly")
}
.GetAggregConvertFun <- function(x, Format) {
AggregConvertFunTable <- rbind(
data.frame(ConvertFun = "mean",
x = c("Prod", "Rout", "Exp", "SnowPack", "ThermalState",
"Gratio", "Temp", "Gthreshold", "Glocalmax", "LayerTempMean", "T"),
stringsAsFactors = FALSE), # R < 4.0 compatibility: avoids mixing numeric and factor into numeric
data.frame(ConvertFun = "sum",
x = c("PotEvap", "Precip", "Pn", "Ps", "AE", "Perc", "PR", "Q9",
"Q1", "Exch", "AExch1", "AExch2", "AExch", "QR", "QRExp",
"QD", "Qsim", "Pliq", "Psol", "PotMelt", "Melt", "PliqAndMelt",
"LayerPrecip", "LayerFracSolidPrecip", "Qmm", "Qls", "E", "P", "Qupstream"),
stringsAsFactors = FALSE) # R < 4.0 compatibility: avoids mixing numeric and factor into numeric
)
res <- sapply(x, function(iX) {
iRes <- AggregConvertFunTable$ConvertFun[AggregConvertFunTable$x == iX]
iRes <- ifelse(test = any(is.na(iRes)), yes = NA, no = iRes) # R < 4.0 compatibility
})
if (Format %in% c("%d", "%m")) {
res <- rep("mean", length(res))
}
return(res)
}