ErrorCrit_KGE2 <- function(InputsCrit, OutputsModel, warnings = TRUE, verbose = TRUE) { ##Arguments_check________________________________ if (!inherits(InputsCrit, "InputsCrit")) { stop("InputsCrit must be of class 'InputsCrit'") } if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) { stop("InputsCrit must be of class 'Single'. Use the ErrorCrit function on objects of class 'Multi' or 'Compo' with KGE'") } if (!inherits(OutputsModel, "OutputsModel")) { stop("OutputsModel must be of class 'OutputsModel'") } ##Initialisation_________________________________ CritName <- NA CritVar <- InputsCrit$varObs if (InputsCrit$transfo == "") { CritName <- "KGE'[CritVar]" } if (InputsCrit$transfo == "sqrt") { CritName <- "KGE'[sqrt(CritVar)]" } if (InputsCrit$transfo == "log") { CritName <- "KGE'[log(CritVar)]" } if (InputsCrit$transfo == "inv") { CritName <- "KGE'[1/CritVar]" } if (InputsCrit$transfo == "sort") { CritName <- "KGE'[sort(CritVar)]" } CritName <- gsub(pattern = "CritVar", replacement = CritVar, x = CritName) CritValue <- NA CritBestValue <- +1 Multiplier <- -1 ### must be equal to -1 or +1 only ##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") } 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") } } 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) } ##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) { OutputsCrit <- list(NA) names(OutputsCrit) <- c("CritValue") return(OutputsCrit) } if (sum(!TS_ignore) == 1) { OutputsCrit <- list(NA) names(OutputsCrit) <- c("CritValue") return(OutputsCrit) } ### to avoid a problem in standard deviation computation if (inherits(OutputsModel, "hourly")) { WarningTS <- 365 } if (inherits(OutputsModel, "daily")) { WarningTS <- 365 } if (inherits(OutputsModel, "monthly")) { WarningTS <- 12 } if (inherits(OutputsModel, "yearly")) { WarningTS <- 3 } if (sum(!TS_ignore) < WarningTS & warnings) { warning("\t criterion computed on less than ", WarningTS, " time-steps") } ##Other_variables_preparation meanVarObs <- mean(VarObs[!TS_ignore]) meanVarSim <- mean(VarSim[!TS_ignore]) iCrit <- 0 SubCritPrint <- NULL SubCritNames <- NULL SubCritValues <- NULL ##SubErrorCrit_____KGE_rPearson__________________ iCrit <- iCrit + 1 SubCritPrint[iCrit] <- paste(CritName, " cor(sim, obs, \"pearson\") =", sep = "") SubCritValues[iCrit] <- NA SubCritNames[iCrit] <- "r" Numer <- sum((VarObs[!TS_ignore] - meanVarObs) * (VarSim[!TS_ignore] - meanVarSim)) Deno1 <- sqrt(sum((VarObs[!TS_ignore] - meanVarObs)^2)) Deno2 <- sqrt(sum((VarSim[!TS_ignore] - meanVarSim)^2)) if (Numer == 0) { if (Deno1 == 0 & Deno2 == 0) { Crit <- 1 } else { Crit <- 0 } } else { Crit <- Numer / (Deno1 * Deno2) } if (is.numeric(Crit) & is.finite(Crit)) { SubCritValues[iCrit] <- Crit } ##SubErrorCrit_____KGE_gamma______________________ iCrit <- iCrit + 1 SubCritPrint[iCrit] <- paste(CritName, " cv(sim)/cv(obs) =", sep = "") SubCritValues[iCrit] <- NA SubCritNames[iCrit] <- "gamma" if (meanVarSim == 0) { if (sd(VarSim[!TS_ignore]) == 0) { CVsim <- 1 } else { CVsim <- 99999 } } else { CVsim <- sd(VarSim[!TS_ignore]) / meanVarSim } if (meanVarObs == 0) { if (sd(VarObs[!TS_ignore]) == 0) { CVobs <- 1 } else { CVobs <- 99999 } } else { CVobs <- sd(VarObs[!TS_ignore]) / meanVarObs } if (CVsim == 0 & CVobs == 0) { Crit <- 1 } else { Crit <- CVsim / CVobs } if (is.numeric(Crit) & is.finite(Crit)) { SubCritValues[iCrit] <- Crit } ##SubErrorCrit_____KGE_beta______________________ iCrit <- iCrit + 1 SubCritPrint[iCrit] <- paste(CritName, " mean(sim)/mean(obs) =", sep = "") SubCritValues[iCrit] <- NA SubCritNames[iCrit] <- "beta" if (meanVarSim == 0 & meanVarObs == 0) { Crit <- 1 } else { Crit <- meanVarSim / meanVarObs } if (is.numeric(Crit) & is.finite(Crit)) { SubCritValues[iCrit] <- Crit } ##ErrorCrit______________________________________ if (sum(is.na(SubCritValues)) == 0) { CritValue <- (1 - sqrt((SubCritValues[1] - 1)^2 + (SubCritValues[2] - 1)^2 + (SubCritValues[3] - 1)^2)) } ##Verbose______________________________________ if (verbose) { message("Crit. ", CritName, " = ", sprintf("%.4f", CritValue)) message(paste("\tSubCrit.", SubCritPrint, sprintf("%.4f", SubCritValues), "\n", sep = " ")) } ##Output_________________________________________ OutputsCrit <- list(CritValue = CritValue, CritName = CritName, SubCritValues = SubCritValues, SubCritNames = SubCritNames, CritBestValue = CritBestValue, Multiplier = Multiplier, Ind_notcomputed = Ind_TS_ignore ) class(OutputsCrit) <- c("KGE2", "ErrorCrit") return(OutputsCrit) }