plot_OutputsModel.R

#*****************************************************************************************************************
#' Function which creates a screen plot giving an overview of the model outputs
#'
#' Dashboard of results including various graphs (depending on the model):
#' (1) time series of total precipitation and simulated flows (and observed flows if provided)
#' (2) interannual median monthly simulated flow (and observed flows if provided)
#' (3) correlation plot between simulated and observed flows (if observed flows provided)
#' (4) cumulative frequency plot for simulated flows (and observed flows if provided)
#*****************************************************************************************************************
#' @title   Default preview of model outputs
#' @author  Laurent Coron (June 2014)
## @example tests/example_plot_OutputsModel.R
#' @encoding UTF-8
#' @export
#_FunctionInputs__________________________________________________________________________________________________
#' @param  OutputsModel      [object of class \emph{OutputsModel}] list of model outputs (which must at least include DatesR, Precip and Qsim) [POSIXlt, mm, mm]
#' @param  Qobs              (optional) [numeric] time series of observed flow (for the same time-steps than simulated) [mm]
#' @param  IndPeriod_Plot    (optional) [numeric] indices of the time-steps to be plotted (among the OutputsModel series)
#' @param  BasinArea         (optional) [numeric] basin area [km2], used to plot flow axes in m3/s
#' @param  quiet             (optional) [boolean] boolean indicating if the function is run in quiet mode or not, default=FALSE
#_FunctionOutputs_________________________________________________________________________________________________
#' @return  screen plot window
#*****************************************************************************************************************'
plot_OutputsModel <- function(OutputsModel,Qobs=NULL,IndPeriod_Plot=NULL,BasinArea=NULL,quiet=FALSE){


  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(!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 ==                     60*60){ BOOL_TS <- TRUE; plotunit <- "[mm/h]";     }
    if(inherits(OutputsModel,"daily"  ) & TimeStep ==                  24*60*60){ BOOL_TS <- TRUE; plotunit <- "[mm/d]";     }
    if(inherits(OutputsModel,"monthly") & TimeStep %in% c(28,29,30,31)*24*60*60){ BOOL_TS <- TRUE; plotunit <- "[mm/month]"; }
    if(inherits(OutputsModel,"yearly" ) & TimeStep %in%     c(365,366)*24*60*60){ BOOL_TS <- TRUE; plotunit <- "[mm/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); }
  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 & !quiet){ warning("\t zeroes detected in Qobs -> some plots in the log space will not be created using all time-steps \n"); }
  if(BOOL_QsimZero & !quiet){ warning("\t zeroes detected in Qsim -> some plots in the log space will not be created using all time-steps \n"); }
  BOOL_FilterZero <- TRUE;

  ##Options
  BLOC <- TRUE; if(BLOC){
    cexaxis <- 1.0; cexlab <- 0.9; cexleg=1.0; lwd=1.8; lineX=2.6; lineY=2.6; bgleg <- rgb(1,1,1,alpha=0.7); bgleg <- NA;