Skip to content
GitLab

Source

Target

Hide whitespace changes
Inline Side-by-side
Showing
with 2260 additions and 392 deletions
+2260 -392
R/Utils.R 0 → 100644
View file @ 4bccfee2
## =================================================================================
## 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.")
# }
# }
# }
## =================================================================================
## 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 <- ""
}
if (GR == "GR1A") {
outGR <- c("PotEvap", "Precip",
"Qsim")
} else if (GR == "GR2M") {
outGR <- c("PotEvap", "Precip", "Prod", "Pn", "Ps",
"AE",
"Perc", "PR",
"Rout",
"AExch",
"Qsim")
} 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", "GR4H")) {
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)
}
## =================================================================================
## functions to extract parts of InputsModel or OutputsModel objects
## =================================================================================
## InputsModel
.ExtractInputsModel <- function(x, i) {
res <- lapply(x, function(x) {
if (is.matrix(x)) {
res0 <- x[i, , drop = FALSE]
}
if (is.vector(x) | inherits(x, "POSIXt")) {
res0 <- x[i]
}
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 (inherits(x, "SD")) {
res$LengthHydro <- x$LengthHydro
res$BasinAreas <- x$BasinAreas
}
class(res) <- class(x)
res
}
'[.InputsModel' <- function(x, i) {
if (!inherits(x, "InputsModel")) {
stop("'x' must be of class 'InputsModel'")
}
if (is.factor(i)) {
i <- as.character(i)
}
if (is.numeric(i)) {
.ExtractInputsModel(x, i)
} else {
NextMethod()
}
}
## 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 (!is.null(x$StateEnd)) {
res$StateEnd <- x$StateEnd
}
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 (is.factor(i)) {
i <- as.character(i)
}
if (is.numeric(i)) {
.ExtractOutputsModel(x, i)
} else {
NextMethod()
}
}
## =================================================================================
## 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)
}
R/plot.OutputsModel.R
View file @ 4bccfee2
plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = NULL, which = "all", log_scale = FALSE,
cex.axis = 1, cex.lab = 0.9, cex.leg = 0.9, lwd = 1, verbose = TRUE, ...) {
plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = NULL, which = "synth", log_scale = FALSE,
cex.axis = 1, cex.lab = 0.9, cex.leg = 0.9, lwd = 1,
AxisTS = function(x) axis.POSIXct(side = 1, x = x$DatesR, ...),
LayoutMat = NULL, LayoutWidths = rep.int(1, ncol(LayoutMat)), LayoutHeights = rep.int(1, nrow(LayoutMat)),
verbose = TRUE, ...) {
## save default graphical & time parameters and resetting on exit
opar <- par(no.readonly = TRUE)
olctime <- Sys.getlocale(category = "LC_TIME")
suppressWarnings(.TrySetLcTimeEN())
on.exit({
par(opar)
Sys.setlocale(category = "LC_TIME", locale = olctime)
})
OutputsModel <- x
## index time series
if (!is.null(IndPeriod_Plot)) {
if (length(IndPeriod_Plot) <= 1 || any(IndPeriod_Plot <= 0)) {
IndPeriod_Plot <- seq_along(OutputsModel$DatesR)
stop('IndPeriod_Plot must be a positive vector of length > 0')
}
# IndPeriod_Plot <- seq_along(IndPeriod_Plot)
OutputsModel <- .IndexOutputsModel(OutputsModel, IndPeriod_Plot)
Qobs <- Qobs[IndPeriod_Plot]
} else {
IndPeriod_Plot <- seq_along(OutputsModel$DatesR)
}
## ---------- check arguments
if (!inherits(OutputsModel, "GR") & !inherits(OutputsModel, "CemaNeige")) {
stop(paste("OutputsModel not in the correct format for default plotting \n", sep = ""))
return(NULL)
}
BOOL_Dates <- FALSE;
if ("DatesR" %in% names(OutputsModel)) { BOOL_Dates <- TRUE; }
BOOL_Pobs <- FALSE;
if ("Precip" %in% names(OutputsModel)) { BOOL_Pobs <- TRUE; }
BOOL_Qsim <- FALSE;
if ("Qsim" %in% names(OutputsModel)) { BOOL_Qsim <- TRUE; }
BOOL_Qobs <- FALSE;
if (BOOL_Qsim & length(Qobs) == length(OutputsModel$Qsim)) {
if (sum(is.na(Qobs)) != length(Qobs)) {
BOOL_Qobs <- TRUE
}
} else if (inherits(OutputsModel, "GR")) {
warning("Incorrect length of Qobs. Time series of observed flow not drawn.")
}
BOOL_Snow <- FALSE;
if ("CemaNeigeLayers" %in% names(OutputsModel)) { if ("SnowPack" %in% names(OutputsModel$CemaNeigeLayers[[1]])) { BOOL_Snow <- TRUE; } }
BOOL_Psol <- FALSE;
if ("CemaNeigeLayers" %in% names(OutputsModel)) { if ("Psol" %in% names(OutputsModel$CemaNeigeLayers[[1]])) { BOOL_Psol <- TRUE; } }
if ( is.null( which)) { stop("which must be a vector of character \n"); return(NULL); }
if (!is.vector( which)) { stop("which must be a vector of character \n"); return(NULL); }
if (!is.character(which)) { stop("which must be a vector of character \n"); return(NULL); }
if (any(!which %in% c("all", "Precip", 'Temp', "SnowPack", "Flows", "Regime", "CumFreq", "CorQQ"))) {
stop("Incorrect element found in argument which:\nit can only contain 'all', 'Precip', 'Temp', 'SnowPack', 'Flows', 'Regime', 'CumFreq' or 'CorQQ'")
return(NULL)
}
if (all(which %in% c("Temp", "SnowPack")) & !inherits(OutputsModel, "CemaNeige")) {
stop("Incorrect element found in argument which:\nwithout CemaNeige it can only contain 'all', 'Precip', 'Flows', 'Regime', 'CumFreq' or 'CorQQ'")
return(NULL)
}
if (length(unique(which %in% c("Temp", "SnowPack"))) == 2 & !inherits(OutputsModel, "CemaNeige")) {
warning("Incorrect element found in argument which:\nit can only contain 'all', 'Precip', 'Flows', 'Regime', 'CumFreq' or 'CorQQ'\nwithout CemaNeige 'Temp' and 'SnowPack' are not available")
}
stop("'OutputsModel' not in the correct format for default plotting")
}
## check 'OutputsModel'
BOOL_Dates <- FALSE
if ("DatesR" %in% names(OutputsModel)) {
BOOL_Dates <- TRUE
}
BOOL_Pobs <- FALSE
if ("Precip" %in% names(OutputsModel)) {
BOOL_Pobs <- TRUE
}
BOOL_EPobs <- FALSE
if ("PotEvap" %in% names(OutputsModel)) {
BOOL_EPobs <- TRUE
}
BOOL_EAobs <- FALSE
if ("AE" %in% names(OutputsModel)) {
BOOL_EAobs <- TRUE
}
BOOL_Qsim <- FALSE
if ("Qsim" %in% names(OutputsModel)) {
BOOL_Qsim <- TRUE
}
BOOL_Qobs <- FALSE
if (BOOL_Qsim & length(Qobs) == length(OutputsModel$Qsim)) {
if (sum(is.na(Qobs)) != length(Qobs)) {
BOOL_Qobs <- TRUE
}
} else if (inherits(OutputsModel, "GR") & !is.null(Qobs)) {
warning("incorrect length of 'Qobs'. Time series of observed flow not drawn")
}
BOOL_Error <- FALSE
if (BOOL_Qsim & BOOL_Qobs) {
BOOL_Error <- TRUE
}
BOOL_Snow <- FALSE
if ("CemaNeigeLayers" %in% names(OutputsModel)) {
if ("SnowPack" %in% names(OutputsModel$CemaNeigeLayers[[1]])) {
BOOL_Snow <- TRUE
}
}
BOOL_Psol <- FALSE
if ("CemaNeigeLayers" %in% names(OutputsModel)) {
if ("Psol" %in% names(OutputsModel$CemaNeigeLayers[[1]])) {
BOOL_Psol <- TRUE
}
}
## check 'which'
whichNeedQobs <- c("Error", "CorQQ")
whichDashboard <- c("all", "synth", "ts", "perf")
whichAll <- c("Precip", "PotEvap", "ActuEvap", "Temp", "SnowPack", "Flows", "Error", "Regime", "CumFreq", "CorQQ")
whichSynth <- c("Precip" , "Temp", "SnowPack", "Flows" , "Regime", "CumFreq", "CorQQ")
whichTS <- c("Precip", "PotEvap", "Temp", "SnowPack", "Flows" )
whichEvap <- c( "PotEvap", "ActuEvap" )
whichPerf <- c( "Error", "Regime", "CumFreq", "CorQQ")
whichCN <- c( "Temp", "SnowPack" )
warnMsgWhich <- "'which' must be a vector of character"
warnMsgNoQobs <- "the %s plot(s) cannot be drawn if there is no 'Qobs'"
warnMsgWhichCN <- sprintf("incorrect element found in argument 'which':\n\twithout CemaNeige, %s are not available \n\tit can only contain %s",
paste0(shQuote(whichCN), collapse = " and "),
paste0(shQuote(c(whichDashboard, whichAll[!whichAll %in% whichCN])), collapse = ", "))
if (is.null(which)) {
stop(warnMsgWhich)
}
if (!is.vector(which)) {
stop(warnMsgWhich)
}
if (!is.character(which)) {
stop(warnMsgWhich)
}
if (any(!which %in% c(whichDashboard, whichAll))) {
stop(sprintf("incorrect element found in argument 'which': %s\nit can only contain %s",
paste0(shQuote(which[!which %in% c(whichDashboard, whichAll)])),
paste0(shQuote(c(whichDashboard, whichAll)), collapse = ", ")))
}
if (all(which %in% whichCN) & !inherits(OutputsModel, "CemaNeige")) {
stop(warnMsgWhichCN)
}
if (length(unique(which %in% whichCN)) == 2 & !inherits(OutputsModel, "CemaNeige")) {
warning(warnMsgWhichCN)
}
if (all(!which %in% c("all", "synth", "ts", whichCN)) & !inherits(OutputsModel, "GR")) {
stop(sprintf("incorrect element found in argument 'which': \nwith CemaNeige alone, only %s are available",
paste0(shQuote(c("all", "synth", "ts", "Temp", "SnowPack")), collapse = ", ")))
}
if (any(!which %in% c("all", "synth", "ts", whichCN)) & !inherits(OutputsModel, "GR")) {
warning(sprintf("incorrect element found in argument 'which': \nwith CemaNeige alone, only %s are available",
paste0(shQuote(c("all", "synth", "ts", "Temp", "SnowPack")), collapse = ", ")))
}
if ("perf" %in% which) {
which <- c(which, whichPerf)
}
if ("ts" %in% which) {
which <- c(which, whichTS)
}
if ("synth" %in% which) {
which <- c(which, whichSynth)
}
if ("all" %in% which) {
which <- c("Precip", "Temp", "SnowPack", "Flows", "Regime", "CumFreq", "CorQQ")
which <- c(which, whichAll)
}
if (is.null(Qobs) & inherits(OutputsModel, "GR")) {
if (length(which) == 1 & (any(which %in% whichNeedQobs))) {
stop(sprintf(warnMsgNoQobs, shQuote(which)))
}
if (length(which) != 1 & any(which %in% whichNeedQobs)) {
BOOL_CorQQ <- FALSE
BOOL_Error <- FALSE
warning(sprintf(warnMsgNoQobs, paste0(shQuote(whichNeedQobs), collapse = " and ")))
}
}
## check dates
if (!BOOL_Dates) {
stop(paste("OutputsModel must contain at least DatesR to allow plotting \n", sep = "")); return(NULL); }
stop("'OutputsModel' must contain at least 'DatesR' to allow plotting")
}
if (inherits(OutputsModel, "GR") & !BOOL_Qsim) {
stop(paste("OutputsModel must contain at least Qsim to allow plotting \n", sep = "")); return(NULL); }
stop("'OutputsModel' must contain at least 'Qsim' to allow plotting")
}
if (BOOL_Dates) {
MyRollMean1 <- function(x, n) {
return(filter(x, rep(1/n, n), sides = 2)); }
MyRollMean2 <- function(x, n) {
return(filter(c(tail(x, n%/%2), x, x[1:(n%/%2)]), rep(1/n, n), sides = 2)[(n%/%2+1):(length(x)+n%/%2)]); }
# MyRollMean1 <- function(x, n) {
# return(filter(x, rep(1 / n, n), sides = 2))
# }
# MyRollMean2 <- function(x, n) {
# return(filter(c(tail(x, n %/% 2), x, x[1:(n %/% 2)]), rep(1 / n, n), sides = 2)[(n %/% 2 + 1):(length(x) + n %/% 2)])
# }
MyRollMean3 <- function(x, n) {
return(filter(x, filter = rep(1/n, n), sides = 2, circular = TRUE))
}
BOOL_TS <- FALSE;
TimeStep <- difftime(tail(OutputsModel$DatesR, 1), tail(OutputsModel$DatesR, 2), units = "secs")[[1]];
if (inherits(OutputsModel, "hourly" ) & TimeStep %in% ( 60*60)) { BOOL_TS <- TRUE; NameTS <- "hour" ; plotunit <- "[mm/h]"; formatAxis <- "%m/%Y"; }
if (inherits(OutputsModel, "daily" ) & TimeStep %in% ( 24*60*60)) { BOOL_TS <- TRUE; NameTS <- "day" ; plotunit <- "[mm/d]"; formatAxis <- "%m/%Y"; }
if (inherits(OutputsModel, "monthly") & TimeStep %in% (c(28, 29, 30, 31)*24*60*60)) { BOOL_TS <- TRUE; NameTS <- "month"; plotunit <- "[mm/month]"; formatAxis <- "%m/%Y"; }
if (inherits(OutputsModel, "yearly" ) & TimeStep %in% ( c(365, 366)*24*60*60)) { BOOL_TS <- TRUE; NameTS <- "year" ; plotunit <- "[mm/y]"; formatAxis <- "%Y" ; }
if (!BOOL_TS) { stop(paste("the time step of the model inputs could not be found \n", sep = "")); return(NULL); }
}
if (length(IndPeriod_Plot) == 0) { IndPeriod_Plot <- 1:length(OutputsModel$DatesR); }
if (inherits(OutputsModel, "CemaNeige")) { NLayers <- length(OutputsModel$CemaNeigeLayers); }
PsolLayerMean <- NULL; if (BOOL_Psol) {
for(iLayer in 1:NLayers) {
if (iLayer == 1) { PsolLayerMean <- OutputsModel$CemaNeigeLayers[[iLayer]]$Psol/NLayers;
} else { PsolLayerMean <- PsolLayerMean + OutputsModel$CemaNeigeLayers[[iLayer]]$Psol/NLayers; } } }
BOOL_QobsZero <- FALSE; if (BOOL_Qobs) { SelectQobsNotZero <- (round(Qobs[IndPeriod_Plot] , 4) != 0); BOOL_QobsZero <- sum(!SelectQobsNotZero, na.rm = TRUE)>0; }
BOOL_QsimZero <- FALSE; if (BOOL_Qsim) { SelectQsimNotZero <- (round(OutputsModel$Qsim[IndPeriod_Plot], 4) != 0); BOOL_QsimZero <- sum(!SelectQsimNotZero, na.rm = TRUE)>0; }
if (BOOL_QobsZero & verbose) { warning("\t zeroes detected in Qobs: some plots in the log space will not be created using all time-steps \n"); }
if (BOOL_QsimZero & verbose) { warning("\t zeroes detected in Qsim: some plots in the log space will not be created using all time-steps \n"); }
BOOL_FilterZero <- TRUE;
## Plots_choices
BOOLPLOT_Precip <- ( "Precip" %in% which & BOOL_Pobs )
BOOLPLOT_Temp <- ( "Temp" %in% which & BOOL_Snow )
BOOLPLOT_SnowPack <- ( "SnowPack" %in% which & BOOL_Snow )
BOOLPLOT_Flows <- ( "Flows" %in% which & (BOOL_Qsim | BOOL_Qobs) )
BOOLPLOT_Regime <- ( "Regime" %in% which & BOOL_TS & BOOL_Qsim & (NameTS %in% c("hour", "day", "month")) )
BOOLPLOT_CumFreq <- ( "CumFreq" %in% which & (BOOL_Qsim | BOOL_Qobs) & BOOL_FilterZero )
BOOLPLOT_CorQQ <- ( "CorQQ" %in% which & (BOOL_Qsim & BOOL_Qobs) & BOOL_FilterZero )
## Options
return(filter(x, filter = rep(1 / n, n), sides = 2, circular = TRUE))
}
BOOL_TS <- FALSE
if (inherits(OutputsModel, "hourly")) {
BOOL_TS <- TRUE
NameTS <- "hour"
plotunit <- "[mm/h]"
} else if (inherits(OutputsModel, "daily")) {
BOOL_TS <- TRUE
NameTS <- "day"
plotunit <- "[mm/d]"
} else if (inherits(OutputsModel, "monthly")) {
BOOL_TS <- TRUE
NameTS <- "month"
plotunit <- "[mm/month]"
} else if (inherits(OutputsModel, "yearly")) {
BOOL_TS <- TRUE
NameTS <- "year"
plotunit <- "[mm/y]"
}
# if (!BOOL_TS) {
# stop("the time step of the model inputs could not be found")
# }
}
if (inherits(OutputsModel, "CemaNeige")) {
NLayers <- length(OutputsModel$CemaNeigeLayers)
}
PsolLayerMean <- NULL
if (BOOL_Psol) {
for (iLayer in 1:NLayers) {
if (iLayer == 1) {
PsolLayerMean <- OutputsModel$CemaNeigeLayers[[iLayer]]$Psol / NLayers
} else {
PsolLayerMean <- PsolLayerMean + OutputsModel$CemaNeigeLayers[[iLayer]]$Psol / NLayers
}
}
}
BOOL_QobsZero <- FALSE
if (BOOL_Qobs) {
SelectQobsNotZero <- round(Qobs, 4) != 0
BOOL_QobsZero <- sum(!SelectQobsNotZero, na.rm = TRUE) > 0
}
BOOL_QsimZero <- FALSE
if (BOOL_Qsim) {
SelectQsimNotZero <- round(OutputsModel$Qsim, 4) != 0
BOOL_QsimZero <- sum(!SelectQsimNotZero, na.rm = TRUE) > 0
}
if ( BOOL_Qobs & !BOOL_Qsim) {
SelectNotZero <- SelectQobsNotZero
}
if (!BOOL_Qobs & BOOL_Qsim) {
SelectNotZero <- SelectQsimNotZero
}
if ( BOOL_Qobs & BOOL_Qsim) {
SelectNotZero <- SelectQobsNotZero & SelectQsimNotZero
}
if (BOOL_QobsZero & verbose) {
warning("zeroes detected in 'Qobs': some plots in the log space will not be created using all time-steps")
}
if (BOOL_QsimZero & verbose) {
warning("zeroes detected in 'Qsim': some plots in the log space will not be created using all time-steps")
}
BOOL_FilterZero <- TRUE
## ---------- plot
## plot choices
BOOLPLOT_Precip <- "Precip" %in% which & BOOL_Pobs
BOOLPLOT_PotEvap <- "PotEvap" %in% which & BOOL_EPobs
BOOLPLOT_ActuEvap <- "ActuEvap" %in% which & BOOL_EAobs
BOOLPLOT_Temp <- "Temp" %in% which & BOOL_Snow
BOOLPLOT_SnowPack <- "SnowPack" %in% which & BOOL_Snow
BOOLPLOT_Flows <- "Flows" %in% which & (BOOL_Qsim | BOOL_Qobs)
BOOLPLOT_Error <- "Error" %in% which & BOOL_Error
BOOLPLOT_Regime <- "Regime" %in% which & BOOL_Qsim & BOOL_TS & (NameTS %in% c("hour", "day", "month"))
BOOLPLOT_CumFreq <- "CumFreq" %in% which & (BOOL_Qsim | BOOL_Qobs) & BOOL_FilterZero
BOOLPLOT_CorQQ <- "CorQQ" %in% which & (BOOL_Qsim & BOOL_Qobs) & BOOL_FilterZero
## options
BLOC <- TRUE
if (BLOC) {
lwdk <- 1.8
line <- 2.6
bg <- NA
matlayout <- NULL; iPlot <- 0;
Sum1 <- sum(c(BOOLPLOT_Precip, BOOLPLOT_SnowPack, BOOLPLOT_Flows))
Sum2 <- sum(c(BOOLPLOT_Regime, BOOLPLOT_CumFreq, BOOLPLOT_CorQQ))
if (BOOLPLOT_Precip) {
matlayout <- rbind(matlayout, c(iPlot+1, iPlot+1, iPlot+1)); iPlot <- iPlot+1; }
if (BOOLPLOT_Temp) {
matlayout <- rbind(matlayout, c(iPlot+1, iPlot+1, iPlot+1), c(iPlot+1, iPlot+1, iPlot+1)); iPlot <- iPlot+1; }
if (BOOLPLOT_SnowPack) {
matlayout <- rbind(matlayout, c(iPlot+1, iPlot+1, iPlot+1), c(iPlot+1, iPlot+1, iPlot+1)); iPlot <- iPlot+1; }
if (BOOLPLOT_Flows) {
matlayout <- rbind(matlayout, c(iPlot+1, iPlot+1, iPlot+1), c(iPlot+1, iPlot+1, iPlot+1)); iPlot <- iPlot+1; }
if ((Sum1 >= 1 & Sum2 != 0) | (Sum1 == 0 & Sum2 == 3)) {
matlayout <- rbind(matlayout, c(iPlot+1, iPlot+2, iPlot+3), c(iPlot+1, iPlot+2, iPlot+3)); iPlot <- iPlot+3; }
if (Sum1 == 0 & Sum2 == 2) {
matlayout <- rbind(matlayout, c(iPlot+1, iPlot+2)); iPlot <- iPlot+2; }
if (Sum1 == 0 & Sum2 == 1) {
matlayout <- rbind(matlayout, iPlot+1); iPlot <- iPlot+1; }
iPlotMax <- iPlot;
# isRStudio <- Sys.getenv("RSTUDIO") == "1";
# if (!isRStudio) {
# if (Sum1 == 1 & Sum2 == 0) { width = 10; height = 05; }
# if (Sum1 == 1 & Sum2 != 0) { width = 10; height = 07; }
# if (Sum1 == 2 & Sum2 == 0) { width = 10; height = 05; }
# if (Sum1 == 2 & Sum2 != 0) { width = 10; height = 07; }
# if (Sum1 == 3 & Sum2 == 0) { width = 10; height = 07; }
# if (Sum1 == 3 & Sum2 != 0) { width = 10; height = 10; }
# if (Sum1 == 0 & Sum2 == 1) { width = 05; height = 05; }
# if (Sum1 == 0 & Sum2 == 2) { width = 10; height = 04; }
# if (Sum1 == 0 & Sum2 == 3) { width = 10; height = 03; }
# dev.new(width = width, height = height)
# }
opar <- par(no.readonly = TRUE)
on.exit(par(opar))
layout(matlayout)
## Set plot arrangement
if (is.null(LayoutMat)) {
matlayout <- NULL
iPlot <- 0
iHght <- NULL
listBOOLPLOT1 <- c(Precip = BOOLPLOT_Precip,
PotEvap = BOOLPLOT_PotEvap | BOOLPLOT_ActuEvap,
Temp = BOOLPLOT_Temp , SnowPack = BOOLPLOT_SnowPack,
Flows = BOOLPLOT_Flows , Error = BOOLPLOT_Error)
listBOOLPLOT2 <- c(Regime = BOOLPLOT_Regime, CumFreq = BOOLPLOT_CumFreq,
CorQQ = BOOLPLOT_CorQQ)
Sum1 <- sum(listBOOLPLOT1)
Sum2 <- sum(listBOOLPLOT2)
Xaxis <- 1:length(IndPeriod_Plot);
if (BOOL_Dates) {
if (NameTS %in% c("hour", "day", "month")) {
Seq1 <- which(OutputsModel$DatesR[IndPeriod_Plot]$mday == 1 & OutputsModel$DatesR[IndPeriod_Plot]$mon %in% c(0, 3, 6, 9))
Seq2 <- which(OutputsModel$DatesR[IndPeriod_Plot]$mday == 1 & OutputsModel$DatesR[IndPeriod_Plot]$mon == 0)
Labels2 <- format(OutputsModel$DatesR[IndPeriod_Plot], format = formatAxis)[Seq2]
for (k in seq_len(Sum1)) {
matlayout <- rbind(matlayout, iPlot + c(1, 1, 1))
iPlot <- iPlot + 1
iHght <- c(iHght, 0.7)
}
## Flows plot is higher than the other TS
listBOOLPLOT1 <- listBOOLPLOT1[listBOOLPLOT1]
listBOOLPLOTF <- (names(listBOOLPLOT1) == "Flows") * BOOLPLOT_Flows
iHght <- iHght + listBOOLPLOTF * listBOOLPLOT1 * 0.3
if (Sum2 >= 1) {
iHght <- c(iHght, 1.0)
}
if (NameTS %in% c("year")) {
Seq1 <- 1:length(OutputsModel$DatesR[IndPeriod_Plot])
Seq2 <- 1:length(OutputsModel$DatesR[IndPeriod_Plot])
Labels2 <- format(OutputsModel$DatesR[IndPeriod_Plot], format = formatAxis)[Seq2]
if ((Sum1 >= 1 & Sum2 != 0) | (Sum1 == 0 & Sum2 == 3)) {
matlayout <- rbind(matlayout, iPlot + c(1, 2, 3))
iPlot <- iPlot + 3
}
if (Sum1 == 0 & Sum2 == 2) {
matlayout <- rbind(matlayout, iPlot + c(1, 2))
iPlot <- iPlot + 2
}
if (Sum1 == 0 & Sum2 == 1) {
matlayout <- rbind(matlayout, iPlot + 1)
iPlot <- iPlot + 1
}
iPlotMax <- iPlot
LayoutWidths <- rep.int(1, ncol(matlayout))
LayoutHeights <- iHght #rep.int(1, nrow(matlayout))
}
if (!is.null(LayoutMat)) {
matlayout <- LayoutMat
}
layout(matlayout, widths = LayoutWidths, heights = LayoutHeights)
Xaxis <- as.POSIXct(OutputsModel$DatesR)
if (!is.null(BasinArea)) {
Factor_MMH_M3S <- 60 * 60
Factor_MMD_M3S <- 60 * 60 * 24
Factor_MMM_M3S <- 60 * 60 * 24 * 365.25 / 12
Factor_MMY_M3S <- 60 * 60 * 24 * 365.25
if (NameTS == "hour" ) Factor_UNIT_M3S <- Factor_MMH_M3S
if (NameTS == "day" ) Factor_UNIT_M3S <- Factor_MMD_M3S
if (NameTS == "month") Factor_UNIT_M3S <- Factor_MMM_M3S
if (NameTS == "year" ) Factor_UNIT_M3S <- Factor_MMY_M3S
Factor_UNIT_M3S <- switch(NameTS,
hour = 60 * 60,
day = 60 * 60 * 24,
month = 60 * 60 * 24 * 365.25 / 12,
year = 60 * 60 * 24 * 365.25)
Factor_UNIT_M3S <- BasinArea / (Factor_UNIT_M3S / 1000)
}
}
kPlot <- 0
## vector of Q values for the y-axis when it is expressed in
Factor <- ifelse(!is.null(BasinArea), Factor_UNIT_M3S, 1)
seqDATA0 <- c(0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000)
......@@ -180,323 +342,514 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
if (BOOLPLOT_Precip) {
kPlot <- kPlot + 1
mar <- c(3, 5, 1, 5)
par(new = FALSE, mar = mar, las = 0)
ylim1 <- range(OutputsModel$Precip[IndPeriod_Plot], na.rm = TRUE)
par(new = FALSE, mar = mar)
ylim1 <- range(OutputsModel$Precip, na.rm = TRUE)
ylim2 <- ylim1 * c(1.0, 1.1)
ylim2 <- rev(ylim2)
lwdP <- lwd * 0.7
if (NameTS %in% c("month", "year")) {
lwdP <- lwd * 2
}
plot(Xaxis, OutputsModel$Precip[IndPeriod_Plot],
plot(Xaxis, OutputsModel$Precip,
type = "h", xaxt = "n", yaxt = "n", yaxs = "i", ylim = ylim2,
col = "royalblue", lwd = lwdP * lwdk, lend = 1,
xlab = "", ylab = "", ...)
axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
par(las = 0)
mtext(side = 2, paste0("precip. ", plotunit), cex = cex.lab, adj = 1, line = line)
par(las = 0)
mtext(side = 2, paste("precip.", plotunit), cex = cex.lab, adj = 1, line = line)
if (BOOL_Psol) {
legend("bottomright", legend = c("solid","liquid"),
col = c("lightblue", "royalblue"), lty = c(1, 1), lwd = c(lwd, lwd),
bty = "o", bg = bg, box.col = bg, cex = cex.leg)
par(new = TRUE)
plot(Xaxis, PsolLayerMean[IndPeriod_Plot],
type = "h", xaxt = "n", yaxt = "n", yaxs = "i", ylim = ylim2,
col = "lightblue", lwd = lwdP * lwdk, lend = 1,
xlab = "", ylab = "", ...)
points(Xaxis, PsolLayerMean,
type = "h", xaxt = "n", yaxt = "n", yaxs = "i", ylim = ylim2,
col = "lightblue", lwd = lwdP * lwdk, lend = 1,
xlab = "", ylab = "", ...)
}
if (BOOL_Dates) {
axis(side = 1, at = Seq1, labels = FALSE, cex.axis = cex.axis, ...)
axis(side = 1, at = Seq2, labels = Labels2, lwd.ticks = 1.5, cex.axis = cex.axis, ...)
AxisTS(OutputsModel)
box()
}
## PotEvap
if (BOOLPLOT_PotEvap) {
kPlot <- kPlot + 1
mar <- c(3, 5, 1, 5)
par(new = FALSE, mar = mar)
if (!BOOLPLOT_ActuEvap) {
ylim1 <- range(OutputsModel$PotEvap, na.rm = TRUE)
xlabE <- "pot. evap."
} else {
axis(side = 1, at = pretty(Xaxis), labels = pretty(Xaxis), cex.axis = cex.axis, ...)
ylim1 <- range(c(OutputsModel$PotEvap,
OutputsModel$AE),
na.rm = TRUE)
xlabE <- "evap."
}
ylim2 <- ylim1 #* c(1.0, 1.1)
plot(Xaxis, OutputsModel$PotEvap,
type = "l", xaxt = "n", yaxt = "n", ylim = ylim2,
col = "green3", lwd = lwd * lwdk,
xlab = "", ylab = "", ...)
if (BOOLPLOT_ActuEvap) {
lines(Xaxis, OutputsModel$AE,
type = "l", xaxt = "n", yaxt = "n", ylim = ylim2,
col = "green4", lwd = lwd * lwdk, lty = 3)
legend("topright", legend = c("pot.", "actu."), col = c("green3", "green4"),
lty = c(1, 3), lwd = c(lwd * 1.0, lwd * 0.8),
bty = "o", bg = bg, box.col = bg, cex = cex.leg)
}
axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
mtext(side = 2, paste(xlabE, plotunit), cex = cex.lab, line = line)
AxisTS(OutputsModel)
box()
}
## ActuEvap
if (BOOLPLOT_ActuEvap & !BOOLPLOT_PotEvap) {
kPlot <- kPlot + 1
mar <- c(3, 5, 1, 5)
par(new = FALSE, mar = mar)
ylim1 <- range(OutputsModel$AE, na.rm = TRUE)
ylim2 <- ylim1 #* c(1.0, 1.1)
plot(Xaxis, OutputsModel$AE,
type = "l", xaxt = "n", yaxt = "n", ylim = ylim2,
col = "green4", lwd = lwd * lwdk,
xlab = "", ylab = "", ...)
axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
mtext(side = 2, paste("actu. evap.", plotunit), cex = cex.lab, line = line)
AxisTS(OutputsModel)
box()
}
## Temp
if (BOOLPLOT_Temp) {
kPlot <- kPlot+1; mar <- c(3, 5, 1, 5);
par(new = FALSE, mar = mar, las = 0)
kPlot <- kPlot + 1
mar <- c(3, 5, 1, 5)
par(new = FALSE, mar = mar)
ylim1 <- c(+99999, -99999)
for(iLayer in 1:NLayers) {
ylim1[1] <- min(ylim1[1], OutputsModel$CemaNeigeLayers[[iLayer]]$Temp);
ylim1[2] <- max(ylim1[2], OutputsModel$CemaNeigeLayers[[iLayer]]$Temp);
if (iLayer == 1) { SnowPackLayerMean <- OutputsModel$CemaNeigeLayers[[iLayer]]$Temp/NLayers;
} else { SnowPackLayerMean <- SnowPackLayerMean + OutputsModel$CemaNeigeLayers[[iLayer]]$Temp/NLayers; }
}
plot(SnowPackLayerMean[IndPeriod_Plot], type = "n", ylim = ylim1, xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
for(iLayer in 1:NLayers) { lines(OutputsModel$CemaNeigeLayers[[iLayer]]$Temp[IndPeriod_Plot], lty = 3, col = "orchid", lwd = lwd * lwdk * 0.8); }
for (iLayer in 1:NLayers) {
ylim1[1] <- min(ylim1[1], OutputsModel$CemaNeigeLayers[[iLayer]]$Temp)
ylim1[2] <- max(ylim1[2], OutputsModel$CemaNeigeLayers[[iLayer]]$Temp)
if (iLayer == 1) {
SnowPackLayerMean <- OutputsModel$CemaNeigeLayers[[iLayer]]$Temp/NLayers
} else {
SnowPackLayerMean <- SnowPackLayerMean + OutputsModel$CemaNeigeLayers[[iLayer]]$Temp/NLayers
}
}
plot(Xaxis, SnowPackLayerMean, type = "n", ylim = ylim1, xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
for (iLayer in 1:NLayers) {
lines(Xaxis, OutputsModel$CemaNeigeLayers[[iLayer]]$Temp, lty = 3, col = "orchid", lwd = lwd * lwdk * 0.8)
}
abline(h = 0, col = "grey", lty = 2)
lines(SnowPackLayerMean[IndPeriod_Plot], type = "l", lwd = lwd * lwdk *1.0, col = "darkorchid4")
lines(Xaxis, SnowPackLayerMean, type = "l", lwd = lwd * lwdk * 1.0, col = "darkorchid4")
axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
par(las = 0); mtext(side = 2, expression(paste("temp. [", degree, "C]", sep = "")), padj = 0.2, line = line, cex = cex.lab); par(las = 0);
legend("topright", legend = c("mean", "layers"), col = c("darkorchid4", "orchid"), lty = c(1, 3), lwd = c(lwd*1.0, lwd*0.8), bty = "o", bg = bg, box.col = bg, cex = cex.leg)
mtext(side = 2, expression(paste("temp. [", degree, "C]"), sep = ""), padj = 0.2, line = line, cex = cex.lab)
legend("topright", legend = c("mean", "layers"), col = c("darkorchid4", "orchid"),
lty = c(1, 3), lwd = c(lwd * 1.0, lwd * 0.8),
bty = "o", bg = bg, box.col = bg, cex = cex.leg)
AxisTS(OutputsModel)
box()
if (BOOL_Dates) {
axis(side = 1, at = Seq1, labels = FALSE, cex.axis = cex.axis, ...);
axis(side = 1, at = Seq2, labels = Labels2, lwd.ticks = 1.5, cex.axis = cex.axis, ...);
} else { axis(side = 1, at = pretty(Xaxis), labels = pretty(Xaxis), cex.axis = cex.axis, ...); }
}
## SnowPack
if (BOOLPLOT_SnowPack) {
kPlot <- kPlot+1; mar <- c(3, 5, 1, 5);
par(new = FALSE, mar = mar, las = 0)
kPlot <- kPlot + 1
mar <- c(3, 5, 1, 5)
par(new = FALSE, mar = mar)
ylim1 <- c(+99999, -99999)
for(iLayer in 1:NLayers) {
ylim1[1] <- min(ylim1[1], OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack);
ylim1[2] <- max(ylim1[2], OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack);
if (iLayer == 1) { SnowPackLayerMean <- OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack/NLayers;
} else { SnowPackLayerMean <- SnowPackLayerMean + OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack/NLayers; }
}
plot(SnowPackLayerMean[IndPeriod_Plot], type = "l", ylim = ylim1, lwd = lwd * lwdk *1.2, col = "royalblue", xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
for(iLayer in 1:NLayers) { lines(OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack[IndPeriod_Plot], lty = 3, col = "royalblue", lwd = lwd * lwdk *0.8); }
for (iLayer in 1:NLayers) {
ylim1[1] <- min(ylim1[1], OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack)
ylim1[2] <- max(ylim1[2], OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack)
if (iLayer == 1) {
SnowPackLayerMean <- OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack/NLayers
} else {
SnowPackLayerMean <- SnowPackLayerMean + OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack/NLayers
}
}
plot(Xaxis, SnowPackLayerMean, type = "l", ylim = ylim1, lwd = lwd * lwdk * 1.2, col = "royalblue", xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
for (iLayer in 1:NLayers) {
lines(Xaxis, OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack, lty = 3, col = "royalblue", lwd = lwd * lwdk * 0.8)
}
axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
par(las = 0); mtext(side = 2, paste0("snow pack ", "[mm]"), line = line, cex = cex.lab); par(las = 0);
legend("topright", legend = c("mean", "layers"), col = c("royalblue", "royalblue"), lty = c(1, 3), lwd = c(lwd*1.2, lwd*0.8), bty = "o", bg = bg, box.col = bg, cex = cex.leg)
mtext(side = 2, paste("snow pack", "[mm]"), line = line, cex = cex.lab)
legend("topright", legend = c("mean", "layers"), col = c("royalblue", "royalblue"),
lty = c(1, 3), lwd = c(lwd * 1.2, lwd * 0.8),
bty = "o", bg = bg, box.col = bg, cex = cex.leg)
AxisTS(OutputsModel)
box()
if (BOOL_Dates) {
axis(side = 1, at = Seq1, labels = FALSE, cex.axis = cex.axis, ...);
axis(side = 1, at = Seq2, labels = Labels2, lwd.ticks = 1.5, cex.axis = cex.axis, ...);
} else { axis(side = 1, at = pretty(Xaxis), labels = pretty(Xaxis), cex.axis = cex.axis, ...); }
}
## Flows
if (BOOLPLOT_Flows & log_scale) {
kPlot <- kPlot+1; mar <- c(3, 5, 1, 5);
par(new = FALSE, mar = mar, las = 0)
DATA2 <- Qobs
DATA2[!SelectQobsNotZero] <- mean(Qobs, na.rm = TRUE) / 10000
DATA2 <- log(DATA2)
kPlot <- kPlot + 1
mar <- c(3, 5, 1, 5)
par(new = FALSE, mar = mar)
if (BOOL_Qobs) {
DATA2 <- Qobs
DATA2[!SelectQobsNotZero] <- mean(Qobs, na.rm = TRUE) / 10000
DATA2 <- log(DATA2)
}
DATA3 <- OutputsModel$Qsim
DATA3[!SelectQsimNotZero] <- mean(OutputsModel$Qsim, na.rm = TRUE) / 10000
DATA3 <- log(DATA3)
ylim1 <- range(DATA3[IndPeriod_Plot], na.rm = TRUE);
if (BOOL_Qobs) { ylim1 <- range(c(ylim1, DATA2[IndPeriod_Plot]), na.rm = TRUE); }
ylim2 <- c(ylim1[1], 1.2*ylim1[2]);
plot(Xaxis, rep(NA, length(Xaxis)), type = "n", ylim = ylim2, xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...);
txtleg <- NULL; colleg <- NULL;
if (BOOL_Qobs) { lines(Xaxis, DATA2[IndPeriod_Plot], lwd = lwd * lwdk , lty = 1, col = par("fg")); txtleg <- c(txtleg, "observed"); colleg <- c(colleg, par("fg")); }
if (BOOL_Qsim) { lines(Xaxis, DATA3[IndPeriod_Plot], lwd = lwd * lwdk , lty = 1, col = "orangered"); txtleg <- c(txtleg, "simulated"); colleg <- c(colleg, "orangered"); }
ylim1 <- range(DATA3, na.rm = TRUE)
if (BOOL_Qobs) {
ylim1 <- range(c(ylim1, DATA2), na.rm = TRUE)
}
ylim2 <- c(ylim1[1], 1.1 * ylim1[2])
plot(Xaxis, rep(NA, length(Xaxis)), type = "n", ylim = ylim2, xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
txtleg <- NULL
colleg <- NULL
if (BOOL_Qobs) {
lines(Xaxis, DATA2, lwd = lwd * lwdk, lty = 1, col = par("fg"))
txtleg <- c(txtleg, "observed")
colleg <- c(colleg, par("fg"))
}
if (BOOL_Qsim) {
lines(Xaxis, DATA3, lwd = lwd * lwdk, lty = 1, col = "orangered")
txtleg <- c(txtleg, "simulated")
colleg <- c(colleg, "orangered")
}
axis(side = 2, at = seqDATA1, labels = seqDATA2, cex.axis = cex.axis, ...)
par(las = 0); mtext(side = 2, paste0("flow ", plotunit), line = line, cex = cex.lab); par(las = 0);
mtext(side = 2, paste("flow", plotunit), line = line, cex = cex.lab)
if (!is.null(BasinArea)) {
Factor <- Factor_UNIT_M3S;
Factor <- Factor_UNIT_M3S
axis(side = 4, at = seqDATA1ba, labels = seqDATA2ba, cex.axis = cex.axis, ...)
par(las = 0); mtext(side = 4, paste0("flow ", "[m3/s]"), line = line, cex = cex.lab); par(las = 0); }
if (BOOL_Dates) {
axis(side = 1, at = Seq1, labels = FALSE, cex.axis = cex.axis, ...);
axis(side = 1, at = Seq2, labels = Labels2, lwd.ticks = 1.5, cex.axis = cex.axis, ...);
} else { axis(side = 1, at = pretty(Xaxis), labels = pretty(Xaxis), cex.axis = cex.axis, ...); }
legend("topright", txtleg, col = colleg, lty = 1, lwd = lwd * lwdk , bty = "o", bg = bg, box.col = bg, cex = cex.leg)
mtext(side = 4, paste("flow", "[m3/s]"), line = line, cex = cex.lab)
}
AxisTS(OutputsModel)
legend("topright", txtleg, col = colleg, lty = 1, lwd = lwd * lwdk, bty = "o", bg = bg, box.col = bg, cex = cex.leg)
legend("bottomright", "log scale", lty = 1, col = NA, bty = "o", bg = bg, box.col = bg, cex = cex.leg)
box()
}
if (BOOLPLOT_Flows & !log_scale) {
kPlot <- kPlot+1; mar <- c(3, 5, 1, 5);
par(new = FALSE, mar = mar, las = 0)
ylim1 <- range(OutputsModel$Qsim[IndPeriod_Plot], na.rm = TRUE);
if (BOOL_Qobs) { ylim1 <- range(c(ylim1, Qobs[IndPeriod_Plot]), na.rm = TRUE); }
ylim2 <- c(ylim1[1], 1.2*ylim1[2]);
plot(Xaxis, rep(NA, length(Xaxis)), type = "n", ylim = ylim2, xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...);
txtleg <- NULL; colleg <- NULL;
if (BOOL_Qobs) { lines(Xaxis, Qobs[IndPeriod_Plot], lwd = lwd * lwdk , lty = 1, col = par("fg")); txtleg <- c(txtleg, "observed"); colleg <- c(colleg, par("fg")); }
if (BOOL_Qsim) { lines(Xaxis, OutputsModel$Qsim[IndPeriod_Plot], lwd = lwd * lwdk , lty = 1, col = "orangered"); txtleg <- c(txtleg, "simulated"); colleg <- c(colleg, "orangered"); }
kPlot <- kPlot + 1
mar <- c(3, 5, 1, 5)
par(new = FALSE, mar = mar)
ylim1 <- range(OutputsModel$Qsim, na.rm = TRUE)
if (BOOL_Qobs) {
ylim1 <- range(c(ylim1, Qobs), na.rm = TRUE)
}
ylim2 <- c(ylim1[1], 1.1 * ylim1[2])
plot(Xaxis, rep(NA, length(Xaxis)), type = "n", ylim = ylim2, xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
txtleg <- NULL
colleg <- NULL
if (BOOL_Qobs) {
lines(Xaxis, Qobs, lwd = lwd * lwdk, lty = 1, col = par("fg"))
txtleg <- c(txtleg, "observed")
colleg <- c(colleg, par("fg"))
}
if (BOOL_Qsim) {
lines(Xaxis, OutputsModel$Qsim, lwd = lwd * lwdk, lty = 1, col = "orangered")
txtleg <- c(txtleg, "simulated")
colleg <- c(colleg, "orangered")
}
axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
par(las = 0); mtext(side = 2, paste("flow", plotunit, sep = " "), line = line, cex = cex.lab); par(las = 0);
mtext(side = 2, paste("flow", plotunit), line = line, cex = cex.lab)
if (!is.null(BasinArea)) {
Factor <- Factor_UNIT_M3S;
axis(side = 4, at = pretty(ylim1*Factor)/Factor, labels = pretty(ylim1*Factor), cex.axis = cex.axis, ...);
par(las = 0); mtext(side = 4, paste("flow", "[m3/s]", sep = " "), line = line, cex = cex.lab); par(las = 0); }
if (BOOL_Dates) {
axis(side = 1, at = Seq1, labels = FALSE, cex.axis = cex.axis, ...);
axis(side = 1, at = Seq2, labels = Labels2, lwd.ticks = 1.5, cex.axis = cex.axis, ...);
} else { axis(side = 1, at = pretty(Xaxis), labels = pretty(Xaxis), cex.axis = cex.axis, ...); }
legend("topright", txtleg, col = colleg, lty = 1, lwd = lwd * lwdk , bty = "o", bg = bg, box.col = bg, cex = cex.leg)
Factor <- Factor_UNIT_M3S
axis(side = 4, at = pretty(ylim1 * Factor)/Factor, labels = pretty(ylim1 * Factor), cex.axis = cex.axis, ...)
mtext(side = 4, paste("flow", "[m3/s]"), line = line, cex = cex.lab)
}
AxisTS(OutputsModel)
legend("topright", txtleg, col = colleg, lty = 1, lwd = lwd * lwdk, bty = "o", bg = bg, box.col = bg, cex = cex.leg)
box()
}
## Error
if (BOOLPLOT_Error) {
kPlot <- kPlot + 1
mar <- c(3, 5, 1, 5)
if (log_scale) {
errorQ <- log(OutputsModel$Qsim) - log(Qobs)
} else {
errorQ <- OutputsModel$Qsim - Qobs
}
par(new = FALSE, mar = mar)
ylim1 <- range(errorQ[SelectNotZero], na.rm = TRUE)
plot(Xaxis, errorQ,
type = "l", xaxt = "n", yaxt = "n", ylim = ylim1,
col = par("fg"), lwd = lwd * lwdk,
xlab = "", ylab = "",
panel.first = abline(h = 0, col = "royalblue"), ...)
axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
mtext(side = 2, paste("flow error", plotunit), cex = cex.lab, line = line)
if (!is.null(BasinArea)) {
Factor <- Factor_UNIT_M3S
axis(side = 4, at = pretty(ylim1 * Factor)/Factor, labels = pretty(ylim1 * Factor), cex.axis = cex.axis, ...)
mtext(side = 4, paste("flow error", "[m3/s]"), line = line, cex = cex.lab)
}
AxisTS(OutputsModel)
if (log_scale) {
legend("bottomright", "log scale", lty = 1, col = NA, bty = "o", bg = bg, box.col = bg, cex = cex.leg)
}
box()
}
## Regime
if (BOOLPLOT_Regime) {
kPlot <- kPlot+1; mar <- c(6, 5, 1, 5); plotunitregime <- "[mm/month]"
par(new = FALSE, mar = mar, las = 0)
kPlot <- kPlot + 1
mar <- c(6, 5, 1, 5)
plotunitregime <- "[mm/month]"
par(new = FALSE, mar = mar)
## Empty plot
if ((NameTS == "hour" & length(IndPeriod_Plot) < 697) |
(NameTS == "day" & length(IndPeriod_Plot) < 30) |
(NameTS == "month" & length(IndPeriod_Plot) < 2) |
(NameTS == "year" & length(IndPeriod_Plot) < 2)) {
plot(0, 0, type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", main = "", ...)
par(las = 0); mtext(side = 1, text = "", line = line, cex = cex.lab); par(las = 0)
text(0, 0, labels = "NO ENOUGH VALUES", col = "grey40"); par(las = 0)
txtlab <- "flow regime"
plot(0, 0, type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ...)
mtext(side = 1, text = "", line = line, cex = cex.lab)
text(0, 0, labels = "NO ENOUGH VALUES", col = "grey40")
txtlab <- "flow"
if (BOOL_Pobs) {
txtlab <- "precip. & flow regime"
txtlab <- "precip. & flow"
}
par(las = 0); mtext(side = 2, paste(txtlab, plotunitregime), line = line, cex = cex.lab); par(las = 0)
mtext(side = 2, paste(txtlab, plotunitregime), line = line, cex = cex.lab)
} else {
## Data_formating_as_table
DataModel <- as.data.frame(matrix(as.numeric(NA), nrow = length(IndPeriod_Plot), ncol = 5));
DataModel[, 1] <- as.numeric(format(OutputsModel$DatesR[IndPeriod_Plot], format = "%Y%m%d%H"));
if (BOOL_Pobs) { DataModel[, 2] <- OutputsModel$Precip[IndPeriod_Plot]; }
if (BOOL_Psol) { DataModel[, 3] <- PsolLayerMean[IndPeriod_Plot]; }
if (BOOL_Qobs) { DataModel[, 4] <- Qobs[IndPeriod_Plot]; }
if (BOOL_Qsim) { DataModel[, 5] <- OutputsModel$Qsim[IndPeriod_Plot]; }
colnames(DataModel) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim");
TxtDatesDataModel <- formatC(DataModel$Dates, format = "d", width = 8, flag = "0");
DataModel <- as.data.frame(matrix(as.numeric(NA), nrow = length(IndPeriod_Plot), ncol = 5))
DataModel[, 1] <- as.numeric(format(OutputsModel$DatesR, format = "%Y%m%d%H"))
if (BOOL_Pobs) {
DataModel[, 2] <- OutputsModel$Precip
}
if (BOOL_Psol) {
DataModel[, 3] <- PsolLayerMean
}
if (BOOL_Qobs) {
DataModel[, 4] <- Qobs
}
if (BOOL_Qsim) {
DataModel[, 5] <- OutputsModel$Qsim
}
colnames(DataModel) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim")
TxtDatesDataModel <- formatC(DataModel$Dates, format = "d", width = 8, flag = "0")
## Building_of_daily_time_series_if_needed
if (NameTS == "month") { DataDaily <- NULL; }
if (NameTS == "day" ) { DataDaily <- DataModel; }
if (NameTS == "hour" ) { DataDaily <- as.data.frame(aggregate(DataModel[, 2:5], by = list(as.numeric(substr(TxtDatesDataModel, 1, 8))), FUN = sum, na.rm = T)); }
if (NameTS == "month") {
DataDaily <- NULL
}
if (NameTS == "day") {
DataDaily <- DataModel
}
if (NameTS == "hour" ) {
DataDaily <- as.data.frame(aggregate(DataModel[, 2:5], by = list(as.numeric(substr(TxtDatesDataModel, 1, 8))), FUN = sum, na.rm = TRUE))
}
if (NameTS %in% c("hour", "day")) {
colnames(DataDaily) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim");
TxtDatesDataDaily <- formatC(DataDaily$Dates, format = "d", width = 8, flag = "0"); }
colnames(DataDaily) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim")
TxtDatesDataDaily <- formatC(DataDaily$Dates, format = "d", width = 8, flag = "0")
}
## Building_of_monthly_time_series_if_needed
if (NameTS == "month") { DataMonthly <- DataModel; }
if (NameTS == "day" ) { DataMonthly <- as.data.frame(aggregate(DataDaily[, 2:5], by = list(as.numeric(substr(TxtDatesDataDaily, 1, 6))), FUN = sum, na.rm = T)); }
if (NameTS == "hour" ) { DataMonthly <- as.data.frame(aggregate(DataDaily[, 2:5], by = list(as.numeric(substr(TxtDatesDataDaily, 1, 6))), FUN = sum, na.rm = T)); }
colnames(DataMonthly) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim");
TxtDatesDataMonthly <- formatC(DataMonthly$Dates, format = "d", width = 6, flag = "0");
if (NameTS == "month") {
DataMonthly <- DataModel
}
if (NameTS == "day" ) {
DataMonthly <- as.data.frame(aggregate(DataDaily[, 2:5], by = list(as.numeric(substr(TxtDatesDataDaily, 1, 6))), FUN = sum, na.rm = TRUE))
}
if (NameTS == "hour" ) {
DataMonthly <- as.data.frame(aggregate(DataDaily[, 2:5], by = list(as.numeric(substr(TxtDatesDataDaily, 1, 6))), FUN = sum, na.rm = TRUE))
}
colnames(DataMonthly) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim")
TxtDatesDataMonthly <- formatC(DataMonthly$Dates, format = "d", width = 6, flag = "0")
## Computation_of_interannual_mean_series
if (!is.null(DataDaily)) {
SeqY <- data.frame(Dates = as.numeric(format(seq(as.Date("1970-01-01", tz = "UTC"),
as.Date("1970-12-31", tz = "UTC"), "day"),
format = "%m%d")))
DataDailyInterAn <- as.data.frame(aggregate(DataDaily[, 2:5], by = list(as.numeric(substr(TxtDatesDataDaily , 5, 8))), FUN = mean, na.rm = T));
DataDailyInterAn <- as.data.frame(aggregate(DataDaily[, 2:5], by = list(as.numeric(substr(TxtDatesDataDaily, 5, 8))), FUN = mean, na.rm = TRUE))
colnames(DataDailyInterAn) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim")
DataDailyInterAn <- merge(SeqY, DataDailyInterAn, by = "Dates", all.x = TRUE, all.y = FALSE)
DataDailyInterAn <- merge(SeqY, DataDailyInterAn, by = "Dates", all.x = TRUE, all.y = FALSE)
}
if (!is.null(DataMonthly)) {
SeqM <- data.frame(Dates = 1:12)
DataMonthlyInterAn <- as.data.frame(aggregate(DataMonthly[, 2:5], by = list(as.numeric(substr(TxtDatesDataMonthly, 5, 6))), FUN = mean, na.rm = T));
DataMonthlyInterAn <- as.data.frame(aggregate(DataMonthly[, 2:5], by = list(as.numeric(substr(TxtDatesDataMonthly, 5, 6))), FUN = mean, na.rm = TRUE))
colnames(DataMonthlyInterAn) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim")
DataMonthlyInterAn <- merge(SeqM, DataMonthlyInterAn, by = "Dates", all.x = TRUE, all.y = FALSE)
DataMonthlyInterAn <- merge(SeqM, DataMonthlyInterAn, by = "Dates", all.x = TRUE, all.y = FALSE)
}
## Smoothing_of_daily_series_and_scale_conversion_to_make_them_become_a_monthly_regime
if (!is.null(DataDaily)) {
## Smoothing
NDaysWindow <- 30;
DataDailyInterAn <- as.data.frame(cbind(DataDailyInterAn$Dates,
MyRollMean3(DataDailyInterAn$Precip, NDaysWindow), MyRollMean3(DataDailyInterAn$Psol, NDaysWindow),
MyRollMean3(DataDailyInterAn$Qobs , NDaysWindow), MyRollMean3(DataDailyInterAn$Qsim, NDaysWindow)));
colnames(DataDailyInterAn) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim");
NDaysWindow <- 30
DataDailyInterAn <- as.data.frame(cbind(DataDailyInterAn$Dates,
MyRollMean3(DataDailyInterAn$Precip, NDaysWindow), MyRollMean3(DataDailyInterAn$Psol, NDaysWindow),
MyRollMean3(DataDailyInterAn$Qobs , NDaysWindow), MyRollMean3(DataDailyInterAn$Qsim, NDaysWindow)))
colnames(DataDailyInterAn) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim")
## Scale_conversion_to_make_them_become_a_monthly_regime
if (plotunitregime != "[mm/month]") { stop(paste("incorrect unit for regime plot \n", sep = "")); return(NULL); }
DataDailyInterAn <- as.data.frame(cbind(DataDailyInterAn[1], DataDailyInterAn[2:5]*30));
if (plotunitregime != "[mm/month]") {
stop("incorrect unit for regime plot")
}
DataDailyInterAn <- as.data.frame(cbind(DataDailyInterAn[1], DataDailyInterAn[2:5]*30))
}
## Plot_preparation
DataPlotP <- DataMonthlyInterAn;
DataPlotP <- DataMonthlyInterAn
if (!is.null(DataDaily)) {
DataPlotQ <- DataDailyInterAn;
SeqX1 <- c( 1, 32, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336, 366);
SeqX2 <- c( 15, 46, 75, 106, 136, 167, 197, 228, 259, 289, 320, 350);
labX <- "30-days rolling mean";
DataPlotQ <- DataDailyInterAn
SeqX1 <- c( 1, 32, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336, 366)
SeqX2 <- c( 15, 46, 75, 106, 136, 167, 197, 228, 259, 289, 320, 350)
labX <- "30-days rolling mean"
} else {
DataPlotQ <- DataMonthlyInterAn;
SeqX1 <- seq(from = 0.5, to = 12.5, by = 1);
SeqX2 <- seq(from = 1 , to = 12 , by = 1);
labX <- "";
DataPlotQ <- DataMonthlyInterAn
SeqX1 <- seq(from = 0.5, to = 12.5, by = 1)
SeqX2 <- seq(from = 1 , to = 12 , by = 1)
labX <- ""
}
xLabels1 <- rep("", 13);
xLabels1 <- rep("", 13)
xLabels2 <- month.abb
ylimQ <- range(c(DataPlotQ$Qobs, DataPlotQ$Qsim), na.rm = TRUE);
if (BOOL_Pobs) { ylimP <- c(max(DataPlotP$Precip, na.rm = TRUE), 0); }
txtleg <- NULL; colleg <- NULL; lwdleg <- NULL; lwdP = 10;
ylimQ <- range(c(DataPlotQ$Qobs, DataPlotQ$Qsim), na.rm = TRUE)
if (BOOL_Pobs) {
ylimP <- c(max(DataPlotP$Precip, na.rm = TRUE), 0)
}
txtleg <- NULL
colleg <- NULL
lwdleg <- NULL
lwdP = 10
## Plot_forcings
if (BOOL_Pobs) {
plot(SeqX2[DataMonthlyInterAn$Dates], DataPlotP$Precip, type = "h",
xlim = range(SeqX1), ylim = c(3*ylimP[1], ylimP[2]), lwd = lwdP, lend = 1, lty = 1, col = "royalblue",
xlim = range(SeqX1), ylim = c(3 * ylimP[1], ylimP[2]), lwd = lwdP, lend = 1, lty = 1, col = "royalblue",
xlab = "", ylab = "", xaxt = "n", yaxt = "n", yaxs = "i", bty = "n", ...)
txtleg <- c(txtleg, "Ptot" ); colleg <- c(colleg, "royalblue"); lwdleg <- c(lwdleg, lwdP/3);
axis(side = 2, at = pretty(0.8*ylimP, n = 3), labels = pretty(0.8*ylimP, n = 3), col.axis = "royalblue", col.ticks = "royalblue", cex.axis = cex.axis, ...);
par(new = TRUE); }
txtleg <- c(txtleg, "Ptot" )
colleg <- c(colleg, "royalblue")
lwdleg <- c(lwdleg, lwdP/3)
axis(side = 2, at = pretty(0.8 * ylimP, n = 3), labels = pretty(0.8 * ylimP, n = 3), col.axis = "royalblue", col.ticks = "royalblue", cex.axis = cex.axis, ...)
par(new = TRUE)
}
if (BOOL_Psol) {
plot(SeqX2, DataPlotP$Psol[DataMonthlyInterAn$Dates], type = "h", xlim = range(SeqX1),
ylim = c(3*ylimP[1], ylimP[2]), lwd = lwdP, lend = 1, lty = 1, col = "lightblue",
xlab = "", ylab = "", xaxt = "n", yaxt = "n", yaxs = "i", bty = "n", ...);
txtleg <- c(txtleg, "Psol" ); colleg <- c(colleg, "lightblue"); lwdleg <- c(lwdleg, lwdP/3);
par(new = TRUE); }
ylim = c(3 * ylimP[1], ylimP[2]), lwd = lwdP, lend = 1, lty = 1, col = "lightblue",
xlab = "", ylab = "", xaxt = "n", yaxt = "n", yaxs = "i", bty = "n", ...)
txtleg <- c(txtleg, "Psol" )
colleg <- c(colleg, "lightblue")
lwdleg <- c(lwdleg, lwdP/3)
par(new = TRUE)
}
## Plot_flows
plot(NULL, type = "n", xlim = range(SeqX1), ylim = c(ylimQ[1], 2*ylimQ[2]), xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
if (BOOL_Qobs) { lines(1:nrow(DataPlotQ), DataPlotQ$Qobs, lwd = lwd * lwdk , lty = 1, col = par("fg") ); txtleg <- c(txtleg, "Qobs" ); colleg <- c(colleg, par("fg") ); lwdleg <- c(lwdleg, lwd); }
if (BOOL_Qsim) { lines(1:nrow(DataPlotQ), DataPlotQ$Qsim, lwd = lwd * lwdk , lty = 1, col = "orangered"); txtleg <- c(txtleg, "Qsim"); colleg <- c(colleg, "orangered"); lwdleg <- c(lwdleg, lwd); }
plot(NULL, type = "n", xlim = range(SeqX1), ylim = c(ylimQ[1], 2 * ylimQ[2]), xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
if (BOOL_Qobs) {
lines(1:nrow(DataPlotQ), DataPlotQ$Qobs, lwd = lwd * lwdk, lty = 1, col = par("fg") )
txtleg <- c(txtleg, "Qobs" )
colleg <- c(colleg, par("fg") )
lwdleg <- c(lwdleg, lwd)
}
if (BOOL_Qsim) {
lines(1:nrow(DataPlotQ), DataPlotQ$Qsim, lwd = lwd * lwdk, lty = 1, col = "orangered")
txtleg <- c(txtleg, "Qsim")
colleg <- c(colleg, "orangered")
lwdleg <- c(lwdleg, lwd)
}
## Axis_and_legend
axis(side = 1, at = SeqX1, tick = TRUE , labels = xLabels1, cex.axis = cex.axis, ...)
axis(side = 1, at = SeqX2, tick = FALSE, labels = xLabels2, cex.axis = cex.axis, ...)
axis(side = 2, at = pretty(ylimQ), labels = pretty(ylimQ), cex.axis = cex.axis, ...)
par(las = 0); mtext(side = 1, labX, line = line, cex = cex.lab); par(las = 0);
posleg <- "topright"; txtlab <- "flow regime";
if (BOOL_Pobs) { posleg <- "right"; txtlab <- "precip. & flow regime"; }
par(las = 0); mtext(side = 2, paste0(txtlab, " ", plotunitregime), line = line, cex = cex.lab); par(las = 0);
mtext(side = 1, labX, line = line, cex = cex.lab)
posleg <- "topright"
txtlab <- "flow"
if (BOOL_Pobs) {
posleg <- "right"
txtlab <- "precip. & flow"
}
mtext(side = 2, paste(txtlab, plotunitregime), line = line, cex = cex.lab)
if (!is.null(BasinArea)) {
Factor <- Factor_UNIT_M3S / (365.25 / 12)
axis(side = 4, at = pretty(ylimQ*Factor)/Factor, labels = pretty(ylimQ*Factor), cex.axis = cex.axis, ...);
par(las = 0); mtext(side = 4, paste0("flow regime ", "[m3/s]"), line = line, cex = cex.lab); par(las = 0); }
### posleg <- "topright"; if (BOOL_Pobs) { posleg <- "right"; }
### legend(posleg, txtleg, col = colleg, lty = 1, lwd = lwdleg, bty = "o", bg = bg, box.col = bg, cex = cex.leg)
axis(side = 4, at = pretty(ylimQ * Factor)/Factor, labels = pretty(ylimQ * Factor), cex.axis = cex.axis, ...)
mtext(side = 4, paste("flow", "[m3/s]"), line = line, cex = cex.lab)
}
box()
}
}
## Cumulative_frequency
if (BOOLPLOT_CumFreq) {
kPlot <- kPlot+1; mar <- c(6, 5, 1, 5);
par(new = FALSE, mar = mar, las = 0)
xlim <- c(0, 1);
if ( BOOL_Qobs & !BOOL_Qsim) { SelectNotZero <- SelectQobsNotZero;
ylim <- range(log(Qobs[IndPeriod_Plot][SelectNotZero]), na.rm = TRUE); }
if (!BOOL_Qobs & BOOL_Qsim) { SelectNotZero <- SelectQsimNotZero;
ylim <- range(log(OutputsModel$Qsim[IndPeriod_Plot][SelectNotZero]), na.rm = TRUE); }
if ( BOOL_Qobs & BOOL_Qsim) { SelectNotZero <- SelectQobsNotZero & SelectQsimNotZero;
ylim <- range(log(c(Qobs[IndPeriod_Plot][SelectNotZero], OutputsModel$Qsim[IndPeriod_Plot][SelectNotZero])), na.rm = TRUE); }
kPlot <- kPlot + 1
mar <- c(6, 5, 1, 5)
par(new = FALSE, mar = mar)
xlim <- c(0, 1)
if ( BOOL_Qobs & !BOOL_Qsim) {
# SelectNotZero <- SelectQobsNotZero
ylim <- range(log(Qobs[SelectNotZero]), na.rm = TRUE)
}
if (!BOOL_Qobs & BOOL_Qsim) {
# SelectNotZero <- SelectQsimNotZero
ylim <- range(log(OutputsModel$Qsim[SelectNotZero]), na.rm = TRUE)
}
if ( BOOL_Qobs & BOOL_Qsim) {
# SelectNotZero <- SelectQobsNotZero & SelectQsimNotZero
ylim <- range(log(c(Qobs[SelectNotZero], OutputsModel$Qsim[SelectNotZero])), na.rm = TRUE)
}
SelectNotZero <- ifelse(is.na(SelectNotZero), FALSE, SelectNotZero)
if (any(SelectNotZero)) {
plot(0, 0, type = "n",
xlim = xlim, ylim = ylim,
xaxt = "n", yaxt = "n",
xlab = "", ylab = "", main = "", ...);
axis(side = 1, at = pretty(xlim), labels = pretty(xlim), cex.axis = cex.axis, ...);
par(las = 0); mtext(side = 1, text = "non-exceedance prob. [-]", line = line, cex = cex.lab); par(las = 0);
axis(side = 2, at = seqDATA1, labels = seqDATA2, cex.axis = cex.axis, ...)
par(las = 0); mtext(side = 2, text = paste("flow ", plotunit, "", sep = ""), line = line, cex = cex.lab); par(las = 0);
txtleg <- NULL; colleg <- NULL
xlab = "", ylab = "", ...)
axis(side = 1, at = pretty(xlim), labels = pretty(xlim), cex.axis = cex.axis, ...)
mtext(side = 1, text = "non-exceedance prob. [-]", line = line, cex = cex.lab)
axis(side = 2, at = seqDATA1, labels = seqDATA2, cex.axis = cex.axis, ...)
mtext(side = 2, text = paste("flow", plotunit), line = line, cex = cex.lab)
txtleg <- NULL
colleg <- NULL
if (BOOL_Qobs) {
DATA2 <- log(Qobs[IndPeriod_Plot][SelectNotZero]);
SeqQuant <- seq(0, 1, by = 1/(length(DATA2))); Quant <- as.numeric(quantile(DATA2, SeqQuant, na.rm = TRUE));
Fn <- ecdf(DATA2); YY <- DATA2; YY <- YY[order( Fn(DATA2) )]; XX <- Fn(DATA2); XX <- XX[order( Fn(DATA2) )];
lines(XX, YY, lwd = lwd, col = par("fg"));
txtleg <- c(txtleg, "observed"); colleg <- c(colleg, par("fg"))
DATA2 <- log(Qobs[SelectNotZero])
Fn <- ecdf(DATA2)
YY <- DATA2
YY <- YY[order(Fn(DATA2))]
XX <- Fn(DATA2)
XX <- XX[order(Fn(DATA2))]
lines(XX, YY, lwd = lwd, col = par("fg"))
txtleg <- c(txtleg, "observed")
colleg <- c(colleg, par("fg"))
}
if (BOOL_Qsim) {
DATA2 <- log(OutputsModel$Qsim[IndPeriod_Plot][SelectNotZero]);
SeqQuant <- seq(0, 1, by = 1/(length(DATA2))); Quant <- as.numeric(quantile(DATA2, SeqQuant, na.rm = TRUE));
Fn <- ecdf(DATA2); YY <- DATA2; YY <- YY[order( Fn(DATA2) )]; XX <- Fn(DATA2); XX <- XX[order( Fn(DATA2) )];
lines(XX, YY, lwd = lwd, col = "orangered");
txtleg <- c(txtleg, "simulated"); colleg <- c(colleg, "orangered")
DATA2 <- log(OutputsModel$Qsim[SelectNotZero])
Fn <- ecdf(DATA2)
YY <- DATA2
YY <- YY[order(Fn(DATA2))]
XX <- Fn(DATA2)
XX <- XX[order(Fn(DATA2))]
lines(XX, YY, lwd = lwd, col = "orangered")
txtleg <- c(txtleg, "simulated")
colleg <- c(colleg, "orangered")
}
if (!is.null(BasinArea)) {
Factor <- Factor_UNIT_M3S;
axis(side = 4, at = seqDATA1, labels = round(seqDATA2*Factor, digits = 2), cex.axis = cex.axis, ...)
par(las = 0); mtext(side = 4, paste0("flow ", "[m3/s]"), line = line, cex = cex.lab); par(las = 0)
Factor <- Factor_UNIT_M3S
axis(side = 4, at = seqDATA1, labels = round(seqDATA2 * Factor, digits = 2), cex.axis = cex.axis, ...)
mtext(side = 4, paste("flow", "[m3/s]"), line = line, cex = cex.lab)
}
legend("topleft", txtleg, col = colleg, lty = 1, lwd = lwd, bty = "o", bg = bg, box.col = bg, cex = cex.leg)
legend("bottomright", "log scale", lty = 1, col = NA, bty = "o", bg = bg, box.col = bg, cex = cex.leg)
} else {
plot(0, 0, type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", main = "", ...);
par(las = 0); mtext(side = 1, text = "non-exceedance prob. [-]", line = line, cex = cex.lab); par(las = 0);
par(las = 0); mtext(side = 2, text = paste("flow ", plotunit, "", sep = ""), line = line, cex = cex.lab); par(las = 0);
plot(0, 0, type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ...)
mtext(side = 1, text = "non-exceedance prob. [-]", line = line, cex = cex.lab)
mtext(side = 2, text = paste("flow", plotunit), line = line, cex = cex.lab)
text(0, 0, labels = "NO COMMON DATA", col = "grey40")
}
box()
......@@ -505,43 +858,43 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
## Correlation_QQ
if (BOOLPLOT_CorQQ) {
kPlot <- kPlot+1; mar <- c(6, 5, 1, 5);
par(new = FALSE, mar = mar, las = 0)
kPlot <- kPlot + 1
mar <- c(6, 5, 1, 5)
par(new = FALSE, mar = mar)
if (any(SelectNotZero)) {
ylim <- log(range(c(Qobs[IndPeriod_Plot][SelectQobsNotZero & SelectQsimNotZero], OutputsModel$Qsim[IndPeriod_Plot][SelectQobsNotZero & SelectQsimNotZero]), na.rm = TRUE))
plot(log(Qobs[IndPeriod_Plot][SelectQobsNotZero & SelectQsimNotZero]),
log(OutputsModel$Qsim[IndPeriod_Plot][SelectQobsNotZero & SelectQsimNotZero]),
ylim <- log(range(c(Qobs[SelectQobsNotZero & SelectQsimNotZero], OutputsModel$Qsim[SelectQobsNotZero & SelectQsimNotZero]), na.rm = TRUE))
plot(log(Qobs[SelectQobsNotZero & SelectQsimNotZero]),
log(OutputsModel$Qsim[SelectQobsNotZero & SelectQsimNotZero]),
type = "p", pch = 1, cex = 0.9, col = par("fg"), lwd = lwd,
xlim = ylim, ylim = ylim, xaxt = "n", yaxt = "n", xlab = "", ylab = "", ...)
abline(a = 0, b = 1, col = "royalblue", lwd = lwd);
axis(side = 1, at = seqDATA1, labels = seqDATA2, cex = cex.leg, cex.axis = cex.axis, ...);
axis(side = 2, at = seqDATA1, labels = seqDATA2, cex = cex.leg, cex.axis = cex.axis, ...);
par(las = 0); mtext(side = 1, paste0("observed flow ", plotunit), line = line, cex = cex.lab); par(las = 0);
par(las = 0); mtext(side = 2, paste0("simulated flow ", plotunit), line = line, cex = cex.lab); par(las = 0);
abline(a = 0, b = 1, col = "royalblue", lwd = lwd)
axis(side = 1, at = seqDATA1, labels = seqDATA2, cex = cex.leg, cex.axis = cex.axis, ...)
axis(side = 2, at = seqDATA1, labels = seqDATA2, cex = cex.leg, cex.axis = cex.axis, ...)
mtext(side = 1, paste("observed flow", plotunit), line = line, cex = cex.lab)
mtext(side = 2, paste("simulated flow", plotunit), line = line, cex = cex.lab)
if (!is.null(BasinArea)) {
Factor <- Factor_UNIT_M3S;
axis(side = 4, at = seqDATA1, labels = round(seqDATA2*Factor, digits = 2), cex.axis = cex.axis, ...);
par(las = 0); mtext(side = 4, paste0("simulated flow ", "[m3/s]"), line = line, cex = cex.lab); par(las = 0)
Factor <- Factor_UNIT_M3S
axis(side = 4, at = seqDATA1, labels = round(seqDATA2 * Factor, digits = 2), cex.axis = cex.axis, ...)
mtext(side = 4, paste("simulated flow", "[m3/s]"), line = line, cex = cex.lab)
}
legend("bottomright", "log scale", lty = 1, col = NA, bty = "o", bg = bg, box.col = bg, cex = cex.leg)
} else {
plot(0, 0, type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ...)
par(las = 0); mtext(side = 1, paste0("observed flow ", plotunit), line = line, cex = cex.lab); par(las = 0);
par(las = 0); mtext(side = 2, paste0("simulated flow ", plotunit), line = line, cex = cex.lab); par(las = 0);
mtext(side = 1, paste("observed flow", plotunit), line = line, cex = cex.lab)
mtext(side = 2, paste("simulated flow", plotunit), line = line, cex = cex.lab)
text(0, 0, labels = "NO COMMON DATA", col = "grey40")
}
box()
}
## Empty_plots
while (kPlot < iPlotMax) {
kPlot <- kPlot+1;
par(new = FALSE)
plot(0, 0, type = "n", xlab = "", ylab = "", axes = FALSE, ...)
if (exists("iPlotMax")) {
while (kPlot < iPlotMax) {
kPlot <- kPlot + 1
par(new = FALSE)
plot(0, 0, type = "n", xlab = "", ylab = "", axes = FALSE, ...)
}
}
## Restoring_layout_options
layout(1);
}
R/plot_OutputsModel.R deleted 100644 → 0
View file @ b34b5d06
plot_OutputsModel <- function(OutputsModel, PlotChoice = "all", ...) {
.Deprecated(new = "plot_OutputsModel", package = NULL,
msg = "Deprecated function. Please, use plot.OutputsModel() or plot() on an object of class OutputsModel.",
old = as.character(sys.call(sys.parent()))[1L])
if (!missing(OutputsModel)) {
warning("Deprecated \"OutputsModel\" argument. Please, use \"x\" instead.")
}
if (!missing(PlotChoice)) {
warning("Deprecated \"PlotChoice\" argument. Please, use \"which\" instead.")
}
plot.OutputsModel(x = OutputsModel, which = PlotChoice, ...)
}
README.md 0 → 100644
View file @ 4bccfee2
# airGR: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
## Overview
This package brings into R the hydrological modelling tools developed at INRAE-Antony ([Catchment Hydrology research group](https://webgr.inrae.fr/home/) of the HYCAR Research Unit, France), including rainfall-runoff models (**GR4H**, **GR5H**, **GR4J**, **GR5J**, **GR6J**, **GR2M**, **GR1A**) that can be applied either on a **lumped** or **semi-distributed** way. A snow accumulation and melt model (**CemaNeige**) and the associated functions for the calibration and evaluation of models are also included. Each model core is coded in Fortran to ensure low computational time. The other package functions (i.e. mainly the calibration algorithm and the computation of the efficiency criteria) are coded in R.
## Installation
``` r
install.packages("airGR")
```
## Functions and objects
The airGR package has been designed to fulfil two major requirements: facilitate the use by non-expert users and allow flexibility regarding the addition of external criteria, models or calibration algorithms. The names of the functions and their arguments were chosen to this end.
The package is mostly based on three families of functions:
- the functions belonging to the `RunModel` family require three arguments: `InputsModel`, `RunOptions` and `Param`; please refer to help pages `CreateInputsModel` and `CreateRunOptions` for further details and examples;
- the functions belonging to the `ErrorCrit` family require two arguments: `InputsCrit` and `OutputsModel`; please refer to help pages `CreateInputsCrit` and `RunModel` for further details and examples;
- the functions belonging to the `Calibration` family require four arguments: `InputsModel`, `RunOptions`, `InputsCrit` and `CalibOptions`; please refer to help pages `CreateInputsModel`, `CreateRunOptions`, `CreateInputsCrit` and `CreateCalibOptions` for further details and examples.
In order to limit the risk of mis-use and increase the flexibility of these main functions, we imposed the structure of their arguments and defined their class. Most users will not need to worry about these imposed structures since functions are provided to prepare these arguments for them: `CreateInputsModel`, `CreateRunOptions`, `CreateInputsCrit`, `CreateCalibOptions`. However, advanced users wishing to supplement the package with their own models will need to comply with these imposed structures and refer to the package source codes to get all the specification requirements.
## Models
Seven hydrological models and one snow melt and accumulation model are implemented in airGR. The hydrological models can be applied either on a lumped way or on a semi-distributed way (on sub-catchments). The snow model can either be used alone or with the daily or hourly hydrological models. Naturally each hydrological model can also be used alone.
These models can be called within airGR using the following functions:
- `RunModel_GR4H`: four-parameter hourly lumped hydrological model (Mathevet, 2005)
- `RunModel_GR5H`: five-parameter hourly lumped hydrological model (Ficchi, 2017; Ficchi *et al.*, 2019)
- `RunModel_GR4J`: four-parameter daily lumped hydrological model (Perrin *et al.*, 2003)
- `RunModel_GR5J`: five-parameter daily lumped hydrological model (Le Moine, 2008)
- `RunModel_GR6J`: six-parameter daily lumped hydrological model (Pushpalatha *et al.*, 2011)
- `RunModel_GR2M`: two-parameter monthly lumped hydrological model (Mouelhi, 2003; Mouelhi *et al.*, 2006a)
- `RunModel_GR1A`: one-parameter yearly lumped hydrological model (Mouelhi, 2003; Mouelhi *et al.*, 2006b)
- `RunModel_CemaNeige`: two-parameter degree-day snow melt and accumulation daily model (Valéry *et al.*, 2014; Riboust *et al.*, 2019)
- `RunModel_CemaNeigeGR4H`: combined use of GR4H and CemaNeige
- `RunModel_CemaNeigeGR5H`: combined use of GR5H and CemaNeige
- `RunModel_CemaNeigeGR4J`: combined use of GR4J and CemaNeige
- `RunModel_CemaNeigeGR5J`: combined use of GR5J and CemaNeige
- `RunModel_CemaNeigeGR6J`: combined use of GR6J and CemaNeige
## How to get started
To learn how to use the functions from the airGR package, it is recommended to follow the five steps described below:
1. refer to the help for `RunModel_GR4J` then run the provided example to assess how to make a simulation;
2. refer to the help for `CreateInputsModel` to understand how the inputs of a model are prepared/organised;
3. refer to the help for `CreateRunOptions` to understand how the run options of a model are parametrised/organised;
4. refer to the help for `ErrorCrit_NSE` and `CreateInputsCrit` to understand how the computation of an error criterion is prepared/made;
5. refer to the help for `Calibration_Michel`, run the provided example and then refer to the help for `CreateCalibOptions` to understand how a model calibration is prepared/made.
To get started with the package, you can refer to the 'get_started' vignette (`vignette("V01_get_started", package = "airGR")`). To know how to use the models on a semi-distributed way, you can refer to the 'sd_model' vignette (`vignette("V05_sd_model", package = "airGR")`). For more information, please visit the [airGR website](https://hydrogr.github.io/airGR/).
## References
- Ficchi, A. (2017). An adaptive hydrological model for multiple time-steps: Diagnostics and improvements based on fluxes consistency. PhD thesis, UPMC - Irstea Antony, Paris, France.
- Ficchi, A., Perrin, C. and Andréassian, V. (2019). Hydrological modelling at multiple sub-daily time steps: model improvement via flux-matching. Journal of Hydrology, 575, 1308-1327, doi: [10.1016/j.jhydrol.2019.05.084](https://www.doi.org/10.1016/j.jhydrol.2019.05.084).
- Le Moine, N. (2008). Le bassin versant de surface vu par le souterrain : une voie d'amélioration des performances et du réalisme des modèles pluie-débit ?, PhD thesis (in French), UPMC - Cemagref Antony, Paris, France, 324 pp.
- Mathevet, T. (2005). Quels modèles pluie-débit globaux pour le pas de temps horaire ? Développement empirique et comparaison de modèles sur un large échantillon de bassins versants, PhD thesis (in French), ENGREF - Cemagref Antony, Paris, France, 463 pp.
- Mouelhi S. (2003). Vers une chaîne cohérente de modèles pluie-débit conceptuels globaux aux pas de temps pluriannuel, annuel, mensuel et journalier, PhD thesis (in French), ENGREF - Cemagref Antony, Paris, France, 323 pp.
- Mouelhi, S., Michel, C., Perrin, C. and Andréassian, V. (2006a). Stepwise development of a two-parameter monthly water balance model, Journal of Hydrology, 318(1-4), 200-214, doi: [10.1016/j.jhydrol.2005.06.014](https://www.doi.org/10.1016/j.jhydrol.2005.06.014).
- Mouelhi, S., Michel, C., Perrin, C. and Andréassian, V. (2006b). Linking stream flow to rainfall at the annual time step: the Manabe bucket model revisited, Journal of Hydrology, 328, 283-296, doi: [10.1016/j.jhydrol.2005.12.022](https://www.doi.org/10.1016/j.jhydrol.2005.12.022).
- Perrin, C., Michel, C. and Andréassian, V. (2003). Improvement of a parsimonious model for streamflow simulation, Journal of Hydrology, 279(1-4), 275-289, doi: [10.1016/S0022-1694(03)00225-7](https://www.doi.org/10.1016/S0022-1694(03)00225-7).
- Pushpalatha, R., Perrin, C., Le Moine, N., Mathevet, T. and Andréassian, V. (2011). A downward structural sensitivity analysis of hydrological models to improve low-flow simulation, Journal of Hydrology, 411(1-2), 66-76, doi: [10.1016/j.jhydrol.2011.09.034](https://www.doi.org/10.1016/j.jhydrol.2011.09.034).
- Riboust, P., Thirel, G., Le Moine, N. and Ribstein, P. (2019). Revisiting a simple degree-day model for integrating satellite data: Implementation of SWE-SCA hystereses. Journal of Hydrology and Hydromechanics, 67(1), 70–81, doi: [10.2478/johh-2018-0004](https://www.doi.org/10.2478/johh-2018-0004).
- Valéry, A., Andréassian, V. and Perrin, C. (2014). "As simple as possible but not simpler": What is useful in a temperature-based snow-accounting routine? Part 2 - Sensitivity analysis of the Cemaneige snow accounting routine on 380 catchments, Journal of Hydrology, 517(0), 1176-1187, doi: [10.1016/j.jhydrol.2014.04.058](https://www.doi.org/10.1016/j.jhydrol.2014.04.058).
airGR.Rproj
View file @ 4bccfee2
......@@ -9,9 +9,12 @@ UseSpacesForTab: Yes
NumSpacesForTab: 2
Encoding: UTF-8
RnwWeave: Sweave
RnwWeave: knitr
LaTeX: pdfLaTeX
AutoAppendNewline: Yes
StripTrailingWhitespace: Yes
BuildType: Package
PackageUseDevtools: Yes
PackageInstallArgs: --no-multiarch --with-keep.source
data-raw/vignettes.R 0 → 100644
View file @ 4bccfee2
## code to prepare datasets used in vignettes
library(airGR)
source("tests/testthat/helper_vignettes.R")
# V02.1_param_optim
RunVignetteChunks("V02.1_param_optim")
save(resGLOB, resPORT,
file = "inst/vignettesData/vignetteParamOptim.rda",
version = 2)
save(algo, Ind_Run, InputsCrit_inv, InputsModel, MOptimGR4J, optMO, RunOptions,
file = "inst/vignettesData/vignetteParamOptimCaramel.rda",
version = 2)
# V02.2_param_mcmc
RunVignetteChunks("V02.2_param_mcmc")
save(gelRub, multDRAM,
file = "inst/vignettesData/vignetteParamMCMC.rda",
version = 2)
# V04_cemaneige_hysteresis
RunVignetteChunks("V04_cemaneige_hysteresis")
save(OutputsCrit_Cal, OutputsCrit_Cal_NoHyst, OutputsCrit_Val, OutputsCrit_Val_NoHyst,
file = "inst/vignettesData/vignetteParamMCMC.rda",
version = 2)
data/L0123001.rda
View file @ 4bccfee2
No preview for this file type
data/X0310010.rda 0 → 100644
View file @ 4bccfee2
File added
data/exampleSimPlot.rda 0 → 100644
View file @ 4bccfee2
File added
inst/CITATION
View file @ 4bccfee2
......@@ -5,36 +5,38 @@ note <- sprintf("R package version %s", version)
citHeader("To cite airGR in publications use these two references:")
bibentry(bibtype="Article",
title = "The Suite of Lumped {GR} Hydrological Models in an {R} package",
author = personList(as.person("L. Coron"), as.person("G. Thirel"), as.person("O. Delaigue"), as.person("C. Perrin"), as.person("V. Andréassian")),
journal = "Environmental Modelling and Software",
year = "2017",
pages = {"166-171"},
volume = {94},
doi = "10.1016/j.envsoft.2017.05.002",
textVersion =
bibentry(bibtype ="Article",
title = "The Suite of Lumped {GR} Hydrological Models in an {R} package",
author = c(as.person("L. Coron"), as.person("G. Thirel"), as.person("O. Delaigue"), as.person("C. Perrin"), as.person("V. Andréassian")),
journal = "Environmental Modelling and Software",
year = "2017",
pages = {"166-171"},
volume = {94},
doi = "10.1016/j.envsoft.2017.05.002",
textVersion =
paste("Coron, L., Thirel, G., Delaigue, O., Perrin, C. and Andréassian, V.",
"(2017).",
"(2017).",
"The Suite of Lumped GR Hydrological Models in an R package.",
"Environmental Modelling and Software,",
"94, 166-171.",
"DOI: 10.1016/j.envsoft.2017.05.002.",
"94, 166-171,",
"doi: 10.1016/j.envsoft.2017.05.002.",
sep = " ")
)
bibentry(bibtype="Manual",
title = "{airGR}: Suite of {GR} Hydrological Models for Precipitation-Runoff Modelling",
author = personList(as.person("L. Coron"), as.person("C. Perrin"), as.person("O. Delaigue"), as.person("G. Thirel"), as.person("C. Michel")),
journal = "R News",
year = year,
note = note,
url = "https://CRAN.R-project.org/package=airGR/",
textVersion =
paste0("Coron, L., Perrin, C., Delaigue, O., Thirel, G. and Michel, C. (",
bibentry(bibtype = "Manual",
title = "{airGR}: Suite of {GR} Hydrological Models for Precipitation-Runoff Modelling",
author = c(as.person("L. Coron"), as.person("O. Delaigue"), as.person("G. Thirel"), as.person("D. Dorchies"), as.person("C. Perrin"), as.person("C. Michel")),
journal = "R News",
year = year,
note = note,
doi = "10.15454/EX11NA",
url = "https://CRAN.R-project.org/package=airGR",
textVersion =
paste0("Coron, L., Delaigue, O., Thirel, G., Dorchies, D., Perrin, C. and Michel, C. (",
year,
"). airGR: Suite of GR Hydrological Models for Precipitation-Runoff Modelling. ",
note,
". URL: https://CRAN.R-project.org/package=airGR/."
note, ", ",
"doi: 10.15454/EX11NA", ", ",
"URL: https://CRAN.R-project.org/package=airGR."
)
)
inst/WORDLIST 0 → 100644
View file @ 4bccfee2
AE
AExch
AgroParisTech
Alpes
Altitudinal
Anctil
Baillon
Banque
BasinInfo
BasinObs
Berthet
BoolCrit
CRIT
CalibOptions
CemaNeige
CemaNeigeGR
CemaNeigeLayers
Cemagref
Cemaneige
Clapière
Coprs
CreateRunOptions
CritBestValue
CritCompo
CritFinal
CritName
CritValue
Customisation
DEoptim
DOI
DatesR
Développement
EI
ENGREF
Ei
Embrun
Eq
ErrorCrit
Exch
ExpStore
FME
Ficchi
Ficchì
FixedParam
Fortin
Franchistéguy
GAPX
GRP
Gelman
Glocalmax
GlocmaxCemaNeigeLayers
Gratio
Gthr
GthrCemaNeigeLayers
Gthreshold
HYCAR
Hautes
Hervieu
HistCrit
HistParamR
Hydrol
Hydrological
Hydrologie
Hydromechanics
HypsoData
INRAE
Imax
IndPeriod
IniResLevels
IniStates
Initialisation
InputsCrit
InputsModel
Interc
Irstea
KGE
Kling
LH
Lavenne
LayerFracSolidPrecip
LayerPrecip
LayerTempMax
LayerTempMean
LayerTempMin
Lerat
Lobligeois
Loumagne
MODIS
Manabe
MatBoolCrit
MeanAnSolidPrecip
Mieux
Modelling
Modélisation
Moine
MultiCrit
MultiCritNames
MultiCritValues
MultiCritWeights
Multiobjective
MyScriptBlocks
Météo
NGSA
NIter
NParam
NRuns
NSE
NSIDC
Nota
ORCID
Otamin
Oudin's
OutputsModel
POSIXct
POSIXlt
POSIXt
Param
ParamFinalR
ParamOut
Paulin
Perc
Pliq
PliqAndMelt
Pn
PotEvap
PotMelt
Precip
Ps
Psol
QD
QRExp
Qd
Qls
Qmm
Qobs
QrExp
Qsim
Quels
Rdata
Regularisation
Riboust
Ribstein
Rmalschains
RunModel
RunOptions
SAFRAN
SCA
SLS
SWE
Safran
SearchRanges
Semidistributed
SnowPack
Soubeyroux
StartParamDistrib
StartParamList
StateEnd
SubCritNames
SubCritValues
Sutcliffe
Syst
ThermalState
Thibault
TransfoParam
Turcotte
UPMC
USACE
VarObs
Vers
Villeneuve
WarmUp
WarmUpQsim
Yilmaz
ZLayers
addeed
additionnal
aggreg
airGRdatassim
airGRteaching
al
améliore
annuel
anthropogenic
appliquée
atomics
authorised
bassin
bassins
behaviour
caRamel
calcullation
cemaneige
centred
chaîne
cohérente
comparaison
conceptuels
connaitre
crues
customise
d'amélioration
d'un
de
degC
des
disaggregated
disaggregation
doi
du
débit
eTG
eaufrance
empirique
et
evapotranspiration
fract
français
fulfil
funcions
genericity
ggmcmc
globaux
grDevices
gradP
gradT
gridded
hess
horaire
http
https
hydroPSO
hydrograph
hydrographs
hydrological
hydrometeorological
hypsometric
hystereses
hysteretic
iLayer
il
illustred
initialisation
initialise
initialised
interannual
intercatchment
intercetion
irrealistically
jhydrol
johh
journalier
labelling
le
linux
macOS
maximise
maximised
mcmc
mco
memetic
mensuel
metres
mimina
minimise
minimised
mis
mistyped
modelling
modèles
modélisation
multiobjective
neige
neutralisation
nivale
nn
notcomputed
nsidc
optim
optimisation
optimised
optimises
organised
param
parametrised
pareto
perf
petits
pluie
pluies
pluriannuel
pluviometry
pre
precip
predefining
priori
précipitations
recomanded
recommanded
regularisation
repectively
representativeness
rescale
rescaled
rescaling
ruraux
réalisme
sd
serie
skilful
snowmelt
snowpack
souterrain
spatiale
streamflow
subcatchment
subcatchments
sur
swnow
th
traceplot
transfo
un
unajusted
uncalibrated
underperforming
une
unitless
verifications
versant
versants
visualisation
visualised
voie
vu
www
échantillon
évaluation
inst/modelsFeatures/FeatModelsGR.csv 0 → 100644
View file @ 4bccfee2
CodeMod;NameMod;NbParam;TimeUnit;Id;Class;Pkg
GR1A;GR1A;1;yearly;NA;GR;airGR
GR2M;GR2M;2;monthly;NA;GR;airGR
GR4J;GR4J;4;daily;NA;GR;airGR
GR5J;GR5J;5;daily;NA;GR;airGR
GR6J;GR6J;6;daily;NA;GR;airGR
GR4H;GR4H;4;hourly;NA;GR;airGR
GR5H;GR5H;5;hourly;NA;GR;airGR
CemaNeige;CemaNeige;2;NA;NA;NA;airGR
CemaNeigeGR4J;CemaNeigeGR4J;6;daily;NA;GR;airGR
CemaNeigeGR5J;CemaNeigeGR5J;7;daily;NA;GR;airGR
CemaNeigeGR6J;CemaNeigeGR6J;8;daily;NA;GR;airGR
CemaNeigeGR4H;CemaNeigeGR4H;6;hourly;NA;GR;airGR
CemaNeigeGR5H;CemaNeigeGR5H;7;hourly;NA;GR;airGR
Lag;Lag;1;NA;NA;SD;airGR
inst/vignettesData/vignetteCNHysteresis.rda 0 → 100644
View file @ 4bccfee2
File added
inst/vignettesData/vignetteParamMCMC.rda
View file @ 4bccfee2
No preview for this file type
inst/vignettesData/vignetteParamOptim.rda
View file @ 4bccfee2
No preview for this file type