Commit ed28715c authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

Merge conflict

refactor: review the x-axis management of the plot.OutputsModel
Refs #122, !6

Merge branch '122-review-the-x-axis-management-of-the-plot-outputsmodel' into dev

# Veuillez entrer un message de validation pour expliquer en quoi cette fusion est
# nécessaire, surtout si cela fusionne une branche amont mise à jour dans une branche de sujet.
#
# Les lignes commençant par '#' seront ignorées, et un message vide
# abandonne la validation.
parents 265725dc 36154a79
Pipeline #24993 passed with stage
in 9 minutes and 4 seconds
......@@ -219,16 +219,41 @@
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()
# }
# }
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)
}
})
}
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 parameters and resetting on exit
## save default graphical & time parameters and resetting on exit
opar <- par(no.readonly = TRUE)
on.exit(par(opar))
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) == 0) {
IndPeriod_Plot <- seq_along(OutputsModel$DatesR)
}
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("'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)) {
......@@ -53,27 +70,27 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
} 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")
......@@ -139,7 +156,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
}
}
## check dates
if (!BOOL_Dates) {
stop("'OutputsModel' must contain at least 'DatesR' to allow plotting")
......@@ -147,7 +164,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
if (inherits(OutputsModel, "GR") & !BOOL_Qsim) {
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))
......@@ -157,42 +174,29 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
# }
MyRollMean3 <- function(x, n) {
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]"
formatAxis <- "%m/%Y"
}
if (inherits(OutputsModel, "daily")) {
} else if (inherits(OutputsModel, "daily")) {
BOOL_TS <- TRUE
NameTS <- "day"
plotunit <- "[mm/d]"
formatAxis <- "%m/%Y"
}
if (inherits(OutputsModel, "monthly")) {
} else if (inherits(OutputsModel, "monthly")) {
BOOL_TS <- TRUE
NameTS <- "month"
plotunit <- "[mm/month]"
formatAxis <- "%m/%Y"
if (format(OutputsModel$DatesR[1L], format = "%d") != "01") {
OutputsModel$DatesR <- as.POSIXlt(format(OutputsModel$DatesR, format = "%Y-%m-01"), tz = "UTC", format = "%Y-%m-%d")
}
}
if (inherits(OutputsModel, "yearly")) {
} else if (inherits(OutputsModel, "yearly")) {
BOOL_TS <- TRUE
NameTS <- "year"
plotunit <- "[mm/y]"
formatAxis <- "%Y"
}
# if (!BOOL_TS) {
# stop("the time step of the model inputs could not be found")
# }
}
if (length(IndPeriod_Plot) == 0) {
IndPeriod_Plot <- 1:length(OutputsModel$DatesR)
}
if (inherits(OutputsModel, "CemaNeige")) {
NLayers <- length(OutputsModel$CemaNeigeLayers)
}
......@@ -208,12 +212,12 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
}
BOOL_QobsZero <- FALSE
if (BOOL_Qobs) {
SelectQobsNotZero <- round(Qobs[IndPeriod_Plot], 4) != 0
SelectQobsNotZero <- round(Qobs, 4) != 0
BOOL_QobsZero <- sum(!SelectQobsNotZero, na.rm = TRUE) > 0
}
BOOL_QsimZero <- FALSE
if (BOOL_Qsim) {
SelectQsimNotZero <- round(OutputsModel$Qsim[IndPeriod_Plot], 4) != 0
SelectQsimNotZero <- round(OutputsModel$Qsim, 4) != 0
BOOL_QsimZero <- sum(!SelectQsimNotZero, na.rm = TRUE) > 0
}
if ( BOOL_Qobs & !BOOL_Qsim) {
......@@ -232,11 +236,11 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
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
......@@ -248,30 +252,30 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
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
## Set plot arrangement
if (is.null(LayoutMat)) {
if (is.null(LayoutMat)) {
matlayout <- NULL
iPlot <- 0
iHght <- NULL
listBOOLPLOT1 <- c(Precip = BOOLPLOT_Precip,
PotEvap = BOOLPLOT_PotEvap | BOOLPLOT_ActuEvap,
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)
for (k in seq_len(Sum1)) {
matlayout <- rbind(matlayout, iPlot + c(1, 1, 1))
iPlot <- iPlot + 1
......@@ -305,24 +309,12 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
matlayout <- LayoutMat
}
layout(matlayout, widths = LayoutWidths, heights = LayoutHeights)
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]
}
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]
}
}
Xaxis <- as.POSIXct(OutputsModel$DatesR)
if (!is.null(BasinArea)) {
Factor_UNIT_M3S <- switch(NameTS,
hour = 60 * 60,
......@@ -332,9 +324,9 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
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)
......@@ -344,115 +336,106 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
seqDATA1ba <- log(seqDATA0 * Factor_UNIT_M3S)
seqDATA2ba <- round(exp(seqDATA1ba), digits = 2)
}
## Precip
if (BOOLPLOT_Precip) {
kPlot <- kPlot + 1
mar <- c(3, 5, 1, 5)
par(new = FALSE, mar = mar)
ylim1 <- range(OutputsModel$Precip[IndPeriod_Plot], na.rm = TRUE)
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, ...)
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 = "", ...)
}
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, ...)
points(Xaxis, PsolLayerMean,
type = "h", xaxt = "n", yaxt = "n", yaxs = "i", ylim = ylim2,
col = "lightblue", lwd = lwdP * lwdk, lend = 1,
xlab = "", ylab = "", ...)
}
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[IndPeriod_Plot], na.rm = TRUE)
ylim1 <- range(OutputsModel$PotEvap, na.rm = TRUE)
xlabE <- "pot. evap."
} else {
ylim1 <- range(c(OutputsModel$PotEvap[IndPeriod_Plot],
OutputsModel$AE[IndPeriod_Plot]),
ylim1 <- range(c(OutputsModel$PotEvap,
OutputsModel$AE),
na.rm = TRUE)
xlabE <- "evap."
}
ylim2 <- ylim1 #* c(1.0, 1.1)
plot(Xaxis, OutputsModel$PotEvap[IndPeriod_Plot],
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[IndPeriod_Plot],
type = "l", xaxt = "n", yaxt = "n", ylim = ylim2,
col = "green4", lwd = lwd * lwdk, lty = 3)
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),
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)
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, ...)
}
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[IndPeriod_Plot], na.rm = TRUE)
ylim1 <- range(OutputsModel$AE, na.rm = TRUE)
ylim2 <- ylim1 #* c(1.0, 1.1)
plot(Xaxis, OutputsModel$AE[IndPeriod_Plot],
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)
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, ...)
}
AxisTS(OutputsModel)
box()
}
## Temp
if (BOOLPLOT_Temp) {
kPlot <- kPlot + 1
......@@ -468,29 +451,26 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
SnowPackLayerMean <- SnowPackLayerMean + OutputsModel$CemaNeigeLayers[[iLayer]]$Temp/NLayers
}
}
plot(SnowPackLayerMean[IndPeriod_Plot], type = "n", ylim = ylim1, xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
plot(Xaxis, SnowPackLayerMean, 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)
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, ...)
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),
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
......@@ -506,33 +486,30 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
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", ...)
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(OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack[IndPeriod_Plot], lty = 3, col = "royalblue", lwd = lwd * lwdk *0.8)
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, ...)
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),
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, ...)
}
}