plot.OutputsModel.R 27.80 KiB
plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = NULL, which = "all", log_scale = FALSE, verbose = TRUE, ...){
  OutputsModel <- x
  if(!inherits(OutputsModel, "GR") & !inherits(OutputsModel, "CemaNeige")){
    stop(paste("OutputsModel not in the correct format for default plotting \n", sep = ""))
    return(NULL)
  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; } }
  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 \n"); return(NULL); } 
  if(!is.vector(   which)){ stop("which must be a vector of character \n"); return(NULL); } 
  if(!is.character(which)){ stop("which must be a vector of character \n"); return(NULL); } 
  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'")
    return(NULL)
  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'")
    return(NULL)
  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")
  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 \n", sep = "")); return(NULL); }
  if(inherits(OutputsModel, "GR") & !BOOL_Qsim){
    stop(paste("OutputsModel must contain at least Qsim to allow plotting \n", sep = "")); return(NULL); }
  if(BOOL_Dates){
    MyRollMean1 <- function(x, n){
      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)]); }
    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 \n", sep = "")); return(NULL); } 
  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 \n"); }
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
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 BOOLPLOT_Precip <- ( "Precip" %in% which & BOOL_Pobs ) BOOLPLOT_Temp <- ( "Temp" %in% which & BOOL_Snow ) BOOLPLOT_SnowPack <- ( "SnowPack" %in% which & BOOL_Snow ) BOOLPLOT_Flows <- ( "Flows" %in% which & (BOOL_Qsim | BOOL_Qobs) ) 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){ cexaxis <- 1.0; cexlab <- 0.9; cexleg = 1.0; lwdLine = 1.8; lineX = 2.6; lineY = 2.6; bgleg <- NA 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; # 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; } # dev.new(width = width, height = height) # } layout(matlayout); 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]; } } if(!is.null(BasinArea)){ Factor_MMH_M3S <- BasinArea/( 60*60/1000); Factor_MMD_M3S <- BasinArea/( 24*60*60/1000); Factor_MMM_M3S <- BasinArea/(365.25/12*24*60*60/1000); Factor_MMY_M3S <- BasinArea/( 365.25*24*60*60/1000);
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
if(NameTS == "hour" ){ Factor_UNIT_M3S <- Factor_MMH_M3S; } if(NameTS == "day" ){ Factor_UNIT_M3S <- Factor_MMD_M3S; } if(NameTS == "month"){ Factor_UNIT_M3S <- Factor_MMM_M3S; } if(NameTS == "year" ){ Factor_UNIT_M3S <- Factor_MMY_M3S; } } } kPlot <- 0 ## vector of Q values for the y-axis when it is expressed in seqDATA1 <- log(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)) seqDATA2 <- exp(seqDATA1) ##Precip if(BOOLPLOT_Precip){ kPlot <- kPlot+1; mar <- c(3, 5, 1, 5); par(new = FALSE, mar = mar, las = 0) ylim1 <- range(OutputsModel$Precip[IndPeriod_Plot], na.rm = TRUE); ylim2 <- ylim1 * c(1.0, 1.1); ylim2 <- rev(ylim2); lwdP <- 0.7; if(NameTS %in% c("month", "year")){ lwdP <- 2; } plot(Xaxis, OutputsModel$Precip[IndPeriod_Plot], type = "h", ylim = ylim2, col = "royalblue", lwd = lwdP, xaxt = "n", yaxt = "n", xlab = "", ylab = "", yaxs = "i", ...); axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cexaxis) par(las = 0); mtext(side = 2, paste("precip.", plotunit, sep = " "), line = lineY, cex = cexlab, adj = 1); par(las = 0); if(BOOL_Psol){ legend("bottomright", c("solid","liquid"), col = c("lightblue", "royalblue"), lty = c(1, 1), lwd = c(lwdLine, lwdLine), bty = "o", bg = bgleg, box.col = bgleg, cex = cexleg) par(new = TRUE); plot(Xaxis, PsolLayerMean[IndPeriod_Plot], type = "h", ylim = ylim2, col = "lightblue", lwd = lwdP, xaxt = "n", yaxt = "n", xlab = "", ylab = "", yaxs = "i", ...); } if(BOOL_Dates){ axis(side = 1, at = Seq1, labels = FALSE, cex.axis = cexaxis); axis(side = 1, at = Seq2, labels = Labels2, lwd.ticks = 1.5, cex.axis = cexaxis); } else { axis(side = 1, at = pretty(Xaxis), labels = pretty(Xaxis), cex.axis = cexaxis); } } ##Temp if(BOOLPLOT_Temp){ 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; } } 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 = lwdLine*0.8); } abline(h = 0, col = "grey", lty = 2) lines(SnowPackLayerMean[IndPeriod_Plot], type = "l", lwd = lwdLine*1.0, col = "darkorchid4") axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cexaxis) par(las = 0); mtext(side = 2, expression(paste("temp. [", degree, "C]", sep = "")), padj = 0.2, line = lineY, cex = cexlab); par(las = 0); legend("topright", c("mean", "layers"), col = c("darkorchid4", "orchid"), lty = c(1, 3), lwd = c(lwdLine*1.0, lwdLine*0.8), bty = "o", bg = bgleg, box.col = bgleg, cex = cexleg) box() if(BOOL_Dates){ axis(side = 1, at = Seq1, labels = FALSE, cex.axis = cexaxis); axis(side = 1, at = Seq2, labels = Labels2, lwd.ticks = 1.5, cex.axis = cexaxis); } else { axis(side = 1, at = pretty(Xaxis), labels = pretty(Xaxis), cex.axis = cexaxis); } } ##SnowPack if(BOOLPLOT_SnowPack){ 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; }