R/VizAnimateMap.R

In esviz: Plotting Functions for Climate Science and Services

#' Animate Maps of Forecast/Observed Values or Scores Over Forecast Time
#' 
#' Create animations of maps in an equi-rectangular or stereographic 
#' projection, showing the anomalies, the climatologies, the mean InterQuartile 
#' Range, Maximum-Mininum, Standard Deviation, Median Absolute Deviation, 
#' the trends, the RMSE, the correlation or the RMSSS, between modelled and 
#' observed data along the forecast time (lead-time) for all input experiments 
#' and input observational datasets.
#' 
#' @param data Matrix of dimensions (nltime, nlat, nlon) or 
#'   (nexp/nmod, nltime, nlat, nlon) or (nexp/nmod, 3/4, nltime, nlat, nlon) or 
#'   (nexp/nmod, nobs, 3/4, nltime, nlat, nlon).
#' @param var Deprecated. Use 'data' instead.
#' @param lon Vector containing longtitudes (degrees).
#' @param lat Vector containing latitudes (degrees).
#' @param toptitle c('','', \dots) array of main title for each animation, 
#'   optional. If RMS, RMSSS, correlations: first exp with successive obs, then 
#'   second exp with successive obs, etc ...
#' @param title_scale Scale factor for the figure top title. Defaults to 1.
#' @param sizetit Deprecated. Use 'title_scale' instead.
#' @param units Units, optional.
#' @param monini Starting month between 1 and 12. Default = 1.
#' @param freq 1 = yearly, 12 = monthly, 4 = seasonal ...
#' @param msk95lev TRUE/FALSE grid points with dots if 95\% significance level 
#'   reached. Default = FALSE.
#' @param brks Limits of colour levels, optional. For example: 
#'   seq(min(data), max(data), (max(data) - min(data)) / 10).
#' @param cols Vector of colours of length(brks) - 1, optional.
#' @param filled.continents Continents filled in grey (TRUE) or represented by 
#'   a black line (FALSE). Default = TRUE. Filling unavailable if crossing 
#'   Greenwich and equi = TRUE. Filling unavailable if square = FALSE and 
#'   equi = TRUE.
#' @param lonmin Westward limit of the domain to plot (> 0 or < 0). 
#'   Default : 0 degrees.
#' @param lonmax Eastward limit of the domain to plot (> 0 or < 0). 
#'   lonmax > lonmin. Default : 360 degrees.
#' @param latmin Southward limit of the domain to plot. Default : -90 degrees.
#' @param latmax Northward limit of the domain to plot. Default : 90 degrees.
#' @param intlat Interval between latitude ticks on y-axis for equi = TRUE or 
#'   between latitude circles for equi = FALSE. Default = 30 degrees.
#' @param intlon Interval between longitude ticks on x-axis. 
#'   Default = 20 degrees.
#' @param drawleg Draw a colorbar. Can be FALSE only if square = FALSE or 
#'   equi = FALSE. Default = TRUE.
#' @param subsampleg Supsampling factor of the interval between ticks on 
#'   colorbar. Default = 1 = every colour level.
#' @param colNA Color used to represent NA. Default = 'white'.
#' @param equi TRUE/FALSE == cylindrical equidistant/stereographic projection. 
#'   Default: TRUE.
#' @param fileout c('', '', \dots) array of output file name for each animation.
#'    If RMS, RMSSS, correlations : first exp with successive obs, then second 
#'   exp with successive obs, etc ...
#' @param ... Arguments to be passed to the method. Only accepts the following 
#'   graphical parameters:\cr
#'   adj ann ask bty cex cex.axis cex.lab cex.main 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 las lheight ljoin lmitre lty 
#'   lwd mai mar mex mfcol mfrow mfg mgp mkh oma omd omi page pch 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}.
#' 
#' @details
#' Examples of input:
#' \enumerate{
#'   \item{
#'   Outputs from clim (exp, obs, memb = FALSE):
#'          (nmod, nltime, nlat, nlon)
#'       or (nobs, nltime, nlat, nlon) 
#'   }
#'   \item{
#'   Model output from load/ano/smoothing:
#'            (nmod, nmemb, sdate, nltime, nlat, nlon)
#'   then passed through spread(data, posdim = 2, narm = TRUE)
#'                     & mean1dim(data, posdim = 3, narm = TRUE)
#'   or through trend(mean1dim(data, 2), posTR = 2):
#'            (nmod, 3, nltime, nlat, nlon)
#'   animates average along start dates of IQR/MaxMin/SD/MAD across members 
#'   or trends of the ensemble-mean computed accross the start dates.
#'   }
#'   \item{
#'   model and observed output from load/ano/smoothing:
#'            (nmod, nmemb, sdate, nltime, nlat, nlon)
#'          & (nobs, nmemb, sdate, nltime, nlat, nlon)
#'   then averaged along members mean1dim(var_exp/var_obs, posdim = 2):
#'            (nmod, sdate, nltime, nlat, nlon)
#'            (nobs, sdate, nltime, nlat, nlon)
#'   then passed through corr(exp, obs, posloop = 1, poscor = 2)
#'                    or RMS(exp, obs, posloop = 1, posRMS = 2):
#'            (nmod, nobs, 3, nltime, nlat, nlon)
#'   animates correlations or RMS between each exp & each obs against leadtime.
#'   }
#' }
#' 
#' @examples
#' clim <- s2dv::Clim(map_temp$exp, map_temp$obs, memb = FALSE, 
#'                    dat_dim = c('dat', 'member'), memb_dim = 'member')
#' lats <- attr(map_temp$exp, "Variables")$common$lat
#' lons <- attr(map_temp$exp, "Variables")$common$lon
#' 
#' \donttest{
#' VizAnimateMap(clim$clim_exp[1, 1, , , ], lon = lons, lat = lats,
#'   toptitle = "climatology of decadal prediction", title_scale = 1, 
#'   units = "K", brks = seq(270, 300, 3), monini = 11, freq = 12, 
#'   msk95lev = FALSE, filled.continents = FALSE, intlon = 10, intlat = 10,
#'   drawleg = FALSE, fileout = "VizAnimateMap.gif")
#' unlink("VizAnimateMap.gif")
#' }
#' 
#' @importFrom grDevices postscript dev.off
#' @importFrom s2dv InsertDim
#' @export
VizAnimateMap <- function(data, lon, lat, toptitle = rep("", 11), title_scale = 1,
                          sizetit = NULL, units = "", monini = 1, freq = 12,
                          msk95lev = FALSE, brks = NULL, cols = NULL,
                          filled.continents = FALSE, lonmin = 0, lonmax = 360,
                          latmin = -90, latmax = 90, intlon = 20, intlat = 30,
                          drawleg = TRUE, subsampleg = 1, colNA = "white",
                          equi = TRUE, fileout = c("output1_animvsltime.gif",
                                                   "output2_animvsltime.gif",
                                                   "output3_animvsltime.gif"),
                          var = NULL, ...) {
  # Process the user graphical parameters that may be passed in the call
  ## Graphical parameters to exclude
  excludedArgs <- c("bg", "col", "fin", "lab", "lend", "new", "pin", "ps")
  userArgs <- .FilterUserGraphicArgs(excludedArgs, ...)

  ## fileout content with extension for consistency between
  ## functions keeping only filename without extension
  ext <- regmatches(fileout, regexpr("\\.[a-zA-Z0-9]*$", fileout))
  if ((length(ext) != 0) && any(ext != ".gif")) {
    warning("some or all extensions of the filenames provided in 'fileout' are not 'gif'. The extensions are being converted to 'gif'.")
  }
  fileout <- sub("\\.[a-zA-Z0-9]*$", "", fileout)

  # Check data
  if (missing(data) || is.null(data)) {
    if (!is.null(var)) {
      data <- var
      warning("The parameter 'var' is deprecated. Use 'data' instead.")
    } else {
      stop("Parameter 'data' cannot be NULL.")
    }
  } else if (!is.null(var)) {
    warning("The parameter 'var' is deprecated. 'data' will be used instead.")
  }
  if (!is.numeric(data) || !is.array(data)) {
    stop("Parameter 'data' must be a numeric array.")
  }
  if (length(dim(data)) < 3 || length(dim(data)) > 6) {
    stop("Parameter 'data' must be an array with 3 to 6 dimensions.")
  }
  if (length(dim(data)) == 3) {
    data <- InsertDim(data, posdim = 1, lendim = 1, name = 'new')
  }
  if (length(dim(data)) == 4) {
    data <- InsertDim(data, posdim = 2, lendim = 3, name = 'new')
  }
  if (length(dim(data)) == 5) {
    data <- InsertDim(data, posdim = 2, lendim = 1, name = 'new')
  }
  
  # Check title_scale
  if (missing(title_scale) && !missing(sizetit)) {
    warning("The parameter 'sizetit' is deprecated. Use 'title_scale' instead.")
    title_scale <- sizetit
  }
  if (!is.numeric(title_scale) || length(title_scale) != 1) {
    stop("Parameter 'title_scale' must be a single numerical value.")
  }

  nleadtime <- dim(data)[4]
  nexp <- dim(data)[1]
  nobs <- dim(data)[2]
  nlat <- dim(data)[5]
  nlon <- dim(data)[6]
  if (length(lon) != nlon | length(lat) != nlat) {
    stop("Inconsistent data dimensions / longitudes +  latitudes")
  }
  colorbar <- ClimPalette()
  if (is.null(brks) == TRUE) {
    ll <- signif(min(data[, , 2, , , ], na.rm = TRUE), 4)
    ul <- signif(max(data[, , 2, , , ], na.rm = TRUE), 4)
    if (is.null(cols) == TRUE) {
      cols <- colorbar(10)
    }
    nlev <- length(cols)
    brks <- signif(seq(ll, ul, (ul - ll)/nlev), 4)
  } else {
    if (is.null(cols) == TRUE) {
      nlev <- length(brks) - 1
      cols <- colorbar(nlev)
    } else {
      nlev <- length(cols)
    }
  }
  lon[which(lon < lonmin)] <- lon[which(lon < lonmin)] + 360
  lon[which(lon > lonmax)] <- lon[which(lon > lonmax)] - 360
  latb <- sort(lat[which(lat >= latmin & lat <= latmax)], index.return = TRUE)
  lonb <- sort(lon[which(lon >= lonmin & lon <= lonmax)], index.return = TRUE)
  
  # Define some plot parameters ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  labind <- 1:nleadtime
  months <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", 
    "Aug", "Sep", "Oct", "Nov", "Dec")
  years <- ((labind - 1) * 12/freq + monini - 1)%/%12
  suffixtit <- months[((labind - 1) * 12/freq + monini - 1)%%12 + 
    1]
  for (jx in 1:nleadtime) {
    y2o3dig <- paste("0", as.character(years[jx]), sep = "")
    suffixtit[jx] <- paste(suffixtit[jx], "-", substr(y2o3dig, 
      nchar(y2o3dig) - 1, nchar(y2o3dig)), sep = "")
  }
  
  # Loop on experimental & observational data
  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  for (jexp in 1:nexp) {
    for (jobs in 1:nobs) {
      postscript(paste(fileout[(jexp - 1) * nobs + jobs], 
        ".png", sep = ""), width = 550, height = 300, 
        bg = "white")
      # Load the user parameters
      oldpar <- par(names(userArgs))
      on.exit(par(oldpar), add = TRUE)
      par(userArgs)
      for (jt in 1:nleadtime) {
        title <- paste(toptitle[(jexp - 1) * nobs + jobs], 
          " Time=", suffixtit[jt], sep = "")
        varbis <- data[jexp, jobs, 2, jt, which(lat >= 
          latmin & lat <= latmax), which(lon >= lonmin & 
          lon <= lonmax)]
        varbis <- varbis[latb$ix, lonb$ix]
        flag <- array(FALSE, dim(varbis))
        if (msk95lev) {
          flag[which(data[jexp, jobs, 1, jt, latb$ix, 
          lonb$ix] > 0 & data[jexp, jobs, 3, jt, latb$ix, 
          lonb$ix] > 0)] <- TRUE
          flag[which(data[jexp, jobs, 1, jt, latb$ix, 
          lonb$ix] < 0 & data[jexp, jobs, 3, jt, latb$ix, 
          lonb$ix] < 0)] <- TRUE
        }
        varbis[which(varbis <= min(brks))] <- min(brks) + 
          (max(brks) - min(brks))/1000
        varbis[which(varbis >= max(brks))] <- max(brks) - 
          (max(brks) - min(brks))/1000
        if (equi) {
          VizEquiMap(t(varbis), lonb$x, latb$x, toptitle = title, 
          title_scale = title_scale, units = units, filled.continents = filled.continents, 
          dots = t(flag), brks = brks, cols = cols, 
          intxlon = intlon, intylat = intlat, drawleg = drawleg, 
          subsampleg = subsampleg, colNA = colNA, ...)
        } else {
          VizStereoMap(t(varbis), lonb$x, latb$x, latlims = c(latmin, 
          latmax), toptitle = title, title_scale = title_scale, 
          units = units, filled.continents = filled.continents, 
          dots = t(flag), brks = brks, cols = cols, 
          intlat = intlat, drawleg = drawleg, subsampleg = subsampleg, 
          colNA = colNA, ...)
        }
      }
      dev.off()
      system(paste("convert -rotate 90 -loop 10 -delay 50 ", 
        fileout[(jexp - 1) * nobs + jobs], ".png ", fileout[(jexp - 
          1) * nobs + jobs], ".gif", sep = ""))
      file.remove(paste0(fileout[(jexp - 1) * nobs + jobs], 
        ".png"))
    }
  }
}