diff --git a/DESCRIPTION b/DESCRIPTION index 1556c1207ab575582c91cb816cc708e41c8437c8..3242392947c1b15c7f9716b6eb08d00d579a440d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: airGR Type: Package Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling -Version: 1.0.9.39 +Version: 1.0.9.40 Date: 2017-09-07 Authors@R: c( person("Laurent", "Coron", role = c("aut", "trl")), diff --git a/R/plot.OutputsModel.R b/R/plot.OutputsModel.R index 13983dd88b8f43e7b028fddffea8a917e0b0bb32..66ef22c15591c34c728f651778e5272f5ea5354a 100644 --- a/R/plot.OutputsModel.R +++ b/R/plot.OutputsModel.R @@ -82,7 +82,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = 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 + ## Plots_choices BOOLPLOT_Precip <- ( "Precip" %in% which & BOOL_Pobs ) BOOLPLOT_Temp <- ( "Temp" %in% which & BOOL_Snow ) BOOLPLOT_SnowPack <- ( "SnowPack" %in% which & BOOL_Snow ) @@ -92,7 +92,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = BOOLPLOT_CorQQ <- ( "CorQQ" %in% which & (BOOL_Qsim & BOOL_Qobs) & BOOL_FilterZero ) - ##Options + ## Options BLOC <- TRUE if (BLOC) { lwdk <- 1.8 @@ -172,7 +172,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = seqDATA2ba <- round(exp(seqDATA1ba), digits = 2) } - ##Precip + ## Precip if (BOOLPLOT_Precip) { kPlot <- kPlot + 1 mar <- c(3, 5, 1, 5) @@ -214,7 +214,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = } - ##Temp + ## Temp if (BOOLPLOT_Temp) { kPlot <- kPlot+1; mar <- c(3, 5, 1, 5); par(new = FALSE, mar = mar, las = 0) @@ -240,7 +240,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = } - ##SnowPack + ## SnowPack if (BOOLPLOT_SnowPack) { kPlot <- kPlot+1; mar <- c(3, 5, 1, 5); par(new = FALSE, mar = mar, las = 0) @@ -264,7 +264,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = } - ##Flows + ## Flows if (BOOLPLOT_Flows & log_scale) { kPlot <- kPlot+1; mar <- c(3, 5, 1, 5); par(new = FALSE, mar = mar, las = 0) @@ -323,115 +323,130 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = } - ##Regime + ## Regime if (BOOLPLOT_Regime) { - kPlot <- kPlot+1; mar <- c(6, 5, 1, 5); plotunitregime <- "[mm/month]"; + kPlot <- kPlot+1; mar <- c(6, 5, 1, 5); plotunitregime <- "[mm/month]" par(new = FALSE, mar = mar, las = 0) - ##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"); - ##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 %in% c("hour", "day")) { - 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"); - ##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)); - colnames(DataDailyInterAn) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim") - 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)); - colnames(DataMonthlyInterAn) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim") - 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"); - ##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)); - } - ##Plot_preparation - 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"; + ## 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" + if (BOOL_Pobs) { + txtlab <- "precip. & flow regime" + } + par(las = 0); mtext(side = 2, paste(txtlab, plotunitregime), line = line, cex = cex.lab); par(las = 0) } else { - DataPlotQ <- DataMonthlyInterAn; - SeqX1 <- seq(from = 0.5, to = 12.5, by = 1); - SeqX2 <- seq(from = 1 , to = 12 , by = 1); - labX <- ""; + ## 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"); + ## 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 %in% c("hour", "day")) { + 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"); + ## 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)); + colnames(DataDailyInterAn) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim") + 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)); + colnames(DataMonthlyInterAn) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim") + 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"); + ## 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)); + } + ## Plot_preparation + 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"; + } else { + DataPlotQ <- DataMonthlyInterAn; + SeqX1 <- seq(from = 0.5, to = 12.5, by = 1); + SeqX2 <- seq(from = 1 , to = 12 , by = 1); + labX <- ""; + } + 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; + ## 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", + 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); } + 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); } + ## 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); } + ## 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); + 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) + box() } - 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; - ##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", - 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); } - 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); } - ##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); } - ##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); - 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) - box() } - ##Cumulative_frequency + ## Cumulative_frequency if (BOOLPLOT_CumFreq) { kPlot <- kPlot+1; mar <- c(6, 5, 1, 5); par(new = FALSE, mar = mar, las = 0) @@ -484,7 +499,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = } - ##Correlation_QQ + ## Correlation_QQ if (BOOLPLOT_CorQQ) { kPlot <- kPlot+1; mar <- c(6, 5, 1, 5); par(new = FALSE, mar = mar, las = 0) @@ -514,14 +529,14 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = box() } - ##Empty_plots + ## Empty_plots while (kPlot < iPlotMax) { kPlot <- kPlot+1; par(new = FALSE) plot(0, 0, type = "n", xlab = "", ylab = "", axes = FALSE, ...) } - ##Restoring_layout_options + ## Restoring_layout_options layout(1);