From 7fea176b8024a565a90724e63278e90b686e2f84 Mon Sep 17 00:00:00 2001
From: ARIADNA BATALLA FERRES <ariadna.batalla@bsc.es>
Date: Mon, 19 May 2025 10:25:17 +0200
Subject: [PATCH 1/4] Update country borders dev in PlotCombinedMap.R

---
 modules/Visualization/R/tmp/PlotCombinedMap.R | 20 ++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/modules/Visualization/R/tmp/PlotCombinedMap.R b/modules/Visualization/R/tmp/PlotCombinedMap.R
index e4d6d7f3..93e04c5f 100644
--- a/modules/Visualization/R/tmp/PlotCombinedMap.R
+++ b/modules/Visualization/R/tmp/PlotCombinedMap.R
@@ -138,6 +138,8 @@ PlotCombinedMap <- function(maps, lon, lat,
                             size_units = 'in', res = 100, drawleg = T, 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 
@@ -443,11 +445,16 @@ PlotCombinedMap <- function(maps, lon, lat,
   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, ...)
   
+  country.borders.PlotEquiMap <- if (!is.null(mask)) FALSE else country.borders
+  
+  do.call(PlotEquiMap, c(list(
+    var = ml_map, lon = lon, lat = lat, brks = tbrks, 
+    cols = tcols, drawleg = FALSE, filled.continents = FALSE,
+    country.borders = country.borders.PlotEquiMap,
+    dots = dots, margin_scale = plot_margin
+  ), args))
+                        
   #----------------------
   # Add overplot on top
   #----------------------
@@ -471,9 +478,9 @@ PlotCombinedMap <- function(maps, lon, lat,
       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).
+      map('world', interior = country.borders, 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).
+      map('world2', interior = country.borders, add = TRUE, lwd = 1, col = coast_color) # Low resolution world map (lon 0 to 360).
     }
     box()
   }
@@ -522,4 +529,3 @@ PlotCombinedMap <- function(maps, lon, lat,
 
 }
 
-
-- 
GitLab


From 93d029dd8ca1d28a3d563e1bb9b9d7ffc78a10d3 Mon Sep 17 00:00:00 2001
From: ARIADNA BATALLA FERRES <ariadna.batalla@bsc.es>
Date: Mon, 19 May 2025 11:08:18 +0200
Subject: [PATCH 2/4] Add VizCombinedMap.R and VizMostLikelyQuantileMap.R and
 update names in plot_most_likely_terciles_map.R

---
 .../R/plot_most_likely_terciles_map.R         |  10 +-
 modules/Visualization/R/tmp/VizCombinedMap.R  | 550 ++++++++++++++++++
 .../R/tmp/VizMostLikelyQuantileMap.R          | 186 ++++++
 3 files changed, 741 insertions(+), 5 deletions(-)
 create mode 100644 modules/Visualization/R/tmp/VizCombinedMap.R
 create mode 100644 modules/Visualization/R/tmp/VizMostLikelyQuantileMap.R

diff --git a/modules/Visualization/R/plot_most_likely_terciles_map.R b/modules/Visualization/R/plot_most_likely_terciles_map.R
index c2d29859..16ca0af0 100644
--- a/modules/Visualization/R/plot_most_likely_terciles_map.R
+++ b/modules/Visualization/R/plot_most_likely_terciles_map.R
@@ -1,7 +1,7 @@
 ## Functions required for normal cat and triangles end.
 ## See https://earth.bsc.es/gitlab/external/cstools/-/issues/125
-source("modules/Visualization/R/tmp/PlotMostLikelyQuantileMap.R")
-source("modules/Visualization/R/tmp/PlotCombinedMap.R")
+source("modules/Visualization/R/tmp/VizMostLikelyQuantileMap.R")
+source("modules/Visualization/R/tmp/VizCombinedMap.R")
 source("modules/Visualization/R/tmp/ColorBar.R")
 source("modules/Visualization/R/tmp/clim.palette.R")
 source("modules/Visualization/R/tmp/Utils.R")
@@ -171,10 +171,10 @@ plot_most_likely_terciles <- function(recipe,
                              i_syear)
         }
         # Plots
-        ## NOTE: PlotLayout() and PlotMostLikelyQuantileMap() are still being worked
+        ## NOTE: PlotLayout() and VizMostLikelyQuantileMap() are still being worked
         ## on. This option does not work with mask or dots for now. 
         suppressWarnings(
-        PlotLayout(PlotMostLikelyQuantileMap, c('bin', 'longitude', 'latitude'),
+        PlotLayout(VizMostLikelyQuantileMap, c('bin', 'longitude', 'latitude'),
                    cat_dim = 'bin',
                    i_var_probs, longitude, latitude,
                    mask = var_mask,
@@ -315,7 +315,7 @@ plot_most_likely_terciles <- function(recipe,
           base_args$probs <- i_var_probs[i, , , ]
           base_args$mask <- mask_i
           base_args$dots <- dots_i
-          cb_info <- do.call(PlotMostLikelyQuantileMap,
+          cb_info <- do.call(VizMostLikelyQuantileMap,
                              args = c(base_args,
                                       list(toptitle = toptitle,
                                            fileout = fileout)))
diff --git a/modules/Visualization/R/tmp/VizCombinedMap.R b/modules/Visualization/R/tmp/VizCombinedMap.R
new file mode 100644
index 00000000..bd97ff9b
--- /dev/null
+++ b/modules/Visualization/R/tmp/VizCombinedMap.R
@@ -0,0 +1,550 @@
+#'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 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 
+#'  VizEquiMap.
+#'@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}.
+#' 
+#'@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{
+#'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')), 
+#'                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{
+#'VizCombinedMap(data, lon = Lon, lat = Lat, map_select_fun = max,
+#'                display_range = range(data), mask = mask,
+#'                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, 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, 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.")
+  }
+  
+  # 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
+  #----------------------
+  #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.
+  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(plot_margin)) {
+    plot_margin <- 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(
+    var = ml_map, lon = lon, lat = lat,
+    brks = tbrks, cols = tcols, drawleg = FALSE,
+    filled.continents = FALSE, country.borders = country.borders.VizEquiMap,
+    dots = dots, margin_scale = plot_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(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 = 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).
+    }
+    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 & !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 = legend_scale * 1.5,
+                         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 = legend_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.
+  }
+  
+}
+
diff --git a/modules/Visualization/R/tmp/VizMostLikelyQuantileMap.R b/modules/Visualization/R/tmp/VizMostLikelyQuantileMap.R
new file mode 100644
index 00000000..079a065c
--- /dev/null
+++ b/modules/Visualization/R/tmp/VizMostLikelyQuantileMap.R
@@ -0,0 +1,186 @@
+#'Plot Maps of Most Likely Quantiles
+#'
+#'@author Veronica Torralba, \email{veronica.torralba@bsc.es}, Nicolau Manubens, 
+#'\email{nicolau.manubens@bsc.es}
+#'@description This function receives as main input (via the parameter 
+#'\code{probs}) a collection of longitude-latitude maps, each containing the 
+#'probabilities (from 0 to 1) of the different grid cells of belonging to a 
+#'category. As many categories as maps provided as inputs are understood to 
+#'exist. The maps of probabilities must be provided on a common rectangular 
+#'regular grid, and a vector with the longitudes and a vector with the latitudes 
+#'of the grid must be provided. The input maps can be provided in two forms, 
+#'either as a list of multiple two-dimensional arrays (one for each category) or 
+#'as a three-dimensional array, where one of the dimensions corresponds to the 
+#'different categories.
+#'
+#'@param probs A list of bi-dimensional arrays with the named dimensions 
+#'  'latitude' (or 'lat') and 'longitude' (or 'lon'), with equal size and in the 
+#'  same order, or a single tri-dimensional array with an additional dimension 
+#'  (e.g. 'bin') for the different categories. The arrays must contain 
+#'  probability values between 0 and 1, and the probabilities for all categories 
+#'  of a grid cell should not exceed 1 when added.
+#'@param lon A numeric vector with the longitudes of the map grid, in the same 
+#'  order as the values along the corresponding dimension in \code{probs}.
+#'@param lat A numeric vector with the latitudes of the map grid, in the same 
+#'  order as the values along the corresponding dimension in \code{probs}.
+#'@param cat_dim The name of the dimension along which the different categories 
+#'  are stored in \code{probs}. This only applies if \code{probs} is provided in 
+#'  the form of 3-dimensional array. The default expected name is 'bin'.
+#'@param bar_titles Vector of character strings with the names to be drawn on 
+#'  top of the color bar for each of the categories. As many titles as 
+#'  categories provided in \code{probs} must be provided.
+#'@param col_unknown_cat Character string with a colour representation of the 
+#'  colour to be used to paint the cells for which no category can be clearly 
+#'  assigned. Takes the value 'white' by default.
+#'@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 sent to \code{VizCombinedMap} and 
+#'  \code{VizEquiMap}.
+#'@seealso \code{VizCombinedMap} and \code{VizEquiMap}
+#'@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{
+#'VizMostLikelyQuantileMap(list(a, b, c), lons, lats, 
+#'                         toptitle = 'Most likely tercile map',
+#'                         bar_titles = paste('% of belonging to', c('a', 'b', 'c')), 
+#'                         brks = 20, width = 10, height = 8)
+#'}
+#'
+#'# More complex example
+#'
+#'# 1. Generation of sample data
+#'lons <- seq(0, 359.5, length = 40)
+#'lats <- seq(-89.5, 89.5, length = 20)
+#'
+#'# Generate sample data (with dimensions time, lat, lon)
+#'sample_data <- sample(1:5, 50 * 20 * 40, replace = TRUE)
+#'dim(sample_data) <- c(time = 50, lat = 20, lon = 40)
+#'
+#'# 2. Binning sample data
+#'n_bins <- 4
+#'prob_thresholds <- 1:n_bins / n_bins
+#'prob_thresholds <- prob_thresholds[1:(n_bins - 1)]
+#'thresholds <- quantile(sample_data, prob_thresholds)
+#'
+#'binning <- function(x, thresholds) {
+#'  n_samples <- length(x)
+#'  n_bins <- length(thresholds) + 1
+#'  thresholds <- c(thresholds, max(x))
+#'  result <- 1:n_bins
+#'  lower_threshold <- min(x) - 1
+#'  for (i in 1:n_bins) {
+#'    result[i] <- sum(x > lower_threshold & x <= thresholds[i]) / n_samples
+#'    lower_threshold <- thresholds[i]
+#'  }
+#'  dim(result) <- c(bin = n_bins)
+#'  result
+#'}
+#'bins <- apply(sample_data, 2:3, binning, thresholds)
+#'names(dim(bins))[1] <- "bin"
+#'
+#'\dontrun{
+#'VizMostLikelyQuantileMap(bins, lons, lats, 
+#'                         toptitle = 'Most likely quantile map',
+#'                         bar_titles = paste('% of belonging to', letters[1:n_bins]),
+#'                         brks = 20, width = 10, height = 10)
+#'}
+#'@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
+VizMostLikelyQuantileMap <- function(probs, lon, lat, cat_dim = 'bin',
+                                     bar_titles = NULL, 
+                                     col_unknown_cat = 'white', drawleg = T,
+                                     ...) {
+  # Check probs
+  error <- FALSE
+  if (is.list(probs)) {
+    if (length(probs) < 1) {
+      stop("Parameter 'probs' must be of length >= 1 if provided as a list.")
+    }
+    check_fun <- function(x) {
+      is.numeric(x) && (length(dim(x)) == 2)
+    }
+    if (!all(sapply(probs, check_fun))) {
+      error <- TRUE
+    }
+    ref_dims <- dim(probs[[1]])
+    equal_dims <- all(sapply(probs, function(x) identical(dim(x), ref_dims)))
+    if (!equal_dims) {
+      stop("All arrays in parameter 'probs' must have the same dimension ",
+           "sizes and names when 'probs' is provided as a list of arrays.")
+    }
+    num_probs <- length(probs)
+    probs <- unlist(probs)
+    dim(probs) <- c(ref_dims, map = num_probs)
+    cat_dim <- 'map'
+  }
+  if (!is.numeric(probs)) {
+    error <- TRUE
+  }
+  if (is.null(dim(probs))) {
+    error <- TRUE
+  }
+  if (length(dim(probs)) != 3) {
+    error <- TRUE
+  }
+  if (error) {
+    stop("Parameter 'probs' 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(probs))
+  
+  # Check cat_dim
+  if (is.character(cat_dim)) {
+    if (is.null(dimnames)) {
+      stop("Specified a dimension name in 'cat_dim' but no dimension names provided ",
+           "in 'probs'.")
+    }
+    cat_dim <- which(dimnames == cat_dim)
+    if (length(cat_dim) < 1) {
+      stop("Dimension 'cat_dim' not found in 'probs'.")
+    }
+    cat_dim <- cat_dim[1]
+  } else if (!is.numeric(cat_dim)) {
+    stop("Parameter 'cat_dim' must be either a numeric value or a ",
+         "dimension name.")
+  }
+  if (length(cat_dim) != 1) {
+    stop("Parameter 'cat_dim' must be of length 1.")
+  }
+  cat_dim <- round(cat_dim)
+  nprobs <- dim(probs)[cat_dim]
+  
+  # Check bar_titles
+  if (is.null(bar_titles)) {
+    if (nprobs == 3) {
+      bar_titles <- c("Below normal (%)", "Normal (%)", "Above normal (%)")
+    } else if (nprobs == 5) {
+      bar_titles <- c("Low (%)", "Below normal (%)", 
+                      "Normal (%)", "Above normal (%)", "High (%)")
+    } else {
+      bar_titles <- paste0("Cat. ", 1:nprobs, " (%)")
+    }
+  }
+  
+  minimum_value <- ceiling(1 / nprobs * 10 * 1.1) * 10
+  
+  # By now, the VizCombinedMap function is included below in this file. 
+  # In the future, VizCombinedMap will be part of s2dverification and will 
+  # be properly imported.
+  VizCombinedMap(probs * 100, lon, lat, map_select_fun = max, 
+                 display_range = c(minimum_value, 100), map_dim = cat_dim, 
+                 bar_titles = bar_titles, col_unknown_map = col_unknown_cat, 
+                 drawleg = drawleg, ...)
+}
\ No newline at end of file
-- 
GitLab


From e79aaa3f6b7082a5e50107be3ad271a0fc213740 Mon Sep 17 00:00:00 2001
From: ARIADNA BATALLA FERRES <ariadna.batalla@bsc.es>
Date: Wed, 28 May 2025 14:56:20 +0200
Subject: [PATCH 3/4] Add shapefile after appluing mask in VizCombinedMap.R

---
 modules/Visualization/R/tmp/VizCombinedMap.R | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/modules/Visualization/R/tmp/VizCombinedMap.R b/modules/Visualization/R/tmp/VizCombinedMap.R
index bd97ff9b..a675fe18 100644
--- a/modules/Visualization/R/tmp/VizCombinedMap.R
+++ b/modules/Visualization/R/tmp/VizCombinedMap.R
@@ -500,6 +500,14 @@ VizCombinedMap <- function(maps, lon, lat,
     } 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)) {
+      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(shapefile, interior = country.borders, #wrap = wrap_vec,
+                fill = filled.continents, add = TRUE, plot = TRUE,
+                lwd = shapefile_lwd, col = shapefile_color)
+    }
     box()
   }
   
-- 
GitLab


From b53ae912416ab2e2a9d47e937039204ba545cf26 Mon Sep 17 00:00:00 2001
From: ARIADNA BATALLA FERRES <ariadna.batalla@bsc.es>
Date: Wed, 28 May 2025 18:02:52 +0200
Subject: [PATCH 4/4] Fix processing of shapefile

---
 modules/Visualization/R/tmp/VizCombinedMap.R | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/modules/Visualization/R/tmp/VizCombinedMap.R b/modules/Visualization/R/tmp/VizCombinedMap.R
index a675fe18..b9707ad6 100644
--- a/modules/Visualization/R/tmp/VizCombinedMap.R
+++ b/modules/Visualization/R/tmp/VizCombinedMap.R
@@ -501,10 +501,15 @@ VizCombinedMap <- function(maps, lon, lat,
       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(shapefile, interior = country.borders, #wrap = wrap_vec,
+      maps::map(shape, interior = country.borders, #wrap = wrap_vec,
                 fill = filled.continents, add = TRUE, plot = TRUE,
                 lwd = shapefile_lwd, col = shapefile_color)
     }
-- 
GitLab