diff --git a/R/plot.OutputsModel.R b/R/plot.OutputsModel.R index c11dc5b8477e77d9baa3c18fcd9f79ec5476acf6..b930f94f48da9f12c1749812446b7c87628b1234 100644 --- a/R/plot.OutputsModel.R +++ b/R/plot.OutputsModel.R @@ -15,12 +15,18 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = }) - 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] - IndPeriod_Plot <- seq_along(IndPeriod_Plot) + } else { + IndPeriod_Plot <- seq_along(OutputsModel$DatesR) } @@ -191,9 +197,6 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = # 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) } @@ -209,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) { @@ -310,7 +313,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = - Xaxis <- as.POSIXct(OutputsModel$DatesR[IndPeriod_Plot]) + Xaxis <- as.POSIXct(OutputsModel$DatesR) if (!is.null(BasinArea)) { Factor_UNIT_M3S <- switch(NameTS, @@ -340,7 +343,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = 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) @@ -348,7 +351,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = 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 = "", ...) @@ -359,7 +362,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = 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) - points(Xaxis, PsolLayerMean[IndPeriod_Plot], + points(Xaxis, PsolLayerMean, type = "h", xaxt = "n", yaxt = "n", yaxs = "i", ylim = ylim2, col = "lightblue", lwd = lwdP * lwdk, lend = 1, xlab = "", ylab = "", ...) @@ -378,23 +381,23 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = 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], + 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"), @@ -416,10 +419,10 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = 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 = "", ...) @@ -448,12 +451,12 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = SnowPackLayerMean <- SnowPackLayerMean + OutputsModel$CemaNeigeLayers[[iLayer]]$Temp/NLayers } } - plot(Xaxis, 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(Xaxis, 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(Xaxis, 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) @@ -483,9 +486,9 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = SnowPackLayerMean <- SnowPackLayerMean + OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack/NLayers } } - plot(Xaxis, 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(Xaxis, 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, ...) @@ -515,21 +518,21 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = DATA3[!SelectQsimNotZero] <- mean(OutputsModel$Qsim, na.rm = TRUE) / 10000 DATA3 <- log(DATA3) - ylim1 <- range(DATA3[IndPeriod_Plot], na.rm = TRUE) + ylim1 <- range(DATA3, na.rm = TRUE) if (BOOL_Qobs) { - ylim1 <- range(c(ylim1, DATA2[IndPeriod_Plot]), na.rm = TRUE) + 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[IndPeriod_Plot], lwd = lwd * lwdk , lty = 1, col = par("fg")) + 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[IndPeriod_Plot], lwd = lwd * lwdk , lty = 1, col = "orangered") + lines(Xaxis, DATA3, lwd = lwd * lwdk , lty = 1, col = "orangered") txtleg <- c(txtleg, "simulated") colleg <- c(colleg, "orangered") } @@ -549,21 +552,21 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = kPlot <- kPlot + 1 mar <- c(3, 5, 1, 5) par(new = FALSE, mar = mar) - ylim1 <- range(OutputsModel$Qsim[IndPeriod_Plot], na.rm = TRUE) + ylim1 <- range(OutputsModel$Qsim, na.rm = TRUE) if (BOOL_Qobs) { - ylim1 <- range(c(ylim1, Qobs[IndPeriod_Plot]), na.rm = TRUE) + 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[IndPeriod_Plot], lwd = lwd * lwdk , lty = 1, col = par("fg")) + 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[IndPeriod_Plot], lwd = lwd * lwdk , lty = 1, col = "orangered") + lines(Xaxis, OutputsModel$Qsim, lwd = lwd * lwdk , lty = 1, col = "orangered") txtleg <- c(txtleg, "simulated") colleg <- c(colleg, "orangered") } @@ -586,9 +589,9 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = mar <- c(3, 5, 1, 5) if (log_scale) { - errorQ <- log(OutputsModel$Qsim[IndPeriod_Plot]) - log(Qobs[IndPeriod_Plot]) + errorQ <- log(OutputsModel$Qsim) - log(Qobs) } else { - errorQ <- OutputsModel$Qsim[IndPeriod_Plot] - Qobs[IndPeriod_Plot] + errorQ <- OutputsModel$Qsim - Qobs } par(new = FALSE, mar = mar) ylim1 <- range(errorQ[SelectNotZero], na.rm = TRUE) @@ -634,18 +637,18 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = } 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")) + DataModel[, 1] <- as.numeric(format(OutputsModel$DatesR, format = "%Y%m%d%H")) if (BOOL_Pobs) { - DataModel[, 2] <- OutputsModel$Precip[IndPeriod_Plot] + DataModel[, 2] <- OutputsModel$Precip } if (BOOL_Psol) { - DataModel[, 3] <- PsolLayerMean[IndPeriod_Plot] + DataModel[, 3] <- PsolLayerMean } if (BOOL_Qobs) { - DataModel[, 4] <- Qobs[IndPeriod_Plot] + DataModel[, 4] <- Qobs } if (BOOL_Qsim) { - DataModel[, 5] <- OutputsModel$Qsim[IndPeriod_Plot] + DataModel[, 5] <- OutputsModel$Qsim } colnames(DataModel) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim") TxtDatesDataModel <- formatC(DataModel$Dates, format = "d", width = 8, flag = "0") @@ -791,15 +794,15 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = xlim <- c(0, 1) if ( BOOL_Qobs & !BOOL_Qsim) { # SelectNotZero <- SelectQobsNotZero - ylim <- range(log(Qobs[IndPeriod_Plot][SelectNotZero]), na.rm = TRUE) + ylim <- range(log(Qobs[SelectNotZero]), na.rm = TRUE) } if (!BOOL_Qobs & BOOL_Qsim) { # SelectNotZero <- SelectQsimNotZero - ylim <- range(log(OutputsModel$Qsim[IndPeriod_Plot][SelectNotZero]), na.rm = TRUE) + ylim <- range(log(OutputsModel$Qsim[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) + ylim <- range(log(c(Qobs[SelectNotZero], OutputsModel$Qsim[SelectNotZero])), na.rm = TRUE) } SelectNotZero <- ifelse(is.na(SelectNotZero), FALSE, SelectNotZero) if (any(SelectNotZero)) { @@ -814,7 +817,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = txtleg <- NULL colleg <- NULL if (BOOL_Qobs) { - DATA2 <- log(Qobs[IndPeriod_Plot][SelectNotZero]) + DATA2 <- log(Qobs[SelectNotZero]) SeqQuant <- seq(0, 1, by = 1/(length(DATA2))) Quant <- as.numeric(quantile(DATA2, SeqQuant, na.rm = TRUE)) Fn <- ecdf(DATA2) @@ -827,7 +830,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = colleg <- c(colleg, par("fg")) } if (BOOL_Qsim) { - DATA2 <- log(OutputsModel$Qsim[IndPeriod_Plot][SelectNotZero]) + DATA2 <- log(OutputsModel$Qsim[SelectNotZero]) SeqQuant <- seq(0, 1, by = 1/(length(DATA2))) Quant <- as.numeric(quantile(DATA2, SeqQuant, na.rm = TRUE)) Fn <- ecdf(DATA2) @@ -862,9 +865,9 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = 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)