R/plot2.R

In hydroGOF: Goodness-of-Fit Functions for Comparison of Simulated and Observed Hydrological Time Series

# File plot2.R
# Part of the hydroGOF R package, https://github.com/hzambran/hydroGOF ; 
#                                 https://cran.r-project.org/package=hydroGOF
#                                 http://www.rforge.net/hydroGOF/
# Copyright 2010-2024 Mauricio Zambrano-Bigiarini
# Distributed under GPL 2 or later

################################################################################
# 'plot2':     Plots 2 time series on the same graph                           #
#              It is a wrapper for the 'plot.zoo'                              #
################################################################################
# Author : Mauricio Zambrano-Bigiarini                                         #
################################################################################
# Started: March 04, 2009                                                      #
# Updates: May 2009                                                            #
#          2010                                                                #
#          21-Jan-2011 ; 15-Apr-2011 ;                                         #
#          25-Aug-2011 ; 31-Aug-2011 ; 14-Sep-2011                             #
#          23-Jan-2012                                                         #
#          15-Apr-2013 ; 15-May-2013                                           #
#          06-Aug-2017                                                         #
#          28-Dec-2022                                                         #
#          22-Mar-2024 ; 23-Mar-204 ; 08-May-2024                              #
################################################################################
                
plot2 <- function (x, y, 
                   plot.type="multiple", 
                   
                   tick.tstep="auto", 
                   lab.tstep="auto", 
                   lab.fmt=NULL,
                   
                   main, 
                   xlab="Time", 
                   ylab,
                   
                   cal.ini=NA, 
                   val.ini=NA, 
                   date.fmt="%Y-%m-%d",                   
                   
                   gof.leg= FALSE, 
                   gof.digits=2, 
                   gofs=c("ME" , "MAE" , "RMSE", "NRMSE", "PBIAS", "RSR", "rSD", 
                         "NSE", "mNSE", "rNSE",     "d",    "md",  "rd",   "r", 
                         "R2",   "bR2", "KGE" ,    "VE"),
                   
                   legend,
                   leg.cex=1,                       
                        
                   col = c("black","blue"),
                   
                   cex=c(0.5, 0.5),
                   cex.axis=1.2,
                   cex.lab=1.2,
                   
                   lwd= c(1,1), 
                   lty= c(1,3), 
                   pch= c(1,9),   
                   
                   pt.style = "ts",
                   add=FALSE,                   
                   
                    ...) {
                 
  # Saving the current graphical parameter settings,
  # and restoring them on exit
  #old.par <- par(no.readonly = TRUE)
  #on.exit(par(old.par))
  
  # Checking that the user provided 'x'
  if ( missing(x) ) stop("Missing argument: 'x'")
         
  # Checking that the user provided 'y'
  if ( missing(y) ) stop("Missing argument: 'y'")

  # Checking 'gofs'.  'rSpearman' and 'pbiasFDC' are not computed
  gofs.all=c(   "ME",    "MAE",    "MSE",  "RMSE", "ubRMSE", 
             "NRMSE",  "PBIAS",   "RSR",    "rSD",    "NSE",  
             "mNSE" ,   "rNSE",  "wNSE",  "wsNSE",      "d",     
                "dr",     "md",    "rd",     "cp",      "r",     
                "R2",    "bR2",    "VE",    "KGE",  "KGElf",  
             "KGEnp",  "KGEkm",  "sKGE",   "APFB",    "HFB")  

  # Removing 'sKGE', 'APFB' and 'HFB' when 'x' and 'y' are not zoo objects
  if ( !( zoo::is.zoo(x) & zoo::is.zoo(y) ) )
    gofs.all <- gofs.all[ -c( (length(gofs.all)-2):(length(gofs.all)) ) ]

  # Checking 'gofs' 
  noNms.index <- which( !(gofs %in% gofs.all) )
  if (length(noNms.index) > 0) {
    noNms <- gofs[noNms.index]
    warning("[Unknown names in 'gofs': ", paste(noNms, collapse = ", "), " (not used) !]")
  } # IF end	   
  gofs.index <- which( (gofs %in% gofs.all) )
  gofs       <- gofs[gofs.index]

  # 'xname' and 'yname' values
  xname <- deparse(substitute(x))
  yname <- deparse(substitute(y))

  # 'legend' value
  if (missing(legend)) legend <- c(xname, yname)

  # 'ylab' value
  if (missing(ylab)) ylab <- c(xname, yname)
  
  # Checking that the user provided a valid argument for 'x' and 'y'
  valid.class <- c("xts", "zoo", "ts", "numeric", "integer")       
  if (length(which(!is.na(match(class(x), valid.class )))) <= 0) 
         stop("Invalid argument: 'class(x)' must be in c('integer', 'numeric', 'ts', 'zoo', 'xts')")
  if (length(which(!is.na(match(class(y), valid.class )))) <= 0)
         stop("Invalid argument: 'class(y)' must be in c('integer', 'numeric', 'ts', 'zoo', 'xts')")
         
  # Checking that the user provided a valid argument for 'plot.type'       
  if (is.na(match(plot.type, c("single", "multiple") ) ) ) 
         stop("Invalid argument: 'plot.type' must be in c('single', 'multiple')")
         
  # If the user wants to draw a legend, it checks that the type of plot is 'single'
  if (gof.leg & (plot.type == "multiple") )
    stop("Invalid argument: For drawing a legend, 'plot.type' must be 'single'")
         
  # Checking that the user provided a valid argument for 'pt.style'       
  if (is.na(match(pt.style, c("ts", "bar") ) ) ) 
         stop("Invalid argument: 'pt.style' must be in c('ts', 'bar')")
         
  # If 'x' is 'ts' or 'zoo' and 'y' is only a vector, y is transformed into 
  # the same class of 'x', with the same times
  if ( (length(which(!is.na(match(class(x), c("ts", "zoo", "xts") )))) <= 0) & 
       (length(which(!is.na(match(class(y), c("integer", "numeric") )))) <= 0) ) {  
  
    # class(time(x))== "Date" for 'daily' and 'monthly' time series
    # class(time(x))== "character" for 'annual' time series
    if ( inherits(time(x), "Date" ) ) {
        y <- vector2zoo(y, dates=time(x)) # hydroTSM::vector2zoo
    } else if ( inherits(time(x), "character") ) {
        y <- vector2zoo(y, dates=time(x), date.fmt="%Y") # hydroTSM::vector2zoo
        time(x) <- time(y) #'annual' time series
    } # ELSE END
    
  } # IF END

         
  # Checking that the user provied the same length for 'sim' and 'obs'      
  #if ( length(x) != length(y) )  
  #       stop("Invalid argument: 'obs' and 'sim' must have the same length")
         
  # If the user didn't provide a title for the plot, the default is used 
  if ( missing(main) ) main <- "Observed vs Simulated"   
  
  # If the user provided a value for 'cal.ini', it is transformed into a Date class
  if ( !missing(cal.ini) ) cal.ini <- as.Date(cal.ini, format=date.fmt)
    
  # If the user provided a value for 'val.ini', it is transformed into a Date class
  if ( !missing(val.ini) ) val.ini <- as.Date(val.ini, format=date.fmt)
  
  # If the legend has to be plotted AND no other plots will be added
  # IF 'add' is TRUE, the layout of the screen is set up by the calling procedure (usually 'ggof')
  if (gof.leg & add==FALSE) {  
            def.par <- par(no.readonly = TRUE) # save default, for resetting...     
            #par(mar=c(5, 2, 4, 0.5) + 0.1)
            # Setting up the screen with 1 rows and 2 columns
            layout( matrix( c(1,1,1,1,1,1,1,1,1,2,2), ncol=11, byrow=TRUE) ) 
            par(mar=c(5, 4, 4, 0) + 0.1)    
            on.exit(par(def.par))      
  } # ELSE end    
  
  # If the legend will not be plotted, the marginns are set to 'almost' the default values
  if (!gof.leg)
    par(mar=c(5, 4, 4, 2) + 0.1) # default values are par(mar=c(5, 4, 4, 4) + 0.1)
  
  # If the plot type is "time series"
  if (pt.style=="ts") {
  
    # If both time series have to be ploted in the same plot area
    if (plot.type == "single") {

      if (length(which(!is.na(match(class(x), c("ts", "zoo", "xts") )))) > 0) x <- zoo::as.zoo(x) 
      if (length(which(!is.na(match(class(y), c("ts", "zoo", "xts") )))) > 0) y <- zoo::as.zoo(y) 

      # Y axis limits
      if (!hasArg(ylim)) ylim <- range( range(x[is.finite(x)]), range(y[is.finite(y)]) )

      # Plotting the Observed Time Series. 
      # It calls to 'plot' or 'plot.zoo'  depending on 'x' class
      if (!hasArg(ylim)) {
         zoo::plot.zoo(x, xaxt = "n", yaxt = "n",
                       type="o", lwd=lwd[1], lty= lty[1], col= col[1], pch= pch[1], 
                       cex = cex[1], cex.axis=cex.axis, cex.lab=cex.lab,
                       main=main, xlab=xlab, ylab= ylab, ylim=ylim, ... )
      } else zoo::plot.zoo(x, xaxt = "n", yaxt = "n", 
                      type="o", lwd=lwd[1], lty= lty[1], col= col[1], pch= pch[1], 
                      cex = cex[1], cex.axis=cex.axis, cex.lab=cex.lab,
                      main=main, xlab=xlab, ylab= ylab, ... )
      axis(2, cex.axis=cex.axis, cex.lab=cex.lab)
     
      lines(y, type="o", lwd=lwd[2], lty= lty[2], col= col[2], pch= pch[2], cex = cex[2])      
               
      # If the user provided a value for 'cal.ini', a vertical line is drawn
      if ( !missing(cal.ini) ) abline(v=as.POSIXct(cal.ini), col="red", lty=1, lwd=2)
      
      # If the user provided a value for 'val.ini', a vertical line is drawn
      if ( !missing(val.ini) ) abline(v=as.POSIXct(val.ini), col="red", lty=1, lwd=2)
               
      # Drawing a legend with 'Obs' vs 'Sim' 
      # y.intersp=0.5, is for the vertical spacin in the legend
      # bty="n" => no box around the legend
      # 'inset=0.03' is usefult when plot.type= "multiple" for having a nice margin to the legend
      legend("topright", legend=legend, y.intersp=0.8, inset=0.03,
             bty="n", cex = leg.cex, col = col, lwd= lwd, lty= lty, pch=pch )  
      
      # Drawing the 'x' axis
      # If the user provided, in some way, valid values for being used as dates, 
      # they will be used, if not, only a numeric index will be used
      if ( (length(which(!is.na(match(class(x), c("ts", "zoo", "xts") )))) > 0) | 
           (length(which(!is.na(match(class(y), c("ts", "zoo", "xts") )))) > 0) ) {
    
        if ( (length(which(!is.na(match(class(x), c("ts", "zoo", "xts") )))) > 0) ) { 
          z <- x
        } else z <- y
    
        # Draws ticks in the X axis, but labels only in years
        drawTimeAxis(z, tick.tstep=tick.tstep, lab.tstep= lab.tstep, lab.fmt=lab.fmt, cex.axis=cex.axis, cex.lab=cex.lab) # hydroTSM::drawTimeAxis

        # manually adding a grid
        grid(nx=NA, ny=NULL)
        abline(v=time(x)[axTicksByTime(x)], col = "lightgray", lty = "dotted")
      
      } else { # When 'numeric' or 'integer' values (not 'zoo' or 'xts') are plotted
             Axis(side = 1, labels = TRUE, cex.axis=cex.axis, cex.lab=cex.lab)
             box()
             # manually adding a grid
             grid()
             }
               
    } else  #plot.type == "multiple"  
          {       
            # all the work (mainly Time axis) is made automatically be the 'plot.zoo' function 
            zoo::plot.zoo(cbind(x, y), plot.type=plot.type, type=c("o","o"), 
                         lwd=lwd, lty= lty, col= col, pch= pch, 
                         cex = cex, 
                         cex.axis=cex.axis, 
                         cex.lab=cex.lab,
                         main=main, xlab=xlab, ylab= ylab,...)
                         
            # If the user provided a value for 'cal.ini', a vertical line is drawn
            if ( !missing(cal.ini) ) abline(v=as.POSIXct(cal.ini), col="red", lty=1, lwd=2)
      
            # If the user provided a value for 'val.ini', a vertical line is drawn
            if ( !missing(val.ini) ) abline(v=as.POSIXct(val.ini), col="red", lty=1, lwd=2)
                         
      } # ELSE end 
      
  } else if (pt.style=="bar") {
    
        # Creation of the table that will be plotted as barplot
        b <- rbind(coredata(x),coredata(y))
        
        # Giving the as name to each bar the YEAR, because the 
        # bar plot is thought for being used ONLY for annual time series
        colnames(b) <- format( time(x), "%Y")
        
        # Barplot  
        barplot(b, beside=TRUE, axis.lty=1, col=col, density=25, angle=c(45,-45), 
                main=main, xlab=xlab, ylab= ylab, legend.text=legend, 
                args.legend=list(x="topleft", y.intersp=0.8, inset=0.00, bty="n", 
                                 cex= cex.axis, col = col), # lwd= lwd, lty= lty, pch=pch),
                cex.axis=cex.axis, cex.lab=cex.lab, cex = cex.lab, ...)  
       
        if ( !missing(cal.ini) ) {
          # Index of the bar corresponding to 'cal.ini'.
          # It is necessary to multiply it by 3 because for each year there are 3 vertical lines
          # It is necessary to substract 2, for shifting the line form the 3 line to the first one
          cal.index <- 3*which(colnames(b) == format( cal.ini, "%Y")) - 2
          # If the user provided a value for 'cal.ini', a vertical line is drawn
          if ( !missing(cal.ini) ) {
           abline(v=cal.index, col="red", lty=1, lwd=2)
          } # IF end
        } # IF end
       
        if ( !missing(val.ini) ) {
          # Index of the bar corresponding to 'val.ini'.
          # It is necessary to multiply it by 3 because for each year there are 3 vertical lines
          # It is necessary to substract 2, for shifting the line form the 3 line to the first one
          val.index <- 3*which(colnames(b) == format( val.ini, "%Y")) - 2
          # If the user provided a value for 'val.ini', a vertical line is drawn
          if ( !missing(val.ini) ) {
            abline(v=val.index, col="red", lty=1, lwd=2)
          } # IF end
        } # IF end
      
    }  # ELSE end        
  
  
  # If the Goodness-of-fit indexes have to be computed and plotted:
  if (gof.leg & plot.type == "single" ) {
  
   # If the user provided 'cal.ini' and 'x' & 'y' are zoo objects, 
   # the warming up period is removed from the computation of the goodness-of-fit
   # measures => all the values before 'cal.ini' are ignored but plotted
   if ( is.zoo(x) & is.zoo(y) ) {
     if ( !is.na(cal.ini) ) {
       x <- window(x, start=cal.ini)
       y <- window(y, start=cal.ini)
     } # IF end 
   } # IF end

   gof.index <- pmatch(gofs, gofs.all)
   gof.index <- gof.index[!is.na(gof.index)]  
   
   #gof.xy <- gof(sim=as.numeric(x), obs=as.numeric(y), do.spearman=FALSE, do.pbfdc=FALSE, digits=gof.digits, ...)
   gof.xy <- gof(sim=x, obs=y, do.spearman=FALSE, do.pbfdc=FALSE, digits=gof.digits, ...)

   gofs.stg  <- gofs.all[gof.index] 
   gofs.num  <- gof.xy[gof.index, 1] 

   leg.text <- paste(gofs.stg, " = ", gofs.num, sep="")
   
   legend.position <- "center"
   par( mar=c(0.5, 0.5, 0.5, 0.5) ) # mar=c(bottom, left, top, right). Default values are: mar=c(5,4,4,2) + 0.1)
   plot.new() 
          
   # 'inset':  The optional 'inset' argument specifies how far the legend is inset from the plot margins.  
   #           If a single value is given, it is used for both margins; 
   #           if two values are given, the first is used for 'x'-distance, the second for 'y'-distance.
	 
   legend(legend.position,  y.intersp=1.2, cex =leg.cex, # bty="n", #inset=0.01,  
          leg.text, 
#          c( paste( "ME =", gof.xy["ME", 1], sep=" "),
#             paste( "MAE =", gof.xy["MAE", 1], sep=" "),
#             #paste( "MSE =", gof.xy["MSE", 1], sep=" "),
#             paste( "RMSE =", gof.xy["RMSE", 1], sep=" "),
#             paste( "NRMSE% =", gof.xy["NRMSE %", 1], sep=" "),
#             paste( "PBIAS% =", gof.xy["PBIAS %", 1], sep=" "),
#             #paste( "pbiasFDC% =", gof.xy["pbiasFDC %", 1], sep=" "),
#             paste( "RSR =", gof.xy["RSR", 1], sep=" "),
#             paste( "rSD =", gof.xy["rSD", 1], sep=" "),             
#             paste( "NSE =", gof.xy["NSE", 1], sep=" "),
#             paste( "mNSE =", gof.xy["mNSE", 1], sep=" "),
#             paste( "rNSE =", gof.xy["rNSE", 1], sep=" "),
#             paste( "d =", gof.xy["d", 1], sep=" "),
#             paste( "md =", gof.xy["md", 1], sep=" "),
#             paste( "rd =", gof.xy["rd", 1], sep=" "),
#             #paste( "cp =", gof.xy["cp", 1], sep=" "),
#             paste( "r =", gof.xy["r", 1], sep=" "),
#             paste( "R2 =", gof.xy["R2", 1], sep=" "), 
#             paste( "bR2 =", gof.xy["bR2", 1], sep=" "),
#             paste( "KGE =", gof.xy["KGE", 1], sep=" "), 
#             paste( "VE =", gof.xy["VE", 1], sep=" ")               
#            ), 
            title="GoF's:", title.col="darkblue",
            bg="azure"
           )
         
  } #IF END
  
} # 'plot2' end