Skip to content

GitLab

Commit bfc6e0fd authored by unknown's avatar unknown
Browse files

v1.0.9.40 plot.OutputsModel now runs when IndPeriod_Run < 30 days with GR4H and GR4J

parent 33986184
Hide whitespace changes
Inline Side-by-side
DESCRIPTION
View file @ bfc6e0fd
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")),
......
R/plot.OutputsModel.R
View file @ bfc6e0fd
......@@ -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);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment