R/plot_out.R

In hydroMOPSO: Multi-Objective Optimisation with Focus on Environmental Models

################################################################################
##                                plot_out                                    ##
################################################################################
# Author : Rodrigo Marinao Rivas                                              ##
################################################################################
# Created: 2024/08/31                                                         ##
################################################################################
# Function: plot_out
# Description: A function that generates plots for hydrological model outputs,
#              including observed data, model uncertainty bands, and the best 
#              compromise solution. This function is used for visualizing the 
#              calibration and validation performance of a model by comparing 
#              model outputs to observed data across various time periods and 
#              scenarios.
################################################################################

plot_out <- function(
     Results,                # (list) Object containing preprocessed 
                             # hydrological results.
     model.out = NULL,       # (list or NULL) Output from the 
                             # hydrological model used for evaluation.
     analysis.period = NULL, # (character or NULL) Time period for analysis
                             # (e.g., "calibration" or "verification").
     model.out.bcs = NULL,   # (list or NULL) Model output representing the 
                             # "Best Compromise Solution" (BCS).
     bcs = NULL,             # (matrix or NULL) Parameters or results 
                             # corresponding to the best compromise solution
                             # (BCS).
     obs.var = NULL,         # (list or NULL) Observed variables to be compared
                             # against the model outputs.
     dimensions = NULL,      # (matrix or NULL) Dimensions of the modeled
                             # problem, typically a vector 
                             # indicating the number of objectives and
                             # variables.
     obj.names = NULL,       # (character or NULL) Names of the objectives.
     dates.cal = NULL,       # (Date or NULL) Dates of the calibration period.
     dates.warmup = NULL,    # (Date or NULL) Dates of the model's warm-up 
                             # period.
     var.names = NULL,       # (character or NULL) Names of the modeled and 
                             # observed variables.
     var.units = NULL,       # (character or NULL) Units of the modeled and 
                             # observed variables.
     xlim = NULL,            # (numeric or NULL) Limits for the x-axis in the 
                             # plots.
     ylim = NULL,            # (numeric or NULL) Limits for the y-axis in the 
                             # plots.
     digits = 4,             # (integer) Number of digits to use for rounding 
                             # values in plots and legends.
     col.band = "skyblue",   # (character) Color used for the model uncertainty 
                             # bands in the plots.
     col.bcs = "mediumblue", # (character) Color used for the line representing 
                             # the best compromise solution (BCS).
     col.obs = "black",      # (character) Color used for the observed variable 
                             # lines in the plots.
     lwd = 0.75,             # (numeric) Line width used in the plots for model 
                             # and observation lines.
     pch.bcs = 15,           # (integer) Symbol type for points in the plot 
                             # representing the best compromise solution (BCS).
     pch.obs = 15,           # (integer) Symbol type for points in the plot 
                             # representing the observations.
     main = "study case #1", # (character) Main title for the plot.
     drty.out = "MOPSO.out", # (character) Output directory where plots will be 
                             # saved if specified to save as PNG files.
     cex.pt = 0.25,          # (numeric) Size of points in the plots.
     cex.main = 1,           # (numeric) Size of the main title text in plot.
     cex.lab = 1,            # (numeric) Size of the axis label text in plots.
     cex.axis = 1,           # (numeric) Size of the axis values text in plots.
     do.png = FALSE,         # (logical) Boolean value indicating whether the 
                             # plots should be saved as PNG files.
     legend.obs = "Observation",
                             # (character) Legend text for observations.
     legend.bcs = "Best compromise solution",
                             # (character) Legend text for the best compromise
                             # solution (BCS).
     legend.band = "Pareto front bands"
                             # (character) Legend text for Pareto front bands.
){

    
    oldpar <- par(no.readonly = TRUE)
    on.exit(par(oldpar))

    if(!is.null(Results[["hydroResults"]])){
      Results <- Results[["hydroResults"]]
    }


    results.names <- c("ModelOut", "AnalysisPeriod", "ModelOutBestCS",
                       "ParticleBestCS", "Obs", "Dimensions", "DatesCal",
                       "WarmUp", "VarNames", "VarUnits")

    if(!is.null(Results)){
        if(all(results.names %in% names(Results))){
            flag.results <- TRUE
        }else{
            flag.results <- FALSE
        }
    }else{
        flag.results <- FALSE
    }

    #

    out.flag <- TRUE
    bcs.flag <- TRUE
    obs.flag <- TRUE

    if(is.null(model.out)){
        if(flag.results){
            model.out <- Results[["ModelOut"]]
        }else{
            out.flag <- FALSE
        }
    }

    if(is.null(model.out.bcs)){
        if(flag.results){
            model.out.bcs <- Results[["ModelOutBestCS"]]        
        }else{
            bcs.flag <- FALSE
        }
    }

    if(is.null(bcs)){
        if(flag.results){
            bcs <- Results[["ParticleBestCS"]]        
        }else{
            bcs.flag <- FALSE
        }
    }


    if(is.null(obs.var)){
        if(flag.results){
            obs.var <- Results[["Obs"]]        
        }else{
            obs.flag <- FALSE
        }
    }


    if(sum(c(out.flag, bcs.flag, obs.flag)) < 2){
        stop("It is necessary to have at least two of the following arguments: 
              1) 'model.out', 2) 'bcs' and 'model.out.bcs', 3) 'obs.var'")
    }



    if(is.null(analysis.period)){
        if(flag.results){
            analysis.period <- Results[["AnalysisPeriod"]]        
        }else{
            stop("The 'analysis.period' argument was not provided or included
                  in the 'Results' argument. Please make sure to include all
                  required arguments, and if using the 'Results' argument,
                  ensure 'AnalysisPeriod' is correctly specified.")
        }
    }

    if(is.null(dimensions)){
        if(flag.results){
            dimensions <- Results[["Dimensions"]]   
        }else{
            stop("The 'dimensions' argument was not provided or included in the
                  'Results' argument. Please make sure to include all required
                  arguments, and if using the 'Results' argument, ensure
                  'Dimensions' is correctly specified.")
        }
    }

    if(is.null(dates.cal)){
        if(flag.results){
            dates.cal <- Results[["DatesCal"]]   
        }else{
            stop("The 'dates.cal' argument was not provided or included in the'
                  ' argument. Please make sure to include all required 
                  arguments, and if using the 'Results' argument, ensure
                  'DatesCal' is correctly specified.")
        }
    }

    if(is.null(dates.warmup)){
        if(flag.results){
            dates.warmup <- Results[["WarmUp"]]   
        }else{
            stop("The 'dates.warmup' argument was not provided or included in
                  the 'Results' argument. Please make sure to include all
                  required arguments, and if using the 'Results' argument, 
                  ensure 'WarmUp' is correctly specified.")
        }
    }

    if(is.null(var.names)){
        if(flag.results){
            var.names <- Results[["VarNames"]]   
        }else{
            stop("The 'var.names' argument was not provided or included in the
                  'Results' argument. Please make sure to include all required
                  arguments, and if using the 'Results' argument, ensure 
                  'VarNames' is correctly specified.")
        }
    }

    if(is.null(var.units)){
        if(flag.results){
            var.units <- Results[["VarUnits"]]   
        }else{
            stop("The 'var.units' argument was not provided or included in the
                  'Results' argument. Please make sure to include all required
                  arguments, and if using the 'Results' argument, ensure 
                  'VarUnits' is correctly specified.")
        }
    }


    if(out.flag & bcs.flag){

        ModelOutBandVsBCS(model.out = model.out,
                          model.out.bcs = model.out.bcs,
                          analysis.period = analysis.period,
                          bcs = bcs,
                          dimensions = dimensions,
                          dates.cal = dates.cal,
                          dates.warmup = dates.warmup,
                          var.names = var.names,
                          var.units = var.units,
                          xlim = xlim,
                          ylim = ylim,
                          col.band = col.band,
                          col.bcs = col.bcs,
                          lwd = lwd,
                          pch.bcs = pch.bcs,
                          main = main,
                          drty.out = drty.out,
                          legend.bcs = legend.bcs, 
                          legend.band = legend.band,
                          cex.pt = cex.pt, 
                          cex.main = cex.main,
                          cex.lab = cex.lab,
                          cex.axis = cex.axis,
                          do.png = do.png)


    }



    if(out.flag & obs.flag){

        ModelOutBandVsObs(model.out = model.out,
                          obs.var = obs.var,
                          analysis.period = analysis.period,
                          dimensions = dimensions,
                          dates.cal = dates.cal,
                          dates.warmup = dates.warmup,
                          var.names = var.names,
                          var.units = var.units,
                          xlim = xlim,
                          ylim = ylim,
                          col.band = col.band,
                          col.obs = col.obs,
                          lwd = lwd,
                          pch.obs = pch.obs,
                          main = main,
                          drty.out = drty.out,
                          legend.band = legend.band,
                          legend.obs = legend.obs,
                          cex.pt = cex.pt, 
                          cex.main = cex.main,
                          cex.lab = cex.lab,
                          cex.axis = cex.axis,
                          do.png = do.png)

    }

    if(out.flag & obs.flag){

        ModelOutBCSVsObs(model.out.bcs = model.out.bcs,
                         bcs = bcs,
                         obs.var = obs.var,
                         analysis.period = analysis.period,
                         dimensions = dimensions,
                         obj.names = obj.names,
                         dates.cal = dates.cal,
                         dates.warmup = dates.warmup,
                         var.names = var.names,
                         var.units = var.units,
                         digits = digits,
                         xlim = xlim,
                         ylim = ylim,
                         col.obs = col.obs,
                         col.bcs = col.bcs,
                         lwd = lwd,
                         pch.bcs = pch.bcs,
                         main = main,
                         drty.out = drty.out,
                         legend.bcs = legend.bcs,
                         legend.obs = legend.obs,
                         cex.pt = cex.pt, 
                         cex.main = cex.main,
                         cex.lab = cex.lab,
                         cex.axis = cex.axis,
                         do.png = do.png)

    }


}








#########

ModelOutBandVsBCS <- function(model.out,
                              analysis.period,
                              model.out.bcs,
                              bcs,
                              dimensions,
                              dates.cal,
                              dates.warmup,
                              var.names,
                              var.units,
                              xlim = NULL,
                              ylim = NULL,
                              col.band = "skyblue",
                              col.bcs = "mediumblue",
                              lwd = 0.75,
                              pch.bcs = 15,
                              main = "study case #1",
                              drty.out = "MOPSO.out",
                              legend.bcs = "Best compromise solution",
                              legend.band = "Pareto front bands",
                              cex.pt = 0.25, 
                              cex.main = 1,
                              cex.lab = 1,
                              cex.axis = 1,
                              do.png = FALSE
                              ){

  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))

  if(do.png){
    if(!dir.exists(paste0(drty.out, "/", analysis.period, "/png.out"))){
      dir.create(paste0(drty.out, "/", analysis.period, "/png.out"),
                 recursive = TRUE)
    }
  }

  nobjs <- dimensions[1,1]
  nvars <- dimensions[1,2]

  dates <- as.Date(time(model.out[[1]]))

  obj.bcs <- bcs[1,c(2:(1+nobjs))]

  if(!is.null(dates.warmup)){
    warmup <- as.Date(dates.warmup[,1])
  }else{
    warmup <- NULL
  }

  dates.cal <- as.Date(dates.cal[,1])


  nplots<- nvars

  if(!is.null(warmup)){
    flag.warmup <- TRUE
  }else{
    flag.warmup <- FALSE
    warning("Warm-up not specified")
  }

  if(analysis.period == "verification"){
    flag.cal <- TRUE
  }else{
    flag.cal <- FALSE
  }


  # Model Out ##################################################################

  if(do.png){
    png(filename = paste0(drty.out, "/", analysis.period,
                          "/png.out/ModelOut_from_Pareto_Optimal_Front.png"),
        width = 3840, height = nplots*1620, res = 280)
  }
  
    
  par(mfrow = c(nplots, 1), oma = c(0,0,0,0))

  
  for(i in 1:length(model.out)){
    
    var.name <- var.names[i]
    var.unit <- var.units[i]
    
    varmin <- apply(coredata(model.out[[i]]), MARGIN = 1, FUN = min)
    varmax <- apply(coredata(model.out[[i]]), MARGIN = 1, FUN = max)
    

    if(is.null(ylim)){
      ylim <- c(min(varmin, na.rm = TRUE),
                max(varmax, na.rm = TRUE))
    }

    if(is.null(xlim)){
      xlim <- c(min(dates), max(dates)) 
    }
    

    plot(zoo(varmax, dates), type = "n", xaxt = "n", xlab = "Date",
         ylab = paste0(var.name, " [",var.unit,"]"),
         cex.main = cex.main, cex.lab = cex.lab, cex.axis = cex.axis,
         main = "", ylim = ylim, xlim = xlim)
    title(main = paste0(var.name, " in ", main), adj = 0)
    if(flag.warmup){
      polygon(x = c(warmup, rev(warmup)),
              y = c(rep(ylim[1]-0.05*(ylim[2]-ylim[1]), length(warmup)),
                    rep(ylim[2]+1.05*(ylim[2]-ylim[1]), length(warmup))),
              border = NA, col = "#DBDBDB")
    }

    if(flag.cal){
      polygon(x = c(dates.cal, rev(dates.cal)),
              y = c(rep(ylim[1]-0.05*(ylim[2]-ylim[1]), length(dates.cal)),
                    rep(ylim[2]+1.05*(ylim[2]-ylim[1]), length(dates.cal))),
              border = NA, col = "#DCEDFF")
    }

    grid()
    if(flag.warmup){
      text(x = xlim[1], y = ylim[2], adj = c(1,1), label = "Warm-Up", srt = 90,
           col = "#666666", font = 2)
    }
    if(flag.cal){
      text(x = dates.cal[1], y = ylim[2], adj = c(0,1), label = "Calibration",
           srt = 0, col = "#0072EA", font = 2)
    }
    box()
    drawTimeAxis(zoo(varmax, dates))
    PlotBandsOnly(lband = zoo(varmin, dates), uband = zoo(varmax, dates),
                  bands.col = col.band)
    lines(zoo(varmin, dates), col = col.band, lwd = 0.5)
    lines(zoo(varmax, dates), col = col.band, lwd = 0.5)
    lines(model.out.bcs[[i]], col = col.bcs, lwd = lwd, type = "o",
          pch = pch.bcs, cex = cex.pt)

    par(xpd = TRUE)

    legend(x = xlim[1] + 0.80*(xlim[2] - xlim[1]),  xjust = 0,
           y = ylim[2] + 0.05*(ylim[2] - ylim[1]), yjust = 0,
           legend = c(legend.bcs, legend.band), 
           col = c(col.bcs, col.band), 
           lty = c(1,0),
           lwd = lwd,
           pch = c(pch.bcs,15),
           bty = "n",
           pt.cex = c(cex.pt,1), 
           cex = 1, 
           text.col = "black", 
           horiz = FALSE, 
           inset = c(0.1, 0.1))

    par(xpd = FALSE)
  }
  
  if(do.png){
    dev.off()
  }
  
}





#############3


ModelOutBandVsObs <- function(model.out,
                              analysis.period,
                              obs.var,
                              dimensions,
                              dates.cal,
                              dates.warmup,
                              var.names,
                              var.units,
                              xlim = NULL,
                              ylim = NULL,
                              col.band = "skyblue",
                              col.obs = "black",
                              lwd = 0.75,
                              pch.obs = 15,
                              main = "study case #1",
                              drty.out = "MOPSO.out",
                              legend.band = "Pareto front bands",
                              legend.obs = "Observation",
                              cex.pt = 0.25, 
                              cex.main = 1,
                              cex.lab = 1,
                              cex.axis = 1,
                              do.png = FALSE
                              ){

  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))

  if(do.png){
    if(!dir.exists(paste0(drty.out, "/", analysis.period, "/png.out"))){
      dir.create(paste0(drty.out, "/", analysis.period, "/png.out"),
                 recursive = TRUE)
    }
  }


  nvars <- dimensions[1,2]
  
  dates <- as.Date(time(model.out[[1]]))

  if(!is.null(dates.warmup)){
    warmup <- as.Date(dates.warmup[,1])
  }else{
    warmup <- NULL
  }

  dates.cal <- as.Date(dates.cal[,1])


  nplots<- nvars

  if(!is.null(warmup)){
    flag.warmup <- TRUE
  }else{
    flag.warmup <- FALSE
    warning("Warm-up not specified")
  }

  if(analysis.period == "verification"){
    flag.cal <- TRUE
  }else{
    flag.cal <- FALSE
  }






  # Obs vs Model Out ===========================================================
  
  if(do.png){
    png(filename=paste0(drty.out, "/", analysis.period,
                      "/png.out/ModelOut_from_Pareto_Optimal_Front_vs_Obs.png"),
                        width = 3840, height = nplots*1620, res = 280)
  }
  

  par(mfrow = c(nplots, 1), oma = c(0, 0, 0, 0))


  
  for(i in 1:length(model.out)){
    
    var.name <- var.names[i]
    var.unit <- var.units[i]
    
    varmin <- apply(coredata(model.out[[i]]), MARGIN = 1, FUN = min)
    varmax <- apply(coredata(model.out[[i]]), MARGIN = 1, FUN = max)



    if(is.null(ylim)){
      ylim <- c(min(c(varmin, coredata(obs.var[[i]])), na.rm = TRUE),
                max(c(varmax, coredata(obs.var[[i]])), na.rm = TRUE))
    }

    if(is.null(xlim)){
      xlim <- c(min(dates), max(dates)) 
    }
    

    
    plot(zoo(varmax, dates), type = "n", xaxt = "n", xlab = "Date",
         ylab = paste0(var.name, " [",var.unit,"]"),
         cex.main = cex.main, cex.lab = cex.lab, cex.axis = cex.axis,
         main = "", ylim = ylim, xlim = xlim)
    title(main = paste0(var.name, " in ", main), adj = 0)
    if(flag.warmup){
      polygon(x = c(warmup, rev(warmup)),
              y = c(rep(ylim[1]-0.05*(ylim[2]-ylim[1]), length(warmup)),
                    rep(ylim[2]+1.05*(ylim[2]-ylim[1]), length(warmup))),
              border = NA, col = "#DBDBDB")
    }
    if(flag.cal){
      polygon(x = c(dates.cal, rev(dates.cal)),
              y = c(rep(ylim[1]-0.05*(ylim[2]-ylim[1]), length(dates.cal)),
                    rep(ylim[2]+1.05*(ylim[2]-ylim[1]), length(dates.cal))),
              border = NA, col = "#DCEDFF")
    }
    grid()
    if(flag.warmup){
      text(x = xlim[1], y = ylim[2], adj = c(1,1), label = "Warm-Up", srt = 90,
           col = "#666666", font = 2)
    }
    if(flag.cal){
      text(x = dates.cal[1], y = ylim[2], adj = c(0,1), label = "Calibration",
           srt = 0, col = "#0072EA", font = 2)
    }
    box()
    drawTimeAxis(zoo(varmax, dates))
    PlotBandsOnly(lband = zoo(varmin, dates), uband = zoo(varmax, dates),
                  bands.col = col.band)
    lines(zoo(varmin, dates), col = col.band, lwd = 0.5)
    lines(zoo(varmax, dates), col = col.band, lwd = 0.5)
    lines(obs.var[[i]], col = col.obs, lwd = lwd, type = "o", pch = pch.obs,
          cex = cex.pt)

    pf <- pfactor(x=obs.var[[i]], lband=zoo(varmin, dates),
                  uband=zoo(varmax, dates), na.rm=TRUE)
    rf <- rfactor(x=obs.var[[i]], lband=zoo(varmin, dates),
                  uband=zoo(varmax, dates), na.rm=TRUE)

    par(xpd = TRUE)
    
    legend(x = xlim[1] + 0.80*(xlim[2] - xlim[1]),  xjust = 0,
           y = ylim[2] + 0.05*(ylim[2] - ylim[1]), yjust = 0,
           legend = c(legend.obs, legend.band), 
           col = c(col.obs, col.band), 
           lty = c(1,0),
           lwd = lwd,
           pch = c(pch.obs,15),
           bty = "n",
           pt.cex = c(cex.pt,1), 
           cex = 1, 
           text.col = "black", 
           horiz = FALSE, 
           inset = c(0.1, 0.1))

    legend("topright", legend=c(paste0("P-factor: ", round(pf, 2)), 
                               paste0("R-factor: ", round(rf, 2))),
           bty="n", cex = 1)

    par(xpd = FALSE)
  }
  
  if(do.png){
    dev.off()
  }


}




#########



ModelOutBCSVsObs <- function(model.out.bcs,
                             bcs,
                             obj.names,
                             obs.var,
                             analysis.period,
                             dimensions,
                             dates.cal,
                             dates.warmup,
                             var.names,
                             var.units,
                             digits = 4,
                             xlim = NULL,
                             ylim = NULL,
                             col.obs = "black",
                             col.bcs = "mediumblue",
                             lwd = 0.75,
                             pch.bcs = 15,
                             main = "study case #1",
                             drty.out = "MOPSO.out",
                             legend.bcs = "Best compromise solution",
                             legend.obs = "Observation",
                             cex.pt = 0.25, 
                             cex.main = 1,
                             cex.lab = 1,
                             cex.axis = 1,
                             do.png = FALSE
                             ){

  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))

  if(do.png){
    if(!dir.exists(paste0(drty.out, "/", analysis.period, "/png.out"))){
      dir.create(paste0(drty.out, "/", analysis.period, "/png.out"),
                 recursive = TRUE)
    }
  }

  nobjs <- dimensions[1,1]
  nvars <- dimensions[1,2]

  dates <- as.Date(time(model.out.bcs[[1]]))

  obj.bcs <- bcs[1,c(2:(1+nobjs))]

  if(is.null(obj.names)){
    obj.names <- colnames(obj.bcs)
  }

  if(!is.null(dates.warmup)){
    warmup <- as.Date(dates.warmup[,1])
  }else{
    warmup <- NULL
  }

  dates.cal <- as.Date(dates.cal[,1])

  nplots<- nvars

  if(!is.null(warmup)){
    flag.warmup <- TRUE
  }else{
    flag.warmup <- FALSE
    warning("Warm-up not specified")
  }

  if(analysis.period == "verification"){
    flag.cal <- TRUE
  }else{
    flag.cal <- FALSE
  }


  # Model Out BCS vs Obs #######################################################

  if(do.png){
    png(filename=paste0(drty.out, "/", analysis.period,
                  "/png.out/ModelOut_BCS_from_Pareto_Optimal_Front_vs_Obs.png"),
                        width = 3840, height = nplots*1620, res = 280)
  }

  SimVsObs(sim = model.out.bcs,
           obs = obs.var,
           obj.names = obj.names,
           obj.values = obj.bcs,
           var.names = var.names,
           var.units = var.units,
           legend.sim = legend.bcs,
           legend.obs = legend.obs,
           xlim = xlim,
           ylim = ylim,
           full.period = dates,
           cal.period = dates.cal,
           warmup.period = warmup,
           main = main,
           analysis.period = analysis.period,
           digits.round = digits,
           cex.main = cex.main,
           cex.lab = cex.lab,
           cex.axis = cex.axis)
  
  if(do.png){
    dev.off()
  }

  
}




###############



SimVsObs <- function(sim,
                     obs,
                     obj.values,
                     obj.names,
                     var.names, 
                     var.units,
                     legend.sim = "Simulated", 
                     legend.obs = "Observed", 
                     cal.period,
                     warmup.period,
                     full.period,
                     xlim = NULL,
                     ylim = NULL,
                     main = "study case #1",
                     analysis.period,
                     digits.round = 8,
                     col.obs = "black",
                     col.sim = "mediumblue",
                     lwd = 0.75,
                     pch.obs = 15,
                     cex.pt = 0.25,
                     cex.main = 1,
                     cex.lab = 1,
                     cex.axis = 1,
                     cex.leg = 1
                     ){

  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))

  nobjs <- length(obj.values)
  nvars <- length(obs)

  par(mfrow = c(nvars, 1), oma = c(0, 0, 0, 0))

  objs.text <- apply(data.frame(obj.names, rep(": ", nobjs), 
                     round(as.numeric(obj.values), digits.round)),
                     MARGIN = 1, FUN = paste, collapse = "")

  objs.text.just0 <- paste(format(objs.text,
                                  width = max(sapply(objs.text, FUN = nchar)),
                                  justify = "right"),
                           collapse = "\n")

  objs.text.just <- paste0(analysis.period, " GoFs\n", 
                           objs.text.just0)

  for(i in 1:nvars){

    var.sim <- sim[[i]]
    var.obs <- obs[[i]] 

    var.name <- var.names[i]
    var.unit <- var.units[i]

    if(!is.null(warmup.period)){
      flag.warmup <- TRUE
    }else{
      flag.warmup <- FALSE
      warning("Warm-up not specified")
    }



    flag.period <- sum(length(full.period)) != 
                   sum(length(cal.period) + length(warmup.period))

    if(analysis.period == "verification" | flag.period){
      flag.cal <- TRUE
    }else{
      flag.cal <- FALSE
    }

    if(is.null(xlim)){
      xlim <- c(min(full.period), max(full.period))
    }

    if(is.null(ylim)){
      ylim <- c(min(c(coredata(var.sim), coredata(var.obs)), na.rm = TRUE),
                max(c(coredata(var.sim), coredata(var.obs)), na.rm = TRUE))
    }


    plot(var.sim, type = "n", xaxt = "n", xlab = "Date",
         ylab = paste0(var.name, " [",var.unit,"]"),
         cex.main = cex.main, cex.lab = cex.lab, cex.axis = cex.axis,
         main = "", ylim = ylim, xlim = xlim)
    title(main = paste0(var.name, " in ", main), adj = 0)
    if(flag.warmup){
      polygon(x = c(warmup.period, rev(warmup.period)),
              y = c(rep(ylim[1]-0.05*(ylim[2]-ylim[1]), length(warmup.period)),
                    rep(ylim[2]+1.05*(ylim[2]-ylim[1]), length(warmup.period))),
              border = NA, col = "#DBDBDB")
    }
    if(flag.cal){
      polygon(x = c(cal.period, rev(cal.period)),
              y = c(rep(ylim[1]-0.05*(ylim[2]-ylim[1]), length(cal.period)),
                    rep(ylim[2]+1.05*(ylim[2]-ylim[1]), length(cal.period))),
              border = NA, col = "#DCEDFF")
    }
    grid()
    if(flag.warmup){
      text(x = xlim[1], y = ylim[2], adj = c(1,1), label = "Warm-Up", srt = 90,
           col = "#666666", font = 2)
    }
    if(flag.cal){
      text(x = cal.period[1], y = ylim[2], adj = c(0,1), label = "Omitted",
           srt = 0, col = "#0072EA", font = 2)
    }

    text(x = xlim[2], y = ylim[2], label = objs.text.just, adj = c(1,1),
         cex = 0.8)
    box()
    drawTimeAxis(var.sim, cex.axis = cex.axis)
    lines(var.obs, col = col.obs, type = "o", lwd = lwd, pch = pch.obs,
          cex = cex.pt)
    lines(var.sim, col = col.sim, lwd = lwd)

    par(xpd = TRUE)

    legend(x = xlim[1] + 0.80*(xlim[2] - xlim[1]),  xjust = 0,
           y = ylim[2] + 0.05*(ylim[2] - ylim[1]), yjust = 0,
           legend = c(legend.obs, legend.sim), 
           col = c(col.obs, col.sim), 
           lty = c(1,1),
           lwd = lwd,
           pch = c(pch.obs,0),
           bty = "n", 
           pt.cex = c(cex.pt,0), 
           cex = cex.leg, 
           text.col = "black", 
           horiz = FALSE, 
           inset = c(0.1, 0.1))

    par(xpd = FALSE)
    
  }

}