R/VizBoxWhisker.R

In esviz: Plotting Functions for Climate Science and Services

#' Box-And-Whisker Plot of Time Series with Ensemble Distribution
#' 
#' Produce time series of box-and-whisker plot showing the distribution of the 
#' members of a forecast vs. the observed evolution. The correlation between
#' forecast and observational data is calculated and displayed. Only works for 
#' n-monthly to n-yearly time series. 
#' 
#' @param exp Forecast array of multi-member time series, e.g., the NAO index 
#'   of one experiment. The expected dimensions are 
#'   c(members, start dates/forecast horizons). A vector with only the time 
#'   dimension can also be provided. Only monthly or lower frequency time 
#'   series are supported. See parameter freq.
#' @param obs Observational vector or array of time series, e.g., the NAO index 
#'   of the observations that correspond the forecast data in \code{exp}.
#'   The expected dimensions are c(start dates/forecast horizons) or 
#'   c(1, start dates/forecast horizons). Only monthly or lower frequency time 
#'   series are supported. See parameter freq.
#' @param toptitle Character string to be drawn as figure title.
#' @param ytitle Character string to be drawn as y-axis title.
#' @param monini Number of the month of the first time step, from 1 to 12.
#' @param yearini Year of the first time step.
#' @param freq Frequency of the provided time series: 1 = yearly, 12 = monthly, 
#   4 = seasonal, ... Default = 12.
#' @param expname Experimental dataset name.
#' @param obsname Name of the observational reference dataset.
#' @param drawleg TRUE/FALSE: whether to draw the legend or not.
#' @param fileout Name of output file. Extensions allowed: eps/ps, jpeg, png, 
#'   pdf, bmp and tiff. \cr
#'   Default = 'output_PlotBox.ps'.
#' @param width File width, in the units specified in the parameter size_units 
#'   (inches by default). Takes 8 by default.
#' @param height File height, in the units specified in the parameter 
#'   size_units (inches by default). Takes 5 by default.
#' @param size_units Units of the size of the device (file or window) to plot 
#'   in. Inches ('in') by default. See ?Devices and the creator function of the 
#'   corresponding device.
#' @param res Resolution of the device (file or window) to plot in. See 
#'   ?Devices and the creator function of the corresponding device.
#' @param ... Arguments to be passed to the method. Only accepts the following
#'   graphical parameters:\cr
#'   ann ask bg cex.lab cex.sub cin col.axis col.lab col.main col.sub cra crt 
#'   csi cxy err family fg fig font font.axis font.lab font.main font.sub lend 
#'   lheight ljoin lmitre mex mfcol mfrow mfg mkh oma omd omi page pin plt pty 
#'   smo srt tck tcl usr xaxp xaxs xaxt xlog xpd yaxp yaxs yaxt ylbias ylog \cr
#'   For more information about the parameters see `par`.
#' 
#' @return A figure in popup window by default, or saved to the specified path
#' via \code{fileout}.
#' 
#' @seealso EOF, ProjectField, NAO
#' 
#' @examples
#' # No example data is available over NAO region, so in this example we will
#' # tweak the longitude and latitude.
#' ano <- s2dv::Ano_CrossValid(map_temp$exp, map_temp$obs, memb = FALSE,
#'                             dat_dim = c('dat', 'member'), memb_dim = 'member')
#' nao <- s2dv::NAO(ano$exp, ano$obs, lat = seq(20, 80, length.out = 11),
#'                  lon = seq(-80, 40, length.out = 16), memb_dim = "member",
#'                  ftime_dim = "time")
#' nao$exp <- drop(aperm(nao$exp, c(2, 1, 3, 4)))
#' nao$obs <- drop(nao$obs)
#' VizBoxWhisker(nao$exp, nao$obs, toptitle = "NAO index", 
#'               ytitle = "NAO index (PC1) TOS", monini = 11, freq = 1,
#'               yearini = 2000, expname = "SEAS5", obsname = "ERA5")
#' 
#' @importFrom grDevices dev.cur dev.new dev.off 
#' @importFrom stats cor
#' @export
VizBoxWhisker <- function(exp, obs, toptitle = '', ytitle = '', monini = 1, 
                          yearini = 0, freq = 1, expname = "exp 1", 
                          obsname = "obs 1", drawleg = TRUE, fileout = NULL, 
                          width = 8, height = 5, size_units = 'in', res = 100, 
                          ...) {

  # Process the user graphical parameters that may be passed in the call
  ## Graphical parameters to exclude
  excludedArgs <- c("adj", "bty", "cex", "cex.axis", "cex.main", "col", "din", 
                    "fin", "lab", "las", "lty", "lwd", "mai", "mar", "mgp", 
                    "new", "pch", "ps")
  userArgs <- .FilterUserGraphicArgs(excludedArgs, ...)

  # If there is any filenames to store the graphics, process them
  # to select the right device 
  if (!is.null(fileout)) {
    deviceInfo <- .SelectDevice(fileout = fileout, width = width, height = height, units = size_units, res = res)
    saveToFile <- deviceInfo$fun
    fileout <- deviceInfo$files
  }

  # Checking exp
  if (is.numeric(exp)) {
    if (is.null(dim(exp)) || length(dim(exp)) == 1) {
      dim(exp) <- c(1, length(exp))
    }
  }
  if (!is.numeric(exp) || length(dim(exp)) != 2) {
    stop("Parameter 'exp' must be a numeric vector or array of dimensions c(forecast horizons/start dates) or c(ensemble members, forecast horizons/start dates)")
  }

  # Checking obs
  if (is.numeric(obs)) {
    if (is.null(dim(obs)) || length(dim(obs)) == 1) {
      dim(obs) <- c(1, length(obs))
    }
  }
  if (!is.numeric(obs) || length(dim(obs)) != 2) {
    stop("Parameter 'obs' must be a numeric vector or array of dimensions c(forecast horizons/start dates) or c(1, forecast horizons/start dates)")
  }

  # Checking consistency in exp and obs
  if (dim(exp)[2] != dim(obs)[2]) {
    stop("'exp' and 'obs' must have data for the same amount of time steps.")
  }

  if (!is.character(toptitle) || !is.character(ytitle)) {
    stop("Parameters 'ytitle' and 'toptitle' must be character strings.")
  }

  if (!is.numeric(monini)) {
    stop("'monini' must be a month number, from 1 to 12.")
  }
  if (monini < 1 || monini > 12) {
    stop("'monini' must be >= 1 and <= 12.")
  }

  if (!is.numeric(yearini)) {
    stop("'yearini' must be a month number, from 1 to 12.")
  }

  if (!is.numeric(freq)) {
    stop("'freq' must be a number <= 12.")
  }

  if (!is.character(expname) || !is.character(obsname)) {
    stop("'expname' and 'obsname' must be character strings.")
  }

  if (!is.logical(drawleg)) {
    stop("Parameter 'drawleg' must be either TRUE or FALSE.")
  }

  if (!is.character(fileout) && !is.null(fileout)) {
    stop("Parameter 'fileout' must be a character string.")
  }

  ntimesteps <- dim(exp)[2]
  lastyear <- (monini + (ntimesteps - 1) * 12 / freq - 1) %/% 12 + yearini
  lastmonth <- (monini + (ntimesteps - 1) * 12 / freq - 1) %% 12 + 1
  #
  #  Define some plot parameters
  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  #
  labind <- seq(1, ntimesteps)
  months <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep",
              "Oct", "Nov", "Dec")
  labyear <- ((labind - 1) * 12 / freq + monini - 1) %/% 12 + yearini
  labmonth <- months[((labind - 1) * 12 / freq + monini - 1) %% 12 + 1]
  for (jx in 1:length(labmonth)) {
    y2o3dig <- paste("0", as.character(labyear[jx]), sep = "")
    labmonth[jx] <- paste(labmonth[jx], "\nYr ", substr(y2o3dig,
                          nchar(y2o3dig) - 1, nchar(y2o3dig)), sep = "")
  }

  # Open connection to graphical device
  if (!is.null(fileout)) {
    saveToFile(fileout)
  } else if (names(dev.cur()) == 'null device') {
    dev.new(units = size_units, res = res, width = width, height = height)
  }

  # Load the user parameters
  oldpar <- par(c(names(userArgs), "mar", "mgp"))
  on.exit(par(oldpar), add = TRUE)
  par(userArgs)

  ## Observed time series.
  #pc.o <- ts(obs[1, ], deltat = 1, start = yr1, end = yr2)
  pc.o <- obs[1, ]
  ## Normalization of obs, forecast members. Fabian
  ## Normalization of forecast should be according to ensemble
  ## mean, to keep info on ensemble spread, no? Lauriane pc.o <-
  ## pc.o/sd(pc.o) sd.fc <- apply(exp,c(1),sd)
  ## exp <- exp/sd.fc mn.fc <-
  ## apply(exp,2, mean) exp <-
  ## exp/sd(mn.fc) Produce plot.
  par(mar = c(5, 6, 4, 2))
  boxplot(exp, add = FALSE, main = toptitle, 
    ylab = "", xlab = "", col = "red", lwd = 2, t = "b", 
    axes = FALSE, cex.main = 2, ylim = c(-max(abs(c(exp, pc.o))), max(abs(c(exp, pc.o)))))
  lines(1:ntimesteps, pc.o, lwd = 3, col = "blue")
  abline(h = 0, lty = 1)
  if (drawleg) {
    legend("bottomleft", c(obsname, expname), lty = c(1, 1), lwd = c(3, 
      3), pch = c(NA, NA), col = c("blue", "red"), horiz = FALSE, 
      bty = "n", inset = 0.05)
  }
  ##mtext(1, line = 3, text = tar, cex = 1.9)
  mtext(3, line = -2, text = paste(" AC =", round(cor(pc.o, 
        apply(exp, c(2), mean)), 2)), cex = 1.9, adj = 0)
  axis(2, cex.axis = 2)
  mtext(2, line = 3, text = ytitle, cex = 1.9)
  par(mgp = c(0, 4, 0))
  ##axis(1, c(1:ntimesteps), NA, cex.axis = 2)
  axis(1, seq(1, ntimesteps, by = 1), labmonth, cex.axis = 2)
  box()

  # If the graphic was saved to file, close the connection with the device
  if(!is.null(fileout)) dev.off()
}