R/VizCombinedMap.R

In esviz: Plotting Functions for Climate Science and Services

#' Plot Multiple Lon-Lat Variables In a Single Map According to a Decision Function
#'  
#' @description Plot a number a two dimensional matrices with (longitude, 
#' latitude) dimensions on a single map with the cylindrical equidistant 
#' latitude and longitude projection.
#'  
#' @author Nicolau Manubens, \email{nicolau.manubens@bsc.es}
#' @author Veronica Torralba, \email{veronica.torralba@bsc.es}
#' 
#' @param maps List of matrices to plot, each with (longitude, latitude) 
#'   dimensions, or 3-dimensional array with the dimensions (longitude, latitude, 
#'   map). Dimension names are required.
#' @param lon Vector of longitudes. Must match the length of the corresponding 
#'   dimension in 'maps'.
#' @param lat Vector of latitudes. Must match the length of the corresponding 
#'   dimension in 'maps'.
#' @param map_select_fun Function that selects, for each grid point, which value 
#'   to take among all the provided maps. This function receives as input a 
#'   vector of values for a same grid point for all the provided maps, and must 
#'   return a single selected value (not its index!) or NA. For example, the 
#'   \code{min} and \code{max} functions are accepted.
#' @param display_range Range of values to be displayed for all the maps. This 
#'   must be a numeric vector c(range min, range max). The values in the 
#'   parameter 'maps' can go beyond the limits specified in this range. If the 
#'   selected value for a given grid point (according to 'map_select_fun') falls 
#'   outside the range, it will be coloured with 'col_unknown_map'.
#' @param map_dim Optional name for the dimension of 'maps' along which the 
#'   multiple maps are arranged. Only applies when 'maps' is provided as a 
#'   3-dimensional array. Takes the value 'map' by default.
#' @param brks Colour levels to be sent to VizEquiMap. This parameter is 
#'   optional and adjusted automatically by the function.
#' @param cols List of vectors of colours to be sent to VizEquiMap for the 
#'   colour bar of each map. This parameter is optional and adjusted 
#'   automatically by the function (up to 5 maps). The colours provided for each 
#'   colour bar will be automatically interpolated to match the number of breaks. 
#'   Each item in this list can be named, and the name will be used as title for 
#'   the corresponding colour bar (equivalent to the parameter 'bar_titles').
#' @param bar_limits A numeric vector of 2 indicating the range of color bar.
#'   The default is NULL, and the function will decide the range automatically.
#' @param triangle_ends A logical vector of two indicating if the lower and upper
#'   triangles of the color bar should be plotted. The default is 
#'   c(FALSE, FALSE).
#' @param col_inf A character string of recognized color name or code indicating
#'   the color of the lower triangle of the color bar. The default is NULL.
#' @param col_sup A character string of recognized color name or code indicating
#'   the color of the upper triangle of the color bar. The default is NULL.
#' @param col_unknown_map Colour to use to paint the grid cells for which a map 
#'   is not possible to be chosen according to 'map_select_fun' or for those 
#'   values that go beyond 'display_range'. Takes the value 'white' by default.
#' @param mask Optional numeric array with dimensions (latitude, longitude), with 
#'   values in the range [0, 1], indicating the opacity of the mask over each 
#'   grid point. Cells with a 0 will result in a totally opaque superimposed 
#'   pixel coloured in 'mask_color', whereas cells with a 1 will have no mask and
#'   remain totally visible.
#' @param mask_color Colour to be used for the superimposed mask (if specified in 
#'   'mask'). Takes the value 'grey' by default.
#' @param col_mask Deprecated. Use 'mask_color' instead.
#' @param dots Array of same dimensions as 'var' or with dimensions 
#'   c(n, dim(var)), where n is the number of dot/symbol layers to add to the 
#'   plot. A value of TRUE at a grid cell will draw a dot/symbol on the 
#'   corresponding square of the plot. By default all layers provided in 'dots' 
#'   are plotted with dots, but a symbol can be specified for each of the 
#'   layers via the parameter 'dot_symbol'.
#' @param bar_titles Optional vector of character strings providing the titles to 
#'   be shown on top of each of the colour bars.
#' @param bar_label_scale Scale factor for the size of the colour bar labels. 
#'   Takes 1 by default.
#' @param legend_scale Deprecated. Use 'bar_label_scale' instead.
#' @param cex_bar_titles Scale factor for the sizes of the bar titles. Takes 1.5 
#'   by default.
#' @param margin_scale Numeric vector of length 4 for the margin sizes in the 
#'   following order: bottom, left, top, and right. If not specified (NULL), the 
#'   default of par("mar"), c(5.1, 4.1, 4.1, 2.1), is used. Default is NULL.
#' @param plot_margin Deprecated. Use 'margin_scale' instead.
#' @param bar_extra_margin A numeric vector of 4 indicating the extra margins to
#'   be added around the color bar, in the format c(y1, x1, y2, x2). The units 
#'   are margin lines. The default values are c(2, 0, 2, 0).
#' @param fileout File where to save the plot. If not specified (default) a 
#'   graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, bmp 
#'   and tiff
#' @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 drawleg Where to draw the common colour bar. Can take values TRUE, 
#'   FALSE or:\cr
#'   'up', 'u', 'U', 'top', 't', 'T', 'north', 'n', 'N'\cr
#'   'down', 'd', 'D', 'bottom', 'b', 'B', 'south', 's', 'S' (default)\cr
#'   'right', 'r', 'R', 'east', 'e', 'E'\cr
#'   'left', 'l', 'L', 'west', 'w', 'W' 
#' @param return_leg A logical value indicating if the color bars information 
#'   should be returned by the function. If TRUE, the function doesn't plot the
#'   color bars but still creates the layout with color bar areas, and the 
#'   arguments for GradientCatsColorBar() or ColorBarContinuous() will be 
#'   returned. It is convenient for users to adjust the color bars manually. The 
#'   default is FALSE, the color bars will be plotted directly.
#' @param ... Additional parameters to be passed on to \code{VizEquiMap}.
#' 
#' @return Invisibly returns \code{NULL} when \code{return_leg = FALSE} 
#' (default). When \code{return_leg = TRUE}, returns a list with color bar 
#' configuration information.
#'  
#' @examples
#' # Simple example
#' x <- array(1:(20 * 10), dim = c(lat = 10, lon = 20)) / 200
#' a <- x * 0.6
#' b <- (1 - x) * 0.6
#' c <- 1 - (a + b)
#' lons <- seq(0, 359.5, length = 20)
#' lats <- seq(-89.5, 89.5, length = 10)
#' VizCombinedMap(list(a, b, c), lons, lats, 
#'                 toptitle = 'Maximum map',
#'                 map_select_fun = max,
#'                 display_range = c(0, 1),
#'                 bar_titles = paste('% of belonging to', c('a', 'b', 'c')),
#'                 drawleg = FALSE, brks = 20, width = 8, height = 5)
#' 
#' # Mask example
#' lons2 <- c(0:40)
#' lats2 <- 51:26
#' data <- rnorm(41 * 26 * 3)
#' dim(data) <- c(map = 3, lon = 41, lat = 26)
#' mask <-  sample(c(0,1), replace = TRUE, size = 41 * 26)
#' dim(mask) <- c(lat = 26, lon = 41)
#' VizCombinedMap(data, lon = lons2, lat = lats2, map_select_fun = max,
#'                 display_range = range(data), mask = mask,
#'                 drawleg = FALSE, width = 14, height = 10)
#' 
#' @seealso \code{VizCombinedMap} and \code{VizEquiMap}
#' 
#' @import utils 
#' @importFrom maps map 
#' @importFrom graphics box image layout mtext par plot.new
#' @importFrom grDevices adjustcolor bmp colorRampPalette dev.cur dev.new dev.off 
#'   hcl jpeg pdf png postscript svg tiff
#' @export
VizCombinedMap <- function(maps, lon, lat, 
                           map_select_fun, display_range, 
                           map_dim = 'map',
                           brks = NULL, cols = NULL,  
                           bar_limits = NULL, triangle_ends = c(FALSE, FALSE), 
                           col_inf = NULL, col_sup = NULL,
                           col_unknown_map = 'white',
                           mask = NULL, mask_color = 'grey', col_mask = NULL,
                           dots = NULL,
                           bar_titles = NULL, bar_label_scale = 1,
                           legend_scale = NULL,
                           cex_bar_titles = 1.5, margin_scale = NULL,
                           plot_margin = NULL, bar_extra_margin = c(2, 0, 2, 0),
                           fileout = NULL, width = 8, height = 5, 
                           size_units = 'in', res = 100, drawleg = TRUE,
                           return_leg = FALSE, ...) {
  args <- list(...)
  country.borders <- if (!is.null(args$country.borders)) args$country.borders else FALSE
  args$country.borders <- NULL # Prevent duplication
  
  # 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
  }
  
  # Check probs
  error <- FALSE
  # Change list into an array
  if (is.list(maps)) {
    if (length(maps) < 1) {
      stop("Parameter 'maps' must be of length >= 1 if provided as a list.")
    }
    check_fun <- function(x) {
      is.numeric(x) && (length(dim(x)) == 2)
    }
    if (!all(sapply(maps, check_fun))) {
      error <- TRUE
    }
    ref_dims <- dim(maps[[1]])
    equal_dims <- all(sapply(maps, function(x) identical(dim(x), ref_dims)))
    if (!equal_dims) {
      stop("All arrays in parameter 'maps' must have the same dimension ",
           "sizes and names when 'maps' is provided as a list of arrays.")
    }
    num_maps <- length(maps)
    maps <- unlist(maps)
    dim(maps) <- c(ref_dims, map = num_maps)
    map_dim <- 'map'
  }
  if (!is.numeric(maps)) {
    error <- TRUE
  }
  if (is.null(dim(maps))) {
    error <- TRUE
  }
  if (length(dim(maps)) != 3) {
    error <- TRUE
  }
  if (error) {
    stop("Parameter 'maps' must be either a numeric array with 3 dimensions ",
         " or a list of numeric arrays of the same size with the 'lon' and ",
         "'lat' dimensions.")
  }
  dimnames <- names(dim(maps))
  
  # Check map_dim
  if (is.character(map_dim)) {
    if (is.null(dimnames)) {
      stop("Specified a dimension name in 'map_dim' but no dimension names provided ",
           "in 'maps'.")
    }
    map_dim <- which(dimnames == map_dim)
    if (length(map_dim) < 1) {
      stop("Dimension 'map_dim' not found in 'maps'.")
    } else {
      map_dim <- map_dim[1]
    }
  } else if (!is.numeric(map_dim)) {
    stop("Parameter 'map_dim' must be either a numeric value or a ",
         "dimension name.")
  }
  if (length(map_dim) != 1) {
    stop("Parameter 'map_dim' must be of length 1.")
  }
  map_dim <- round(map_dim)
  
  # Work out lon_dim and lat_dim
  lon_dim <- NULL
  if (!is.null(dimnames)) {
    lon_dim <- which(dimnames %in% c('lon', 'longitude'))[1]
  }
  if (length(lon_dim) < 1) {
    lon_dim <- (1:3)[-map_dim][1]
  }
  lon_dim <- round(lon_dim)
  
  lat_dim <- NULL
  if (!is.null(dimnames)) {
    lat_dim <- which(dimnames %in% c('lat', 'latitude'))[1]
  }
  if (length(lat_dim) < 1) {
    lat_dim <- (1:3)[-map_dim][2]
  }
  lat_dim <- round(lat_dim)
  
  # Check lon
  if (!is.numeric(lon)) {
    stop("Parameter 'lon' must be a numeric vector.")
  }
  if (length(lon) != dim(maps)[lon_dim]) {
    stop("Parameter 'lon' does not match the longitude dimension in 'maps'.")
  }
  
  # Check lat
  if (!is.numeric(lat)) {
    stop("Parameter 'lat' must be a numeric vector.")
  }
  if (length(lat) != dim(maps)[lat_dim]) {
    stop("Parameter 'lat' does not match the longitude dimension in 'maps'.")
  }
  
  # Check map_select_fun
  if (is.numeric(map_select_fun)) {
    if (length(dim(map_select_fun)) != 2) {
      stop("Parameter 'map_select_fun' must be an array with dimensions ",
           "'lon' and 'lat' if provided as an array.")
    }
    if (!identical(dim(map_select_fun), dim(maps)[-map_dim])) {
      stop("The dimensions 'lon' and 'lat' in the 'map_select_fun' array must ",
           "have the same size, name and order as in the 'maps' parameter.")
    }
  }
  if (!is.function(map_select_fun)) {
    stop("The parameter 'map_select_fun' must be a function or a numeric array.")
  }
  
  # Check bar_label_scale
  if (missing(bar_label_scale) && !missing(legend_scale)) {
    warning("The parameter 'legend_scale' is deprecated. Use 'bar_label_scale' instead.")
    bar_label_scale <- legend_scale
  }
  if (!is.numeric(bar_label_scale)) {
    stop("Parameter 'bar_label_scale' must be numeric.")
  }
  
  # Generate the desired brks and cols. Only nmap, brks, cols, bar_limits, and 
  # bar_titles matter here because plot = F.
  var_limits_maps <- range(maps, na.rm = TRUE)
  if (is.null(bar_limits)) bar_limits <- display_range
  nmap <- dim(maps)[map_dim]
  colorbar <- GradientCatsColorBar(nmap = nmap, 
                                   brks = brks, cols = cols, vertical = FALSE, subsampleg = NULL, 
                                   bar_limits = bar_limits, var_limits = var_limits_maps, 
                                   triangle_ends = triangle_ends, col_inf = col_inf, col_sup = col_sup,
                                   plot = FALSE, draw_separators = TRUE, bar_titles = bar_titles, 
                                   title_scale = cex_bar_titles, label_scale = bar_label_scale * 1.5,
                                   bar_extra_margin = bar_extra_margin)
  
  # Check col_unknown_map
  if (!is.character(col_unknown_map)) {
    stop("Parameter 'col_unknown_map' must be a character string.")
  }
  
  # Check mask_color
  if (missing(mask_color) && !missing(col_mask)) {
    warning("The parameter 'col_mask' is deprecated. Use 'mask_color' instead.")
    mask_color <- col_mask
  }
  if (!is.character(mask_color)) {
    stop("Parameter 'mask_color' must be a character string.")
  }
  
  # Check mask
  if (!is.null(mask)) {
    if (!is.numeric(mask)) {
      stop("Parameter 'mask' must be numeric.")
    }
    if (length(dim(mask)) != 2) {
      stop("Parameter 'mask' must have two dimensions.")
    }
    if ((dim(mask)[1] != dim(maps)[lat_dim]) || 
        (dim(mask)[2] != dim(maps)[lon_dim])) {
      stop("Parameter 'mask' must have dimensions c(lat, lon).")
    }
  }
  
  # Check dots
  if (!is.null(dots)) {
    if (length(dim(dots)) != 2) {
      stop("Parameter 'mask' must have two dimensions.")
    }
    if ((dim(dots)[1] != dim(maps)[lat_dim]) ||
        (dim(dots)[2] != dim(maps)[lon_dim])) {
      stop("Parameter 'mask' must have dimensions c(lat, lon).")
    }
  }
  
  # Check margin_scale
  if (missing(margin_scale) && !missing(plot_margin)) {
    warning("The parameter 'plot_margin' is deprecated. Use 'margin_scale' instead.")
    margin_scale <- plot_margin
  }
  if (!is.null(margin_scale)) {
    if (!is.numeric(margin_scale) || length(margin_scale) != 4) {
      stop("Parameter 'margin_scale' must be a numeric vector of length 4: c(bottom, left, top, right).")
    }
  }
 
  #----------------------
  # Identify the most likely map
  #----------------------
  #TODO: Consider col_inf
  if (!is.null(colorbar$col_inf[[1]])) {
    warning("Lower triangle is not supported now. Please contact maintainer if you have this need.")
  }
  if (!is.null(colorbar$col_sup[[1]])) {

    brks_norm <- vector('list', length = nmap)
    range_width <- vector('list', length = nmap)
    slightly_tune_val <- vector('list', length = nmap)
    for (ii in 1:nmap) {
      brks_norm[[ii]] <- seq(0, 1, length.out = length(colorbar$brks[[ii]]) + 1)  # add one break for col_sup
      slightly_tune_val[[ii]] <- brks_norm[[ii]][2] / (length(brks_norm[[ii]]) * 2)
      range_width[[ii]] <- diff(range(colorbar$brks[[ii]]))
    }
    ml_map <- apply(maps, c(lat_dim, lon_dim), function(x) {
      if (any(is.na(x))) {
        res <- NA
      } else {
        res <- which(x == map_select_fun(x))
        if (length(res) > 0) {
          res <- res_ind <- res[1]
          if (map_select_fun(x) < display_range[1] || map_select_fun(x) > display_range[2]) {
            res <- -0.5
          } else {
            if (map_select_fun(x) > tail(colorbar$brks[[res_ind]], 1)) { # col_sup
              res <- res + 1 - slightly_tune_val[[res_ind]]
            } else {
              res <- res + ((map_select_fun(x) - colorbar$brks[[res_ind]][1]) /
                range_width[[res_ind]] * ((length(brks_norm[[res_ind]]) - 2) / (length(brks_norm[[res_ind]]) - 1)))
              if (map_select_fun(x) == colorbar$brks[[res_ind]][1]) {
                res <- res + slightly_tune_val[[res_ind]]
              }
            }
          }
        } else {
          res <- -0.5
        }
      }
      res
    })

  } else {

    brks_norm <- vector('list', length = nmap)
    range_width <- vector('list', length = nmap)
    slightly_tune_val <- vector('list', length = nmap)
    for (ii in 1:nmap) {
      brks_norm[[ii]] <- seq(0, 1, length.out = length(colorbar$brks[[ii]]))
      slightly_tune_val[[ii]] <- brks_norm[[ii]][2] / (length(brks_norm[[ii]]) * 2)
      range_width[[ii]] <- diff(range(colorbar$brks[[ii]]))
    }
    ml_map <- apply(maps, c(lat_dim, lon_dim), function(x) {
      if (any(is.na(x))) {
        res <- NA
      } else {
        res <- which(x == map_select_fun(x))
        if (length(res) > 0) {
          res <- res_ind <- res[1]
          if (map_select_fun(x) < display_range[1] ||
            map_select_fun(x) > display_range[2]) {
            res <- -0.5
          } else {
              res <- res + ((map_select_fun(x) - colorbar$brks[[res_ind]][1]) /
                range_width[[res_ind]]) 
              if (map_select_fun(x) == colorbar$brks[[res_ind]][1]) {
                res <- res + slightly_tune_val[[res_ind]]
              }
          }
        } else {
          res <- -0.5
        }
      }
      res
    })
  }

  nlat <- length(lat)
  nlon <- length(lon)
  
  #----------------------
  # Set latitudes from minimum to maximum
  #----------------------
  if (lat[1] > lat[nlat]) {
    lat <- lat[nlat:1]
    indices <- list(nlat:1, TRUE)
    ml_map <- do.call("[", c(list(x = ml_map), indices))
    if (!is.null(mask)){
      mask <- mask[nlat:1, ]
    }
    if (!is.null(dots)){
      dots <- dots[nlat:1,]
    }
  }
  
  #----------------------
  # Set layout and parameters
  #----------------------
  # 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)
  }
  #NOTE: I think plot.new() is not necessary in any case.
#  plot.new()
  #TODO: Don't hardcoded. Let users decide.
  oldpar <- par(c("font.main", "mar"))
  on.exit(par(oldpar), add = TRUE)
  par(font.main = 1)
  # If colorbars need to be plotted, re-define layout.
  if (drawleg) {
    layout(matrix(c(rep(1, nmap),2:(nmap + 1)), 2, nmap, byrow = TRUE), heights = c(6, 1.5))
  }

  #----------------------
  # Set colors and breaks and then VizEquiMap
  #----------------------
 if (!is.null(colorbar$col_sup[[1]])) {
    tcols <- c(col_unknown_map, colorbar$cols[[1]], colorbar$col_sup[[1]])
    tbrks <- c(-1, brks_norm[[1]] + rep(1, each = length(brks_norm[[1]])))
    for (k in 2:nmap) {
      tcols <- append(tcols, c(col_unknown_map, colorbar$cols[[k]], colorbar$col_sup[[k]]))
      tbrks <- c(tbrks, brks_norm[[k]] + rep(k, each = length(brks_norm[[k]])))
    }
  } else {  # original code
    tcols <- c(col_unknown_map, colorbar$cols[[1]])
    tbrks <- c(-1, brks_norm[[1]] + rep(1, each = length(brks_norm[[1]])))
    for (k in 2:nmap) {
      tcols <- append(tcols, c(col_unknown_map, colorbar$cols[[k]]))
      tbrks <- c(tbrks, brks_norm[[k]] + rep(k, each = length(brks_norm[[k]])))
    }
  }

  if (is.null(margin_scale)) {
    margin_scale <- c(5, 4, 4, 2) + 0.1  # default of par()$mar
  }
  
  country.borders.VizEquiMap <- if (!is.null(mask)) FALSE else country.borders
  
  do.call(VizEquiMap, c(list(
    data = ml_map, lon = lon, lat = lat,
    brks = tbrks, cols = tcols, drawleg = FALSE,
    filled.continents = FALSE, country.borders = country.borders.VizEquiMap,
    dots = dots, margin_scale = margin_scale, bar_extra_margin = bar_extra_margin 
  ), args))
  
  #----------------------
  # Add overplot on top
  #----------------------
  if (!is.null(mask)) {
    dims_mask <- dim(mask)
    latb <- sort(lat, index.return = TRUE)
    dlon <- lon[2:dims_mask[2]] - lon[1:(dims_mask[2] - 1)]
    wher <- which(dlon > (mean(dlon) + 1))
    if (length(wher) > 0) {
      lon[(wher + 1):dims_mask[2]] <- lon[(wher + 1):dims_mask[2]] - 360
    }
    lonb <- sort(lon, index.return = TRUE)

    cols_mask <- sapply(seq(from = 0, to = 1, length.out = 10), 
                        function(x) adjustcolor(mask_color, alpha.f = 1 - x))
    image(lonb$x, latb$x, t(mask)[lonb$ix, latb$ix], 
          axes = FALSE, col = cols_mask, 
          breaks = seq(from = 0, to = 1, by = 0.1), 
          xlab = '', ylab = '', add = TRUE, xpd = TRUE)
    if (!exists('coast_color')) {
      coast_color <- 'black'
    }
    if (min(lon) < 0) {
      map('world', interior = country.borders, add = TRUE, lwd = 1, col = coast_color) # Low resolution world map (lon -180 to 180).
    } else {
      map('world2', interior = country.borders, add = TRUE, lwd = 1, col = coast_color) # Low resolution world map (lon 0 to 360).
    }
    if (!is.null(args$shapefile)) {
      if (is.list(args$shapefile)) {
        shape <- args$shapefile
      } else {
        shape <- readRDS(file = args$shapefile)
      }
      shapefile_color <- if (!is.null(args$shapefile_color)) args$shapefile_color else "black"
      shapefile_lwd <- if (!is.null(args$shapefile_lwd)) args$shapefile_lwd else 1
      filled.continents <- if (!is.null(args$filled.continents)) args$filled.continents else FALSE
      maps::map(shape, interior = country.borders, #wrap = wrap_vec,
                fill = filled.continents, add = TRUE, plot = TRUE,
                lwd = shapefile_lwd, col = shapefile_color)
    }
    box()
  }
  
  #----------------------
  # Add colorbars 
  #----------------------
  
  ## TODO: review placement of top margin manipulation
  # if ('toptitle' %in% names(args)) {
  #   size_title <- 1
  #   if ('title_scale' %in% names(args)) {
  #     size_title <- args[['title_scale']]
  #   }
  #   old_mar <- par('mar')
  #   old_mar[3] <- old_mar[3] - (2 * size_title + 1)
  #   par(mar = old_mar)
  # }

  if (drawleg & !return_leg) {
    GradientCatsColorBar(nmap = dim(maps)[map_dim],
                         brks = colorbar$brks, cols = colorbar$cols, vertical = FALSE,
                         subsampleg = NULL, bar_limits = bar_limits, var_limits = var_limits_maps,
                         triangle_ends = triangle_ends, col_inf = colorbar$col_inf, 
                         col_sup = colorbar$col_sup, plot = TRUE, draw_separators = TRUE,
                         bar_titles = bar_titles, title_scale = cex_bar_titles, 
                         label_scale = bar_label_scale * 1.5, bar_extra_margin = bar_extra_margin)
  }

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

  if (return_leg) {
    tmp <- list(nmap = dim(maps)[map_dim],
                brks = colorbar$brks, cols = colorbar$cols, vertical = FALSE,
                subsampleg = NULL, bar_limits = bar_limits,
                var_limits = var_limits_maps,
                triangle_ends = triangle_ends, col_inf = colorbar$col_inf, col_sup = colorbar$col_sup,
                plot = TRUE, draw_separators = TRUE,
                bar_titles = bar_titles, title_scale = cex_bar_titles, label_scale = bar_label_scale * 1.5,
                extra_margin = bar_extra_margin)
    warning("The device is not off yet. Use dev.off() after plotting the color bars.")
    return(tmp)
    #NOTE: The device is not off! Can keep plotting the color bars.
  }

}