Commit 2d46dc4c authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

v1.2.9.3 CLEAN: syntax revision of plot.OutputsModel

parent a1933bc5
Package: airGR
Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
Version: 1.2.9.2
Date: 2019-03-12
Version: 1.2.9.3
Date: 2019-03-13
Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
person("Charles", "Perrin", role = c("aut", "ths"), comment = c(ORCID = "0000-0001-8552-1881")),
......
......@@ -13,7 +13,7 @@ output:
### 1.2.9.2 Release Notes (2019-03-12)
### 1.2.9.3 Release Notes (2019-03-13)
......
......@@ -2,83 +2,158 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
cex.axis = 1, cex.lab = 0.9, cex.leg = 0.9, lwd = 1, verbose = TRUE, ...) {
OutputsModel <- x
if (!inherits(OutputsModel, "GR") & !inherits(OutputsModel, "CemaNeige")) {
stop(paste("OutputsModel not in the correct format for default plotting", sep = ""))
stop("'OutputsModel' not in the correct format for default plotting")
}
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") }
if (!is.vector( which)) { stop("which must be a vector of character") }
if (!is.character(which)) { stop("which must be a vector of character") }
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")
}
if (!is.vector(which)) {
stop("'which' must be a vector of character")
}
if (!is.character(which)) {
stop("'which' must be a vector of character")
}
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'")
stop("Incorrect element found in argument 'which':\nit can only contain 'all', 'Precip', 'Temp', 'SnowPack', 'Flows', 'Regime', 'CumFreq' or 'CorQQ'")
}
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'")
stop("Incorrect element found in argument 'which':\nwithout CemaNeige it can only contain 'all', 'Precip', 'Flows', 'Regime', 'CumFreq' or 'CorQQ'")
}
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")
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")
}
if ("all" %in% which) {
which <- c("Precip", "Temp", "SnowPack", "Flows", "Regime", "CumFreq", "CorQQ")
}
if (!BOOL_Dates) {
stop(paste("OutputsModel must contain at least DatesR to allow plotting", sep = "")) }
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", sep = "")) }
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)); }
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)]); }
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))
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", sep = "")) }
}
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"); }
if (BOOL_QsimZero & verbose) { warning("\t zeroes detected in Qsim: some plots in the log space will not be created using all time-steps"); }
BOOL_FilterZero <- 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("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)
}
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")
}
if (BOOL_QsimZero & verbose) {
warning( "\t zeroes detected in 'Qsim': some plots in the log space will not be created using all time-steps")
}
BOOL_FilterZero <- TRUE
## Plots_choices
BOOLPLOT_Precip <- ( "Precip" %in% which & BOOL_Pobs )
BOOLPLOT_Temp <- ( "Temp" %in% which & BOOL_Snow )
......@@ -87,67 +162,83 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
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
BLOC <- TRUE
if (BLOC) {
lwdk <- 1.8
line <- 2.6
bg <- NA
matlayout <- NULL; iPlot <- 0;
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;
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; }
# 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)
Xaxis <- 1:length(IndPeriod_Plot);
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]
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]
Seq1 <- 1:length(OutputsModel$DatesR[IndPeriod_Plot])
Seq2 <- 1:length(OutputsModel$DatesR[IndPeriod_Plot])
Labels2 <- format(OutputsModel$DatesR[IndPeriod_Plot], format = formatAxis)[Seq2]
}
}
if (!is.null(BasinArea)) {
Factor_MMH_M3S <- 60 * 60
Factor_MMD_M3S <- 60 * 60 * 24
......@@ -160,7 +251,7 @@ 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
......@@ -172,7 +263,7 @@ 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
......@@ -193,7 +284,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
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)
mtext(side = 2, paste("precip.", plotunit), cex = cex.lab, adj = 1, line = line)
par(las = 0)
if (BOOL_Psol) {
......@@ -217,116 +308,180 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
## Temp
if (BOOLPLOT_Temp) {
kPlot <- kPlot+1; mar <- c(3, 5, 1, 5);
kPlot <- kPlot + 1
mar <- c(3, 5, 1, 5)
par(new = FALSE, mar = mar, las = 0)
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; }
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) {
lines(OutputsModel$CemaNeigeLayers[[iLayer]]$Temp[IndPeriod_Plot], 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")
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)
par(las = 0)
mtext(side = 2, expression(paste0("temp. [", degree, "C]")), 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)
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, ...); }
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);
kPlot <- kPlot + 1
mar <- c(3, 5, 1, 5)
par(new = FALSE, mar = mar, las = 0)
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; }
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) {
lines(OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack[IndPeriod_Plot], 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)
par(las = 0)
mtext(side = 2, paste("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)
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, ...); }
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);
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)
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[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")
}
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);
par(las = 0)
mtext(side = 2, paste("flow", plotunit), line = line, cex = cex.lab)
par(las = 0)
if (!is.null(BasinArea)) {
Factor <- Factor_UNIT_M3S;
Factor <- Factor_UNIT_M3S
axis(side = 4, at = seqDATA1ba, labels = seqDATA2ba, cex.axis = cex.axis, ...)