Utils.R 7.65 KB
Newer Older
1

2
3
4
5
6
7
8
9
10
11
12
13
## =================================================================================
## function to check
## =================================================================================

# .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.")
#     }
#   }
# }

14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

## =================================================================================
## function to manage Fortran outputs
## =================================================================================

.FortranOutputs <- function(GR = NULL, isCN = FALSE) {
  
  outGR <- NULL
  outCN <- NULL
  
  if (is.null(GR)) {
    GR <- ""
  }
  if (GR == "GR1A") {
    outGR <- c("PotEvap", "Precip",
               "Qsim")
  } else if (GR == "GR2M") {
    outGR <- c("PotEvap", "Precip", "Prod", "Pn",
               "AE",
               "Perc", "PR",
               "Rout", "Exch",
               "Qsim")
  } else if (GR == "GR4H") {
    outGR <- c("PotEvap", "Precip", "Prod",
               "AE",
               "Perc", "PR",
               "Q9", "Q1",
               "Rout", "Exch",
               "AExch", "QR",
43
44
45
46
               "QD",
               "Qsim")
  } else if (GR == "GR5H") {
    outGR <- c("PotEvap", "Precip", "Interc", "Prod", "Pn", "Ps",
47
               "AE", "EI", "ES",
48
49
50
51
52
               "Perc", "PR",
               "Q9", "Q1",
               "Rout", "Exch", 
               "AExch1", "AExch2",
               "AExch", "QR",
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
               "QD",
               "Qsim")
  } else if (GR %in% c("GR4J", "GR5J")) {
    outGR <- c("PotEvap", "Precip", "Prod", "Pn", "Ps",
               "AE",
               "Perc", "PR",
               "Q9", "Q1",
               "Rout", "Exch", 
               "AExch1", "AExch2",
               "AExch", "QR",
               "QD",
               "Qsim")
  } else if (GR == "GR6J") {
    outGR <- c("PotEvap", "Precip", "Prod", "Pn", "Ps",
               "AE",
               "Perc", "PR",
               "Q9", "Q1",
               "Rout", "Exch",
               "AExch1", "AExch2",
               "AExch", "QR",
               "QRExp", "Exp",
               "QD",
               "Qsim")
  }
  if (isCN) {
    outCN <- c("Pliq", "Psol", 
               "SnowPack", "ThermalState", "Gratio", 
               "PotMelt", "Melt", "PliqAndMelt", "Temp", 
               "Gthreshold", "Glocalmax")
  }
  
  res <- list(GR = outGR, CN = outCN)
  
}



90
91
92
## =================================================================================
## function to manage inputs of specific ErrorCrit_*() functions
## =================================================================================
93

94
.ErrorCrit <- function(InputsCrit, Crit, OutputsModel, warnings) {
95
96
97
98
99
100
  
  ## Arguments check
  if (!inherits(InputsCrit, "InputsCrit")) {
    stop("'InputsCrit' must be of class 'InputsCrit'", call. = FALSE)
  }
  if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) {
101
    if (Crit == "RMSE") {
102
103
      stop("'InputsCrit' must be of class 'Single'. Use the 'ErrorCrit' function on objects of class 'Multi' with RMSE", call. = FALSE)
    } else {
104
      stop(paste0("'InputsCrit' must be of class 'Single'. Use the 'ErrorCrit' function on objects of class 'Multi' or 'Compo' with ", Crit), call. = FALSE)
105
106
107
108
109
110
111
112
    }
  }
  
  
  ## Initialisation
  CritName <- NA
  CritVar  <- InputsCrit$VarObs
  if (InputsCrit$transfo == "") {
113
    CritName <- paste0(Crit, "[CritVar]")
114
  }
115
  if (InputsCrit$transfo %in% c("sqrt", "log", "sort", "boxcox")) {
116
    CritName <- paste0(Crit, "[", InputsCrit$transfo, "(CritVar)]")
117
118
  }
  if (InputsCrit$transfo == "inv") {
119
    CritName <- paste0(Crit, "[1/CritVar]")
120
  }
121
122
123
  if (grepl("\\^", InputsCrit$transfo)) {
    transfoPow <- suppressWarnings(as.numeric(gsub("\\^", "", InputsCrit$transfo)))
    CritName <- paste0(Crit, "[CritVar^", transfoPow, "]")
124
125
126
  }
  CritName  <- gsub(pattern = "CritVar", replacement = CritVar, x = CritName)
  CritValue <- NA
127
  if (Crit %in% c("RMSE")) {
128
129
130
    CritBestValue <- +1
    Multiplier    <- +1
  }
131
  if (Crit %in% c("NSE", "KGE", "KGE2")) {
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
    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 (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)
    }  
  }
161
  if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon) & !(InputsCrit$transfo == "boxcox")) {
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
    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)
  }
185
  if (InputsCrit$transfo == "boxcox") {
186
187
188
    muTransfoVarObs <- (0.01 * mean(VarObs, na.rm = TRUE))^0.25
    VarSim <- (VarSim^0.25 - muTransfoVarObs) / 0.25
    VarObs <- (VarObs^0.25 - muTransfoVarObs) / 0.25
189
  }
190
  if (grepl("\\^", InputsCrit$transfo)) {
191
192
193
194
195
196
197
198
199
200
201
    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
  }
202
  if (sum(!TS_ignore) == 0 | (sum(!TS_ignore) == 1 & Crit %in% c("KGE", "KGE2"))) { 
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
    CritCompute <- FALSE
  } else {
    CritCompute <- TRUE
  }
  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", call. = FALSE)
  }
  
  
  ## Outputs
  OutputsCritCheck <- list(WarningTS = WarningTS,
                           VarObs = VarObs, 
                           VarSim = VarSim, 
                           CritBestValue = CritBestValue, 
                           Multiplier = Multiplier, 
                           CritName = CritName, 
                           CritVar = CritVar, 
232
233
234
                           CritCompute = CritCompute,
                           TS_ignore = TS_ignore,
                           Ind_TS_ignore = Ind_TS_ignore)
235
}
236
237