Forked from HYCAR-Hydro / airGR
Source project has a limited visibility.
ErrorCrit_KGE.R 6.52 KiB
ErrorCrit_KGE <- 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)
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
} 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))