#'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 PlotEquiMap. This parameter is 
#'  optional and adjusted automatically by the function.
#'@param cols List of vectors of colours to be sent to PlotEquiMap 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 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 no mask, whereas cells with a 1 
#'  will result in a totally opaque superimposed pixel coloured in 'col_mask'.
#'@param col_mask Colour to be used for the superimposed mask (if specified in 
#'  'mask'). Takes the value 'grey' by default.
#'@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 legend_scale Scale factor for the size of the colour bar labels. Takes 
#'  1 by default.
#'@param cex_bar_titles Scale factor for the sizes of the bar titles. Takes 1.5 
#'  by default.
#'@param plot_margin Numeric vector of length 4 for the margin sizes in the 
#'  following order: bottom, left, top, and right. If not specified, use the 
#'  default of par("mar"), c(5.1, 4.1, 4.1, 2.1). Used as 'margin_scale' in 
#'  s2dv::PlotEquiMap.
#'@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 ... Additional parameters to be passed on to \code{PlotEquiMap}.
#' 
#'@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)
#'\dontrun{
#'PlotCombinedMap(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')), 
#'                brks = 20, width = 12, height = 10)
#'}
#'
#'Lon <- c(0:40, 350:359)
#'Lat <- 51:26
#'data <- rnorm(51 * 26 * 3)
#'dim(data) <- c(map = 3, lon = 51, lat = 26)
#'mask <-  sample(c(0,1), replace = TRUE, size = 51 * 26)
#'dim(mask) <- c(lat = 26, lon = 51)
#'\dontrun{
#'PlotCombinedMap(data, lon = Lon, lat = Lat, map_select_fun = max,
#'                display_range = range(data), mask = mask,
#'                width = 14, height = 10) 
#'}
#'
#'@seealso \code{PlotCombinedMap} and \code{PlotEquiMap}
#' 
#'@importFrom s2dv PlotEquiMap ColorBar
#'@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
PlotCombinedMap <- function(maps, lon, lat, 
                            map_select_fun, display_range, 
                            map_dim = 'map',
                            brks = NULL, cols = NULL,  
                            bar_limits = NULL, triangle_ends = c(F, F), col_inf = NULL, col_sup = NULL,
                            col_unknown_map = 'white',
                            mask = NULL, col_mask = 'grey',
                            dots = NULL,
                            bar_titles = NULL, legend_scale = 1,
                            cex_bar_titles = 1.5,
                            plot_margin = NULL, bar_extra_margin = c(2, 0, 2, 0),
                            fileout = NULL, width = 8, height = 5, 
                            size_units = 'in', res = 100, drawleg = T,
                            ...) {
  args <- list(...)

  # If there is any filenames to store the graphics, process them
  # to select the right device 
  if (!is.null(fileout)) {
    .SelectDevice <- utils::getFromNamespace(".SelectDevice", "s2dv")
    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.")
  }
  
  # 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 = legend_scale * 1.5,
                                   extra_margin = bar_extra_margin)
  
  # Check legend_scale
  if (!is.numeric(legend_scale)) {
    stop("Parameter 'legend_scale' must be numeric.")
  }
  
  # Check col_unknown_map
  if (!is.character(col_unknown_map)) {
    stop("Parameter 'col_unknown_map' must be a character string.")
  }
  
  # Check col_mask
  if (!is.character(col_mask)) {
    stop("Parameter 'col_mask' 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).")
    }
  }
 
  #----------------------
  # Identify the most likely map
  #----------------------
  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 <- display_range[2] - display_range[1] #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)
    }
    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[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) - display_range[1]) / range_width
            if (map_select_fun(x) == display_range[1]) {
              res <- res + slightly_tune_val
            }
          }
        } 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.
  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 PlotEquiMap
  #----------------------
 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(plot_margin)) {
    plot_margin <- c(5, 4, 4, 2) + 0.1  # default of par()$mar
  }

  PlotEquiMap(var = ml_map, lon = lon, lat = lat, 
              brks = tbrks, cols = tcols, drawleg = FALSE, 
              filled.continents = FALSE, dots = dots, margin_scale = plot_margin, ...)
  
  #----------------------
  # 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(col_mask, alpha.f = 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 = FALSE, add = TRUE, lwd = 1, col = coast_color) # Low resolution world map (lon -180 to 180).
    } else {
      map('world2', interior = FALSE, add = TRUE, lwd = 1, col = coast_color) # Low resolution world map (lon 0 to 360).
    }
    box()
  }
  
  #----------------------
  # Add colorbars 
  #----------------------
  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) {
    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 = legend_scale * 1.5,
                         extra_margin = bar_extra_margin)
  }
  
  # If the graphic was saved to file, close the connection with the device
  if (!is.null(fileout)) dev.off()
}