From e827d257093bac9228c8107bcb8e4b7dd98058a0 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Wed, 8 Jan 2025 10:59:23 +0100
Subject: [PATCH 01/24] (WIP) Extend plot_ensemble_mean() to other forecast
 metrics

---
 modules/Visualization/R/plot_ensemble_mean.R | 72 ++++++++++-----
 modules/Visualization/Visualization.R        | 93 +++++++++++---------
 2 files changed, 102 insertions(+), 63 deletions(-)

diff --git a/modules/Visualization/R/plot_ensemble_mean.R b/modules/Visualization/R/plot_ensemble_mean.R
index 31a20e55..0657b25d 100644
--- a/modules/Visualization/R/plot_ensemble_mean.R
+++ b/modules/Visualization/R/plot_ensemble_mean.R
@@ -1,4 +1,16 @@
-plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, output_conf) {
+# recipe: the SUNSET recipe
+# fcst: an s2dv_cube object containing the forecast to be plotted. The forecast
+#   should have the following dimensions: 
+#   'dat', 'var', 'syear', 'time', 'latitude', 'longitude', 'ensemble'
+# mask: 
+# dots: 
+# outdir: output directory
+# output_conf: list as returned by read_yaml("output_conf.yml")
+# method: Computation to be applied to the forecast. Options:
+#   'median' (default), 'iqr', 'mean'.
+plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL, 
+                              outdir, output_conf,
+                              method = 'median') {
   ## TODO: Add 'anomaly' to plot title
   # Abort if frequency is daily
   if (recipe$Analysis$Variables$freq %in% c("daily", "daily_mean")) {
@@ -34,23 +46,37 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
   } else {
     projection <- "cylindrical_equidistant"
   }
-  # Compute ensemble mean
-  ensemble_mean <- s2dv::MeanDims(fcst$data, 'ensemble')
+  # Compute ensemble mean or other
+  if (!is.null(method)) {
+    method <- tolower(method)
+  }
+  if (method == 'mean') {
+    forecast_product <- s2dv::MeanDims(fcst$data, 'ensemble')
+  } else if (method == 'iqr') {
+    forecast_produc <- Apply(fcst$data, target_dim = 'ensemble',
+                             fun = function(x) {
+                             IQR(x)})$output1
+  } else {
+    ## return median by default?
+    forecast_product <- Apply(fcst$data, target_dim = 'ensemble',
+                              fun = function(x) {
+                              median(x)})$output1
+  }
   # Loop over variable dimension
   for (var in 1:fcst$dims[['var']]) {
     variable <- fcst$attrs$Variable$varName[[var]]
     var_long_name <- fcst$attrs$Variable$metadata[[variable]]$long_name
     units <- fcst$attrs$Variable$metadata[[variable]]$units
     # Subset ensemble mean by variable
-    var_ens_mean <- ClimProjDiags::Subset(ensemble_mean,
-                                          along = c("dat", "var",
-                                                    "sday", "sweek"),
-                                          indices = list(1, var, 1, 1),
-                                          drop = 'selected')
-    var_ens_mean <- Reorder(var_ens_mean, c("syear",
-                                            "time",
-                                            "longitude",
-                                            "latitude"))
+    var_fcst <- ClimProjDiags::Subset(forecast_product,
+                                      along = c("dat", "var",
+                                                "sday", "sweek"),
+                                      indices = list(1, var, 1, 1),
+                                      drop = 'selected')
+    var_fcst <- Reorder(var_fcst, c("syear",
+                                    "time",
+                                    "longitude",
+                                    "latitude"))
     ## TODO: Redefine column colors, possibly depending on variable
     if (variable == 'prlr') {
       palette = "BrBG"
@@ -62,7 +88,7 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
     # Define brks, centered around zero in the case of anomalies
     if (grepl("anomaly", var_long_name)) {
       variable <- paste(variable, "anomaly")
-      max_value <- max(abs(var_ens_mean), na.rm = TRUE)
+      max_value <- max(abs(var_fcst), na.rm = TRUE)
       ugly_intervals <- seq(-max_value, max_value, max_value/20)
       brks <- pretty(ugly_intervals, n = 12, min.n = 8)
       cols <- grDevices::hcl.colors(length(brks) - 1, palette, rev = rev)
@@ -74,10 +100,13 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
         cols <- c("#33BFD1", "white", "#FF764D")
         col_fun <- colorRampPalette(cols)
       }
-      if (all(is.na(var_ens_mean))) {
+      if (all(is.na(var_fcst))) {
         brks <- NULL
       } else {
-        brks <- pretty(range(var_ens_mean, na.rm = T), n = 15, min.n = 8)
+        max_value <- max(abs(var_fcst), na.rm = T)
+        brks <- pretty(c(-1 * max_value, max_value), n = 15, min.n = 8)
+        ## old:
+        ## brks <- pretty(range(var_fcst, na.rm = T), n = 15, min.n = 8)
       }
       if (is.null(brks)) {
         cols <- NULL
@@ -87,11 +116,12 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
     }
     
     for (i_syear in start_date) {
-      i_var_ens_mean <- ClimProjDiags::Subset(var_ens_mean,
+      i_var_fcst <- ClimProjDiags::Subset(var_fcst,
                                               along = c("syear"),
                                               indices = which(start_date == i_syear),
                                               drop = 'selected')
-      outfile <- paste0(outdir[[var]], "forecast_ensemble_mean-", i_syear)
+      outfile <- paste0(outdir[[var]], "forecast_ensemble_", method,
+                        "-", i_syear)
       # Mask
       if (!is.null(mask)) {
         outfile <- paste0(outfile, "_enscormask")
@@ -144,8 +174,8 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
         output_configuration <- output_conf$multipanel$forecast_ensemble_mean
         base_args <- list(fun = "PlotEquiMap",
                           plot_dims = c('longitude', 'latitude'),
-                          var = i_var_ens_mean, lon = longitude,
-                          lat = latitude, mask = mask, dots = dots,
+                          var = i_var_fcst, lon = longitude,
+                          lat = latitude, mask = var_mask, dots = var_dots,
                           filled.continents = FALSE, toptitle = toptitle,
                           title_scale = 0.7, subtitle_scale = 0.8,
                           subtitle_margin_scale = 2, titles = titles,
@@ -235,7 +265,7 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
                      "Units: ", units)
           }
           # Modify base arguments
-          base_args[[1]] <- i_var_ens_mean[i, , ]
+          base_args[[1]] <- i_var_fcst[i, , ]
           if (is.null(attributes(fcst$attrs$time_bounds))) {
             fileout <- paste0(outfile, "_ft", sprintf("%02d", i), ".pdf")
           } else {
@@ -259,5 +289,5 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
     }
   }
   info(recipe$Run$logger,
-       "##### FORECAST ENSEMBLE MEAN PLOTS SAVED TO OUTPUT DIRECTORY #####")
+       "##### FORECAST MAPS SAVED TO OUTPUT DIRECTORY #####")
 }
diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R
index 749fe3b6..091b55f1 100644
--- a/modules/Visualization/Visualization.R
+++ b/modules/Visualization/Visualization.R
@@ -139,51 +139,60 @@ Visualization <- function(recipe,
       if (is.null(recipe$Analysis$Workflow$Visualization$dots)) {
         recipe$Analysis$Workflow$Visualization$dots <- FALSE
       }
-      # Plot without mask or dots
-      if ((recipe$Analysis$Workflow$Visualization$mask_ens
-           %in% c('both', FALSE)) ||
-          (recipe$Analysis$Workflow$Visualization$dots
-           %in% c('both', FALSE))) {
-        plot_ensemble_mean(recipe, data$fcst, outdir, 
-                           mask = NULL, dots = NULL,
-                           output_conf = output_conf)
+      if (is.null(recipe$Analysis$Workflow$Visualization$dots)) {
+        recipe$Analysis$Workflow$Visualization$forecast_method <- "median"
       }
-      # Plots with masked
-      if (recipe$Analysis$Workflow$Visualization$mask_ens %in%
-          c('both', TRUE)) {
-        if (is.null(skill_metrics)) {
-          error(recipe$Run$logger,
-                paste0("For the forecast ensemble mean  plot, skill_metrics ",
-                       "need to be provided to be masked."))
-        } else if (!('enscorr' %in% names(skill_metrics))) {
-          error(recipe$Run$logger,
-                paste0("For the forecast ensemble mean plot, enscor metric ",
-                       "need to be provided to be masked"))
-        } else {
-          plot_ensemble_mean(recipe, data$fcst,
-                             mask = skill_metrics$enscorr,
-                             dots = NULL,
-                             outdir, output_conf = output_conf)
+      
+      # Loop over the methods requested
+      for (method in strsplit(recipe$Analysis$Workflow$Visualization$forecast_method, ", | |,")[[1]]) {
+        # Plot without mask or dots
+        if ((recipe$Analysis$Workflow$Visualization$mask_ens
+             %in% c('both', FALSE)) ||
+            (recipe$Analysis$Workflow$Visualization$dots
+             %in% c('both', FALSE))) {
+          plot_forecast_map(recipe, data$fcst,
+                            mask = NULL, dots = NULL,
+                            outdir = outdir, method = method,
+                            output_conf = output_conf)
         }
-      }
-      # Plots with dotted negative correlated in ens-mean-fcst
-      if (recipe$Analysis$Workflow$Visualization$dots %in% c('both', TRUE)) {
-        if (is.null(skill_metrics)) {
-          error(recipe$Run$logger,
-                paste0("For the forecast ensemble mean plot, skill_metrics ",
-                       "need to be provided for the dots"))
-        } else if (!('enscorr' %in% names(skill_metrics))) {
-          error(recipe$Run$logger,
-                paste0("For the forecast ensemble mean plot, enscor metric ",
-                       "needs to be provided for the dots"))
-        } else {
-          plot_ensemble_mean(recipe, data$fcst,
-                             mask = NULL,
-                             dots = skill_metrics$enscorr,
-                             outdir, output_conf = output_conf)
+        # Plots with masked
+        if (recipe$Analysis$Workflow$Visualization$mask_ens %in%
+            c('both', TRUE)) {
+          if (is.null(skill_metrics)) {
+            error(recipe$Run$logger,
+                  paste0("For the forecast ensemble mean  plot, skill_metrics ",
+                         "need to be provided to be masked."))
+          } else if (!('enscorr' %in% names(skill_metrics))) {
+            error(recipe$Run$logger,
+                  paste0("For the forecast ensemble mean plot, enscor metric ",
+                         "need to be provided to be masked"))
+          } else {
+            plot_forecast_map(recipe, data$fcst,
+                              mask = skill_metrics$enscorr,
+                              dots = NULL,
+                              outdir = outdir, method = method,
+                              output_conf = output_conf)
+          }
         }
-      }
-
+        # Plots with dotted negative correlated in ens-mean-fcst
+        if (recipe$Analysis$Workflow$Visualization$dots %in% c('both', TRUE)) {
+          if (is.null(skill_metrics)) {
+            error(recipe$Run$logger,
+                  paste0("For the forecast ensemble mean plot, skill_metrics ",
+                         "need to be provided for the dots"))
+          } else if (!('enscorr' %in% names(skill_metrics))) {
+            error(recipe$Run$logger,
+                  paste0("For the forecast ensemble mean plot, enscor metric ",
+                         "needs to be provided for the dots"))
+          } else {
+            plot_forecast_map(recipe, data$fcst,
+                              mask = NULL,
+                              dots = skill_metrics$enscorr,
+                              outdir = outdir, method = method,
+                              output_conf = output_conf)
+          }
+        }
+      } ## End loop over methods
     } else {
       error(recipe$Run$logger,
             paste0("The forecast ensemble mean plot has been requested, but ",
-- 
GitLab


From ec1530c71967402d15415007ae1178c1933f1882 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Wed, 8 Jan 2025 11:25:16 +0100
Subject: [PATCH 02/24] Fix bug, replace 'Ensemble Mean' with method-specific
 title string

---
 modules/Visualization/R/plot_ensemble_mean.R | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/modules/Visualization/R/plot_ensemble_mean.R b/modules/Visualization/R/plot_ensemble_mean.R
index 0657b25d..c4c375ae 100644
--- a/modules/Visualization/R/plot_ensemble_mean.R
+++ b/modules/Visualization/R/plot_ensemble_mean.R
@@ -53,9 +53,9 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
   if (method == 'mean') {
     forecast_product <- s2dv::MeanDims(fcst$data, 'ensemble')
   } else if (method == 'iqr') {
-    forecast_produc <- Apply(fcst$data, target_dim = 'ensemble',
-                             fun = function(x) {
-                             IQR(x)})$output1
+    forecast_product <- Apply(fcst$data, target_dim = 'ensemble',
+                              fun = function(x) {
+                              IQR(x)})$output1
   } else {
     ## return median by default?
     forecast_product <- Apply(fcst$data, target_dim = 'ensemble',
@@ -164,11 +164,11 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
         }
         if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
           toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
-                             "\n", "Forecast Ensemble Mean / ", "Issued on ", 
+                             "\n", "Forecast Ensemble", method, " / ", "Issued on ",
                              format(ymd(start_date), "%d-%m-%Y"))
         } else {
           toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
-                             "\n", "Forecast Ensemble Mean / ", "Init.: ", i_syear)
+                             "\n", "Forecast Ensemble ", method, "/ ", "Init.: ", i_syear)
         }
         # Plots
         output_configuration <- output_conf$multipanel$forecast_ensemble_mean
@@ -231,7 +231,7 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
           if (tolower(recipe$Analysis$Horizon) == 'seasonal') {
             toptitle <- paste0(system_name, " / ",
                                str_to_title(var_long_name),
-                               "\n", "Ensemble Mean / ",
+                               "\n", "Ensemble ", method, " / ",
                                time_labels[i], " ", years[i],
                                " / Start date: ",
                                format(as.Date(i_syear, format = "%Y%m%d"),
@@ -240,7 +240,7 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
           } else if (tolower(recipe$Analysis$Horizon) == 'subseasonal') {
             toptitle <- paste0(system_name, " / ",
                                str_to_title(var_long_name),
-                               "\n", "Ensemble Mean / ",
+                               "\n", "Ensemble ", method, " / ",
                                "Issued on ", 
                                format(ymd(start_date), "%d-%m-%Y"),
                                "\n", time_labels[i], years[i])
@@ -248,7 +248,7 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
           } else {
             toptitle <- paste0(system_name, " / ",
                                str_to_title(var_long_name),
-                               "\n", "Ensemble Mean / ",
+                               "\n", "Ensemble ", method, " / ",
                                time_labels[i], " ", years[i],
                                " / Start date: ",
                                i_syear)
-- 
GitLab


From 7e0da484fa5a814db9a5a38d808e9efd926c4848 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Wed, 8 Jan 2025 12:54:34 +0100
Subject: [PATCH 03/24] Add 'method' parameter for forecast map

---
 recipes/atomic_recipes/recipe_test_multivar.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/recipes/atomic_recipes/recipe_test_multivar.yml b/recipes/atomic_recipes/recipe_test_multivar.yml
index 7de6bc19..1e476f4e 100644
--- a/recipes/atomic_recipes/recipe_test_multivar.yml
+++ b/recipes/atomic_recipes/recipe_test_multivar.yml
@@ -47,6 +47,7 @@ Analysis:
       save: 'all'
     Visualization:
       plots: most_likely_terciles, skill_metrics, forecast_ensemble_mean
+      forecast_method: mean, median
       multi_panel: no
       projection: cylindrical_equidistant
       dots: both
-- 
GitLab


From 052fce1a272b257ecc2432aef6f5963c4e4e8989 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Fri, 10 Jan 2025 16:40:07 +0100
Subject: [PATCH 04/24] Add new function 'get_variable_longname

---
 .../Visualization/R/get_variable_longname.R   | 28 +++++++++++++++++++
 modules/Visualization/R/plot_ensemble_mean.R  |  4 ++-
 modules/Visualization/R/plot_metrics.R        |  7 +++--
 .../R/plot_most_likely_terciles_map.R         |  4 ++-
 modules/Visualization/Visualization.R         |  4 +--
 5 files changed, 40 insertions(+), 7 deletions(-)
 create mode 100644 modules/Visualization/R/get_variable_longname.R

diff --git a/modules/Visualization/R/get_variable_longname.R b/modules/Visualization/R/get_variable_longname.R
new file mode 100644
index 00000000..e655f10e
--- /dev/null
+++ b/modules/Visualization/R/get_variable_longname.R
@@ -0,0 +1,28 @@
+#' Return the longname of the variable.
+#' The function attempts to retrieve the longname from the 
+#' modules/Visualization/var_long_names.yml configuration file,
+#' using the entry specified in 'format'.
+#' If the file is not found or does not contain the expected data,
+#' the longname from the s2dv_cube is used instead.
+#'
+#'@param data_cube the s2dv_cube object or its metadata, containing the 
+#'  necessary information.
+#'@param index The index for the variable
+#'@param format Character string indicating the name of the format key.
+#'
+#'@return Character string of the variable longname
+.get_variable_longname <- function(data_cube, index, format = NULL) {
+  var_name <- data_cube$attrs$Variable$varName[[index]]
+  var_long_name <- tryCatch(
+    expr = {
+      var_long_name <- read_yaml("modules/Visualization/var_long_names.yml")[[format]][[var_name]]
+      return(var_long_name)
+    },
+    error = function(e) {
+      var_long_name <- str_to_title(data_cube$attrs$Variable$metadata[[var_name]]$long_name)
+      warning("No configuration file found for this variable long name. Using s2dv_cube metadata.")
+      return(var_long_name)
+    }
+  )
+  return(var_long_name)
+}
diff --git a/modules/Visualization/R/plot_ensemble_mean.R b/modules/Visualization/R/plot_ensemble_mean.R
index c4c375ae..77421651 100644
--- a/modules/Visualization/R/plot_ensemble_mean.R
+++ b/modules/Visualization/R/plot_ensemble_mean.R
@@ -65,7 +65,9 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
   # Loop over variable dimension
   for (var in 1:fcst$dims[['var']]) {
     variable <- fcst$attrs$Variable$varName[[var]]
-    var_long_name <- fcst$attrs$Variable$metadata[[variable]]$long_name
+    var_long_name <- .get_variable_longname(data_cube = fcst,
+                                            index = var,
+                                            format = NULL) ## TODO: change
     units <- fcst$attrs$Variable$metadata[[variable]]$units
     # Subset ensemble mean by variable
     var_fcst <- ClimProjDiags::Subset(forecast_product,
diff --git a/modules/Visualization/R/plot_metrics.R b/modules/Visualization/R/plot_metrics.R
index 9850b2c7..875366ba 100644
--- a/modules/Visualization/R/plot_metrics.R
+++ b/modules/Visualization/R/plot_metrics.R
@@ -2,7 +2,7 @@ library(stringr)
 library(lubridate)
 
 plot_metrics <- function(recipe, data_cube, metrics,
-                          outdir, significance = F, output_conf) {
+                         outdir, significance = F, output_conf) {
   # recipe: Auto-S2S recipe
   # archive: Auto-S2S archive
   # data_cube: s2dv_cube object with the corresponding hindcast data
@@ -223,8 +223,9 @@ plot_metrics <- function(recipe, data_cube, metrics,
         
         # Get variable name and long name
         var_name <- data_cube$attrs$Variable$varName[[var]]
-        var_long_name <- data_cube$attrs$Variable$metadata[[var_name]]$long_name
-        
+        var_long_name <- .get_variable_longname(data_cube = data_cube,
+                                                index = var,
+                                                format = NULL) ## TODO: change
         # Multi-panel or single-panel plots
         if (recipe$Analysis$Workflow$Visualization$multi_panel) {
           # Define titles
diff --git a/modules/Visualization/R/plot_most_likely_terciles_map.R b/modules/Visualization/R/plot_most_likely_terciles_map.R
index c97da758..67236328 100644
--- a/modules/Visualization/R/plot_most_likely_terciles_map.R
+++ b/modules/Visualization/R/plot_most_likely_terciles_map.R
@@ -68,7 +68,9 @@ plot_most_likely_terciles <- function(recipe,
   for (var in 1:fcst$dims[['var']]) {
     ## NOTE: Variables starting with var_* represent the variable counter.
     variable <- fcst$attrs$Variable$varName[[var]]
-    var_long_name <- fcst$attrs$Variable$metadata[[variable]]$long_name
+    var_long_name <- .get_variable_longname(data_cube = fcst,
+                                            index = var,
+                                            format = NULL) ## TODO: change
     # Choose colors depending on the variable
     if (variable %in% c('prlr')) { ## add others
       cols <- list(c("#FFC473", "#FFAB38"), 
diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R
index 091b55f1..181c65ef 100644
--- a/modules/Visualization/Visualization.R
+++ b/modules/Visualization/Visualization.R
@@ -7,6 +7,7 @@
 source("modules/Visualization/R/plot_metrics.R")
 source("modules/Visualization/R/get_proj_code.R")
 source("modules/Visualization/R/get_plot_time_labels.R")
+source("modules/Visualization/R/get_variable_longname.R")
 ## TODO: Remove after the next s2dv release
 source("modules/Visualization/R/tmp/PlotRobinson.R")
 source("modules/Visualization/R/plot_most_likely_terciles_map.R")
@@ -141,8 +142,7 @@ Visualization <- function(recipe,
       }
       if (is.null(recipe$Analysis$Workflow$Visualization$dots)) {
         recipe$Analysis$Workflow$Visualization$forecast_method <- "median"
-      }
-      
+      } 
       # Loop over the methods requested
       for (method in strsplit(recipe$Analysis$Workflow$Visualization$forecast_method, ", | |,")[[1]]) {
         # Plot without mask or dots
-- 
GitLab


From 13fd930ab188b7c9ed3180e786275450aac8caf6 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Mon, 13 Jan 2025 14:45:56 +0100
Subject: [PATCH 05/24] Use VizEquiMap() and VizRobinson()

---
 example_scripts/test_seasonal.R               |    3 +-
 modules/Visualization/R/plot_ensemble_mean.R  |   20 +-
 modules/Visualization/R/plot_metrics.R        |   25 +-
 modules/Visualization/R/tmp/ClimPalette.R     |   79 +
 .../Visualization/R/tmp/ColorBarContinuous.R  |  594 +++++++
 modules/Visualization/R/tmp/VizEquiMap.R      | 1445 +++++++++++++++++
 modules/Visualization/R/tmp/VizRobinson.R     |  567 +++++++
 modules/Visualization/R/tmp/zzz.R             |  228 +--
 modules/Visualization/Visualization.R         |    9 +-
 9 files changed, 2836 insertions(+), 134 deletions(-)
 create mode 100644 modules/Visualization/R/tmp/ClimPalette.R
 create mode 100644 modules/Visualization/R/tmp/ColorBarContinuous.R
 create mode 100644 modules/Visualization/R/tmp/VizEquiMap.R
 create mode 100644 modules/Visualization/R/tmp/VizRobinson.R

diff --git a/example_scripts/test_seasonal.R b/example_scripts/test_seasonal.R
index 2c7d673a..bd5e9f03 100644
--- a/example_scripts/test_seasonal.R
+++ b/example_scripts/test_seasonal.R
@@ -30,4 +30,5 @@ skill_metrics <- Skill(recipe, data)
 # Compute percentiles and probability bins
 probabilities <- Probabilities(recipe, data)
 # Plot data
-Visualization(recipe, data, skill_metrics, probabilities, significance = T)
+Visualization(recipe, data, skill_metrics = skill_metrics, probabilities = probabilities,
+              significance = T)
diff --git a/modules/Visualization/R/plot_ensemble_mean.R b/modules/Visualization/R/plot_ensemble_mean.R
index 77421651..de2de50b 100644
--- a/modules/Visualization/R/plot_ensemble_mean.R
+++ b/modules/Visualization/R/plot_ensemble_mean.R
@@ -174,7 +174,7 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
         }
         # Plots
         output_configuration <- output_conf$multipanel$forecast_ensemble_mean
-        base_args <- list(fun = "PlotEquiMap",
+        base_args <- list(fun = "VizEquiMap",
                           plot_dims = c('longitude', 'latitude'),
                           var = i_var_fcst, lon = longitude,
                           lat = latitude, mask = var_mask, dots = var_dots,
@@ -190,7 +190,7 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
       } else {
         # Define function and parameters depending on projection
         if (projection == 'cylindrical_equidistant') {
-          fun <- PlotEquiMap
+          fun <- VizEquiMap
           base_args <- list(var = NULL, dots = NULL, mask = NULL,
                             lon = longitude, lat = latitude,
                             dot_symbol = 20, title_scale = 0.6,
@@ -257,15 +257,13 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
             forecast_time_caption <- paste0("Forecast month: ", sprintf("%02d", i))
           }
           # Define caption
-          if (identical(fun, PlotRobinson)) {
-            ## TODO: Customize technical details
-            base_args[['caption']] <-
-              paste0("Nominal start date: ", format(ymd(start_date), 
-                                                    "%d-%m-%Y"), "\n",
-                     forecast_time_caption, "\n",
-                     "Reference: ", recipe$Analysis$Datasets$Reference, "\n",
-                     "Units: ", units)
-          }
+          ## TODO: Customize technical details
+          base_args[['caption']] <-
+            paste0("Nominal start date: ", format(ymd(start_date), 
+                                                  "%d-%m-%Y"), "\n",
+                   forecast_time_caption, "\n",
+                   "Reference: ", recipe$Analysis$Datasets$Reference, "\n",
+                   "Units: ", units)
           # Modify base arguments
           base_args[[1]] <- i_var_fcst[i, , ]
           if (is.null(attributes(fcst$attrs$time_bounds))) {
diff --git a/modules/Visualization/R/plot_metrics.R b/modules/Visualization/R/plot_metrics.R
index 875366ba..4c39b6d9 100644
--- a/modules/Visualization/R/plot_metrics.R
+++ b/modules/Visualization/R/plot_metrics.R
@@ -262,7 +262,7 @@ plot_metrics <- function(recipe, data_cube, metrics,
           }
           ## TODO: Combine PlotLayout with PlotRobinson?
           output_configuration <- output_conf$multipanel$skill_metrics
-          base_args <- list(fun = "PlotEquiMap",
+          base_args <- list(fun = "VizEquiMap",
                             plot_dims = c('longitude', 'latitude'),
                             var = asplit(metric, MARGIN = 1),
                             lon = longitude, lat = latitude,
@@ -281,7 +281,7 @@ plot_metrics <- function(recipe, data_cube, metrics,
         } else {
           # Define function and parameters depending on projection
           if (projection == 'cylindrical_equidistant') {
-            fun <- PlotEquiMap
+            fun <- VizEquiMap
             base_args <- list(var = NULL, dots = NULL,
                               lon = longitude, lat = latitude,
                               dot_symbol = 20, dot_size = 1,
@@ -292,7 +292,7 @@ plot_metrics <- function(recipe, data_cube, metrics,
                               bar_label_digits = 3, bar_label_scale = 1.5,
                               axes_label_scale = 1, width = 8, height = 5)
           } else {
-            fun <- PlotRobinson
+            fun <- VizRobinson
             common_projections <- c("robinson", "stereographic", "lambert_europe")
             if (projection %in% common_projections) {
               target_proj <- get_proj_code(projection)
@@ -430,7 +430,7 @@ plot_metrics <- function(recipe, data_cube, metrics,
                     base_args[[10]] <- metric_significance[[i]][[1]]
                   } else {
                     # The position of arguments in base_args requires this cond
-                    # so PlotEquiMap plots dots when requested
+                    # so VizEquiMap plots dots when requested
                     base_args[[2]] <- metric_significance[[i]][[1]]
                   }
                   sign_file_label <- '_dots'
@@ -444,15 +444,14 @@ plot_metrics <- function(recipe, data_cube, metrics,
               significance_caption <- NULL
               sign_file_label <- NULL
             }
-            if (identical(fun, PlotRobinson)) {
-              ## TODO: Customize alpha and other technical details depending on the metric
-              base_args[['caption']] <- 
-                paste0("Nominal start date: ", nominal_startdate_caption, "\n",
-                       forecast_time_caption, "\n",
-                       "Reference: ", recipe$Analysis$Datasets$Reference, "\n",
-                       "Units: ", data_cube$attrs$Variable$metadata[[var_name]]$units, "\n",
-                       significance_caption)
-            }
+            ## TODO: Customize alpha and other technical details depending on the metric
+            base_args[['caption']] <- 
+              paste0("Nominal start date: ", nominal_startdate_caption, "\n",
+                     forecast_time_caption, "\n",
+                     "Reference: ", recipe$Analysis$Datasets$Reference, "\n",
+                     "Units: ", data_cube$attrs$Variable$metadata[[var_name]]$units, "\n",
+                     significance_caption)
+            
             fileout <- paste0(outfile, "_ft", forecast_time,
                               sign_file_label, ".pdf")
             # Plot
diff --git a/modules/Visualization/R/tmp/ClimPalette.R b/modules/Visualization/R/tmp/ClimPalette.R
new file mode 100644
index 00000000..4b943f48
--- /dev/null
+++ b/modules/Visualization/R/tmp/ClimPalette.R
@@ -0,0 +1,79 @@
+#'Generate Climate Color Palettes
+#'
+#'Generates a colorblind friendly color palette with color ranges useful in 
+#'climate temperature variable plotting.
+#'
+#'@param palette A character string of palette. The current choices:
+#'  \itemize{
+#'    \item{'bluered': from blue through white to red.} 
+#'    \item{'redblue': from red through white to blue.}
+#'    \item{'yellowred': from yellow through orange to red.}
+#'    \item{'redyellow': from red through orange to yellow.}
+#'    \item{'purpleorange': from purple through white to orange.}
+#'    \item{'orangepurple': from orange through white to purple.}
+#'  }
+#'@param n A number indicating how many colors to generate.
+#'
+#'@return
+#'ClimPalette() returns the function that generates the color palette and the 
+#'attribute 'na_color'.\cr
+#'ClimColors() returns a vector of the colors.
+#'
+#'@examples
+#'lims <- seq(-1, 1, length.out = 21)
+#'
+#'cb <- ColorBarContinuous(lims, color_fun = ClimPalette('redyellow'), plot = FALSE)
+#'
+#'cols <- ClimColors(20)
+#'cb <- ColorBarContinuous(lims, cols, plot = FALSE)
+#'
+#'@importFrom grDevices colorRampPalette
+#'@export
+ClimPalette <- function(palette = "bluered") {
+  if (palette == "bluered") {
+    colorbar <- colorRampPalette(rev(c("#67001f", "#b2182b", "#d6604d",
+                                       "#f4a582", "#fddbc7", "#f7f7f7", 
+                                       "#d1e5f0", "#92c5de", "#4393c3", 
+                                       "#2166ac", "#053061")))
+    attr(colorbar, 'na_color') <- 'pink'
+  } else if (palette == "redblue") {
+    colorbar <- colorRampPalette(c("#67001f", "#b2182b", "#d6604d",
+                                   "#f4a582", "#fddbc7", "#f7f7f7", 
+                                   "#d1e5f0", "#92c5de", "#4393c3", 
+                                   "#2166ac", "#053061"))
+    attr(colorbar, 'na_color') <- 'pink'
+  } else if (palette == "yellowred") {
+    colorbar <- colorRampPalette(c("#ffffcc", "#ffeda0", "#fed976",
+                                   "#feb24c", "#fd8d3c", "#fc4e2a",
+                                   "#e31a1c", "#bd0026", "#800026"))
+    attr(colorbar, 'na_color') <- 'pink'
+  } else if (palette == "redyellow") {
+    colorbar <- colorRampPalette(rev(c("#ffffcc", "#ffeda0", "#fed976",
+                                       "#feb24c", "#fd8d3c", "#fc4e2a",
+                                       "#e31a1c", "#bd0026", "#800026")))
+    attr(colorbar, 'na_color') <- 'pink'
+  } else if (palette == "purpleorange") {
+    colorbar <- colorRampPalette(c("#2d004b", "#542789", "#8073ac",
+				   "#b2abd2", "#d8daeb", "#f7f7f7",
+				   "#fee0b6", "#fdb863", "#e08214", 
+				   "#b35806", "#7f3b08"))
+    attr(colorbar, 'na_color') <- 'pink'
+  } else if (palette == "orangepurple") {
+    colorbar <- colorRampPalette(rev(c("#2d004b", "#542789", "#8073ac",
+				       "#b2abd2", "#d8daeb", "#f7f7f7",
+				       "#fee0b6", "#fdb863", "#e08214",
+				       "#b35806", "#7f3b08")))
+    attr(colorbar, 'na_color') <- 'pink'
+  } else {
+    stop("Parameter 'palette' must be one of 'bluered', 'redblue', 'yellowred'",
+	 "'redyellow', 'purpleorange' or 'orangepurple'.")
+  }
+  colorbar
+}
+
+#'@rdname ClimPalette
+#'@export
+ClimColors <- function(n, palette = "bluered") {
+  ClimPalette(palette)(n)
+}
+
diff --git a/modules/Visualization/R/tmp/ColorBarContinuous.R b/modules/Visualization/R/tmp/ColorBarContinuous.R
new file mode 100644
index 00000000..b7727ba0
--- /dev/null
+++ b/modules/Visualization/R/tmp/ColorBarContinuous.R
@@ -0,0 +1,594 @@
+#'Draws a Continuous Color Bar
+#'
+#'Generates a color bar to use as colouring function for map plots and 
+#'optionally draws it (horizontally or vertically) to be added to map 
+#'multipanels or plots. It is possible to draw triangles at the ends of the 
+#'colour bar to represent values that go beyond the range of interest. A 
+#'number of options is provided to adjust the colours and the position and 
+#'size of the components. The drawn colour bar spans a whole figure region 
+#'and is compatible with figure layouts.\cr\cr 
+#'The generated colour bar consists of a set of breaks that define the 
+#'length(brks) - 1 intervals to classify each of the values in each of the 
+#'grid cells of a two-dimensional field. The corresponding grid cell of a 
+#'given value of the field will be coloured in function of the interval it 
+#'belongs to.\cr\cr
+#'The only mandatory parameters are 'var_limits' or 'brks' (in its second 
+#'format, see below).
+#'
+#'@param brks Can be provided in two formats:
+#'\itemize{
+#'  \item{A single value with the number of breaks to be generated 
+#'  automatically, between the minimum and maximum specified in 'var_limits' 
+#'  (both inclusive). Hence the parameter 'var_limits' is mandatory if 'brks' 
+#'  is provided with this format. If 'bar_limits' is additionally provided, 
+#'  values only between 'bar_limits' will be generated. The higher the value 
+#'  of 'brks', the smoother the plot will look.}
+#'  \item{A vector with the actual values of the desired breaks. Values will 
+#'  be reordered by force to ascending order. If provided in this format, no 
+#'  other parameters are required to generate/plot the colour bar.}
+#'}
+#'  This parameter is optional if 'var_limits' is specified. If 'brks' not 
+#'  specified but 'cols' is specified, it will take as value length(cols) + 1. 
+#'  If 'cols' is not specified either, 'brks' will take 21 as value.
+#'@param cols Vector of length(brks) - 1 valid colour identifiers, for each 
+#'  interval defined by the breaks. This parameter is optional and will be 
+#'  filled in with a vector of length(brks) - 1 colours generated with the 
+#'  function provided in 'color_fun' (\code{clim.colors} by default).\cr 'cols' 
+#'  can have one additional colour at the beginning and/or at the end with the 
+#'  aim to colour field values beyond the range of interest represented in the 
+#'  colour bar. If any of these extra colours is provided, parameter 
+#'  'triangle_ends' becomes mandatory in order to disambiguate which of the 
+#'  ends the colours have been provided for.
+#'@param vertical TRUE/FALSE for vertical/horizontal colour bar 
+#'  (disregarded if plot = FALSE).
+#'@param subsampleg The first of each subsampleg breaks will be ticked on the 
+#'  colorbar. Takes by default an approximation of a value that yields a 
+#'  readable tick arrangement (extreme breaks always ticked). If set to 0 or 
+#'  lower, no labels are drawn. See the code of the function for details or 
+#'  use 'extra_labels' for customized tick arrangements.
+#'@param bar_limits Vector of two numeric values with the extremes of the 
+#'  range of values represented in the colour bar. If 'var_limits' go beyond 
+#'  this interval, the drawing of triangle extremes is triggered at the 
+#'  corresponding sides, painted in 'col_inf' and 'col_sup'. Either of them 
+#'  can be set as NA and will then take as value the corresponding extreme in 
+#'  'var_limits' (hence a triangle end won't be triggered for these sides). 
+#'  Takes as default the extremes of 'brks' if available, else the same values 
+#'  as 'var_limits'.
+#'@param var_limits Vector of two numeric values with the minimum and maximum 
+#'  values of the field to represent. These are used to know whether to draw 
+#'  triangle ends at the extremes of the colour bar and what colour to fill 
+#'  them in with. If not specified, take the same value as the extremes of 
+#'  'brks'. Hence the parameter 'brks' is mandatory if 'var_limits' is not 
+#'  specified.
+#'@param triangle_ends Vector of two logical elements, indicating whether to 
+#'  force the drawing of triangle ends at each of the extremes of the colour 
+#'  bar. This choice is automatically made from the provided 'brks', 
+#'  'bar_limits', 'var_limits', 'col_inf' and 'col_sup', but the behaviour 
+#'  can be manually forced to draw or not to draw the triangle ends with this 
+#'  parameter. If 'cols' is provided, 'col_inf' and 'col_sup' will take 
+#'  priority over 'triangle_ends' when deciding whether to draw the triangle 
+#'  ends or not.
+#'@param col_inf Colour to fill the inferior triangle end with. Useful if 
+#'  specifying colours manually with parameter 'cols', to specify the colour 
+#'  and to trigger the drawing of the lower extreme triangle, or if 'cols' is 
+#'  not specified, to replace the colour automatically generated by ColorBar().
+#'@param col_sup Colour to fill the superior triangle end with. Useful if 
+#'  specifying colours manually with parameter 'cols', to specify the colour 
+#'  and to trigger the drawing of the upper extreme triangle, or if 'cols' is 
+#'  not specified, to replace the colour automatically generated by ColorBar().
+#'@param color_fun Function to generate the colours of the color bar. Must 
+#'  take an integer and must return as many colours. The returned colour vector 
+#'  can have the attribute 'na_color', with a colour to draw NA values. This 
+#'  parameter is set by default to ClimPalette().
+#'@param plot Logical value indicating whether to only compute its breaks and 
+#'  colours (FALSE) or to also draw it on the current device (TRUE).
+#'@param draw_ticks Whether to draw ticks for the labels along the colour bar 
+#' (TRUE) or not (FALSE). TRUE by default. Disregarded if 'plot = FALSE'.
+#'@param draw_separators Whether to draw black lines in the borders of each of 
+#'  the colour rectancles of the colour bar (TRUE) or not (FALSE). FALSE by 
+#'  default. Disregarded if 'plot = FALSE'.
+#'@param triangle_ends_scale Scale factor for the drawn triangle ends of the 
+#'  colour bar, if drawn at all. Takes 1 by default (rectangle triangle 
+#'  proportional to the thickness of the colour bar). Disregarded if 
+#'  'plot = FALSE'.
+#'@param extra_labels Numeric vector of extra labels to draw along axis of 
+#'  the colour bar. The number of provided decimals will be conserved. 
+#'  Disregarded if 'plot = FALSE'.
+#'@param title Title to draw on top of the colour bar, most commonly with the 
+#'  units of the represented field in the neighbour figures. Empty by default.
+#'@param title_scale Scale factor for the 'title' of the colour bar. 
+#'  Takes 1 by default.
+#'@param label_scale Scale factor for the labels of the colour bar. 
+#'  Takes 1 by default.
+#'@param tick_scale Scale factor for the length of the ticks of the labels 
+#'  along the colour bar. Takes 1 by default.
+#'@param extra_margin Extra margins to be added around the colour bar, 
+#'  in the format c(y1, x1, y2, x2). The units are margin lines. Takes 
+#'  rep(0, 4) by default.
+#'@param label_digits Number of significant digits to be displayed in the 
+#'  labels of the colour bar, usually to avoid too many decimal digits 
+#'  overflowing the figure region. This does not have effect over the labels 
+#'  provided in 'extra_labels'. Takes 4 by default.
+#'@param ... Arguments to be passed to the method. Only accepts the following 
+#'  graphical parameters:\cr adj ann ask bg bty cex.lab cex.main cex.sub cin 
+#'  col.axis col.lab col.main col.sub cra crt csi cxy err family fg fig fin 
+#'  font font.axis font.lab font.main font.sub lend lheight ljoin lmitre lty 
+#'  lwd mai mex mfcol mfrow mfg mkh oma omd omi page pch pin plt pty smo srt 
+#'  tck tcl usr xaxp xaxs xaxt xlog xpd yaxp yaxs yaxt ylbias ylog.\cr For more 
+#' information about the parameters see `par`.
+#'
+#'@return 
+#'\item{brks}{
+#'  Breaks used for splitting the range in intervals.
+#'}
+#'\item{cols}{
+#'  Colours generated for each of the length(brks) - 1 intervals. 
+#'  Always of length length(brks) - 1.
+#'}
+#'\item{col_inf}{
+#'  Colour used to draw the lower triangle end in the colour 
+#'  bar (NULL if not drawn at all).
+#'}
+#'\item{col_sup}{
+#'  Colour used to draw the upper triangle end in the colour 
+#'  bar (NULL if not drawn at all).
+#'}
+#'
+#'@examples
+#'cols <- c("dodgerblue4", "dodgerblue1", "forestgreen", "yellowgreen", "white",
+#'          "white", "yellow", "orange", "red", "saddlebrown")
+#'lims <- seq(-1, 1, 0.2)
+#'cb <- ColorBarContinuous(lims, cols, plot = FALSE)
+#'
+#'@importFrom grDevices col2rgb rgb
+#'@import utils
+#'@export
+ColorBarContinuous <- function(brks = NULL, cols = NULL, vertical = TRUE, 
+                     subsampleg = NULL, bar_limits = NULL, var_limits = NULL, 
+                     triangle_ends = NULL, col_inf = NULL, col_sup = NULL, 
+                     color_fun = ClimPalette(), plot = TRUE, 
+                     draw_ticks = TRUE, draw_separators = FALSE, 
+                     triangle_ends_scale = 1, extra_labels = NULL, 
+                     title = NULL, title_scale = 1, 
+                     label_scale = 1, tick_scale = 1, 
+                     extra_margin = rep(0, 4), label_digits = 4, ...) {
+  # Required checks
+  if ((is.null(brks) || length(brks) < 2) && is.null(bar_limits) && is.null(var_limits)) {
+    stop("At least one of 'brks' with the desired breaks, 'bar_limits' or ",
+         "'var_limits' must be provided to generate the colour bar.")
+  }
+
+  # Check brks
+  if (!is.null(brks)) {
+    if (!is.numeric(brks)) {
+      stop("Parameter 'brks' must be numeric if specified.")
+    } else if (length(brks) > 1) {
+      reorder <- sort(brks, index.return = TRUE)
+      if (!is.null(cols)) {
+        cols <- cols[reorder$ix[which(reorder$ix <= length(cols))]]
+      }
+      brks <- reorder$x
+    } 
+  }
+
+  # Check bar_limits
+  if (!is.null(bar_limits)) {
+    if (!(all(is.na(bar_limits) | is.numeric(bar_limits)) && (length(bar_limits) == 2))) {
+      stop("Parameter 'bar_limits' must be a vector of two numeric elements or NAs.")
+    }
+  }
+
+  # Check var_limits
+  if (!is.null(var_limits)) {
+    if (!(is.numeric(var_limits) && (length(var_limits) == 2))) {
+      stop("Parameter 'var_limits' must be a numeric vector of length 2.")
+    } else if (anyNA(var_limits)) {
+      stop("Parameter 'var_limits' must not contain NA values.")
+    } else if (any(is.infinite(var_limits))) {
+      stop("Parameter 'var_limits' must not contain infinite values.")
+    }
+  }
+
+  # Check cols
+  if (!is.null(cols)) {
+    if (!is.character(cols)) {
+      stop("Parameter 'cols' must be a vector of character strings.")
+    } else if (any(!sapply(cols, .IsColor))) {
+      stop("Parameter 'cols' must contain valid colour identifiers.")
+    }
+  }
+
+  # Check color_fun
+  if (!is.function(color_fun)) {
+    stop("Parameter 'color_fun' must be a colour-generator function.")
+  }
+
+  # Check integrity among brks, bar_limits and var_limits
+  if (is.null(brks) || (length(brks) < 2)) {
+    if (is.null(brks)) {
+      if (is.null(cols)) {
+        brks <- 21
+      } else {
+        brks <- length(cols) + 1
+      }
+    }
+    if (is.null(bar_limits) || anyNA(bar_limits)) {
+      # var_limits is defined
+      if (is.null(bar_limits)) {
+        bar_limits <- c(NA, NA)
+      }
+      half_width <- 0.5 * (var_limits[2] - var_limits[1]) / (brks - 1)  
+      bar_limits[which(is.na(bar_limits))] <- c(var_limits[1] - half_width, var_limits[2] + half_width)[which(is.na(bar_limits))]
+      brks <- seq(bar_limits[1], bar_limits[2], length.out = brks)
+    } else if (is.null(var_limits)) {
+      # bar_limits is defined
+      var_limits <- bar_limits
+      half_width <- 0.5 * (var_limits[2] - var_limits[1]) / (brks - 1)  
+      brks <- seq(bar_limits[1], bar_limits[2], length.out = brks)
+      var_limits[1] <- var_limits[1] + half_width / 50
+    } else {
+      # both bar_limits and var_limits are defined
+      brks <- seq(bar_limits[1], bar_limits[2], length.out = brks)
+    }
+  } else if (is.null(bar_limits)) {
+    if (is.null(var_limits)) {
+      # brks is defined
+      bar_limits <- c(head(brks, 1), tail(brks, 1))
+      var_limits <- bar_limits
+      half_width <- 0.5 * (var_limits[2] - var_limits[1]) / (length(brks) - 1)  
+      var_limits[1] <- var_limits[1] + half_width / 50
+    } else {
+      # brks and var_limits are defined
+      bar_limits <- c(head(brks, 1), tail(brks, 1))
+    }
+  } else {
+    # brks and bar_limits are defined
+    # or
+    # brks, bar_limits and var_limits are defined
+    if (head(brks, 1) != bar_limits[1] || tail(brks, 1) != bar_limits[2]) {
+      stop("Parameters 'brks' and 'bar_limits' are inconsistent.")
+    }
+  }
+
+  # Check col_inf
+  if (!is.null(col_inf)) {
+    if (!.IsColor(col_inf)) {
+      stop("Parameter 'col_inf' must be a valid colour identifier.")
+    }
+  }  
+
+  # Check col_sup
+  if (!is.null(col_sup)) {
+    if (!.IsColor(col_sup)) {
+      stop("Parameter 'col_sup' must be a valid colour identifier.")
+    }
+  }
+
+  # Check triangle_ends
+  if (!is.null(triangle_ends) && (!is.logical(triangle_ends) || length(triangle_ends) != 2)) {
+    stop("Parameter 'triangle_ends' must be a logical vector with two elements.")
+  }
+  teflc <- triangle_ends_from_limit_cols <- c(!is.null(col_inf), !is.null(col_sup))
+  if (is.null(triangle_ends) && is.null(col_inf) && is.null(col_sup)) {
+    triangle_ends <- c(FALSE, FALSE)
+    if (bar_limits[1] >= var_limits[1]) {
+      triangle_ends[1] <- TRUE
+    }
+    if (bar_limits[2] < var_limits[2]) {
+      triangle_ends[2] <- TRUE
+    }
+  } else if (!is.null(triangle_ends) && is.null(col_inf) && is.null(col_sup)) {
+    triangle_ends <- triangle_ends
+  } else if (is.null(triangle_ends) && (!is.null(col_inf) || !is.null(col_sup))) {
+    triangle_ends <- teflc
+  } else if (any(teflc != triangle_ends)) {
+    if (!is.null(brks) && length(brks) > 1 && !is.null(cols) && length(cols) >= length(brks)) {
+      triangle_ends <- teflc
+    } else if (!is.null(cols)) {
+      triangle_ends <- teflc
+    } else {
+      triangle_ends <- triangle_ends
+    }
+  }
+  if (plot && !is.null(var_limits)) {
+    if ((bar_limits[1] > var_limits[1]) && !triangle_ends[1]) {
+      warning("There are variable values smaller than the lower limit ",
+               "of the colour bar and the lower triangle end has been ",
+               "disabled. These will be painted in the colour for NA values.")
+    }
+    if ((bar_limits[2] < var_limits[2]) && !triangle_ends[2]) {
+      warning("There are variable values greater than the higher limit ",
+               "of the colour bar and the higher triangle end has been ",
+               "disabled. These will be painted in the colour for NA values.")
+    }
+  }
+
+  # Generate colours if needed
+  if (is.null(cols)) {
+    cols <- color_fun(length(brks) - 1 + sum(triangle_ends))
+    attr_bk <- attributes(cols)
+    if (triangle_ends[1]) {
+      if (is.null(col_inf)) col_inf <- head(cols, 1)
+      cols <- cols[-1]
+    }
+    if (triangle_ends[2]) {
+      if (is.null(col_sup)) col_sup <- tail(cols, 1)
+      cols <- cols[-length(cols)]
+    }
+    attributes(cols) <- attr_bk
+  } else if ((length(cols) != (length(brks) - 1))) {
+    stop("Incorrect number of 'brks' and 'cols'. There must be one more break than the number of colours.")
+  }
+
+  # Check vertical
+  if (!is.logical(vertical)) {
+    stop("Parameter 'vertical' must be TRUE or FALSE.")
+  }
+
+  # Check extra_labels
+  if (is.null(extra_labels)) {
+    extra_labels <- numeric(0)
+  }
+  if (!is.numeric(extra_labels)) {
+    stop("Parameter 'extra_labels' must be numeric.")
+  } else {
+    if (any(extra_labels > bar_limits[2]) || any(extra_labels < bar_limits[1])) {
+      stop("Parameter 'extra_labels' must not contain ticks beyond the color bar limits.")
+    }
+  }
+  extra_labels <- sort(extra_labels)
+
+  # Check subsampleg
+  primes <- function(x) {
+    # Courtesy of Chase. See http://stackoverflow.com/questions/6424856/r-function-for-returning-all-factors
+    x <- as.integer(x)
+    div <- seq_len(abs(x))
+    factors <- div[x %% div == 0L]
+    factors <- list(neg = -factors, pos = factors)
+    return(factors)
+  }
+  remove_final_tick <- FALSE
+  added_final_tick <- TRUE
+  if (is.null(subsampleg)) {
+    subsampleg <- 1
+    while (length(brks) / subsampleg > 15 - 1) {
+      next_factor <- primes((length(brks) - 1) / subsampleg)$pos
+      next_factor <- next_factor[length(next_factor) - ifelse(length(next_factor) > 2, 1, 0)]
+      subsampleg <- subsampleg * next_factor
+    }
+    if (subsampleg > (length(brks) - 1) / 4) {
+      subsampleg <- max(1, round(length(brks) / 4))
+      extra_labels <- c(extra_labels, bar_limits[2])
+      added_final_tick <- TRUE
+      if ((length(brks) - 1) %% subsampleg < (length(brks) - 1) / 4 / 2) {
+        remove_final_tick <- TRUE
+      }
+    }
+  } else if (!is.numeric(subsampleg)) {
+    stop("Parameter 'subsampleg' must be numeric.")
+  }
+  subsampleg <- round(subsampleg)
+  draw_labels <- TRUE
+  if ((subsampleg) < 1) {
+    draw_labels <- FALSE
+  }
+
+  # Check plot
+  if (!is.logical(plot)) {
+    stop("Parameter 'plot' must be logical.")
+  }
+
+  # Check draw_separators
+  if (!is.logical(draw_separators)) {
+    stop("Parameter 'draw_separators' must be logical.")
+  }
+
+  # Check triangle_ends_scale
+  if (!is.numeric(triangle_ends_scale)) {
+    stop("Parameter 'triangle_ends_scale' must be numeric.")
+  }
+
+  # Check draw_ticks
+  if (!is.logical(draw_ticks)) {
+    stop("Parameter 'draw_ticks' must be logical.")
+  }
+
+  # Check title
+  if (is.null(title)) {
+    title <- ''
+  }
+  if (!is.character(title)) {
+    stop("Parameter 'title' must be a character string.")
+  }
+
+  # Check title_scale
+  if (!is.numeric(title_scale)) {
+    stop("Parameter 'title_scale' must be numeric.")
+  }
+
+  # Check label_scale
+  if (!is.numeric(label_scale)) {
+    stop("Parameter 'label_scale' must be numeric.")
+  }
+
+  # Check tick_scale
+  if (!is.numeric(tick_scale)) {
+    stop("Parameter 'tick_scale' must be numeric.")
+  }
+
+  # Check extra_margin
+  if (!is.numeric(extra_margin) || length(extra_margin) != 4) {
+    stop("Parameter 'extra_margin' must be a numeric vector of length 4.")
+  }
+
+  # Check label_digits
+  if (!is.numeric(label_digits)) {
+    stop("Parameter 'label_digits' must be numeric.")
+  }
+  label_digits <- round(label_digits)
+
+  # Process the user graphical parameters that may be passed in the call
+  ## Graphical parameters to exclude
+  excludedArgs <- c("cex", "cex.axis", "col", "lab", "las", "mar", "mgp", "new", "ps")
+  userArgs <- .FilterUserGraphicArgs(excludedArgs, ...)
+
+  #
+  #  Plotting colorbar
+  # ~~~~~~~~~~~~~~~~~~~
+  #
+  if (plot) {
+    pars_to_save <- c('mar', 'cex', names(userArgs), 'mai', 'mgp', 'las', 'xpd')
+    saved_pars <- par(pars_to_save) 
+    par(mar = c(0, 0, 0, 0), cex = 1)
+    image(1, 1, t(t(1)), col = rgb(0, 0, 0, 0), axes = FALSE, xlab = '', ylab = '')
+    # Get the availale space
+    figure_size <- par('fin')
+    cs <- par('csi')
+    # This allows us to assume we always want to plot horizontally
+    if (vertical) {
+      figure_size <- rev(figure_size)
+    }
+#    pannel_to_redraw <- par('mfg')
+#    .SwitchToFigure(pannel_to_redraw[1], pannel_to_redraw[2])
+    # Load the user parameters
+    par(new = TRUE)
+    par(userArgs)
+    # Set up color bar plot region
+    margins <- c(0.0, 0, 0.0, 0)
+    cex_title <- 1 * title_scale
+    cex_labels <- 0.9 * label_scale
+    cex_ticks <- -0.3 * tick_scale
+    spaceticklab <- max(-cex_ticks, 0)
+    if (vertical) {
+      margins[1] <- margins[1] + (1.2 * cex_labels * 3 + spaceticklab) * cs
+      margins <- margins + extra_margin[c(4, 1:3)] * cs
+    } else {
+      margins[1] <- margins[1] + (1.2 * cex_labels * 1 + spaceticklab) * cs
+      margins <- margins + extra_margin * cs
+    }
+    if (title != '') {
+      margins[3] <- margins[3] + (1.0 * cex_title) * cs
+    }
+    margins[3] <- margins[3] + sqrt(figure_size[2] / (margins[1] + margins[3])) * 
+                               figure_size[2] / 6 * ifelse(title != '', 0.5, 0.8)
+    # Set side margins
+    margins[2] <- margins[2] + figure_size[1] / 16
+    margins[4] <- margins[4] + figure_size[1] / 16
+    triangle_ends_prop <- 1 / 32 * triangle_ends_scale 
+    triangle_ends_cex <- triangle_ends_prop * figure_size[2]
+    if (triangle_ends[1]) {
+      margins[2] <- margins[2] + triangle_ends_cex
+    }
+    if (triangle_ends[2]) {
+      margins[4] <- margins[4] + triangle_ends_cex
+    }
+    ncols <- length(cols)
+    # Set up the points of triangles
+    # Compute the proportion of horiz. space occupied by one plot unit
+    prop_unit <- (1 - (margins[2] + margins[4]) / figure_size[1]) / ncols
+    # Convert triangle height to plot inits
+    triangle_height <- triangle_ends_prop / prop_unit
+    left_triangle <- list(x = c(1, 1 - triangle_height, 1) - 0.5,
+                          y = c(1.4, 1, 0.6))
+    right_triangle <- list(x = c(ncols, ncols + triangle_height, ncols) + 0.5,
+                           y = c(1.4, 1, 0.6))
+    # Draw the color squares and title
+    if (vertical) {
+      par(mai = c(margins[2:4], margins[1]), 
+          mgp = c(0, spaceticklab + 0.2, 0), las = 1)
+      d <- 4
+      image(1, 1:ncols, t(1:ncols), axes = FALSE, col = cols, 
+            xlab = '', ylab = '')
+      title(ylab = title, line = cex_title * (0.2 + 0.1), cex.lab = cex_title)
+      # Draw top and bottom border lines    
+      lines(c(0.6, 0.6), c(1 - 0.5, ncols + 0.5))
+      lines(c(1.4, 1.4), c(1 - 0.5, ncols + 0.5))
+      # Rotate triangles
+      names(left_triangle) <- rev(names(left_triangle))
+      names(right_triangle) <- rev(names(right_triangle))
+    } else {
+      # The term - cex_labels / 4 * (3 / cex_labels - 1) was found by
+      # try and error
+      par(mai = margins, 
+          mgp = c(0, cex_labels / 2 + spaceticklab 
+                     - cex_labels / 4 * (3 / cex_labels - 1), 0),
+          las = 1)
+      d <- 1
+      image(1:ncols, 1, t(t(1:ncols)), axes = FALSE, col = cols, 
+            xlab = '', ylab = '')
+      title(title, line = cex_title * (0.2 + 0.1), cex.main = cex_title)
+      # Draw top and bottom border lines    
+      lines(c(1 - 0.5, ncols + 0.5), c(0.6, 0.6))
+      lines(c(1 - 0.5, ncols + 0.5), c(1.4, 1.4))
+      tick_length <- -0.4
+    }
+    # Draw the triangles
+    par(xpd = TRUE)
+    if (triangle_ends[1]) {
+      # Draw left triangle
+      polygon(left_triangle$x, left_triangle$y, col = col_inf, border = NA)
+      lines(left_triangle$x, left_triangle$y)       
+    }
+    if (triangle_ends[2]) {
+      # Draw right triangle
+      polygon(right_triangle$x, right_triangle$y, col = col_sup, border = NA)
+      lines(right_triangle$x, right_triangle$y)
+    }
+    par(xpd = FALSE)
+
+    # Put the separators
+    if (vertical) {
+      if (draw_separators) {
+        for (i in 1:(ncols - 1)) {
+          lines(c(0.6, 1.4), c(i, i) + 0.5)
+        }
+      }
+      if (draw_separators || is.null(col_inf)) {
+        lines(c(0.6, 1.4), c(0.5, 0.5))
+      }
+      if (draw_separators || is.null(col_sup)) {
+        lines(c(0.6, 1.4), c(ncols + 0.5, ncols + 0.5))
+      }
+    } else {
+      if (draw_separators) {
+        for (i in 1:(ncols - 1)) {
+          lines(c(i, i) + 0.5, c(0.6, 1.4))
+        }
+      }
+      if (draw_separators || is.null(col_inf)) {
+        lines(c(0.5, 0.5), c(0.6, 1.4))
+      }
+      if (draw_separators || is.null(col_sup)) {
+        lines(c(ncols + 0.5, ncols + 0.5), c(0.6, 1.4))
+      }
+    }
+    # Put the ticks
+    plot_range <- length(brks) - 1
+    var_range <- tail(brks, 1) - head(brks, 1)
+    extra_labels_at <- ((extra_labels - head(brks, 1)) / var_range) * plot_range + 0.5
+    at <- seq(1, length(brks), subsampleg)
+    labels <- brks[at]
+    # Getting rid of next-to-last tick if too close to last one
+    if (remove_final_tick) {
+      at <- at[-length(at)]
+      labels <- labels[-length(labels)]
+    }
+    labels <- signif(labels, label_digits)
+    if (added_final_tick) {
+      extra_labels[length(extra_labels)] <- signif(tail(extra_labels, 1), label_digits)
+    }
+    at <- at - 0.5
+    at <- c(at, extra_labels_at)
+    labels <- c(labels, extra_labels)
+    tick_reorder <- sort(at, index.return = TRUE)
+    at <- tick_reorder$x
+    if (draw_labels) {
+      labels <- labels[tick_reorder$ix]
+    } else {
+      labels <- FALSE
+    }
+    axis(d, at = at, tick = draw_ticks, labels = labels, cex.axis = cex_labels, tcl = cex_ticks)
+    par(saved_pars)
+  }
+  invisible(list(brks = brks, cols = cols, col_inf = col_inf, col_sup = col_sup))
+}
diff --git a/modules/Visualization/R/tmp/VizEquiMap.R b/modules/Visualization/R/tmp/VizEquiMap.R
new file mode 100644
index 00000000..7751203e
--- /dev/null
+++ b/modules/Visualization/R/tmp/VizEquiMap.R
@@ -0,0 +1,1445 @@
+#'Maps A Two-Dimensional Variable On A Cylindrical Equidistant Projection
+#'
+#'Map longitude-latitude array (on a regular rectangular or gaussian grid) 
+#'on a cylindrical equidistant latitude and longitude projection with coloured 
+#'grid cells. Only the region for which data has been provided is displayed. 
+#'A colour bar (legend) can be plotted and adjusted. It is possible to draw 
+#'superimposed arrows, dots, symbols, contour lines and boxes. A number of 
+#'options is provided to adjust the position, size and colour of the 
+#'components. Some parameters are provided to add and adjust the masks that
+#'include continents, oceans, and lakes. This plot function is compatible with
+#'figure layouts if colour bar is disabled.
+#'
+#'@param var Array with the values at each cell of a grid on a regular 
+#'  rectangular or gaussian grid. The array is expected to have two 
+#'  dimensions: c(latitude, longitude). Longitudes can be in ascending or 
+#'  descending order and latitudes in any order. It can contain NA values 
+#'  (coloured with 'colNA'). Arrays with dimensions c(longitude, latitude) 
+#'  will also be accepted but 'lon' and 'lat' will be used to disambiguate so 
+#'  this alternative is not appropriate for square arrays. It is allowed that
+#'  the positions of the longitudinal and latitudinal coordinate dimensions 
+#'  are interchanged.
+#'@param lon Numeric vector of longitude locations of the cell centers of the 
+#'  grid of 'var', in ascending or descending order (same as 'var'). Expected 
+#'  to be regularly spaced, within either of the ranges [-180, 180] or 
+#'  [0, 360]. Data for two adjacent regions split by the limits of the 
+#'  longitude range can also be provided, e.g. \code{lon = c(0:50, 300:360)} 
+#'  ('var' must be provided consitently).
+#'@param lat Numeric vector of latitude locations of the cell centers of the 
+#'  grid of 'var', in any order (same as 'var'). Expected to be from a regular 
+#'  rectangular or gaussian grid, within the range [-90, 90].
+#'@param varu Array of the zonal component of wind/current/other field with 
+#'  the same dimensions as 'var'. It is allowed that the positions of the 
+#'  longitudinal and latitudinal coordinate dimensions are interchanged.
+#'@param varv Array of the meridional component of wind/current/other field 
+#'  with the same dimensions as 'var'. It is allowed that the positions of the 
+#'  longitudinal and latitudinal coordinate dimensions are interchanged.
+#'@param toptitle Top title of the figure, scalable with parameter 
+#'  'title_scale'.
+#'@param sizetit Scale factor for the figure top title provided in parameter 
+#'  'toptitle'. Deprecated. Use 'title_scale' instead.
+#'@param caption A character string of the caption located at the left-bottom of
+#'  the plot. Captions with multiple lines can be constructed using string 
+#'  manipulation functions like \code{paste()} or \code{paste0()}, using 
+#'  \code{"\n"} to indicate line breaks. 
+#'@param units Title at the top of the colour bar, most commonly the units of 
+#'  the variable provided in parameter 'var'.
+#'@param brks,cols,bar_limits,triangle_ends Usually only providing 'brks' is 
+#'  enough to generate the desired colour bar. These parameters allow to 
+#'  define n breaks that define n - 1 intervals to classify each of the values 
+#'  in 'var'. The corresponding grid cell of a given value in 'var' will be 
+#'  coloured in function of the interval it belongs to. These parameters are 
+#'  sent to \code{ColorBar()} to generate the breaks and colours. Additional 
+#'  colours for values beyond the limits of the colour bar are also generated 
+#'  and applied to the plot if 'bar_limits' or 'brks' and 'triangle_ends' are 
+#'  properly provided to do so. See ?ColorBar for a full explanation.
+#'@param col_inf,col_sup,colNA Colour identifiers to colour the values in 
+#'  'var' that go beyond the extremes of the colour bar and to colour NA 
+#'  values, respectively. 'colNA' takes attr(cols, 'na_color') if available by 
+#'  default, where cols is the parameter 'cols' if provided or the vector of 
+#'  colors returned by 'color_fun'. If not available, it takes 'pink' by 
+#'  default. 'col_inf' and 'col_sup' will take the value of 'colNA' if not 
+#'  specified. See ?ColorBar for a full explanation on 'col_inf' and 'col_sup'.
+#'@param color_fun,subsampleg,bar_extra_labels,draw_bar_ticks Set of 
+#'  parameters to control the visual aspect of the drawn colour bar 
+#'  (1/3). See ?ColorBar for a full explanation.
+#'@param draw_separators,triangle_ends_scale,bar_label_digits Set of 
+#'  parameters to control the visual aspect of the drawn colour bar 
+#'  (2/3). See ?ColorBar for a full explanation.
+#'@param bar_label_scale,units_scale,bar_tick_scale,bar_extra_margin Set of  
+#'  parameters to control the visual aspect of the drawn colour bar (3/3). 
+#'  See ?ColorBar for a full explanation.
+#'@param square Logical value to choose either to draw a coloured square for 
+#'  each grid cell in 'var' (TRUE; default) or to draw contour lines and fill 
+#'  the spaces in between with colours (FALSE). In the latter case, 
+#'  'filled.continents' will take the value FALSE if not specified.
+#'@param filled.continents Colour to fill in drawn projected continents. 
+#'  If 'square = FALSE', it is set as FALSE.
+#'  If set to FALSE (default), the continents are not filled.
+#'@param filled.oceans A logical value or the color name to fill in drawn 
+#'  projected oceans. The default value is FALSE. If it is TRUE, the default
+#'  colour is "light blue".
+#'@param country.borders A logical value indicating if the country borders 
+#'  should be plotted (TRUE) or not (FALSE). It only works when 
+#'  'filled.continents' is FALSE. The default value is FALSE.
+#'@param coast_color Colour of the coast line of the drawn projected continents.
+#'   Takes the value gray(0.5) by default.
+#'@param coast_width Line width of the coast line of the drawn projected 
+#'  continents. Takes the value 1 by default.
+#'@param lake_color Colour of the lake or other water body inside continents.
+#'  The default value is NULL.
+#'@param shapefile A character string of the path to a .rds file or a list 
+#'  object containinig shape file data. If it is a .rds file, it should contain
+#'  a list. The list should contains 'x' and 'y' at least, which indicate the 
+#'  location of the shape. The default value is NULL.
+#'@param shapefile_color Line color of the shapefile.
+#'@param shapefile_lwd Line width of the shapefile. The default value is 1. 
+#'@param contours Array of same dimensions as 'var' to be added to the plot 
+#'  and displayed with contours. Parameter 'brks2' is required to define the 
+#'  magnitude breaks for each contour curve. Disregarded if 'square = FALSE'.
+#'  It is allowed that the positions of the longitudinal and latitudinal 
+#'  coordinate dimensions are interchanged.
+#'@param brks2 Vector of magnitude breaks where to draw contour curves for the 
+#'  array provided in 'contours' or if 'square = FALSE'.
+#'@param contour_lwd Line width of the contour curves provided via 'contours' 
+#'  and 'brks2', or if 'square = FALSE'.
+#'@param contour_color Line color of the contour curves provided via 'contours' 
+#'  and 'brks2', or if 'square = FALSE'.
+#'@param contour_lty Line type of the contour curves. Takes 1 (solid) by 
+#'  default. See help on 'lty' in par() for other accepted values.
+#'@param contour_draw_label A logical value indicating whether to draw the 
+#'  contour labels or not. The default value is TRUE.
+#'@param contour_label_scale Scale factor for the superimposed labels when 
+#'  drawing contour levels.
+#'@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'. It is allowed that the positions of
+#'  the longitudinal and latitudinal coordinate dimensions are interchanged.
+#'@param dot_symbol Single character/number or vector of characters/numbers 
+#'  that correspond to each of the symbol layers specified in parameter 'dots'. 
+#'  If a single value is specified, it will be applied to all the layers in 
+#'  'dots'. Takes 4 (cross) by default. See 'pch' in par() for 
+#'  additional accepted options.
+#'@param dot_size Scale factor for the dots/symbols to be plotted, specified 
+#'  in 'dots'. If a single value is specified, it will be applied to all 
+#'  layers in 'dots'. Takes 1 by default.
+#'@param mask An array with the same dimensions as 'var' of [0, 1] or logical
+#'  indicating the grids to not plot data. The value 0 or FALSE is the point not
+#'  to be plotted.
+#'@param mask_color Color of the mask. The default value is 'white'.
+#'@param arr_subsamp Subsampling factor to select a subset of arrows in 
+#'  'varu' and 'varv' to be drawn. Only one out of arr_subsamp arrows will 
+#'  be drawn. Takes 1 by default.
+#'@param arr_scale Scale factor for drawn arrows from 'varu' and 'varv'. 
+#'  Takes 1 by default.
+#'@param arr_ref_len Length of the refence arrow to be drawn as legend at the 
+#'  bottom of the figure (in same units as 'varu' and 'varv', only affects the
+#'  legend for the wind or variable in these arrays). Defaults to 15.
+#'@param arr_units Units of 'varu' and 'varv', to be drawn in the legend. 
+#'  Takes 'm/s' by default.
+#'@param arr_scale_shaft Parameter for the scale of the shaft of the arrows 
+#'  (which also depend on the number of figures and the arr_scale parameter). 
+#'  Defaults to 1.
+#'@param arr_scale_shaft_angle Parameter for the scale of the angle of the 
+#'  shaft of the arrows (which also depend on the number of figure and the 
+#'  arr_scale parameter). Defaults to 1.
+#'@param axelab Whether to draw longitude and latitude axes or not. 
+#'  TRUE by default.
+#'@param labW Whether to label the longitude axis with a 'W' instead of minus 
+#'  for negative values. Defaults to FALSE.
+#'@param lab_dist_x A numeric of the distance of the longitude labels to the 
+#'  box borders. The default value is NULL and is automatically adjusted by 
+#'  the function.
+#'@param lab_dist_y A numeric of the distance of the latitude labels to the 
+#'  box borders. The default value is NULL and is automatically adjusted by 
+#'  the function.
+#'@param degree_sym A logical indicating whether to include degree symbol 
+#'  (30° N) or not (30N; default).
+#'@param intylat Interval between latitude ticks on y-axis, in degrees. 
+#'  Defaults to 20.
+#'@param intxlon Interval between latitude ticks on x-axis, in degrees. 
+#'  Defaults to 20.
+#'@param xlonshft A numeric of the degrees to shift the latitude ticks. The 
+#'  default value is 0.
+#'@param ylatshft A numeric of the degrees to shift the longitude ticks. The
+#'  default value is 0.
+#'@param xlabels A vector of character string of the custumized x-axis labels.
+#'  The values should correspond to each tick, which is decided by the longitude
+#'  and parameter 'intxlon'. The default value is NULL and the labels will be
+#'  automatically generated.
+#'@param ylabels A vector of character string of the custumized y-axis labels.
+#'  The values should correspond to each tick, which is decided by the latitude
+#'  and parameter 'intylat'. The default value is NULL and the labels will be
+#'  automatically generated.
+#'@param axes_tick_scale Scale factor for the tick lines along the longitude 
+#'  and latitude axes.
+#'@param axes_label_scale Scale factor for the labels along the longitude 
+#'  and latitude axes.
+#'@param drawleg Whether to plot a color bar (legend, key) or not. Defaults to 
+#'  TRUE. It is not possible to plot the colour bar if 'add = TRUE'. Use 
+#'  ColorBar() and the return values of PlotEquiMap() instead.
+#'@param vertical TRUE/FALSE for vertical/horizontal colour bar. Default is 
+#'  FALSE. Parameters 'width' and 'height' might need to be modified to
+#'  accommodate the vertical colour bar.
+#'@param boxlim Limits of a box to be added to the plot, in degrees: 
+#'  c(x1, y1, x2, y2). A list with multiple box specifications can also be 
+#'  provided.
+#'@param boxcol Colour of the box lines. A vector with a colour for each of 
+#'  the boxes is also accepted. Defaults to 'purple2'.
+#'@param boxlwd Line width of the box lines. A vector with a line width for 
+#'  each of the boxes is also accepted. Defaults to 5.
+#'@param margin_scale Scale factor for the margins around the map plot, with 
+#'  the format c(y1, x1, y2, x2). Defaults to rep(1, 4). If drawleg = TRUE, 
+#'  then margin_scale[1] is subtracted 1 unit.
+#'@param title_scale Scale factor for the figure top title. Defaults to 1.
+#'@param caption_size Scale factor for the figure caption. Default is 0.8 (1 if
+#'  vertical = TRUE).
+#'@param numbfig Number of figures in the layout the plot will be put into. 
+#'  A higher numbfig will result in narrower margins and smaller labels, 
+#'  axe labels, ticks, thinner lines, ... Defaults to 1.
+#'@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 include_lower_boundary Logical value indicating whether to include
+#'  the minimum value of the field. Takes TRUE by default.
+#'@param include_upper_boundary Logical value indicating whether to include
+#'  the maximum value of the field. Takes TRUE by default.
+#'@param \dots Arguments to be passed to the method. Only accepts the following 
+#'  graphical parameters:\cr
+#'  adj ann ask bg bty cex.sub cin col.axis col.lab col.main col.sub cra crt 
+#'  csi cxy err family fg font font.axis font.lab font.main font.sub lend 
+#'  lheight ljoin lmitre mex mfcol mfrow mfg mkh omd omi page pch pin plt 
+#'  pty smo srt tcl usr xaxp xaxs xaxt xlog xpd yaxp yaxs yaxt ylbias ylog \cr
+#'  For more information about the parameters see `par`.
+#'
+#'@return 
+#'\item{brks}{
+#'  Breaks used for colouring the map (and legend if drawleg = TRUE).
+#'}
+#'\item{cols}{
+#'  Colours used for colouring the map (and legend if drawleg = TRUE). Always 
+#'    of length length(brks) - 1.
+#'}
+#'\item{col_inf}{
+#'  Colour used to draw the lower triangle end in the colour bar (NULL if not 
+#'   drawn at all).
+#' }
+#'\item{col_sup}{
+#'  Colour used to draw the upper triangle end in the colour bar (NULL if not 
+#'   drawn at all).
+#'}
+#'
+#'@examples
+#'  \dontrun{
+#'ano <- s2dv::Ano_CrossValid(map_temp$exp, map_temp$obs, memb = FALSE,
+#'                            dat_dim = c('dat', 'member'), memb_dim = 'member')
+#'var <- s2dv::MeanDims(ano$exp, "member")
+#'lats <- attr(map_temp$exp, "Variables")$common$lat
+#'lons <- attr(map_temp$exp, "Variables")$common$lon
+#'
+#'VizEquiMap(var[1, 1, 1, 1, , ], lon = lons, lat = lats,
+#'           toptitle = 'Near-surface temperature anomaly, Nov. 2000',
+#'           filled.continents = FALSE, title_scale = 0.7,
+#'           caption = paste0("This is a test caption.", "\n",
+#'                            "This is a new line.")
+#'  }
+#'@import graphics maps utils
+#'@importFrom grDevices dev.cur dev.new dev.off gray
+#'@importFrom stats cor
+#' @importFrom s2dv InsertDim
+#'@export
+VizEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
+                        toptitle = NULL, sizetit = NULL, caption = NULL,
+                        units = NULL, brks = NULL, cols = NULL, bar_limits = NULL, 
+                        triangle_ends = NULL, col_inf = NULL, col_sup = NULL, 
+                        colNA = NULL, color_fun = ClimPalette(),
+                        square = TRUE, filled.continents = FALSE,
+                        filled.oceans = FALSE, country.borders = FALSE,
+                        coast_color = NULL, coast_width = 1, lake_color = NULL,
+                        shapefile = NULL, shapefile_color = NULL, shapefile_lwd = 1,
+                        contours = NULL, brks2 = NULL, contour_lwd = 0.5,
+                        contour_color = 'black', contour_lty = 1,
+                        contour_draw_label = TRUE, contour_label_scale = 1,
+                        dots = NULL, dot_symbol = 4, dot_size = 1,
+                        mask = NULL, mask_color = 'white',
+                        arr_subsamp = floor(length(lon) / 30), arr_scale = 1, 
+                        arr_ref_len = 15, arr_units = "m/s", 
+                        arr_scale_shaft = 1, arr_scale_shaft_angle = 1,
+                        axelab = TRUE, labW = FALSE, 
+                        lab_dist_x = NULL, lab_dist_y = NULL, degree_sym = FALSE,
+                        intylat = 20, intxlon = 20,
+                        xlonshft = 0, ylatshft = 0, xlabels = NULL, ylabels = NULL,
+                        axes_tick_scale = 1, axes_label_scale = 1,
+                        drawleg = TRUE, vertical = FALSE, subsampleg = NULL, 
+                        bar_extra_labels = NULL, draw_bar_ticks = TRUE, 
+                        draw_separators = FALSE, triangle_ends_scale = 1, 
+                        bar_label_digits = 4, bar_label_scale = 1, 
+                        units_scale = 1, bar_tick_scale = 1, 
+                        bar_extra_margin = rep(0, 4),
+                        boxlim = NULL, boxcol = 'purple2', boxlwd = 5, 
+                        margin_scale = rep(1, 4), title_scale = 1, 
+                        caption_size = 0.8, numbfig = NULL, fileout = NULL, 
+                        width = 8, height = 5, size_units = 'in', 
+                        res = 100, include_lower_boundary = TRUE,
+                        include_upper_boundary = TRUE, ...) {
+  # Process the user graphical parameters that may be passed in the call
+  ## Graphical parameters to exclude
+  excludedArgs <- c("cex", "cex.axis", "cex.lab", "cex.main", "col", "din", "fig", "fin", "lab", "las", "lty", "lwd", "mai", "mar", "mgp", "new", "oma", "ps", "tck")
+  userArgs <- .FilterUserGraphicArgs(excludedArgs, ...)
+
+  # If there is any filenames to store the graphics, process them
+  # to select the right device 
+  if (!is.null(fileout)) {
+    deviceInfo <- .SelectDevice(fileout = fileout, width = width, height = height, units = size_units, res = res)
+    saveToFile <- deviceInfo$fun
+    fileout <- deviceInfo$files
+  }
+
+  # Check lon, lat
+  if (!is.numeric(lon) || !is.numeric(lat)) {
+    stop("Parameters 'lon' and 'lat' must be numeric vectors.")
+  }
+
+  # Check var
+  if (is.null(var)) {
+    stop("Parameter 'var' cannot be NULL.")
+  }
+  if (!is.array(var)) {
+    stop("Parameter 'var' must be a numeric array.")
+  }
+
+  transpose <- FALSE
+  if (!is.null(names(dim(var)))) {
+    if (any(names(dim(var)) %in% .KnownLonNames()) &&
+        any(names(dim(var)) %in% .KnownLatNames())) {
+      lon_dim <- names(dim(var))[names(dim(var)) %in% .KnownLonNames()]
+      lat_dim <- names(dim(var))[names(dim(var)) %in% .KnownLatNames()]
+    } else {
+      names(dim(var)) <- NULL
+      lat_dim <- NULL
+      lon_dim <- NULL
+      warning("Dimension names of 'var' doesn't correspond to any coordinates names supported by s2dv package.")
+    }
+  } else {
+    lon_dim <- NULL
+    lat_dim <- NULL
+    warning("Parameter 'var' should have dimension names. Coordinates 'lon' and 'lat' have been assigned into the corresponding coordinates dimensions.")
+  }
+
+  if (length(dim(var)) > 2) {
+    if (!is.null(lon_dim) & !is.null(lat_dim)) {
+      dimnames <- names(dim(var))
+      dim(var) <- dim(var)[which((dimnames == lon_dim | dimnames == lat_dim | dim(var) != 1))]
+    } else {
+      if (all(dim(var) == 1)) {
+        dim(var) <- c(1, 1)
+      } else if (length(dim(var)[which(dim(var) > 1)]) == 2) {
+        var <- drop(var)
+      } else if (length(dim(var)[which(dim(var) > 1)]) == 1) {
+        dim(var) <- c(dim(var)[which(dim(var) > 1)], 1)
+      }
+    }
+  }
+
+  if (length(dim(var)) != 2) {
+    stop("Parameter 'var' must be a numeric array with two dimensions.")
+  }
+
+  if ((dim(var)[1] == length(lon) && dim(var)[2] == length(lat)) ||
+      (dim(var)[2] == length(lon) && dim(var)[1] == length(lat))) {
+    if (dim(var)[2] == length(lon) && dim(var)[1] == length(lat)) {
+      if (length(lon) == length(lat)) {
+        if (is.null(names(dim(var)))) {
+          warning("Parameter 'var' should have dimension names. Coordinates 'lon' and 'lat' have been assigned into the first and second dimensions.")
+        } else {
+          if (names(dim(var)[1]) == lat_dim) {
+            transpose <- TRUE
+          }
+        }
+      } else {
+        transpose <- TRUE
+      }   
+    }
+  } else {
+    stop("Parameters 'lon' and 'lat' must have as many elements as the number of cells along longitudes and latitudes in the input array 'var'.")
+  }
+
+  if (!is.null(names(dim(var)))) {
+    if (names(dim(var)[1]) == lon_dim) {
+      if (transpose) {
+        stop("Coordinates dimensions of 'var' doesn't correspond to lat or lon.")
+      }
+    } else if (names(dim(var)[2]) == lon_dim) {
+      if (!transpose) {
+        stop("Coordinates dimensions of 'var' doesn't correspond to lat or lon.")
+      }
+    }
+  }
+  
+  # Transpose the input matrices because the base plot functions work directly 
+  # with dimensions c(lon, lat).
+
+  if (transpose) {
+    var <- t(var)
+  }
+
+  transpose <- FALSE
+
+  names(dim(var)) <- c(lon_dim, lat_dim)
+  dims <- dim(var)
+
+  # Check varu and varv
+  if (!is.null(varu) && !is.null(varv)) {
+    if (!is.array(varu) || !(length(dim(varu)) == 2)) {
+      stop("Parameter 'varu' must be a numerical array with two dimensions.")
+    }
+    if (!is.array(varv) || !(length(dim(varv)) == 2)) {
+      stop("Parameter 'varv' must be a numerical array with two dimensions.")
+    }
+  } else if (!is.null(varu) || !is.null(varv)) {
+    stop("Only one of the components 'varu' or 'varv' has been provided. Both must be provided.")
+  }
+
+  if (!is.null(varu) && !is.null(varv)) {
+    if (!all(dim(varu) %in% dim(varv)) || !all(names(dim(varv)) %in% names(dim(varu)))) {
+      stop("Parameter 'varu' and 'varv' must have equal dimensions and dimension names.")
+    } else if (any(dim(varu) != dim(varv)) || any(names(dim(varv)) != names(dim(varu)))) {
+      varv <- t(varv)
+      names(dim(varv)) <- names(dim(varu))
+    }
+
+    if (is.null(lon_dim)) {
+      names(dim(varu)) <- NULL
+      names(dim(varv)) <- NULL
+    } else {
+      if (!is.null(names(dim(varu)))) {
+        if (!(lon_dim %in% names(dim(varu)) && lat_dim %in% names(dim(varu)))) {
+          stop("Parameters 'varu' and 'varv' must have same dimension names as 'var'.")
+        } else if (dim(varu)[lon_dim] != dim(var)[lon_dim] || dim(varu)[lat_dim] != dim(var)[lat_dim]) {
+          stop("Parameters 'varu' and 'varv' must have same dimensions as 'var'.")
+        }
+      } else {
+        warning("Parameters 'varu' and 'varv' should have dimension names. Coordinates 'lon' and 'lat' have been assigned into the corresponding coordinates dimensions.")
+      }
+    }
+    
+
+    if ((dim(varu)[1] == dims[1] && dim(varu)[2] == dims[2]) ||
+        (dim(varu)[2] ==  dims[1] && dim(varu)[1] == dims[2])) {
+      if (dim(varu)[2] == dims[1] && dim(varu)[1] == dims[2]) {
+        if (length(lon) == length(lat)) {
+          if (is.null(names(dim(varu)))) {
+            warning("Parameters 'varu' and 'varv' should have dimension names. Coordinates 'lon' and 'lat' have been assigned into the first and second dimensions.")
+          } else {
+            if (names(dim(varu)[1]) == lat_dim) {
+              transpose <- TRUE
+            }
+          }
+        } else {
+          transpose <- TRUE
+        }   
+      } 
+    } else {
+      stop("Parameters 'lon' and 'lat' must have as many elements as the number of cells along longitudes and latitudes in the input array 'varu' and 'varv'.")
+    }
+    
+    if (transpose) {
+      varu <- t(varu)
+      varv <- t(varv)
+    }
+
+    transpose <- FALSE
+
+  }
+
+  # Check contours
+  if (!is.null(contours)) {
+    if (!is.array(contours) || !(length(dim(contours)) == 2)) {
+      stop("Parameter 'contours' must be a numerical array with two dimensions.")
+    }
+  }
+
+
+  if (!is.null(contours)) {
+
+    if (is.null(lon_dim)) {
+      names(dim(contours)) <- NULL
+    } else {
+      if (!is.null(names(dim(contours)))) {
+        if (!(lon_dim %in% names(dim(contours)) && lat_dim %in% names(dim(contours)))) {
+          stop("Parameters 'contours' must have same dimension names as 'var'.")
+        } else if (dim(contours)[lon_dim] != dim(var)[lon_dim] || dim(contours)[lat_dim] != dim(var)[lat_dim]) {
+          stop("Parameters 'contours' must have same dimensions as 'var'.")
+        }
+      } else {
+        warning("Parameters 'contours' should have dimension names. Coordinates 'lon' and 'lat' have been assigned into the corresponding coordinates dimensions.")
+      }
+    }
+
+
+    transpose <- FALSE
+    if ((dim(contours)[1] == dims[1] && dim(contours)[2] == dims[2]) ||
+        (dim(contours)[2] ==  dims[1] && dim(contours)[1] == dims[2])) {
+      if (dim(contours)[2] == dims[1] && dim(contours)[1] == dims[2]) {
+        if (length(lon) == length(lat)) {
+          if (is.null(names(dim(contours)))) {
+            warning("Parameter 'contours' should have dimension names. Coordinates 'lon' and 'lat' have been assigned into the first and second dimensions.")
+          } else {
+            if (names(dim(contours)[1]) == lat_dim) {
+              transpose <- TRUE
+            }
+          }
+        } else {
+          transpose <- TRUE
+        }   
+      }
+    } else {
+      stop("Parameters 'lon' and 'lat' must have as many elements as the number of cells along longitudes and latitudes in the input array 'contours'.")
+    }
+
+    if (transpose) {
+      contours <- t(contours)
+    }
+
+    transpose <- FALSE
+
+  }
+
+  # Check toptitle
+  if (is.null(toptitle) || is.na(toptitle)) {
+    toptitle <- ''
+  }
+  if (!is.character(toptitle)) {
+    stop("Parameter 'toptitle' must be a character string.")
+  }
+
+  # Check sizetit
+  if (!is.null(sizetit)) {
+    warning("Parameter 'sizetit' is obsolete. Use 'title_scale' instead.")
+    if (!is.numeric(sizetit) || length(sizetit) != 1) {
+      stop("Parameter 'sizetit' must be a single numeric value.")
+    }
+    title_scale <- sizetit
+  }
+
+  # Check caption
+  if (!is.null(caption)) {
+    if (!is.character(caption)) {
+      stop("Parameter 'caption' must be a character string.")
+    } else {
+      num_lines <- length(strsplit(caption, "\n")[[1]])
+    }
+  }
+  
+  # Check include_lower_boundary and include_upper_boundary
+  if (!is.null(include_lower_boundary) && (!is.logical(include_lower_boundary) || length(include_lower_boundary) != 1)) {
+    stop("Parameter 'include_lower_boundary' must be a logical element.")
+  }
+  if (!is.null(include_upper_boundary) && (!is.logical(include_upper_boundary) || length(include_upper_boundary) != 1)) {
+    stop("Parameter 'include_upper_boundary' must be a logical element.")
+  }
+  
+  # Check vertical
+  if (!is.logical(vertical)) {
+    stop("Parameter 'vertical' must be TRUE or FALSE.")
+  }
+
+  tmp <- .create_var_limits(data = var, brks = brks, 
+                            bar_limits = bar_limits, drawleg = drawleg)
+  var_limits <- tmp$var_limits
+  drawleg <- tmp$drawleg
+
+  # Check: brks, cols, subsampleg, bar_limits, color_fun, bar_extra_labels, draw_bar_ticks
+  #        draw_separators, triangle_ends_scale, label_scale, units, units_scale, 
+  #        bar_label_digits
+  # Build: brks, cols, bar_limits, col_inf, col_sup
+  colorbar <- ColorBarContinuous(brks, cols, vertical = vertical, subsampleg, bar_limits, 
+                       var_limits, triangle_ends, col_inf, col_sup, color_fun, FALSE, 
+                       extra_labels = bar_extra_labels, draw_ticks = draw_bar_ticks,
+                       draw_separators = draw_separators, 
+                       triangle_ends_scale = triangle_ends_scale,
+                       label_scale = bar_label_scale, title = units, 
+                       title_scale = units_scale, tick_scale = bar_tick_scale,
+                       extra_margin = bar_extra_margin, label_digits = bar_label_digits)
+  brks <- colorbar$brks
+  cols <- colorbar$cols
+  col_inf <- colorbar$col_inf
+  col_sup <- colorbar$col_sup
+  bar_limits <- c(head(brks, 1), tail(brks, 1))
+
+  # Adjust 'var' values according to 'include_lower_boundary' and 'include_upper_boundary'.
+  # This adjustment ensures that, by default, values at the lower limit of the color bars ('brks[1]') are included. 
+  # Refer to issue #15 in the esviz GitLab for more details.
+  if (include_lower_boundary) {
+    var[var == head(brks, 1)] <- head(brks, 1) + head(diff(brks), 1)/10
+  }
+  if (!include_upper_boundary) {
+    var[var == tail(brks, 1)] <- tail(brks, 1) + tail(diff(brks), 1)/10
+  }  
+
+  # Check colNA
+  if (is.null(colNA)) {
+    if ('na_color' %in% names(attributes(cols))) {
+      colNA <- attr(cols, 'na_color')
+      if (!.IsColor(colNA)) {
+        stop("The 'na_color' provided as attribute of the colour vector must be a valid colour identifier.")
+      }
+    } else {
+      colNA <- 'pink'
+    }
+  } else if (!.IsColor(colNA)) {
+    stop("Parameter 'colNA' must be a valid colour identifier.")
+  }
+
+  # Check square
+  if (!is.logical(square)) {
+    stop("Parameter 'square' must be logical.")
+  }
+
+  # Check filled.continents
+  if (is.null(filled.continents)) {
+    if (!square) {
+      filled.continents <- FALSE
+    } else {
+      filled.continents <- TRUE
+    }
+  }
+  if (!.IsColor(filled.continents) && !is.logical(filled.continents)) {
+    stop("Parameter 'filled.continents' must be logical or a colour identifier.")
+  } else if (!is.logical(filled.continents)) {
+    continent_color <- filled.continents
+    filled.continents <- TRUE
+  } else {
+    continent_color <- gray(0.5)
+  }
+
+  # Check filled.oceans
+  if (!.IsColor(filled.oceans) & !is.logical(filled.oceans)) {
+    stop("Parameter 'filled.oceans' must be logical or a colour identifier.")
+  } else if (!is.logical(filled.oceans)) {
+    ocean_color <- filled.oceans
+    filled.oceans <- TRUE
+  } else if (filled.oceans) {
+    ocean_color <- "light blue"
+  }
+
+  # Check country.borders
+  if (!is.logical(country.borders)) {
+    stop("Parameter 'country.borders' must be logical.")
+  }
+
+  # Check coast_color
+  if (is.null(coast_color)) {
+    if (filled.continents) {
+      coast_color <- continent_color
+    } else {
+      coast_color <- 'black'
+    }
+  }
+  if (!.IsColor(coast_color)) {
+    stop("Parameter 'coast_color' must be a valid colour identifier.")
+  }
+
+  # Check coast_width
+  if (!is.numeric(coast_width)) {
+    stop("Parameter 'coast_width' must be numeric.")
+  }
+
+  # Check lake_color
+  if (!is.null(lake_color)) {
+    if (!.IsColor(lake_color)) {
+      stop("Parameter 'lake_color' must be a valid colour identifier.")
+    }
+  }
+
+  # Check shapefile
+  if (!is.null(shapefile)) {
+    if (is.list(shapefile)) {
+      shape <- shapefile
+      if (any(!c('x', 'y') %in% names(shape))) {
+        stop("The list names of the object in 'shapefile' .rds file should ", 
+             "have at least 'x' and 'y'.")
+      }
+      if (length(shape$x) != length(shape$y)) {
+        stop("The length of x and y in 'shapefile' list should be equal.")
+      }
+    } else if (!is.character(shapefile)) {
+      stop("Parameter 'shapefile' must be a .rds file or a list.")
+    } else {  # .rds file
+      if (!file.exists(shapefile)) {
+        stop("Parameter 'shapefile' is not a valid file.")
+      }
+      if (!grepl("\\.rds$", shapefile)) {
+        stop("Parameter 'shapefile' must be a .rds file or a list.")
+      }
+      shape <- readRDS(file = shapefile)
+      if (!is.list(shape)) {
+        stop("Parameter 'shapefile' should be a .rds file of a list object.")
+      }
+      if (any(!c('x', 'y') %in% names(shape))) {
+        stop("The list names of the object in 'shapefile' .rds file should ",
+             "have at least 'x' and 'y'.")
+      }
+      if (length(shape$x) != length(shape$y)) {
+        stop("The length of x and y in 'shapefile' list should be equal.")
+      }
+    }
+  }
+
+  # Check shapefile_col
+  if (is.null(shapefile_color)) {
+    if (filled.continents) {
+      shapefile_color <- continent_color
+    } else {
+      shapefile_color <- 'black'
+    }
+  }
+  if (!.IsColor(shapefile_color)) {
+    stop("Parameter 'shapefile_color' must be a valid colour identifier.")
+  }
+
+  # Check brks2
+  if (is.null(brks2)) {
+    if (is.null(contours)) { 
+      if (!square) {
+        brks2 <- brks
+        contours <- var 
+      }
+    } else {
+      ll <- signif(min(contours, na.rm = TRUE), 2)
+      ul <- signif(max(contours, na.rm = TRUE), 2)
+      brks2 <- signif(seq(ll, ul, length.out = length(brks)), 2)
+    }
+  }
+
+  # Check contour_lwd
+  if (!is.numeric(contour_lwd)) {
+    stop("Parameter 'contour_lwd' must be numeric.")
+  }
+
+  # Check contour_color
+  if (!.IsColor(contour_color)) {
+    stop("Parameter 'contour_color' must be a valid colour identifier.")
+  }
+
+  # Check contour_lty
+  if (!is.numeric(contour_lty) && !is.character(contour_lty)) {
+    stop("Parameter 'contour_lty' must be either a number or a character string.")
+  }
+
+  # Check contour_draw_label
+  if (!is.logical(contour_draw_label)) {
+    stop("Parameter 'contour_draw_label' must be logical.")
+  }
+
+  # Check contour_label_scale
+  if (!is.numeric(contour_label_scale)) {
+    stop("Parameter 'contour_label_scale' must be numeric.")
+  }
+
+  # Check dots
+  if (!is.null(dots)) {
+    if (!is.array(dots) || !(length(dim(dots)) %in% c(2, 3))) {
+      stop("Parameter 'dots' must be a logical array with two or three dimensions.")
+    }
+    if (length(dim(dots)) == 2) {
+      dim(dots) <- c(1, dim(dots))
+    }
+
+    if (is.null(lon_dim)) {
+      names(dim(dots)) <- NULL
+    } else {
+      if (!is.null(names(dim(dots)))) {
+        if (!(lon_dim %in% names(dim(dots)) && lat_dim %in% names(dim(dots)))) {
+          stop("Parameters 'dots' must have same dimension names as 'var'.")
+        } else if (dim(dots)[lon_dim] != dim(var)[lon_dim] || dim(dots)[lat_dim] != dim(var)[lat_dim]) {
+          stop("Parameters 'dots' must have same dimensions as 'var'.")
+        }
+      } else {
+        warning("Parameters 'dots' should have dimension names. Coordinates 'lon' and 'lat' have been assigned into the corresponding coordinates dimensions.")
+      }
+    }
+
+    transpose <- FALSE
+    if ((dim(dots)[2] == dims[1] && dim(dots)[3] == dims[2]) ||
+        (dim(dots)[3] ==  dims[1] && dim(dots)[2] == dims[2])) {
+      if (dim(dots)[3] == dims[1] && dim(dots)[2] == dims[2]) {
+        if (length(lon) == length(lat)) {
+          if (is.null(names(dim(dots)))) {
+            warning("Parameter 'dots' should have dimension names. Coordinates 'lon' and 'lat' have been assigned into the first and second dimensions.")
+          } else {
+            if (names(dim(dots)[2]) == lat_dim) {
+              transpose <- TRUE
+            }
+          }
+        } else {
+          transpose <- TRUE
+        }
+      }
+    } else {
+      stop("Parameter 'dots' must have same number of longitudes and latitudes as 'var'.")
+    }
+
+    if (transpose) {
+      dots <- aperm(dots, c(1, 3, 2))
+    }
+
+    transpose <- FALSE
+
+  }
+
+  # Check dot_symbol and dot_size
+  if (!is.null(dots)) {
+    if (!is.numeric(dot_symbol) && !is.character(dot_symbol)) {
+      stop("Parameter 'dot_symbol' must be a numeric or character string vector.")
+    }
+    if (length(dot_symbol) == 1) {
+      dot_symbol <- rep(dot_symbol, dim(dots)[1])
+    } else if (length(dot_symbol) < dim(dots)[1]) {
+      stop("Parameter 'dot_symbol' does not contain enough symbols.")
+    }
+    if (!is.numeric(dot_size)) {
+      stop("Parameter 'dot_size' must be numeric.")
+    }
+    if (length(dot_size) == 1) {
+      dot_size <- rep(dot_size, dim(dots)[1])
+    } else if (length(dot_size) < dim(dots)[1]) {
+      stop("Parameter 'dot_size' does not contain enough sizes.")
+    }
+  }
+  
+  # Check mask
+  if (!is.null(mask)) {
+    mask <- drop(mask)
+    if (!is.array(mask) || any(!names(dim(mask)) %in% c(lon_dim, lat_dim))) {
+      stop("Parameter 'mask' must have two dimensions named as the longitude and latitude dimensions in 'var'.")
+    } else {
+      if (!identical(names(dim(mask)), names(dim(var)))) {
+        mask <- aperm(mask, match(names(dim(mask)), names(dim(var))))
+      }
+    }
+    if (!identical(dim(mask), dim(var))) {
+      stop("Parameter 'mask' must have the same dimensions as 'var'.")
+    } else if (is.numeric(mask)) {
+      if (all(mask %in% c(0, 1))) {
+        mask <- array(as.logical(mask), dim = dim(mask))
+      } else {
+        stop("Parameter 'mask' must have only TRUE/FALSE or 0/1.")
+      }
+    } else if (is.logical(mask)) {
+      if (!all(mask %in% c(T, F))) {
+        stop("Parameter 'mask' must have only TRUE/FALSE or 0/1.")
+      }
+    } else {
+      stop("Parameter 'mask' must be a logical or numerical array.")
+    }
+  }
+  
+  # Check mask_color
+  if (!is.null(mask_color)) {
+    if (!.IsColor(mask_color)) {
+      stop("Parameter 'mask_color' must be a valid colour identifier.")
+    }
+  }
+
+  # Check arrow parameters
+  if (!is.numeric(arr_subsamp)) {
+    stop("Parameter 'arr_subsamp' must be numeric.")
+  }
+  if (!is.numeric(arr_scale)) {
+    stop("Parameter 'arr_scale' must be numeric.")
+  }
+  if (!is.numeric(arr_ref_len)) {
+    stop("Parameter 'arr_ref_len' must be numeric.")
+  }
+  if (!is.character(arr_units)) {
+    stop("Parameter 'arr_units' must be character.")
+  }
+  if (!is.numeric(arr_scale_shaft)) {
+    stop("Parameter 'arr_scale_shaft' must be numeric.")
+  }
+  if (!is.numeric(arr_scale_shaft_angle)) {
+    stop("Parameter 'arr_scale_shaft_angle' must be numeric.")
+  }
+
+  # Check axis parameters
+  if (!is.logical(axelab)) {
+    stop("Parameter 'axelab' must be logical.")
+  }
+  if (!is.logical(labW)) {
+    stop("Parameter 'labW' must be logical.")
+  }
+  if (!is.null(lab_dist_x)) {
+    if (!is.numeric(lab_dist_x)) {
+      stop("Parameter 'lab_dist_x' must be numeric.")
+    }
+  }
+  if (!is.null(lab_dist_y)) {
+    if (!is.numeric(lab_dist_y)) {
+      stop("Parameter 'lab_dist_y' must be numeric.")
+    }
+  }
+  if (!is.numeric(intylat)) {
+    stop("Parameter 'intylat' must be numeric.")
+  } else {
+    intylat <- round(intylat)
+  }
+  if (!is.numeric(intxlon)) {
+    stop("Parameter 'intxlon' must be numeric.")
+  } else {
+    intxlon <- round(intxlon)
+  }
+  if (!is.numeric(xlonshft) | length(xlonshft) != 1) {
+    stop("Parameter 'xlonshft' must be a number.")
+  }
+  if (!is.numeric(ylatshft) | length(ylatshft) != 1) {
+    stop("Parameter 'ylatshft' must be a number.")
+  }
+  if (!is.null(xlabels)) {
+    if (!is.character(xlabels) | !is.vector(xlabels)) {
+      stop("Parameter 'xlabels' must be a vector of character string.")
+    }
+  }
+  if (!is.null(ylabels)) {
+    if (!is.character(ylabels) | !is.vector(ylabels)) {
+      stop("Parameter 'ylabels' must be a vector of character string.")
+    }
+  }
+
+  # Check legend parameters
+  if (!is.logical(drawleg)) {
+    stop("Parameter 'drawleg' must be logical.")
+  }
+
+  # Check box parameters
+  if (!is.null(boxlim)) {
+    if (!is.list(boxlim)) {
+      boxlim <- list(boxlim)
+    }
+    for (i in 1:length(boxlim)) {
+      if (!is.numeric(boxlim[[i]]) || length(boxlim[[i]]) != 4) {
+        stop("Parameter 'boxlim' must be a a numeric vector or a list of numeric vectors of length 4 (with W, S, E, N box limits).")
+      }
+    }
+    if (!is.character(boxcol)) {
+      stop("Parameter 'boxcol' must be a character string or a vector of character strings.")
+    } else {
+      if (length(boxlim) != length(boxcol)) {
+        if (length(boxcol) == 1) {
+          boxcol <- rep(boxcol, length(boxlim))
+        } else {
+          stop("Parameter 'boxcol' must have a colour for each box in 'boxlim' or a single colour for all boxes.")
+        }
+      }
+    }
+    if (!is.numeric(boxlwd)) {
+      stop("Parameter 'boxlwd' must be numeric.")
+    } else {
+      if (length(boxlim) != length(boxlwd)) {
+        if (length(boxlwd) == 1) {
+          boxlwd <- rep(boxlwd, length(boxlim))
+        } else {
+          stop("Parameter 'boxlwd' must have a line width for each box in 'boxlim' or a single line width for all boxes.")
+        }
+      }
+    }
+  }
+
+  # Check margin_scale
+  if (!is.numeric(margin_scale) || length(margin_scale) != 4) {
+    stop("Parameter 'margin_scale' must be a numeric vector of length 4.")
+  }
+
+  # Check title_scale
+  if (!is.numeric(title_scale)) {
+    stop("Parameter 'title_scale' must be numeric.")
+  }
+  
+  # Check caption_size
+  if (!is.numeric(caption_size)) {
+    stop("Parameter 'caption_size' must be numeric.")
+  }
+  if (vertical) {
+    if (missing(caption_size)) {
+      caption_size <- 1
+    }
+  }
+
+  # Check axes_tick_scale
+  if (!is.numeric(axes_tick_scale)) {
+    stop("Parameter 'axes_tick_scale' must be numeric.")
+  }
+
+  # Check axes_label_scale
+  if (!is.numeric(axes_label_scale)) {
+    stop("Parameter 'axes_label_scale' must be numeric.")
+  }
+
+  # Check numbfig
+  if (!is.null(numbfig)) {
+    if (!is.numeric(numbfig)) {
+      stop("Parameter 'numbfig' must be numeric.")
+    } else {
+      numbfig <- round(numbfig)
+      scale <- 1 / numbfig ** 0.3
+      axes_tick_scale <- axes_tick_scale * scale
+      axes_label_scale <- axes_label_scale * scale
+      title_scale <- title_scale * scale
+      margin_scale <- margin_scale * scale
+      arr_scale <- arr_scale * scale
+      dot_size <- dot_size * scale
+      contour_label_scale <- contour_label_scale * scale
+      contour_lwd <- contour_lwd * scale
+    }
+  }
+
+  #
+  #  Input arguments 
+  # ~~~~~~~~~~~~~~~~~
+  #
+  latb <- sort(lat, index.return = TRUE)
+  dlon <- diff(lon)
+  wher <- which(dlon > (mean(dlon) + 1))
+  if (length(wher) > 0) {
+    warning("Detect gap in 'lon' vector, which is considered as crossing the border.")
+    lon[(wher + 1):dims[1]] <- lon[(wher + 1):dims[1]] - 360
+  }
+  lonb <- sort(lon, index.return = TRUE)
+  latmin <- floor(min(lat) / 10) * 10
+  latmax <- ceiling(max(lat) / 10) * 10
+  lonmin <- floor(min(lon) / 10) * 10
+  lonmax <- ceiling(max(lon) / 10) * 10
+
+  #
+  #  Plotting the map
+  # ~~~~~~~~~~~~~~~~~~
+  #
+
+  # 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)
+  }
+
+  #
+  #  Defining the layout
+  # ~~~~~~~~~~~~~~~~~~~~~
+  # 
+  if (drawleg) {
+    margin_scale[1] <- margin_scale[1] - 1 
+    
+  }
+  margins <- rep(0.4, 4) * margin_scale
+  margins[4] <- margins[4] + 1
+  cex_title <- 2 * title_scale
+  cex_axes_labels <- 1.3 * axes_label_scale
+  cex_axes_ticks <- -0.5 * axes_tick_scale
+  spaceticklab <- 0
+  if (axelab) {
+    # Y axis label
+    if (!is.null(ylabels)) {
+      ypos <- seq(latmin, latmax, intylat) + ylatshft
+      if (length(ypos) != length(ylabels)) {
+        stop(paste0("Parameter 'ylabels' must have the same length as the latitude ",
+                    "vector spaced by 'intylat' (length = ", length(ypos), ")."))
+      }
+      ylabs <- ylabels
+    } else {
+      ypos <- seq(latmin, latmax, intylat) + ylatshft
+      letters <- array('', length(ypos))
+      if (degree_sym == FALSE) {
+        letters[ypos < 0] <- 'S'
+        letters[ypos > 0] <- 'N'
+      } else {
+        letters[ypos < 0] <- paste(intToUtf8(176), 'S')
+        letters[ypos > 0] <- paste(intToUtf8(176), 'N')
+      }
+      ylabs <- paste(as.character(abs(ypos)), letters, sep = '')
+    }
+
+    # X axis label
+    if (!is.null(xlabels)) {
+      xpos <- seq(lonmin, lonmax, intxlon) + xlonshft
+      if (length(xpos) != length(xlabels)) {
+        stop(paste0("Parameter 'xlabels' must have the same length as the longitude ",
+                    "vector spaced by 'intxlon' (length = ", length(xpos), ")."))
+      }
+      xlabs <- xlabels
+    } else {
+      xpos <- seq(lonmin, lonmax, intxlon) + xlonshft
+      letters <- array('', length(xpos))
+      if (labW) {
+        xpos2 <- xpos  
+        xpos2[xpos2 > 180] <- 360 - xpos2[xpos2 > 180]
+      }
+      if (degree_sym == FALSE) {
+        letters[xpos < 0] <- 'W'
+        letters[xpos > 0] <- 'E'
+      } else {
+        letters[xpos < 0] <- paste(intToUtf8(176), 'W')
+        letters[xpos > 0] <- paste(intToUtf8(176), 'E')
+      }
+      if (labW) {
+        letters[xpos == 0] <- ' '
+        letters[xpos == 180] <- ' '
+        if (degree_sym == FALSE) {
+          letters[xpos > 180] <- 'W'
+        } else {
+          letters[xpos > 180] <- paste(intToUtf8(176), 'W')
+        }  
+        xlabs <- paste(as.character(abs(xpos2)), letters, sep = '')
+      } else {
+        xlabs <- paste(as.character(abs(xpos)), letters, sep = '')
+      }
+    }
+    spaceticklab <- max(-cex_axes_ticks, 0)
+    margins[1] <- margins[1] + 1.2 * cex_axes_labels + spaceticklab
+    margins[2] <- margins[2] + 1.2 * cex_axes_labels + spaceticklab
+  }
+  bar_extra_margin[2] <- bar_extra_margin[2] + margins[2]
+  bar_extra_margin[4] <- bar_extra_margin[4] + margins[4]
+  if (toptitle != '') {
+    margins[3] <- margins[3] + cex_title + 1
+  }
+  if (!is.null(varu)) {
+    margins[1] <- margins[1] + 2.2 * units_scale
+  }
+
+  if (drawleg) {
+    if (!is.null(caption)) {
+      margins[2] <- margins[2] + num_lines*0.5
+      margins[4] <- margins[4] + num_lines*0.5
+      if (vertical) { # vertical bar, caption
+        layout(matrix(c(1, 2, 3, 3), ncol = 2, nrow = 2, byrow = TRUE),
+               widths = c(5, 1.3),
+               heights = c(5, 0.2 + num_lines*caption_size/6))
+        
+      } else { # horizontal bar, caption
+        layout(matrix(c(1, 2, 3), ncol = 1, nrow = 3),
+               heights = c(5, 1, 0.2 + num_lines*caption_size/4))
+      }
+    } else { 
+      if (vertical) { # vertical bar, no caption
+        layout(matrix(1:2, ncol = 2, nrow = 1), widths = c(5, 1.3))
+      } else { # horizontal bar, no caption
+        layout(matrix(1:2, ncol = 1, nrow = 2), heights = c(5, 1))
+      }
+    }
+  } else {
+    if (!is.null(caption)) {
+      margins[2] <- margins[2] + num_lines*0.4
+      margins[4] <- margins[4] + num_lines*0.4
+      layout(matrix(1:2, ncol = 1, nrow = 2), heights = c(5, 0.1 + num_lines*caption_size/4))
+    }
+  }
+  plot.new()
+  # Load the user parameters
+  par(userArgs)
+  par(mar = margins, cex.main = cex_title, cex.axis = cex_axes_labels,
+      mgp = c(0, spaceticklab, 0), las = 0)
+
+  #NOTE: Here creates the window for later plot. If 'usr' for par() is not specified,
+  #      use the lat/lon as the borders. If 'usr' is specified, use the assigned values.
+  if (is.null(userArgs$usr)) { 
+    #NOTE: The grids are assumed to be equally spaced
+    xlim_cal <- c(lonb$x[1] - (lonb$x[2] - lonb$x[1]) / 2, 
+                  lonb$x[length(lonb$x)] + (lonb$x[2] - lonb$x[1]) / 2)
+    ylim_cal <- c(latb$x[1] - (latb$x[2] - latb$x[1]) / 2,
+                  latb$x[length(latb$x)] + (latb$x[2] - latb$x[1]) / 2)
+    plot.window(xlim =  xlim_cal, ylim = ylim_cal, xaxs = 'i', yaxs = 'i')
+# Below is Old code. The border grids are only half plotted.
+#    plot.window(xlim = range(lonb$x, finite = TRUE), ylim = range(latb$x, finite = TRUE),
+#                xaxs = 'i', yaxs = 'i')
+  } else {
+    plot.window(xlim = par("usr")[1:2], ylim  = par("usr")[3:4], xaxs = 'i', yaxs = 'i')
+  }
+
+  if (axelab) {
+    lab_distance_y <- ifelse(is.null(lab_dist_y), spaceticklab + 0.2, lab_dist_y)
+    lab_distance_x <- ifelse(is.null(lab_dist_x), spaceticklab + cex_axes_labels / 2 - 0.3, lab_dist_x)
+
+    axis(2, at = ypos, labels = ylabs, cex.axis = cex_axes_labels, tcl = cex_axes_ticks,
+         mgp = c(0, lab_distance_y, 0))
+    axis(1, at = xpos, labels = xlabs, cex.axis = cex_axes_labels, tcl = cex_axes_ticks,
+         mgp = c(0, lab_distance_x, 0))
+  }
+  title(toptitle, cex.main = cex_title)
+  rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col = colNA)
+  col_inf_image <- ifelse(is.null(col_inf), colNA, col_inf)
+  col_sup_image <- ifelse(is.null(col_sup), colNA, col_sup)
+  if (square) {
+    # If lat and lon are both regular-spaced, "useRaster = TRUE" can avoid
+    # artifact white lines on the figure. If not, useRaster has to be FALSE (default)
+    tryCatch({
+      image(lonb$x, latb$x, var[lonb$ix, latb$ix], 
+            col = c(col_inf_image, cols, col_sup_image), 
+            breaks = c(-.Machine$double.xmax, brks, .Machine$double.xmax),
+            axes = FALSE, xlab = "", ylab = "", add = TRUE, useRaster = TRUE)
+    }, error = function(x) {
+      image(lonb$x, latb$x, var[lonb$ix, latb$ix], 
+            col = c(col_inf_image, cols, col_sup_image), 
+            breaks = c(-.Machine$double.xmax, brks, .Machine$double.xmax),
+            axes = FALSE, xlab = "", ylab = "", add = TRUE)
+    })
+  } else {
+    .filled.contour(lonb$x, latb$x, var[lonb$ix, latb$ix], 
+                    levels = c(.Machine$double.xmin, brks, .Machine$double.xmax), 
+                    col = c(col_inf_image, cols, col_sup_image))
+  }
+  if (!is.null(contours)) {
+#NOTE: 'labcex' is the absolute size of contour labels. Parameter 'contour_label_scale'
+#      is provided in PlotEquiMap() but it was not used. Here, 'cex_axes_labels' was used
+#      and it was calculated from 'axes_label_scale', the size of lat/lon axis label.
+#      It is changed to use contour_label_scale*par('cex').
+    contour(lonb$x, latb$x, contours[lonb$ix, latb$ix], levels = brks2,
+            method = "edge", add = TRUE,
+#            labcex = cex_axes_labels,
+            labcex = contour_label_scale * par("cex"),
+            lwd = contour_lwd, lty = contour_lty,
+            col = contour_color, drawlabels = contour_draw_label)
+  }
+
+
+  #
+  #  Adding black dots or symbols
+  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  #
+  if (!is.null(dots)) {
+    data_avail <- !is.na(var)
+    for (counter in 1:(dim(dots)[1])) {
+      points <- which(dots[counter, , ] & data_avail, arr.ind = TRUE)
+      points(lon[points[, 1]], lat[points[, 2]], 
+             pch = dot_symbol[counter], 
+             cex = dot_size[counter] * 3 / sqrt(sqrt(length(var))),
+             lwd = dot_size[counter] * 3 / sqrt(sqrt(length(var))))
+    }
+  }
+
+  #
+  #  Adding a mask
+  # ~~~~~~~~~~~~~~~
+  #
+  if (!is.null(mask)){
+    
+    mask_logical <- as.logical(mask)
+    
+    for (i in 1:length(mask_logical)) {
+      if (!mask_logical[i]) {
+        
+        # Calculate the longitude and latitude indices
+        lon_idx <- (i - 1) %% length(lon) + 1
+        lat_idx <- ceiling(i / length(lon))
+        
+        # Get the longitude and latitude for this point
+        lon_current <- lon[lon_idx]
+        lat_current <- lat[lat_idx]
+        
+        # Coordinates for the corners of the rectangle (entire grid cell)
+        lon_min <- ifelse(lon_idx == 1,
+                          lon_current - (lon[2] - lon[1])/2,
+                          (lon[lon_idx] + lon[lon_idx - 1]) / 2)
+        lon_max <- ifelse(lon_idx == length(lon),
+                          lon_current + (lon[lon_idx] - lon[lon_idx - 1])/2,
+                          (lon[lon_idx + 1] + lon[lon_idx]) / 2)
+        lat_min <- ifelse(lat_idx == 1,
+                          lat_current - (lat[2] - lat[1])/2,
+                          (lat[lat_idx] + lat[lat_idx - 1]) / 2)
+        lat_max <- ifelse(lat_idx == length(lat),
+                          lat_current + (lat[lat_idx] - lat[lat_idx - 1])/2,
+                          (lat[lat_idx + 1] + lat[lat_idx]) / 2)
+        
+        # Draw a rectangle over the masked area
+        rect(lon_min, lat_min, lon_max, lat_max, col = mask_color, border = NA)
+      }
+    }
+  }
+  
+  #
+  #  Plotting continents
+  # ~~~~~~~~~~~~~~~~~~~~~
+  #
+  wrap_vec <- c(lonb$x[1], lonb$x[1] + 360)
+  old_lwd <- par('lwd')
+  par(lwd = coast_width)
+  # If [0, 360], use GEOmap; if [-180, 180], use maps::map
+  # UPDATE: Use maps::map for both cases. The difference between GEOmap and
+  #         maps is trivial. The only thing we can see for now is that 
+  #         GEOmap has better lakes.
+  coast <- maps::map(interior = country.borders, wrap = wrap_vec,
+                     fill = filled.continents, add = TRUE, plot = FALSE)
+
+  if (filled.continents) {
+    polygon(coast, col = continent_color, border = coast_color, lwd = coast_width)
+  } else {
+    lines(coast, col = coast_color, lwd = coast_width)
+  }
+  if (!is.null(lake_color)) {
+    maps::map('lakes', add = TRUE, wrap = wrap_vec, fill = filled.continents, col = lake_color) 
+  }
+  par(lwd = old_lwd)
+
+  # filled.oceans
+  if (filled.oceans) {
+      old_lwd <- par('lwd')
+      par(lwd = coast_width)
+
+      outline <- maps::map(wrap = wrap_vec, fill = T, plot = FALSE)  # must be fill = T
+      xbox <- wrap_vec + c(-2, 2)
+      ybox <- c(-92, 92) 
+      outline$x <- c(outline$x, NA, c(xbox, rev(xbox), xbox[1]))
+      outline$y <- c(outline$y, NA, rep(ybox, each = 2), ybox[1])
+      polypath(outline, col = ocean_color, rule = 'evenodd', border = NA)
+
+      par(lwd = old_lwd)
+  }
+
+  # Plot shapefile
+  #NOTE: the longitude range cannot cut shapefile range, or not all the shapefile will be plotted.
+  if (!is.null(shapefile)) {
+    maps::map(shape, interior = country.borders, #wrap = wrap_vec,
+              fill = filled.continents, add = TRUE, plot = TRUE, 
+              lwd = shapefile_lwd, col = shapefile_color)
+  }
+
+  box()
+  # Draw rectangle on the map
+  if (!is.null(boxlim)) {
+    counter <- 1
+    for (box in boxlim) {
+      if (box[1] > box[3]) {
+        box[1] <- box[1] - 360
+      }
+      if (length(box) != 4) {
+        stop(paste("The", counter, "st box defined in the parameter 'boxlim' is ill defined."))
+      } else if (box[2] < latmin || box[4] > latmax || 
+                 box[1] < lonmin || box[3] > lonmax) {
+        stop(paste("The limits of the", counter, "st box defined in the parameter 'boxlim' are invalid."))
+      } else if (box[1] < 0 && box[3] > 0) {
+        #segments south
+        segments(box[1], box[2], 0, box[2], col = boxcol[counter], lwd = boxlwd[counter])
+        segments(0, box[2], box[3], box[2], col = boxcol[counter], lwd = boxlwd[counter]) 
+        #segments north
+        segments(box[1], box[4], 0, box[4], col = boxcol[counter], lwd = boxlwd[counter])
+        segments(0, box[4], box[3], box[4], col = boxcol[counter], lwd = boxlwd[counter]) 
+        #segments west
+        segments(box[1], box[2], box[1], box[4], col = boxcol[counter], 
+                 lwd = boxlwd[counter])  
+        #segments est
+        segments(box[3], box[2], box[3],box[4], col = boxcol[counter], 
+                 lwd = boxlwd[counter])          
+      } else {
+        rect(box[1], box[2], box[3], box[4], border = boxcol[counter], col = NULL, 
+             lwd = boxlwd[counter], lty = 'solid')
+      }
+      counter <- counter + 1
+    }
+  }
+  #
+  #  PlotWind
+  # ~~~~~~~~~~
+  #
+  if (!is.null(varu) && !is.null(varv)) {
+    # Create a two dimention array of longitude and latitude
+    lontab <- InsertDim(lonb$x, 2, length(latb$x), name = 'lat')
+    lattab <- InsertDim(latb$x, 1, length(lonb$x), name = 'lon')
+    varplotu <- varu[lonb$ix, latb$ix]
+    varplotv <- varv[lonb$ix, latb$ix]
+
+    # Select a subsample af the points to an arrow
+    #for each "subsample" grid point
+    sublon <- seq(1,length(lon), arr_subsamp)
+    sublat <- seq(1,length(lat), arr_subsamp)
+
+    uaux <- lontab[sublon, sublat] + varplotu[sublon, sublat] * 0.5 * arr_scale
+    vaux <- lattab[sublon, sublat] + varplotv[sublon, sublat] * 0.5 * arr_scale
+
+    lenshaft <- 0.18 * arr_scale * arr_scale_shaft
+    angleshaft <- 12 * arr_scale_shaft_angle
+    # Plot Wind
+    arrows(lontab[sublon, sublat], lattab[sublon, sublat],
+           uaux, vaux,
+           angle = angleshaft,
+           length = lenshaft)
+    
+    # Plotting an arrow at the bottom of the plot for the legend
+    posarlon <- lonb$x[1] + (lonmax - lonmin) * 0.1
+    posarlat <- latmin - ((latmax - latmin) + 1) / par('pin')[2] * 
+                         (spaceticklab + 0.2 + cex_axes_labels + 0.6 * units_scale) * par('csi')
+
+    arrows(posarlon, posarlat,
+           posarlon + 0.5 * arr_scale * arr_ref_len, posarlat,
+           length = lenshaft, angle = angleshaft,
+           xpd = TRUE)
+    #save the parameter value
+    xpdsave <- par('xpd')
+    #desactivate xpd to be able to plot in margen
+    par(xpd = NA)
+    #plot text
+    mtext(paste(as.character(arr_ref_len), arr_units, sep = ""),
+          line = spaceticklab + 0.2 + cex_axes_labels + 1.2 * units_scale, side = 1,
+          at = posarlon + (0.5 * arr_scale * arr_ref_len) / 2,
+          cex = units_scale)
+    #come back to the previous xpd value
+    par(xpd = xpdsave)
+  }
+  
+  #
+  #  Adding a caption
+  # ~~~~~~~~~~~~~~~~~
+  #
+  if (!is.null(caption)) {
+    if (drawleg) {
+      if (vertical) {
+        par(mfg = c(2, 1))
+        base_line <- 1
+      } else {
+        par(mfg = c(3, 1))
+        base_line <- 1
+      }
+    } else {
+      par(mfg = c(2, 1))
+      base_line <- 1
+    }
+    mtext(caption, side = 1, line = base_line, adj = 0,
+          cex = caption_size, col = "black")
+  }
+  
+  #
+  #  Colorbar
+  # ~~~~~~~~~~
+  #
+  if (drawleg) {
+    if (vertical) {
+      par(mfg = c(1, 2))
+    } else {
+      if (!is.null(caption)) {
+        par(mfg = c(2, 1))
+      }
+    }
+    ColorBarContinuous(brks, cols, vertical = vertical, subsampleg, bar_limits,
+             var_limits, triangle_ends, col_inf = col_inf, col_sup = col_sup, 
+             extra_labels = bar_extra_labels, draw_ticks = draw_bar_ticks,
+             draw_separators = draw_separators, title = units, 
+             title_scale = units_scale, triangle_ends_scale = triangle_ends_scale, 
+             label_scale = bar_label_scale, tick_scale = bar_tick_scale,
+             extra_margin = bar_extra_margin, label_digits = bar_label_digits)
+  }
+
+  # If the graphic was saved to file, close the connection with the device
+  if (!is.null(fileout)) dev.off()
+
+  invisible(list(brks = brks, cols = cols, col_inf = col_inf, col_sup = col_sup))
+}
diff --git a/modules/Visualization/R/tmp/VizRobinson.R b/modules/Visualization/R/tmp/VizRobinson.R
new file mode 100644
index 00000000..63f5c7fa
--- /dev/null
+++ b/modules/Visualization/R/tmp/VizRobinson.R
@@ -0,0 +1,567 @@
+#'Plot map in Robinson or other projections
+#'
+#'Transform a regular grid longitude-latitude data to a different projection and 
+#'plot the map. The target projection must be a valid CRS string, preferrably be 
+#'EPSG or ESRI code; check \link[sf]{st_crs} for more explanation. This function
+#'is mainly tested for Robinson projection (ESRI:54030), but it can work with 
+#'other projection types in theory.\cr
+#'The map can be plotted by points or polygon. A legend can be plotted as either
+#'a color bar or a discrete ggplot legend. Dots can be drawn on top of the data,
+#'which can be used for significance test. A mask can be added to not plot the 
+#'data specified. A number of options is provided to adjust aesthetics, like 
+#'position, size, colors, etc. 
+#'
+#'@param data A numeric array with longitude and latitude dimensions. The grid
+#'  should be regular grid. It can contain NA values. 
+#'@param lon A numeric vector of longitude locations of the cell centers of the 
+#'  grid of 'data'. Expected to be regularly spaced, within the range of either
+#'  [-180, 180] or [0, 360]. 
+#'@param lat A numeric vector of latitude locations of the cell centers of the 
+#'  grid of 'data'. Expected to be regularly spaced, within the range [-90, 90]
+#'  of ascending or descending order.
+#'@param lon_dim A character string indicating the longitude dimension name in
+#'  'data'. If it is NULL, the function tries to find the name in 
+#'  \code{esviz:::.KnownLonNames}. The default value is NULL.
+#'@param lat_dim A character string indicating the latitude dimension name in
+#'  'data'. If it is NULL, the function tries to find the name in 
+#'  \code{esviz:::.KnownLatNames}. The default value is NULL.
+#'@param target_proj A character string indicating the target projection. It 
+#'  should be a valid crs string. The default projection is Robinson: 
+#'  "ESRI:54030". Note that the character string may work differently depending
+#'  on PROJ and GDAL module version. If package version 'sf' is lower than 
+#'  "1.0.10" and an error appears regarding the target crs, you can try with 
+#'  numeric crs (e.g. target_proj = 54030).
+#'@param drawleg A character string indicating the legend style. It can be 
+#'  'bar' (color bar by \code{ColorBarContinuous()}), 'ggplot2' (discrete legend 
+#'  by ggplot2), or FALSE (no legend). The default value is 'bar'.
+#'@param style A character string indicating the plotting style. It can be 
+#'  'point' or 'polygon'. The default value is 'point'. Note that 'polygon' may
+#'  be time- and memory-consuming for global or high-resolution data.
+#'@param dots An array with the same dimensions as 'data' of [0, 1] or logical
+#'  indicating the grids to plot dots. The value 0 or FALSE is the point to be
+#'  dotted.
+#'@param mask An array with the same dimensions as 'data' of [0, 1] or logical
+#'  indicating the grids to not plot data. The value 0 or FALSE is the point not
+#'  to be plotted.
+#'@param brks,cols,bar_limits,triangle_ends Usually only providing 'brks' is 
+#'  enough to generate the desired color bar. These parameters allow to 
+#'  define n breaks that define n - 1 intervals to classify each of the values 
+#'  in 'data'. The corresponding grid cell of a given value in 'data' will be 
+#'  colored in function of the interval it belongs to. These parameters are 
+#'  sent to \code{ColorBarContinuous()} to generate the breaks and colours. 
+#'  Additional colors for values beyond the limits of the colour bar are also 
+#'  generated and applied to the plot if 'bar_limits' or 'brks' and 
+#'  'triangle_ends' are properly provided to do so. See ?ColorBarContinuous for
+#'  a full explanation.
+#'@param col_inf,col_sup,colNA Colour identifiers to color the values that 
+#'  excess the extremes of the color bar and to color NAs, respectively. 'colNA'
+#'  takes attr(cols, 'na_color') if available by default, where cols is the 
+#'  parameter 'cols' if provided or the vector of colors returned by 
+#'  'color_fun'. 'col_inf' and 'col_sup' will take the value of 'colNA' if not 
+#'  specified. See ?ColorBarContinuous for a full explanation.
+#'@param color_fun,bar_extra_margin Set of 
+#'  parameters to control the visual aspect of the drawn colour bar 
+#'  (1/3). See ?ColorBarContinuous for a full explanation.
+#'@param vertical A logical value indicating the direction of colorbar if 
+#' parameter 'drawleg' is 'bar'. The default value is TRUE.
+#'@param toptitle A character string of the top title of the figure, scalable 
+#'  with parameter 'title_size'.
+#'@param caption A character string of the caption located at left-bottom of the
+#'  plot.
+#'@param units A character string of the data units, which is the title of the 
+#'  legend.
+#'@param crop_coastlines A named numeric vector [lonmin, lonmax, latmin, latmax] 
+#'  indicating the region to plot coastlines. Note that the longitude range 
+#'  cannot exceed 180 degrees.
+#'@param point_size A number of the size of the data points if "style = 'point'".
+#'  The default is 'auto' and the function tries to find the appropriate size.
+#'@param title_size A number of the size of the top title. The default is 16.
+#'@param dots_size A number of the size of the dots. The default is 0.5.
+#'@param dots_shape A number indicating the dot shape recognized by parameter 
+#'  'shape' in \code{geom_point()}.
+#'@param coastlines_width A number indicating the width of the coastlines.
+#'@param fileout A character string of the path to save the plot. If not 
+#'  specified (default), a graphic device will pop up. The extension should be  
+#'  accepted by \code{ggsave()}.
+#'@param width A number of the plot width, in the units specified in parameter
+#'  'size_units'. The default is 8.
+#'@param height A number of the plot height, in the units specified in parameter 
+#'  'size_units'. The default is 4.
+#'@param size_units A character string of the units of the size of the device
+#'  (file or window) to plot in. The default is 'in' (inches). See ?ggsave and
+#'  ?Devices for details of the corresponding device.
+#'@param res Resolution of the device (file or window) to plot in. The default 
+#'  value is 300. See ?ggsave 'dpi' and ?Devices for details of the 
+#'  corresponding device.
+#'
+#'@return A map plot with speficied projection, either in pop-up window or a 
+#'  saved file.
+#'
+#'@examples
+#'data <- array(rep(seq(-10, 10, length.out = 181), 360) + rnorm(360),
+#'              dim = c(lat = 181, lon = 360))
+#'dots <- data
+#'dots[which(dots < 4 & dots > -4)] <- 0
+#'dots[which(dots != 0)] <- 1
+#'  \dontrun{
+#'VizRobinson(data, lon = 0:359, lat = -90:90, dots = dots,
+#'            brks = seq(-10, 10, length.out = 11),
+#'            toptitle = 'synthetic example', vertical = FALSE,
+#'            caption = 'Robinson Projection',
+#'            bar_extra_margin = c(0, 1, 0, 1), width = 8, height = 6)
+#'VizRobinson(data, lon = 0:359, lat = -90:90, mask = dots, drawleg = 'ggplot2',
+#'           target_proj = "+proj=moll", brks = seq(-10, 10, length.out = 11),
+#'           color_fun = ClimPalette("purpleorange"), colNA = 'green',
+#'           toptitle = 'synthetic example', caption = 'Mollweide Projection',
+#'           width = 8, height = 6)
+#'  }
+#'@import sf ggplot2 rnaturalearth cowplot utils
+#'@importFrom dplyr mutate group_by summarise
+#'@importFrom ClimProjDiags Subset
+#' @importFrom rlang .data
+#'@export
+VizRobinson <- function(data, lon, lat, lon_dim = NULL, lat_dim = NULL,
+                        target_proj = NULL, drawleg = 'bar', style = 'point', 
+                        dots = NULL, mask = NULL, brks = NULL, cols = NULL, bar_limits = NULL,
+                        triangle_ends = NULL, col_inf = NULL, col_sup = NULL, colNA = NULL, 
+                        color_fun = ClimPalette(), bar_extra_margin = rep(0, 4), vertical = TRUE,
+                        toptitle = NULL, caption = NULL, units = NULL, crop_coastlines = NULL,
+                        point_size = "auto", title_size = 16, dots_size = 0.5,
+                        dots_shape = 47, coastlines_width = 0.3,
+                        fileout = NULL, width = 8, height = 4, size_units = "in", 
+                        res = 300) {
+
+  # Sanity check
+  # data
+  data <- drop(data)
+  if (length(dim(data)) != 2) {
+    stop("Parameter 'data' must have two dimensions.")
+  }
+  dims <- dim(data)
+  # lon, lon_dim
+  if (is.null(lon_dim)) {
+    lon_dim <- names(dims)[names(dims) %in% .KnownLonNames()]
+    if (identical(lon_dim, character(0))) {
+      stop("Cannot find known longitude name in data dimension. Please define parameter 'lon_dim'.")
+    }
+  }
+  if (is.unsorted(lon)) {
+    warning("Parameter 'lon' should be sorted to guarantee the correct result.")
+  }
+  # lat, lat_dim
+  if (is.null(lat_dim)) {
+    lat_dim <- names(dims)[names(dims) %in% .KnownLatNames()]
+    if (identical(lat_dim, character(0))) {
+      stop("Cannot find known latitude name in data dimension. Please define parameter 'lat_dim'.")
+    }
+  }
+  if (!all(names(dims) %in% c(lat_dim, lon_dim))) {
+    stop("Dimensions names in paramter 'data' should match 'lat_dim' and 'lon_dim.")
+  }
+  if (length(lon) != dims[lon_dim]) {
+    stop("Length of parameter 'lon' should match longitude dimension in 'data'.")
+  }
+  if (length(lat) != dims[lat_dim]) {
+    stop("Length of parameter 'lat' should match latitude dimension in 'data'.")
+  }
+
+  # Reorder data
+  data <- aperm(data, match(names(dim(data)), c(lon_dim, lat_dim)))
+
+  # Make lat always from 90 to -90
+  sort_lat <- FALSE
+  if (!is.unsorted(lat)) {
+    lat <- rev(lat)
+    data <- ClimProjDiags::Subset(data, along = lat_dim, indices = seq(length(lat), 1, -1))
+    sort_lat <- TRUE
+  }
+
+  # original_proj: it can only be regular grid now
+  original_proj <- st_crs(4326)
+  # target_proj
+  if (is.null(target_proj)) {
+    if (packageVersion("sf") < "1.0.10") {
+      target_proj <- 54030
+    } else {
+      target_proj <- "ESRI:54030"
+    }
+  }
+  target_proj_tmp <- st_crs(target_proj)
+  if (is.na(target_proj_tmp)) {
+    warning(paste0("Try ESRI code: ESRI:", target_proj))
+    target_proj <- st_crs(paste0("ESRI:", target_proj))
+  } else {
+    target_proj <- target_proj_tmp
+  }
+
+  # drawleg
+  if (!drawleg %in% c('bar', 'ggplot2', FALSE)) {
+    stop("Parameter 'drawleg' must be FALSE, 'ggplot2' or 'bar'.")
+  }
+  # style
+  if (!style %in% c('point', 'polygon') || length(style) != 1) {
+    stop("Parameter 'style' must be 'point' or 'polygon'.")
+  }
+  if (style == 'polygon') {
+    # polygon is slow for global map (and may be wrong) Confirm if users want to proceed
+    if ((abs(diff(range(lon))) > 350 & abs(diff(range(lat))) > 175) |
+        (prod(dim(data)) >= (180 * 360))) {
+      if (!isTRUE(utils::askYesNo("The region seems to be global and style 'polygon' is chosen. It may be time- and memory-consuming to plot the map. Are you sure that you want to continue?"))) {
+       return(invisible())
+      }
+    }
+  }
+  # dots
+  if (!is.null(dots)) {
+    dots <- drop(dots)
+    if (!is.array(dots) || any(!names(dim(dots)) %in% c(lon_dim, lat_dim))) {
+      stop("Parameter 'dots' must have two dimensions named as longitude and latitude dimensions in 'data'.")
+    } else {
+      dots <- aperm(dots, match(names(dim(dots)), c(lon_dim, lat_dim)))
+    }
+    if (!identical(dim(dots), dim(data))) {
+      stop("Parameter 'dots' must have the same dimensions as 'data'.")
+    } else if (is.numeric(dots)) {
+      if (all(dots %in% c(0, 1))) {
+        dots <- array(as.logical(dots), dim = dim(dots))
+      } else {
+        stop("Parameter 'dots' must have only TRUE/FALSE or 0/1.")
+      }
+    } else if (is.logical(dots)) {
+      if (!all(dots %in% c(T, F))) {
+        stop("Parameter 'dots' must have only TRUE/FALSE or 0/1.")
+      }
+    } else {
+      stop("Parameter 'dots' must be a logical or numerical array.")
+    }
+  }
+  # mask
+  if (!is.null(mask)) {
+    mask <- drop(mask)
+    if (!is.array(mask) || any(!names(dim(mask)) %in% c(lon_dim, lat_dim))) {
+      stop("Parameter 'mask' must have two dimensions named as longitude and latitude dimensions in 'data'.")
+    } else {
+      mask <- aperm(mask, match(names(dim(mask)), c(lon_dim, lat_dim)))
+    }
+    if (!identical(dim(mask), dim(data))) {
+      stop("Parameter 'mask' must have the same dimensions as 'data'.")
+    } else if (is.numeric(mask)) {
+      mask[which(is.na(mask))] <- 0
+      if (all(mask %in% c(0, 1))) {
+        mask <- array(as.logical(mask), dim = dim(mask))
+      } else {
+        stop("Parameter 'mask' must have only TRUE/FALSE or 0/1.")
+      }
+    } else if (is.logical(mask)) {
+      mask[which(is.na(mask))] <- F
+      if (!all(mask %in% c(T, F))) {
+        stop("Parameter 'mask' must have only TRUE/FALSE or 0/1.")
+      }
+    } else {
+      stop("Parameter 'mask' must be a logical or numerical array.")
+    }
+  }
+
+  tmp <- .create_var_limits(data = data, brks = brks, 
+                            bar_limits = bar_limits, drawleg = drawleg)
+  var_limits <- tmp$var_limits
+  drawleg <- tmp$drawleg
+
+  # Color bar
+  ## Check: brks, cols, bar_limits, color_fun, bar_extra_margin, units 
+  ## Build: brks, cols, bar_limits, col_inf, col_sup
+  colorbar <- ColorBarContinuous(brks = brks, cols = cols, vertical = vertical, subsampleg = NULL,
+                       bar_limits = bar_limits, var_limits = var_limits, triangle_ends = triangle_ends, 
+                       col_inf = col_inf, col_sup = col_sup, color_fun = color_fun, 
+                       plot = FALSE, draw_ticks = TRUE, draw_separators = FALSE,
+                       triangle_ends_scale = 1, extra_labels = NULL,
+                       title = units, title_scale = 1, # units_scale
+                       label_scale = 1, tick_scale = 1, #bar_tick_scale
+                       extra_margin = bar_extra_margin, label_digits = 4)
+  brks <- colorbar$brks
+  cols <- colorbar$cols
+  col_inf <- colorbar$col_inf
+  col_sup <- colorbar$col_sup
+  bar_limits <- c(head(brks, 1), tail(brks, 1))
+  # colNA
+  if (is.null(colNA)) {
+    if ('na_color' %in% names(attributes(cols))) {
+      colNA <- attr(cols, 'na_color')
+      if (!.IsColor(colNA)) {
+        stop("The 'na_color' provided as attribute of the colour vector must be a valid colour identifier.")
+      }
+    } else {
+      colNA <- 'pink'
+    }
+  } else if (!.IsColor(colNA)) {
+    stop("Parameter 'colNA' must be a valid colour identifier.")
+  }
+  # toptitle
+  if (!is.null(toptitle) && !is.character(toptitle)) {
+    stop("Parameter 'toptitle' must be a character string.")
+  }
+  # caption
+  if (!is.null(caption) && !is.character(caption)) {
+    stop("Parameter 'caption' must be a character string.")
+  }
+  # crop_coastlines
+  if (!is.null(crop_coastlines)) {
+    # if crop_coastlines doesn't have name, [lonmin, lonmax, latmin, latmax]
+    if (is.null(names(crop_coastlines))) {
+      names(crop_coastlines) <- c("lonmin", "lonmax", "latmin", "latmax")
+    } else if (!identical(sort(names(crop_coastlines)), sort(c("latmax", "latmin", "lonmax", "lonmin")))) {
+      stop("Parameter 'crop_coastlines' needs to have names 'latmax', 'latmin', 'lonmax', 'lonmin'.")
+    }
+  }
+
+  # point_size
+  if (point_size == 'auto') {
+    # 360x181 with default setting, 0.05
+    point_size <- round(0.05 * (360 * 181) / (length(lon) * length(lat)), 2)
+  } else if (!is.numeric(point_size)) {
+    stop("Parameter 'point_size' must be a number.")
+  }
+  #
+
+#=================================================================
+
+  # Adapt ColorBarContinuous parameters to ggplot plot
+  # If drawleg is FALSE, still tune with bar legend way
+  if (isFALSE(drawleg) || drawleg == 'bar') {
+    # the colorbar triangle color. If it is NULL (no triangle plotted), use colNA
+    col_inf_image <- ifelse(is.null(col_inf), colNA, col_inf)
+    col_sup_image <- ifelse(is.null(col_sup), colNA, col_sup)
+    cols_ggplot <- c(col_inf_image, cols, col_sup_image)
+
+    # Add triangles to brks
+    brks_ggplot <- brks
+    if (var_limits[2] > tail(brks, 1)) {
+      brks_ggplot <- c(brks_ggplot, max(data, na.rm = T))
+    } else {
+      brks_ggplot <- c(brks_ggplot, tail(brks, 1) + diff(tail(brks, 2)))
+    }
+    if (var_limits[1] < brks[1]) {
+      brks_ggplot <- c(min(data, na.rm = T), brks_ggplot)
+    } else {
+      brks_ggplot <- c(brks[1] - diff(brks[1:2]), brks_ggplot)
+    }
+
+  } else {  # ggplot2 legend
+    brks_ggplot <- brks
+    cols_ggplot <- cols
+  }
+
+  # Build data dataframe
+  lonlat_df <- data.frame(lon = rep(as.vector(lon), length(lat)),
+                          lat = sort(rep(as.vector(lat), length(lon)), decreasing = TRUE))
+  data_df <- lonlat_df %>%
+    dplyr::mutate(dat = as.vector(data))
+
+  lonlat_df_ori <- NULL
+
+  # Remove the points where mask = FALSE
+  if (!is.null(mask)) {
+    # Save original lonlat_df to plot with expected region
+    lonlat_df_ori <- st_as_sf(lonlat_df, coords = c("lon", "lat"), crs = original_proj)
+    lonlat_df_ori <- st_transform(lonlat_df_ori, crs = target_proj)
+    lonlat_df_ori <- as.data.frame(st_coordinates(lonlat_df_ori))
+    names(lonlat_df_ori) <- c('long', 'lat')
+
+    if (sort_lat) {
+      mask <- ClimProjDiags::Subset(mask, along = lat_dim, indices = seq(length(lat), 1, -1))
+    }
+    mask_df <- data.frame(lon = rep(as.vector(lon), length(lat)),
+                          lat = sort(rep(as.vector(lat), length(lon)), decreasing = TRUE),
+                          mask = as.vector(mask))
+    data_df <- data_df[mask_df$mask == TRUE, ]
+    lonlat_df <- data_df[, 1:2]
+  }
+
+  #NOTE: if target_proj = "ESRI:54030", Nord3v2 has different behavior from hub and ws!!
+  data_df <- st_as_sf(data_df, coords = c("lon", "lat"), crs = original_proj)
+  data_df <- st_transform(data_df, crs = target_proj)
+  data_df <- data_df %>%
+    dplyr::mutate(long = st_coordinates(data_df)[, 1],
+                  lat = st_coordinates(data_df)[, 2])
+
+  # Re-project dots
+  if (!is.null(dots)) {
+    if (sort_lat) {
+      dots <- ClimProjDiags::Subset(dots, along = lat_dim, indices = seq(length(lat), 1, -1))
+    }
+    dots_df <- data.frame(lon = rep(as.vector(lon), length(lat)),
+                          lat = sort(rep(as.vector(lat), length(lon)), decreasing = TRUE),
+                          dot = as.vector(dots))
+
+    dots_df <- st_as_sf(dots_df, coords = c("lon", "lat"), crs = original_proj)
+    dots_df <- st_transform(dots_df, crs = target_proj)
+    dots_df <- dots_df %>%
+      dplyr::mutate(long = st_coordinates(dots_df)[, 1],
+                    lat = st_coordinates(dots_df)[, 2])
+    dots_df <- dplyr::filter(dots_df, .data$dot == FALSE)
+  }
+ 
+  # coastlines
+  coastlines <- rnaturalearth::ne_coastline(scale = "medium", returnclass = "sf")
+  ## crop the coastlines to the desired range
+  if (!is.null(crop_coastlines)) {
+    suppressWarnings({
+    coastlines <- st_crop(coastlines, 
+                          xmin = as.numeric(crop_coastlines['lonmin']), 
+                          xmax = as.numeric(crop_coastlines['lonmax']), 
+                          ymin = as.numeric(crop_coastlines['latmin']), 
+                          ymax = as.numeric(crop_coastlines['latmax']))
+    })
+  }
+  coastlines <- st_transform(coastlines, crs = target_proj) 
+
+  if (style == 'polygon') {
+    # Calculate polygon points from regular lat/lon
+    #NOTE: The original grid must be regular grid with same space
+    d_lon <- abs(lon[2] - lon[1]) / 2
+    d_lat <- abs(lat[2] - lat[1]) / 2
+    lon_poly <- lat_poly <- rep(NA, 4 * dim(lonlat_df)[1])
+    for (ii in 1:dim(lonlat_df)[1]) {
+      lon_poly[(ii*4-3):(ii*4)] <- c(lonlat_df$lon[ii] - d_lon, lonlat_df$lon[ii] + d_lon,
+                                     lonlat_df$lon[ii] + d_lon, lonlat_df$lon[ii] - d_lon)
+      lat_poly[(ii*4-3):(ii*4)] <- c(lonlat_df$lat[ii] - d_lat, lonlat_df$lat[ii] - d_lat,
+                                     lonlat_df$lat[ii] + d_lat, lonlat_df$lat[ii] + d_lat)
+    }
+    #  # To prevent out-of-global lon
+    #  lon_poly[which(lon_poly >= 180)] <- 179.9
+    #  lon_poly[which(lon_poly < -180)] <- -180
+    # To prevent out-of-global lat
+    lat_poly[which(lat_poly > 90)] <- 90
+    lat_poly[which(lat_poly < -90)] <- -90
+
+    lonlat_df <- data.frame(lon = lon_poly, lat = lat_poly)
+    # Transfer lon/lat to projection
+    proj_lonlat <- st_as_sf(lonlat_df, coords = c("lon", "lat"), crs = original_proj)
+    #NOTE: if target_proj = "ESRI:54030", on Nord3v2, st_transform has lon and lat swapped!
+    proj_lonlat <- st_transform(proj_lonlat, crs = target_proj)
+    lonlat_df_proj <- st_coordinates(proj_lonlat)
+
+    # Use id to create groups for each polygon
+    ids <- factor(paste0("id_", 1:dim(data_df)[1]))
+    values <- data.frame(id = ids,
+                         value = data_df$dat)
+    positions <- data.frame(id = rep(ids, each = 4),
+                            x = lonlat_df_proj[, 1],
+                            y = lonlat_df_proj[, 2])
+    datapoly <- merge(values, positions, by = "id")
+    datapoly <- st_as_sf(datapoly, coords = c("x", "y"), crs = target_proj)
+    datapoly <- datapoly %>%
+                dplyr::group_by(.data$id) %>%
+                dplyr::summarise() %>%  #NOTE: VERY SLOW if plot global
+                dplyr::mutate(value = values[order(values$id), ]$value) %>%
+                st_cast("POLYGON") %>%
+                st_convex_hull()  # maintain outer polygen (no bowtie shape)
+  }
+
+  # Plots
+  if (style == 'polygon') {
+    res_p <- ggplot(data = data_df) +  #NOTE: must be data_df?
+             geom_sf(data = datapoly,
+                     aes(col = cut(.data$value, breaks = brks_ggplot, include.lowest = T),
+                         fill = cut(.data$value, breaks = brks_ggplot, include.lowest = T)))
+  } else if (style == 'point') {
+    res_p <- ggplot(data = data_df) +
+             geom_point(aes(x = .data$long, y = .data$lat,
+                            col = cut(.data$dat, breaks = brks_ggplot, include.lowest = T)),
+                        #NOTE: These two lines make point size vary with lat
+                        #size = point_size /  (data_df$lat / min(data_df$lat))) +
+                        #size = (sort(rep(as.vector(lat), length(lon))) / max(lat)) * point_size) +
+                        size = point_size)
+  }
+
+  if (is.null(lonlat_df_ori)) {
+    coord_sf_lim <- c(range(data_df$long), range(data_df$lat))
+  } else {
+    coord_sf_lim <- c(range(lonlat_df_ori$long), range(lonlat_df_ori$lat))
+  }
+  res_p <- res_p +
+           geom_sf(data = coastlines, colour ='black', size = coastlines_width) +
+           # Remove background grid and lat/lon label; add white background
+           theme_void() + theme(plot.background = element_rect(fill = 'white', colour = 'white')) +
+           # crop the projection
+           coord_sf(xlim = coord_sf_lim[1:2], ylim = coord_sf_lim[3:4],
+                    expand = TRUE, datum = target_proj)
+
+  if (!is.null(dots)) {
+    res_p <- res_p + geom_point(data = dots_df, aes(x = .data$long, y = .data$lat),
+                                shape = dots_shape, size = dots_size)
+                                #NOTE: This line makes point size vary with lat
+                                #size = dots_size / (dots_df$lat / min(dots_df$lat)))
+  }
+
+  if (identical(drawleg, 'ggplot2')) {
+    if (style == 'polygon') {
+    res_p <- res_p + scale_colour_manual(values = cols_ggplot,
+                                         aesthetics = c("colour", "fill"),
+                                         drop = FALSE, na.value = colNA) +
+             guides(fill = guide_legend(title = units, override.aes = list(size = 1)),
+                    color = "none")
+   } else if (style == 'point') {
+    res_p <- res_p + scale_colour_manual(values = cols_ggplot,
+                                         drop = FALSE, na.value = colNA) +
+             guides(colour = guide_legend(title = units, override.aes = list(size = 1)))
+   }
+
+  } else {  # bar or NULL
+    if (style == 'polygon') {
+      res_p <- res_p + scale_colour_manual(values = cols_ggplot,
+                                           aesthetics = c("colour", "fill"),
+                                           drop = FALSE, na.value = colNA)
+    } else if (style == 'point') {
+      res_p <- res_p + scale_colour_manual(values = cols_ggplot,
+                                           drop = FALSE, na.value = colNA)
+    }
+    # Remove ggplot legend
+    res_p <- res_p + theme(legend.position = "none", plot.margin = margin(0.5, 0, 0, 0, 'cm'))
+  }
+
+  if (!is.null(toptitle)) {
+    res_p <- res_p + ggtitle(toptitle) +
+      theme(plot.title = element_text(size = title_size, hjust = 0.5, vjust = 3))
+  }
+  if (!is.null(caption)) {
+    res_p <- res_p + labs(caption = caption) +
+      theme(plot.caption = element_text(hjust = 0, vjust = 1, margin = margin(0, 0, 0, 0, 'cm')))
+  }
+
+  # bar legend fun to put in cowplot::plot_grid
+  if (identical(drawleg, 'bar')) {
+    fun_legend <- function() {
+      if (vertical) {
+        par(mar = c(7.1, 2.2, 7.1, 3.1), mgp = c(3, 1, 0))
+      } else {
+        par(mar = c(1.1, 1.2, 0.1, 1.1), mgp = c(3, 1, 0))
+      }
+      ColorBarContinuous(brks = brks, cols = cols, vertical = vertical, subsampleg = NULL,
+               bar_limits = bar_limits, var_limits = var_limits, triangle_ends = triangle_ends,
+               col_inf = col_inf, col_sup = col_sup, color_fun = color_fun,
+               plot = TRUE, draw_ticks = TRUE, draw_separators = FALSE,
+               triangle_ends_scale = 1, extra_labels = NULL,
+               title = units, title_scale = 1, # units_scale
+               label_scale = 1, tick_scale = 1, #bar_tick_scale
+               extra_margin = bar_extra_margin, label_digits = 4)
+    }
+    if (vertical) {
+      res_p <- cowplot::plot_grid(res_p, fun_legend, rel_widths = c(6, 1))
+    } else {
+      res_p <- cowplot::plot_grid(res_p, fun_legend, rel_heights = c(5, 1), ncol = 1)
+    }
+    res_p <- res_p +  theme(plot.background = element_rect(fill = "white", colour = "white"))
+  }
+
+  if (!is.null(fileout)) {
+    ext <- regmatches(fileout, regexpr("[a-zA-Z0-9]*$", fileout))
+    ggsave(fileout, res_p, width = width, height = height, dpi = res, units = size_units, 
+           device = ext)
+  } else { # pop-up window
+    dev.new(units = size_units, res = res, width = width, height = height)
+    res_p 
+  } 
+
+}
+
diff --git a/modules/Visualization/R/tmp/zzz.R b/modules/Visualization/R/tmp/zzz.R
index f2871dfd..cf2eb1c7 100644
--- a/modules/Visualization/R/tmp/zzz.R
+++ b/modules/Visualization/R/tmp/zzz.R
@@ -1,86 +1,111 @@
-# Function to permute arrays of non-atomic elements (e.g. POSIXct)
-.aperm2 <- function(x, new_order) {
-  old_dims <- dim(x)
-  attr_bk <- attributes(x)
-  if ('dim' %in% names(attr_bk)) {
-    attr_bk[['dim']] <- NULL
-  }
-  if (is.numeric(x)) {
-    x <- aperm(x, new_order)
-  } else {
-    y <- array(1:length(x), dim = dim(x))
-    y <- aperm(y, new_order)
-    x <- x[as.vector(y)]
-  }
-  dim(x) <- old_dims[new_order]
-  attributes(x) <- c(attributes(x), attr_bk)
-  x
+.KnownLonNames <- function() {
+  known_lon_names <- c('lon', 'longitude', 'x', 'i', 'nav_lon')
 }
 
-# verbose-only printing function
-.printv <- function(value, verbosity = TRUE) {
-  if (verbosity) {
-    print(value)
-  }
+.KnownLatNames <- function() {
+  known_lat_names <- c('lat', 'latitude', 'y', 'j', 'nav_lat')
 }
 
-# normalize a time series
-.standardize <- function(timeseries) {
-  out <- (timeseries - mean(timeseries, na.rm = T)) / sd(timeseries, na.rm = T)
-  return(out)
+.IsColor <- function(x) {
+  res <- try(col2rgb(x), silent = TRUE)
+  return(!"try-error" %in% class(res))
 }
 
-.selbox <- function(lon, lat, xlim = NULL, ylim = NULL) {
-  if (!is.null(xlim)) {
-    # This transforms c(-20, -10) to c(340, 350) but c(-20, 10) is unchanged
-    # Bring them all to the same units in the 0:360 range
-    xlim1 <- xlim[1] %% 360
-    xlim2 <- xlim[2] %% 360
-    lonm <- lon %% 360
-    if (lonm[1] > tail(lonm, 1)) {
-      lonm <- lon
-    }
-    if (xlim1 > xlim2) {
-      # If box crosses 0
-      ilonsel <- (lonm >= xlim1) | (lonm <= xlim2)
-    } else {
-      ilonsel <- (lonm >= xlim1) & (lonm <= xlim2)
-    }
-    if (!any(ilonsel)) {
-      stop("No intersection between longitude bounds and data domain.")
-    }
-  } else {
-    ilonsel <- 1:length(lon)
-  }
-  if (!is.null(ylim)) {
-    ilatsel <- (lat >= ylim[1]) & (lat <= ylim[2])
-  } else {
-    ilatsel <- 1:length(lat)
+.FilterUserGraphicArgs <- function(excludedArgs, ...) {
+  # This function filter the extra graphical parameters passed by the user in
+  # a plot function, excluding the ones that the plot function uses by default.
+  # Each plot function has a different set of arguments that are not allowed to
+  # be modified.
+  args <- list(...)
+  userArgs <- list()
+  for (name in names(args)) {
+      if ((name != "") & !is.element(name, excludedArgs)) {
+          # If the argument has a name and it is not in the list of excluded
+          # arguments, then it is added to the list that will be used
+          userArgs[[name]] <- args[[name]]
+      } else {
+        warning(paste0("the argument '", name, "' can not be 
+        modified and the new value will be ignored"))
+      }
   }
-  return(list(ilon = ilonsel, ilat = ilatsel))
+  userArgs
 }
 
-# produce a 2d matrix of area weights
-.area.weight <- function(ics, ipsilon, root = T) {
-  field <- array(NA, dim = c(length(ics), length(ipsilon)))
-  if (root == T) {
-    for (j in 1:length(ipsilon)) {
-      field[, j] <- sqrt(cos(pi / 180 * ipsilon[j]))
-    }
-  }
+.SelectDevice <- function(fileout, width, height, units, res) {
+  # This function is used in the plot functions to check the extension of the 
+  # files where the graphics will be stored and select the right R device to 
+  # save them.
+  # If the vector of filenames ('fileout') has files with different 
+  # extensions, then it will only accept the first one, changing all the rest 
+  # of the filenames to use that extension.
+
+  # We extract the extension of the filenames: '.png', '.pdf', ...
+  ext <- regmatches(fileout, regexpr("\\.[a-zA-Z0-9]*$", fileout))
 
-  if (root == F) {
-    for (j in 1:length(ipsilon)) {
-      field[, j] <- cos(pi / 180 * ipsilon[j])
+  if (length(ext) != 0) {
+    # If there is an extension specified, select the correct device
+    ## units of width and height set to accept inches
+    if (ext[1] == ".png") {
+      saveToFile <- function(fileout) {
+        png(filename = fileout, width = width, height = height, res = res, units = units)
+      }
+    } else if (ext[1] == ".jpeg") {
+      saveToFile <- function(fileout) {
+        jpeg(filename = fileout, width = width, height = height, res = res, units = units)
+      }
+    } else if (ext[1] %in% c(".eps", ".ps")) {
+      saveToFile <- function(fileout) {
+        postscript(file = fileout, width = width, height = height)
+      }
+    } else if (ext[1] == ".pdf") {
+      saveToFile <- function(fileout) {
+        pdf(file = fileout, width = width, height = height)
+      }
+    } else if (ext[1] == ".svg") {
+      saveToFile <- function(fileout) {
+        svg(filename = fileout, width = width, height = height)
+      }
+    } else if (ext[1] == ".bmp") {
+      saveToFile <- function(fileout) {
+        bmp(filename = fileout, width = width, height = height, res = res, units = units)
+      }
+    } else if (ext[1] == ".tiff") {
+      saveToFile <- function(fileout) {
+        tiff(filename = fileout, width = width, height = height, res = res, units = units)
+      }
+    } else {
+      warning("file extension not supported, it will be used '.eps' by default.")
+      ## In case there is only one filename
+      fileout[1] <- sub("\\.[a-zA-Z0-9]*$", ".eps", fileout[1])
+      ext[1] <- ".eps"
+      saveToFile <- function(fileout) {
+        postscript(file = fileout, width = width, height = height)
+      }
     }
+    # Change filenames when necessary
+    if (any(ext != ext[1])) {
+      warning(paste0("some extensions of the filenames provided in 'fileout' are not ", ext[1],". The extensions are being converted to ", ext[1], "."))
+      fileout <- sub("\\.[a-zA-Z0-9]*$", ext[1], fileout)
+    }
+  } else {
+    # Default filenames when there is no specification
+    warning("there are no extensions specified in the filenames, default to '.eps'")
+    fileout <- paste0(fileout, ".eps")
+    saveToFile <- postscript
   }
 
-  return(field)
+  # return the correct function with the graphical device, and the correct 
+  # filenames
+  list(fun = saveToFile, files = fileout)
 }
 
 #Draws Color Bars for Categories
-#A wrapper of s2dv::ColorBar to generate multiple color bars for different 
+#A wrapper of ColorBar to generate multiple color bars for different 
+#categories, and each category has different color set.
+#Draws Color Bars for Categories
+#A wrapper of ColorBarContinuous to generate multiple color bars for different 
 #categories, and each category has different color set.
+#'@import utils
 GradientCatsColorBar <- function(nmap, brks = NULL, cols = NULL, vertical = TRUE, subsampleg = NULL,
                                  bar_limits, var_limits = NULL,
                                  triangle_ends = NULL, col_inf = NULL, col_sup = NULL, plot = TRUE,
@@ -206,13 +231,12 @@ GradientCatsColorBar <- function(nmap, brks = NULL, cols = NULL, vertical = TRUE
 
   if (plot) {
     for (k in 1:nmap) {
-#TODO: Add s2dv:: 
-      ColorBar(brks = brks[[k]], cols = cols[[k]], vertical = FALSE, subsampleg = subsampleg,
-                     bar_limits = bar_limits[[k]], #var_limits = var_limits,
-                     triangle_ends = triangle_ends, col_inf = col_inf[[k]], col_sup = col_sup[[k]], plot = TRUE,
-                     draw_separators = draw_separators,
-                     title = bar_titles[[k]], title_scale = title_scale,
-                     label_scale = label_scale, extra_margin = extra_margin)
+      ColorBarContinuous(brks = brks[[k]], cols = cols[[k]], vertical = FALSE, subsampleg = subsampleg,
+                         bar_limits = bar_limits[[k]], #var_limits = var_limits,
+                         triangle_ends = triangle_ends, col_inf = col_inf[[k]], col_sup = col_sup[[k]],  plot = TRUE,
+                         draw_separators = draw_separators,
+                         title = bar_titles[[k]], title_scale = title_scale,
+                         label_scale = label_scale, extra_margin = extra_margin)
     }
   } else {
     return(list(brks = brks, cols = cols, col_inf = col_inf, col_sup = col_sup))
@@ -220,37 +244,27 @@ GradientCatsColorBar <- function(nmap, brks = NULL, cols = NULL, vertical = TRUE
 
 }
 
-.KnownLonNames <- function() {
-  known_lon_names <- c('lon', 'lons', 'longitude', 'x', 'i', 'nav_lon')
-}
-
-.KnownLatNames <- function() {
-  known_lat_names <- c('lat', 'lats', 'latitude', 'y', 'j', 'nav_lat')
-}
-
-.KnownTimeNames <- function() {
-  known_time_names <- c('time', 'ftime', 'sdate', 'sdates', 'syear', 'sweek', 'sday', 'leadtimes')
-}
-
-.KnownForecastTimeNames <- function() {
-  known_time_names <- c('time', 'ftime', 'ltime', 'leadtimes')
-}
-
-.KnownStartDateNames <- function() {
-  known_time_names <- c('sdate', 'sdates', 'syear', 'sweek', 'sday')
-}
-
-.KnownMemberNames <- function() {
-  known_time_names <- c('memb', 'member', 'members', 'ensemble', 'ensembles')
-}
-
-.isNullOb <- function(x) is.null(x) | all(sapply(x, is.null))
-
-.rmNullObs <- function(x) {
-   x <- base::Filter(Negate(.isNullOb), x)
-   lapply(x, function(x) if (is.list(x)) .rmNullObs(x) else x)
-}
-
-# Definition of a global variable to store the warning message used in Calibration
-warning_shown <- FALSE
-
+# Decide var_limits for ColorBarContinuous()
+.create_var_limits <- function(data, brks, bar_limits, drawleg) {
+  if (!all(is.na(data))) {
+    var_limits <- c(min(data[!is.infinite(data)], na.rm = TRUE),
+                    max(data[!is.infinite(data)], na.rm = TRUE))
+  } else {
+    warning("All the data are NAs. The map will be filled with colNA.")
+    if (!is.null(brks) && length(brks) > 1) {
+      #NOTE: var_limits be like this to avoid warnings from ColorBar
+      var_limits <- c(min(brks, na.rm = TRUE) + diff(brks)[1], 
+                      max(brks, na.rm = TRUE))
+    } else if (!is.null(bar_limits)) {
+      var_limits <- c(bar_limits[1] + 0.01, bar_limits[2])
+    } else {
+      var_limits <- c(-0.5, 0.5) # random range since colorbar is not going to be plotted
+      if (!isFALSE(drawleg)) {
+        drawleg <- FALSE
+        warning("All data are NAs. Color bar won't be drawn. If you want to have ",
+                 "color bar still, define parameter 'brks' or 'bar_limits'.")
+      }
+    }
+  }
+  return(list(var_limits = var_limits, drawleg = drawleg))
+}
\ No newline at end of file
diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R
index 181c65ef..ee7270ae 100644
--- a/modules/Visualization/Visualization.R
+++ b/modules/Visualization/Visualization.R
@@ -8,10 +8,15 @@ source("modules/Visualization/R/plot_metrics.R")
 source("modules/Visualization/R/get_proj_code.R")
 source("modules/Visualization/R/get_plot_time_labels.R")
 source("modules/Visualization/R/get_variable_longname.R")
-## TODO: Remove after the next s2dv release
-source("modules/Visualization/R/tmp/PlotRobinson.R")
 source("modules/Visualization/R/plot_most_likely_terciles_map.R")
 source("modules/Visualization/R/plot_ensemble_mean.R")
+## TODO: Remove after esviz release
+source("modules/Visualization/R/tmp/VizEquiMap.R")
+source("modules/Visualization/R/tmp/VizRobinson.R")
+source("modules/Visualization/R/tmp/PlotRobinson.R")
+source("modules/Visualization/R/tmp/zzz.R")
+source("modules/Visualization/R/tmp/ColorBarContinuous.R")
+source("modules/Visualization/R/tmp/ClimPalette.R")
 
 Visualization <- function(recipe,
                           data,
-- 
GitLab


From 612f0f9f965f8791260140262484c5cf1f4adf3f Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Mon, 13 Jan 2025 14:47:12 +0100
Subject: [PATCH 06/24] Do not display variable units in metric captions

---
 modules/Visualization/R/plot_metrics.R | 1 -
 1 file changed, 1 deletion(-)

diff --git a/modules/Visualization/R/plot_metrics.R b/modules/Visualization/R/plot_metrics.R
index 4c39b6d9..c263cd10 100644
--- a/modules/Visualization/R/plot_metrics.R
+++ b/modules/Visualization/R/plot_metrics.R
@@ -449,7 +449,6 @@ plot_metrics <- function(recipe, data_cube, metrics,
               paste0("Nominal start date: ", nominal_startdate_caption, "\n",
                      forecast_time_caption, "\n",
                      "Reference: ", recipe$Analysis$Datasets$Reference, "\n",
-                     "Units: ", data_cube$attrs$Variable$metadata[[var_name]]$units, "\n",
                      significance_caption)
             
             fileout <- paste0(outfile, "_ft", forecast_time,
-- 
GitLab


From 1481acf02adc9bc772bda0c7bc5b2a48557922f9 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Mon, 13 Jan 2025 15:30:59 +0100
Subject: [PATCH 07/24] Change plot_ensemble_mean.R to plot_forecast_map.R

---
 .../R/{plot_ensemble_mean.R => plot_forecast_map.R}             | 0
 modules/Visualization/Visualization.R                           | 2 +-
 2 files changed, 1 insertion(+), 1 deletion(-)
 rename modules/Visualization/R/{plot_ensemble_mean.R => plot_forecast_map.R} (100%)

diff --git a/modules/Visualization/R/plot_ensemble_mean.R b/modules/Visualization/R/plot_forecast_map.R
similarity index 100%
rename from modules/Visualization/R/plot_ensemble_mean.R
rename to modules/Visualization/R/plot_forecast_map.R
diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R
index ee7270ae..c14c76de 100644
--- a/modules/Visualization/Visualization.R
+++ b/modules/Visualization/Visualization.R
@@ -9,7 +9,7 @@ source("modules/Visualization/R/get_proj_code.R")
 source("modules/Visualization/R/get_plot_time_labels.R")
 source("modules/Visualization/R/get_variable_longname.R")
 source("modules/Visualization/R/plot_most_likely_terciles_map.R")
-source("modules/Visualization/R/plot_ensemble_mean.R")
+source("modules/Visualization/R/plot_forecast_map.R")
 ## TODO: Remove after esviz release
 source("modules/Visualization/R/tmp/VizEquiMap.R")
 source("modules/Visualization/R/tmp/VizRobinson.R")
-- 
GitLab


From 65e1913e54cd8b5463f517fd34380c837bc05c7f Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Tue, 14 Jan 2025 10:33:38 +0100
Subject: [PATCH 08/24] Improve default values for colorbar and margin sizes

---
 modules/Visualization/R/plot_metrics.R | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/modules/Visualization/R/plot_metrics.R b/modules/Visualization/R/plot_metrics.R
index c263cd10..d52aad1d 100644
--- a/modules/Visualization/R/plot_metrics.R
+++ b/modules/Visualization/R/plot_metrics.R
@@ -192,7 +192,7 @@ plot_metrics <- function(recipe, data_cube, metrics,
         metric <- Reorder(metric, c("time", "longitude", "latitude"))
         # If the significance has been requested and the variable has it,
         # retrieve it and reorder its dimensions.
-        if (significance != FALSE) { # Both, dots, mask or TRUE	
+        if (significance != FALSE) { # Both, dots, mask or TRUE 
           significance_name <- paste0(name, "_significance")
           if ((significance_name %in% names(metrics))) {
             metric_significance <- var_metric[[significance_name]]
@@ -256,7 +256,7 @@ plot_metrics <- function(recipe, data_cube, metrics,
                            " / ", hcst_period)
             titles <- paste("Valid from", format(week_valid_ini, "%d-%m"),
                             "to", format(week_valid_end, "%d-%m"))
-	  } else {
+          } else {
             toptitle <- "Unknown"
             titles <- "Unknown"
           }
@@ -285,12 +285,14 @@ plot_metrics <- function(recipe, data_cube, metrics,
             base_args <- list(var = NULL, dots = NULL,
                               lon = longitude, lat = latitude,
                               dot_symbol = 20, dot_size = 1,
-                              title_scale = 0.6, margin_scale = c(1, 5, 5, 5),
+                              title_scale = 0.6,
                               filled.continents = F, brks = brks, cols = cols, 
                               col_inf = col_inf, col_sup = col_sup, 
                               units = units, font.main = 2,
-                              bar_label_digits = 3, bar_label_scale = 1.5,
-                              axes_label_scale = 1, width = 8, height = 5)
+                              bar_label_digits = 3, bar_label_scale = 1,
+                              axes_label_scale = 1, width = 8, height = 5,
+                              bar_extra_margin = rep(0.9, 4),
+                              margin_scale = c(3, 5, 5, 5)) 
           } else {
             fun <- VizRobinson
             common_projections <- c("robinson", "stereographic", "lambert_europe")
@@ -439,7 +441,7 @@ plot_metrics <- function(recipe, data_cube, metrics,
                   sign_file_label <- '_mask'
                 }
               }
-              significance_caption <- "alpha = 0.05"
+              significance_caption <- "\n alpha = 0.05"
             } else {
               significance_caption <- NULL
               sign_file_label <- NULL
@@ -448,7 +450,7 @@ plot_metrics <- function(recipe, data_cube, metrics,
             base_args[['caption']] <- 
               paste0("Nominal start date: ", nominal_startdate_caption, "\n",
                      forecast_time_caption, "\n",
-                     "Reference: ", recipe$Analysis$Datasets$Reference, "\n",
+                     "Reference: ", recipe$Analysis$Datasets$Reference,
                      significance_caption)
             
             fileout <- paste0(outfile, "_ft", forecast_time,
-- 
GitLab


From 616706d705110b532dc8770a8e183c17ed0c95a3 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Tue, 14 Jan 2025 13:25:29 +0100
Subject: [PATCH 09/24] Bugfixes for plotting variables

---
 modules/Visualization/R/plot_forecast_map.R | 6 ++++++
 modules/Visualization/Visualization.R       | 4 ++--
 2 files changed, 8 insertions(+), 2 deletions(-)

diff --git a/modules/Visualization/R/plot_forecast_map.R b/modules/Visualization/R/plot_forecast_map.R
index de2de50b..ffc683bb 100644
--- a/modules/Visualization/R/plot_forecast_map.R
+++ b/modules/Visualization/R/plot_forecast_map.R
@@ -53,6 +53,8 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
   if (method == 'mean') {
     forecast_product <- s2dv::MeanDims(fcst$data, 'ensemble')
   } else if (method == 'iqr') {
+    ## Change to uppercase to display correctly in plot title?
+    method <- "IQR"
     forecast_product <- Apply(fcst$data, target_dim = 'ensemble',
                               fun = function(x) {
                               IQR(x)})$output1
@@ -134,6 +136,8 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
         dim_mask <- dim(var_mask)
         var_mask <- as.numeric(var_mask <= 0)
         dim(var_mask) <- dim_mask
+      } else {
+        var_mask <- NULL
       }
       # Dots
       if (!is.null(dots)) {
@@ -145,6 +149,8 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
         dim_dots <- dim(var_dots)
         var_dots <- as.numeric(var_dots <= 0)
         dim(var_dots) <- dim_dots
+      } else {
+        var_dots <- NULL
       }
 
       # Plot title labels look different depending on horizon and aggregation
diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R
index c14c76de..3b0d5f22 100644
--- a/modules/Visualization/Visualization.R
+++ b/modules/Visualization/Visualization.R
@@ -137,7 +137,7 @@ Visualization <- function(recipe,
   }
 
   # Plot forecast ensemble mean
-  if ("forecast_ensemble_mean" %in% plots) {
+  if ("forecast_map" %in% plots) {
     if (!is.null(data$fcst)) {
       if (is.null(recipe$Analysis$Workflow$Visualization$mask_ens)) {
         recipe$Analysis$Workflow$Visualization$mask_ens <- FALSE
@@ -145,7 +145,7 @@ Visualization <- function(recipe,
       if (is.null(recipe$Analysis$Workflow$Visualization$dots)) {
         recipe$Analysis$Workflow$Visualization$dots <- FALSE
       }
-      if (is.null(recipe$Analysis$Workflow$Visualization$dots)) {
+      if (is.null(recipe$Analysis$Workflow$Visualization$forecast_method)) {
         recipe$Analysis$Workflow$Visualization$forecast_method <- "median"
       } 
       # Loop over the methods requested
-- 
GitLab


From 846c33ce5bd48b096a0adbf377d49f02529e88c7 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Tue, 14 Jan 2025 13:25:50 +0100
Subject: [PATCH 10/24] Add check to replace 'forecast_ensemble_mean' with
 'forecast_map'

---
 tools/check_recipe.R | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/tools/check_recipe.R b/tools/check_recipe.R
index e6ab4b2d..3be5e081 100644
--- a/tools/check_recipe.R
+++ b/tools/check_recipe.R
@@ -677,8 +677,13 @@ check_recipe <- function(recipe) {
   
   # Visualization
   if ("Visualization" %in% names(recipe$Analysis$Workflow)) {
-    PLOT_OPTIONS <- c("skill_metrics", "forecast_ensemble_mean",
+    PLOT_OPTIONS <- c("skill_metrics", "forecast_map",
                       "most_likely_terciles", "statistics")
+    ## TODO: Remove this?
+    # Replace 'forecast_ensemble_mean' with new 'forecast_map' option:
+    recipe$Analysis$Workflow$Visualization$plots <- gsub(pattern = "forecast_ensemble_mean",
+                                                         replacement = "forecast_map",
+                                                         x = recipe$Analysis$Workflow$Visualization$plots)
     # Separate plots parameter and check if all elements are in PLOT_OPTIONS
     if (is.null(recipe$Analysis$Workflow$Visualization$plots)) {
       error(recipe$Run$logger,
-- 
GitLab


From e193b821bc2c5141a4eb709c802788eb77c331e1 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Tue, 14 Jan 2025 13:25:59 +0100
Subject: [PATCH 11/24] Update template

---
 recipe_template.yml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/recipe_template.yml b/recipe_template.yml
index b9cd5300..c1287f11 100644
--- a/recipe_template.yml
+++ b/recipe_template.yml
@@ -138,7 +138,8 @@ Analysis:
       # enclosed within brackets. For now, they are INDEPENDENT from skill metrics. (Optional)
       save: 'percentiles_only' # Options: 'all', 'none', 'bins_only', 'percentiles_only' (Mandatory, str)
     Visualization:
-      plots: skill_metrics, most_likely_terciles, forecast_ensemble_mean # Types of plots to generate (Optional, str)
+      plots: skill_metrics, most_likely_terciles, forecast_map # Types of plots to generate (Optional, str)
+      forecast_method: mean, median, iqr # One or more methods for the 'forecast_map' plot. Available methods are: 'mean', 'median', 'IQR'. Default is 'median'. (Optional, str)
       multi_panel: yes # Multi-panel plot or single-panel plots. Default is 'no/false'. (Optional, bool)
       projection: 'cylindrical_equidistant' # Options: 'cylindrical_equidistant', 'robinson', 'lambert_europe'. Default is cylindrical equidistant. (Optional, str)
       significance: 'dots' # Type of mask for statistical significance. Options are 'dots', and yes/no. 'dots'. 'mask' and 'both' options are not available for projections other than cylindrical_equidistant.
-- 
GitLab


From 448ce8dae7af0c2d4dd12cad3fd8513573fbeeef Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Tue, 14 Jan 2025 13:26:15 +0100
Subject: [PATCH 12/24] Update unit tests

---
 tests/testthat/test-decadal_monthly_1.R        | 2 +-
 tests/testthat/test-decadal_monthly_2.R        | 2 +-
 tests/testthat/test-seasonal_monthly.R         | 2 +-
 tests/testthat/test-seasonal_monthly_timeagg.R | 2 +-
 tests/testthat/test-subseasonal_weekly.R       | 8 ++++----
 5 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/tests/testthat/test-decadal_monthly_1.R b/tests/testthat/test-decadal_monthly_1.R
index 2b8c93c9..41084e8b 100644
--- a/tests/testthat/test-decadal_monthly_1.R
+++ b/tests/testthat/test-decadal_monthly_1.R
@@ -273,7 +273,7 @@ test_that("5. Visualization", {
 plots <- paste0(recipe$Run$output_dir, "/plots/")
 expect_equal(
 all(basename(list.files(plots, recursive = T)) %in%
-c("forecast_ensemble_mean-2021.pdf", "forecast_most_likely_tercile-2021.pdf",
+c("forecast_ensemble_median-2021.pdf", "forecast_most_likely_tercile-2021.pdf",
   "rpss.pdf")),
 TRUE
 )
diff --git a/tests/testthat/test-decadal_monthly_2.R b/tests/testthat/test-decadal_monthly_2.R
index 4139e5a2..44de32b8 100644
--- a/tests/testthat/test-decadal_monthly_2.R
+++ b/tests/testthat/test-decadal_monthly_2.R
@@ -262,7 +262,7 @@ test_that("5. Visualization", {
 plots <- paste0(recipe$Run$output_dir, "/plots/")
 expect_equal(
 all(basename(list.files(plots, recursive = T)) %in%
-c("bss10_specs.pdf", "enscorr_specs.pdf", "forecast_ensemble_mean-2020.pdf", "forecast_ensemble_mean-2021.pdf", "forecast_most_likely_tercile-2020.pdf", "forecast_most_likely_tercile-2021.pdf", "frps_specs.pdf", "frps.pdf", "rpss_specs.pdf")
+c("bss10_specs.pdf", "enscorr_specs.pdf", "forecast_ensemble_median-2020.pdf", "forecast_ensemble_median-2021.pdf", "forecast_most_likely_tercile-2020.pdf", "forecast_most_likely_tercile-2021.pdf", "frps_specs.pdf", "frps.pdf", "rpss_specs.pdf")
 ),
 TRUE
 )
diff --git a/tests/testthat/test-seasonal_monthly.R b/tests/testthat/test-seasonal_monthly.R
index 198c947a..df70526c 100644
--- a/tests/testthat/test-seasonal_monthly.R
+++ b/tests/testthat/test-seasonal_monthly.R
@@ -235,7 +235,7 @@ plots <- paste0(recipe$Run$output_dir, "/plots/")
 expect_equal(
 all(basename(list.files(plots, recursive = T)) %in%
 c("crpss-november.pdf", "enscorr_specs-november.pdf", "enscorr-november.pdf",
-  "forecast_ensemble_mean-20201101.pdf", "forecast_most_likely_tercile-20201101.pdf",
+  "forecast_ensemble_median-20201101.pdf", "forecast_most_likely_tercile-20201101.pdf",
   "rpss-november.pdf")),
 TRUE
 )
diff --git a/tests/testthat/test-seasonal_monthly_timeagg.R b/tests/testthat/test-seasonal_monthly_timeagg.R
index ffa389c4..c05ed2e7 100644
--- a/tests/testthat/test-seasonal_monthly_timeagg.R
+++ b/tests/testthat/test-seasonal_monthly_timeagg.R
@@ -193,7 +193,7 @@ plots <- paste0(recipe$Run$output_dir, "/plots/")
 expect_equal(
 all(basename(list.files(plots, recursive = T)) %in%
 c("crpss-november.pdf", "enscorr_specs-november.pdf", "enscorr-november.pdf",
-  "forecast_ensemble_mean-20201101.pdf", "forecast_most_likely_tercile-20201101.pdf",
+  "forecast_ensemble_median-20201101.pdf", "forecast_most_likely_tercile-20201101.pdf",
   "rpss-november.pdf")),
 TRUE
 )
diff --git a/tests/testthat/test-subseasonal_weekly.R b/tests/testthat/test-subseasonal_weekly.R
index 72a1cb26..6e565e76 100644
--- a/tests/testthat/test-subseasonal_weekly.R
+++ b/tests/testthat/test-subseasonal_weekly.R
@@ -165,10 +165,10 @@ Visualization(recipe = recipe, data = data)
 plots <- paste0(recipe$Run$output_dir, "/plots/")
 expect_equal(
 all(basename(list.files(plots, recursive = T)) %in%
-c("forecast_ensemble_mean-20240104_ft01.png",
-  "forecast_ensemble_mean-20240104_ft02.png",
-  "forecast_ensemble_mean-20240104_ft03.png",
-  "forecast_ensemble_mean-20240104_ft04.png")),
+c("forecast_ensemble_median-20240104_ft01.png",
+  "forecast_ensemble_median-20240104_ft02.png",
+  "forecast_ensemble_median-20240104_ft03.png",
+  "forecast_ensemble_median-20240104_ft04.png")),
 TRUE
 )
 expect_equal(
-- 
GitLab


From 05260b3b6ea0078aec91c0d9296293eae0094676 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Fri, 17 Jan 2025 09:34:12 +0100
Subject: [PATCH 13/24] Use .drop_dims() before plotting probability array

---
 .../R/plot_most_likely_terciles_map.R         | 35 +++++++++++--------
 1 file changed, 20 insertions(+), 15 deletions(-)

diff --git a/modules/Visualization/R/plot_most_likely_terciles_map.R b/modules/Visualization/R/plot_most_likely_terciles_map.R
index 67236328..963a195f 100644
--- a/modules/Visualization/R/plot_most_likely_terciles_map.R
+++ b/modules/Visualization/R/plot_most_likely_terciles_map.R
@@ -49,22 +49,27 @@ plot_most_likely_terciles <- function(recipe,
     init_date <- 1
   }
   # Retrieve and rearrange probability bins for the forecast
-  if (is.null(probabilities$probs_fcst$prob_b33) ||
-      is.null(probabilities$probs_fcst$prob_33_to_66) ||
-      is.null(probabilities$probs_fcst$prob_a66)) {
-    stop("The forecast tercile probability bins are not present inside ",
-         "'probabilities', the most likely tercile map cannot be plotted.")
+  if (is.array(probabilities$probs_fcst)) {
+    ## TODO: Adjust array dimensions if needed. What if there is more than one 
+    ## set of categories?
+    probs_fcst <- probabilities$probs_fcst
+    probs_fcst <- .drop_dims(probs_fcst)
+  } else {
+    if (is.null(probabilities$probs_fcst$prob_b33) ||
+        is.null(probabilities$probs_fcst$prob_33_to_66) ||
+        is.null(probabilities$probs_fcst$prob_a66)) {
+      stop("The forecast tercile probability bins are not present inside ",
+           "'probabilities', the most likely tercile map cannot be plotted.")
+    }
+    
+    probs_fcst <- abind(probabilities$probs_fcst$prob_b33,
+                        probabilities$probs_fcst$prob_33_to_66,
+                        probabilities$probs_fcst$prob_a66,
+                        along = 0)
+    names(dim(probs_fcst)) <- c("bin",
+                                names(dim(probabilities$probs_fcst$prob_b33)))
   }
-  
-  probs_fcst <- abind(probabilities$probs_fcst$prob_b33,
-                      probabilities$probs_fcst$prob_33_to_66,
-                      probabilities$probs_fcst$prob_a66,
-                      along = 0)
-  names(dim(probs_fcst)) <- c("bin",
-                              names(dim(probabilities$probs_fcst$prob_b33)))
-  
-  ## TODO: Improve this section
-  # Drop extra dims, add time dim if missing:
+
   for (var in 1:fcst$dims[['var']]) {
     ## NOTE: Variables starting with var_* represent the variable counter.
     variable <- fcst$attrs$Variable$varName[[var]]
-- 
GitLab


From b0acbbf605ae4d978dec7e1fac01942c25156302 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Fri, 17 Jan 2025 10:36:44 +0100
Subject: [PATCH 14/24] Add shapefile option for VizEquiMap plots

---
 modules/Visualization/R/plot_forecast_map.R     | 16 ++++++++++++++--
 modules/Visualization/R/plot_metrics.R          | 17 ++++++++++++++---
 .../R/plot_most_likely_terciles_map.R           | 17 ++++++++++++++---
 recipe_template.yml                             |  1 +
 4 files changed, 43 insertions(+), 8 deletions(-)

diff --git a/modules/Visualization/R/plot_forecast_map.R b/modules/Visualization/R/plot_forecast_map.R
index ffc683bb..3bac0757 100644
--- a/modules/Visualization/R/plot_forecast_map.R
+++ b/modules/Visualization/R/plot_forecast_map.R
@@ -16,6 +16,16 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
   if (recipe$Analysis$Variables$freq %in% c("daily", "daily_mean")) {
     stop("Visualization functions not yet implemented for daily data.")
   }
+  # Read shapefile
+  if (!is.null(recipe$Analysis$Workflow$Visualization$shapefile)) {
+    library(rgdal)
+    shp <- readOGR(recipe$Analysis$Workflow$Visualization$shapefile)
+    s1 <- spTransform(shp, CRS("+proj=longlat"))
+    shapefile <- SpatialPolygons2map(s1)
+  } else {
+    shapefile <- NULL
+  }
+  # Coordinates and names
   latitude <- fcst$coords$lat
   longitude <- fcst$coords$lon
   archive <- get_archive(recipe)
@@ -190,7 +200,8 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
                           units = units, cols = cols, brks = brks,
                           fileout = paste0(outfile, ".pdf"),
                           bar_label_digits = 4, extra_margin = rep(1, 4),
-                          bar_label_scale = 1.5, axes_label_scale = 1.1)
+                          bar_label_scale = 1.5, axes_label_scale = 1.1,
+                          shapefile = shapefile)
         base_args[names(output_configuration)] <- output_configuration
         do.call(PlotLayout, base_args)
       } else {
@@ -203,7 +214,8 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
                             font.main = 2, margin_scale = c(1, 5, 5, 5),
                             filled.continents = F, brks = brks, cols = cols,
                             bar_label_digits = 4, bar_label_scale = 1.5,
-                            axes_label_scale = 1, units = units)
+                            axes_label_scale = 1, units = units,
+                            shapefile = shapefile)
         } else {
           fun <- PlotRobinson
           common_projections <- c("robinson", "stereographic", "lambert_europe")
diff --git a/modules/Visualization/R/plot_metrics.R b/modules/Visualization/R/plot_metrics.R
index d52aad1d..e7c892d6 100644
--- a/modules/Visualization/R/plot_metrics.R
+++ b/modules/Visualization/R/plot_metrics.R
@@ -19,7 +19,16 @@ plot_metrics <- function(recipe, data_cube, metrics,
   if (!is.list(metrics) || is.null(names(metrics))) {
     stop("The element 'metrics' must be a list of named arrays.")
   }
-  
+  # Read shapefile
+  if (!is.null(recipe$Analysis$Workflow$Visualization$shapefile)) {
+    library(rgdal)
+    shp <- readOGR(recipe$Analysis$Workflow$Visualization$shapefile)
+    s1 <- spTransform(shp, CRS("+proj=longlat"))
+    shapefile <- SpatialPolygons2map(s1)
+  } else {
+    shapefile <- NULL
+  }
+  # Coordinates and names
   latitude <- data_cube$coords$latitude
   longitude <- data_cube$coords$longitude
   archive <- get_archive(recipe)
@@ -275,7 +284,8 @@ plot_metrics <- function(recipe, data_cube, metrics,
                             bar_label_digits = 3,
                             bar_extra_margin = rep(0.9, 4),
                             extra_margin = rep(1, 4), bar_label_scale = 1.5,
-                            axes_label_scale = 1.3, width = 11, height = 11)
+                            axes_label_scale = 1.3, width = 11, height = 11,
+                            shapefile = shapefile)
           base_args[names(output_configuration)] <- output_configuration
           do.call(PlotLayout, base_args)
         } else {
@@ -292,7 +302,8 @@ plot_metrics <- function(recipe, data_cube, metrics,
                               bar_label_digits = 3, bar_label_scale = 1,
                               axes_label_scale = 1, width = 8, height = 5,
                               bar_extra_margin = rep(0.9, 4),
-                              margin_scale = c(3, 5, 5, 5)) 
+                              margin_scale = c(3, 5, 5, 5),
+                              shapefile = shapefile) 
           } else {
             fun <- VizRobinson
             common_projections <- c("robinson", "stereographic", "lambert_europe")
diff --git a/modules/Visualization/R/plot_most_likely_terciles_map.R b/modules/Visualization/R/plot_most_likely_terciles_map.R
index 963a195f..790e989d 100644
--- a/modules/Visualization/R/plot_most_likely_terciles_map.R
+++ b/modules/Visualization/R/plot_most_likely_terciles_map.R
@@ -23,7 +23,16 @@ plot_most_likely_terciles <- function(recipe,
   if (recipe$Analysis$Variables$freq %in% c("daily", "daily_mean")) {
     stop("Visualization functions not yet implemented for daily data.")
   }
-  
+ 
+  if (!is.null(recipe$Analysis$Workflow$Visualization$shapefile)) {
+    library(rgdal)
+    shp <- readOGR(recipe$Analysis$Workflow$Visualization$shapefile)
+    s1 <- spTransform(shp, CRS("+proj=longlat"))
+    shapefile <- SpatialPolygons2map(s1)
+  } else {
+    shapefile <- NULL
+  }
+
   latitude <- fcst$coords$lat
   longitude <- fcst$coords$lon
   archive <- get_archive(recipe)
@@ -175,7 +184,8 @@ plot_most_likely_terciles <- function(recipe,
                    extra_margin = rep(1, 4),
                    bar_label_scale = 1.2,
                    axes_label_scale = 1.1,
-                   triangle_ends = c(F, F)) # , width = 11, height = 8)
+                   triangle_ends = c(F, F),
+                   shapefile = shapefile) # , width = 11, height = 8)
         )
       } else {
         output_configuration <- output_conf$cylindrical_equidistant$most_likely_terciles
@@ -201,7 +211,8 @@ plot_most_likely_terciles <- function(recipe,
                           plot_margin = c(1, 5, 5, 5),
                           # plot_margin = c(5.1, 4.1, 4.1, 2.1),
                           return_leg = T,
-                          triangle_ends = c(F, T), width = 10, height = 8)
+                          triangle_ends = c(F, T), width = 10, height = 8,
+                          shapefile = shapefile)
         base_args[names(output_configuration)] <- output_configuration
         for (i in 1:length(time_labels)) {
           ## variables ending in *_i represent each forecast time
diff --git a/recipe_template.yml b/recipe_template.yml
index c1287f11..8adc02dd 100644
--- a/recipe_template.yml
+++ b/recipe_template.yml
@@ -146,6 +146,7 @@ Analysis:
       mask_terciles: no # Whether to dot or mask the non-significant points by rpss in the most likely tercile plot. yes/true, no/false or 'both'. Default is no/false. (Optional, str)
       dots_terciles: yes # Whether to dot the non-significant by rpss in the most likely tercile plot. yes/true, no/false or 'both'. Default is no/false. (Optional, str)
       mask_ens: no # Whether to mask the non-significant points by rpss in the forecast ensemble mean plot. yes/true, no/false or 'both'. Default is no/false. (Optional, str)
+      shapefile: # path to a shapefile (*.shp) to include in the plots. Only available for the cylindrical equidistant projection. (Optional, str)
       file_format: 'PNG' # Final file format of the plots. Formats available: PNG, JPG, JPEG, EPS. Defaults to PDF.
     Scorecards:
       execute: yes # yes/no
-- 
GitLab


From b77f3745433699b912a49d007b3541d88397e43b Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Fri, 17 Jan 2025 11:23:58 +0100
Subject: [PATCH 15/24] recipes for testing

---
 .../atomic_recipes/recipe_test_crossval.yml   | 67 +++++++++++++++++++
 .../atomic_recipes/recipe_test_multivar.yml   | 11 +--
 2 files changed, 73 insertions(+), 5 deletions(-)
 create mode 100644 recipes/atomic_recipes/recipe_test_crossval.yml

diff --git a/recipes/atomic_recipes/recipe_test_crossval.yml b/recipes/atomic_recipes/recipe_test_crossval.yml
new file mode 100644
index 00000000..f1314576
--- /dev/null
+++ b/recipes/atomic_recipes/recipe_test_crossval.yml
@@ -0,0 +1,67 @@
+Description:
+  Author: V. Agudetse
+  Description: Multiple variables in the same atomic recipe
+Analysis:
+  Horizon: Seasonal
+  Variables:
+    name: tas
+    freq: monthly_mean
+    units: C
+  Datasets:
+    System:
+      name: ECMWF-SEAS5.1
+    Multimodel: False
+    Reference:
+      name: ERA5
+  Time:
+    sdate: '1101'
+    fcst_year: '2020'
+    hcst_start: '1999'
+    hcst_end: '2010'
+    ftime_min: 1
+    ftime_max: 1
+  Region:
+    name: "Global"
+    latmin: -10
+    latmax: 10
+    lonmin: 0
+    lonmax: 20
+  Regrid:
+    method: bilinear
+    type: to_system
+  Workflow:
+    Calibration:
+      method: raw
+      save: 'exp_only'
+    Anomalies:
+      compute: yes
+      cross_validation: yes
+      save: 'none'
+    Skill:
+      metric: RPS RPSS CRPS CRPSS BSS10 BSS90 EnsCorr mean_bias mean_bias_SS
+      save: 'all'
+      cross_validation: yes
+    Statistics:
+      metric: cov std var n_eff
+      save: 'all'
+    Probabilities:
+      percentiles: [[1/3, 2/3], [1/10, 9/10]]
+      save: 'all'
+    Visualization:
+      plots: most_likely_terciles, skill_metrics, forecast_ensemble_mean
+      forecast_method: mean, median
+      multi_panel: no
+      projection: cylindrical_equidistant
+      dots: both
+      mask_terciles: both
+      file_format: PNG
+    Indicators:
+      index: no
+  ncores: 14
+  remove_NAs: yes
+  Output_format: S2S4E
+Run:
+  Loglevel: INFO
+  Terminal: yes
+  output_dir: /esarchive/scratch/vagudets/auto-s2s-outputs/
+  code_dir: /esarchive/scratch/vagudets/repos/auto-s2s/
diff --git a/recipes/atomic_recipes/recipe_test_multivar.yml b/recipes/atomic_recipes/recipe_test_multivar.yml
index 1e476f4e..feafa8e1 100644
--- a/recipes/atomic_recipes/recipe_test_multivar.yml
+++ b/recipes/atomic_recipes/recipe_test_multivar.yml
@@ -21,10 +21,11 @@ Analysis:
     ftime_min: 1
     ftime_max: 1
   Region:
-    latmin: -10
-    latmax: 10
-    lonmin: 0
-    lonmax: 20
+    name: "Global"
+    latmin: -90 # -10
+    latmax: 90 # 10
+    lonmin: 0 # 0
+    lonmax: 359.9 # 20
   Regrid:
     method: bilinear
     type: to_system
@@ -55,7 +56,7 @@ Analysis:
       file_format: PNG
     Indicators:
       index: no
-  ncores: 7
+  ncores: 14
   remove_NAs: yes
   Output_format: S2S4E
 Run:
-- 
GitLab


From cf7339fbe74a4ed38fe2d07881dca5c9e075c97e Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Fri, 17 Jan 2025 14:54:21 +0100
Subject: [PATCH 16/24] Use sf instead of rgdal for shapefiles

---
 modules/Visualization/R/plot_forecast_map.R             | 5 ++---
 modules/Visualization/R/plot_metrics.R                  | 5 ++---
 modules/Visualization/R/plot_most_likely_terciles_map.R | 5 ++---
 tools/libs.R                                            | 2 ++
 4 files changed, 8 insertions(+), 9 deletions(-)

diff --git a/modules/Visualization/R/plot_forecast_map.R b/modules/Visualization/R/plot_forecast_map.R
index 3bac0757..8db28416 100644
--- a/modules/Visualization/R/plot_forecast_map.R
+++ b/modules/Visualization/R/plot_forecast_map.R
@@ -18,9 +18,8 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
   }
   # Read shapefile
   if (!is.null(recipe$Analysis$Workflow$Visualization$shapefile)) {
-    library(rgdal)
-    shp <- readOGR(recipe$Analysis$Workflow$Visualization$shapefile)
-    s1 <- spTransform(shp, CRS("+proj=longlat"))
+    shp <- read_sf(recipe$Analysis$Workflow$Visualization$shapefile)
+    s1 <- spTransform(as(shp, "Spatial"), CRS("+proj=longlat"))
     shapefile <- SpatialPolygons2map(s1)
   } else {
     shapefile <- NULL
diff --git a/modules/Visualization/R/plot_metrics.R b/modules/Visualization/R/plot_metrics.R
index e7c892d6..bde8b099 100644
--- a/modules/Visualization/R/plot_metrics.R
+++ b/modules/Visualization/R/plot_metrics.R
@@ -21,9 +21,8 @@ plot_metrics <- function(recipe, data_cube, metrics,
   }
   # Read shapefile
   if (!is.null(recipe$Analysis$Workflow$Visualization$shapefile)) {
-    library(rgdal)
-    shp <- readOGR(recipe$Analysis$Workflow$Visualization$shapefile)
-    s1 <- spTransform(shp, CRS("+proj=longlat"))
+    shp <- read_sf(recipe$Analysis$Workflow$Visualization$shapefile)
+    s1 <- spTransform(as(shp, "Spatial"), CRS("+proj=longlat"))
     shapefile <- SpatialPolygons2map(s1)
   } else {
     shapefile <- NULL
diff --git a/modules/Visualization/R/plot_most_likely_terciles_map.R b/modules/Visualization/R/plot_most_likely_terciles_map.R
index 790e989d..f4a361c1 100644
--- a/modules/Visualization/R/plot_most_likely_terciles_map.R
+++ b/modules/Visualization/R/plot_most_likely_terciles_map.R
@@ -25,9 +25,8 @@ plot_most_likely_terciles <- function(recipe,
   }
  
   if (!is.null(recipe$Analysis$Workflow$Visualization$shapefile)) {
-    library(rgdal)
-    shp <- readOGR(recipe$Analysis$Workflow$Visualization$shapefile)
-    s1 <- spTransform(shp, CRS("+proj=longlat"))
+    shp <- read_sf(recipe$Analysis$Workflow$Visualization$shapefile)
+    s1 <- spTransform(as(shp, "Spatial"), CRS("+proj=longlat"))
     shapefile <- SpatialPolygons2map(s1)
   } else {
     shapefile <- NULL
diff --git a/tools/libs.R b/tools/libs.R
index 40146786..481895fb 100644
--- a/tools/libs.R
+++ b/tools/libs.R
@@ -15,6 +15,8 @@ library(clock)
 library(RColorBrewer)
 library(configr)
 library(sf)
+library(sp)
+library(maps)
 library(ggplot2)
 library(rnaturalearth)
 library(cowplot)
-- 
GitLab


From 3d79a0798dfca3d7a2513652e53bd33c392053c8 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Fri, 17 Jan 2025 16:36:50 +0100
Subject: [PATCH 17/24] Fix typo

---
 modules/Visualization/R/plot_metrics.R | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/modules/Visualization/R/plot_metrics.R b/modules/Visualization/R/plot_metrics.R
index bde8b099..ed94dfb9 100644
--- a/modules/Visualization/R/plot_metrics.R
+++ b/modules/Visualization/R/plot_metrics.R
@@ -21,9 +21,10 @@ plot_metrics <- function(recipe, data_cube, metrics,
   }
   # Read shapefile
   if (!is.null(recipe$Analysis$Workflow$Visualization$shapefile)) {
-    shp <- read_sf(recipe$Analysis$Workflow$Visualization$shapefile)
-    s1 <- spTransform(as(shp, "Spatial"), CRS("+proj=longlat"))
-    shapefile <- SpatialPolygons2map(s1)
+    shapefile <- read_sf(recipe$Analysis$Workflow$Visualization$shapefile)
+    shapefile <- spTransform(as(shapefile, "Spatial"),
+                             CRSobj = CRS("+proj=longlat"))
+    shapefile <- SpatialPolygons2map(shapefile)
   } else {
     shapefile <- NULL
   }
@@ -101,7 +102,7 @@ plot_metrics <- function(recipe, data_cube, metrics,
         units <- NULL
         # Define plot characteristics and metric name to display in plot
         if (name %in% c("rpss", "bss90", "bss10", "frpss", "crpss",
-                        "rpss_specs", "bss90:_specs", "bss10_specs",
+                        "rpss_specs", "bss90_specs", "bss10_specs",
                         "rmsss", "msss")) {
           display_name <- toupper(strsplit(name, "_")[[1]][1])
           metric <- var_metric[[name]]
-- 
GitLab


From acdb30034fa1803238cc97fce23e19e02f162980 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Fri, 17 Jan 2025 16:37:19 +0100
Subject: [PATCH 18/24] Avoid duplicating objects

---
 modules/Visualization/R/plot_forecast_map.R             | 7 ++++---
 modules/Visualization/R/plot_most_likely_terciles_map.R | 7 ++++---
 2 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/modules/Visualization/R/plot_forecast_map.R b/modules/Visualization/R/plot_forecast_map.R
index 8db28416..1a590084 100644
--- a/modules/Visualization/R/plot_forecast_map.R
+++ b/modules/Visualization/R/plot_forecast_map.R
@@ -18,9 +18,10 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
   }
   # Read shapefile
   if (!is.null(recipe$Analysis$Workflow$Visualization$shapefile)) {
-    shp <- read_sf(recipe$Analysis$Workflow$Visualization$shapefile)
-    s1 <- spTransform(as(shp, "Spatial"), CRS("+proj=longlat"))
-    shapefile <- SpatialPolygons2map(s1)
+    shapefile <- read_sf(recipe$Analysis$Workflow$Visualization$shapefile)
+    shapefile <- spTransform(as(shapefile, "Spatial"),
+                             CRSobj = CRS("+proj=longlat"))
+    shapefile <- SpatialPolygons2map(shapefile)
   } else {
     shapefile <- NULL
   }
diff --git a/modules/Visualization/R/plot_most_likely_terciles_map.R b/modules/Visualization/R/plot_most_likely_terciles_map.R
index f4a361c1..599e104b 100644
--- a/modules/Visualization/R/plot_most_likely_terciles_map.R
+++ b/modules/Visualization/R/plot_most_likely_terciles_map.R
@@ -25,9 +25,10 @@ plot_most_likely_terciles <- function(recipe,
   }
  
   if (!is.null(recipe$Analysis$Workflow$Visualization$shapefile)) {
-    shp <- read_sf(recipe$Analysis$Workflow$Visualization$shapefile)
-    s1 <- spTransform(as(shp, "Spatial"), CRS("+proj=longlat"))
-    shapefile <- SpatialPolygons2map(s1)
+    shapefile <- read_sf(recipe$Analysis$Workflow$Visualization$shapefile)
+    shapefile <- spTransform(as(shapefile, "Spatial"),
+                             CRSobj = CRS("+proj=longlat"))
+    shapefile <- SpatialPolygons2map(shapefile)
   } else {
     shapefile <- NULL
   }
-- 
GitLab


From 97175d26f4d0287ab878df4f3e29d6a04d09e9f8 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Tue, 21 Jan 2025 09:15:34 +0100
Subject: [PATCH 19/24] Testing

---
 recipes/atomic_recipes/recipe_test_crossval.yml |  4 ++--
 recipes/atomic_recipes/recipe_test_multivar.yml | 12 +++++++-----
 2 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/recipes/atomic_recipes/recipe_test_crossval.yml b/recipes/atomic_recipes/recipe_test_crossval.yml
index f1314576..f4dea7ab 100644
--- a/recipes/atomic_recipes/recipe_test_crossval.yml
+++ b/recipes/atomic_recipes/recipe_test_crossval.yml
@@ -31,7 +31,7 @@ Analysis:
     type: to_system
   Workflow:
     Calibration:
-      method: raw
+      method: mse_min
       save: 'exp_only'
     Anomalies:
       compute: yes
@@ -45,7 +45,7 @@ Analysis:
       metric: cov std var n_eff
       save: 'all'
     Probabilities:
-      percentiles: [[1/3, 2/3], [1/10, 9/10]]
+      percentiles: [[1/3, 2/3], [1/10], [9/10]] # , [9/10]]
       save: 'all'
     Visualization:
       plots: most_likely_terciles, skill_metrics, forecast_ensemble_mean
diff --git a/recipes/atomic_recipes/recipe_test_multivar.yml b/recipes/atomic_recipes/recipe_test_multivar.yml
index feafa8e1..3e26ff0b 100644
--- a/recipes/atomic_recipes/recipe_test_multivar.yml
+++ b/recipes/atomic_recipes/recipe_test_multivar.yml
@@ -21,11 +21,12 @@ Analysis:
     ftime_min: 1
     ftime_max: 1
   Region:
-    name: "Global"
-    latmin: -90 # -10
-    latmax: 90 # 10
-    lonmin: 0 # 0
-    lonmax: 359.9 # 20
+    - {name: "Iberia", latmin: 36, latmax: 44, lonmin: -10, lonmax: 5}
+    # name: "Global"
+    # latmin: -90 # -10
+    # latmax: 90 # 10
+    # lonmin: 0 # 0
+    # lonmax: 359.9 # 20
   Regrid:
     method: bilinear
     type: to_system
@@ -53,6 +54,7 @@ Analysis:
       projection: cylindrical_equidistant
       dots: both
       mask_terciles: both
+      shapefile: /esarchive/scratch/cdelgado/aspect_outputs/casestudy-wine/shp_spain/recintos_provinciales_inspire_peninbal_etrs89/recintos_provinciales_inspire_peninbal_etrs89.shp
       file_format: PNG
     Indicators:
       index: no
-- 
GitLab


From b1fa6cbbe318fd9f6522e97dc059e26c2a55fd3a Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Thu, 23 Jan 2025 13:01:24 +0100
Subject: [PATCH 20/24] Refactor file_format conversion and logo so that they
 are run individually for each plot

---
 modules/Visualization/R/plot_forecast_map.R   | 18 ++++-
 modules/Visualization/R/plot_metrics.R        | 20 +++++-
 .../R/plot_most_likely_terciles_map.R         | 18 ++++-
 modules/Visualization/Visualization.R         | 67 +++++++++----------
 tools/add_logo.R                              | 37 ++++++----
 tools/libs.R                                  |  1 +
 6 files changed, 108 insertions(+), 53 deletions(-)

diff --git a/modules/Visualization/R/plot_forecast_map.R b/modules/Visualization/R/plot_forecast_map.R
index 1a590084..53685586 100644
--- a/modules/Visualization/R/plot_forecast_map.R
+++ b/modules/Visualization/R/plot_forecast_map.R
@@ -10,7 +10,8 @@
 #   'median' (default), 'iqr', 'mean'.
 plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL, 
                               outdir, output_conf,
-                              method = 'median') {
+                              method = 'median',
+                              logo = NULL) {
   ## TODO: Add 'anomaly' to plot title
   # Abort if frequency is daily
   if (recipe$Analysis$Variables$freq %in% c("daily", "daily_mean")) {
@@ -302,6 +303,21 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
                   args = c(base_args,
                            list(toptitle = toptitle,
                                 fileout = fileout)))
+          if (!is.null(recipe$Analysis$Workflow$Visualization$file_format) &&
+              tolower(recipe$Analysis$Workflow$Visualization$file_format) != "pdf") {
+            extension <- tolower(recipe$Analysis$Workflow$Visualization$file_format)
+            system_command <- paste("convert -density 300", fileout,
+                                    "-resize 40% -alpha remove",
+                                     paste0(tools::file_path_sans_ext(fileout), ".",
+                                            extension))
+            system(system_command)
+            unlink(fileout)
+            fileout <- paste0(tools::file_path_sans_ext(fileout), ".",
+                              extension)
+          }
+          if (!is.null(logo)) {
+            add_logo(file = fileout, logo = logo)
+          }
         }        
       }
     }
diff --git a/modules/Visualization/R/plot_metrics.R b/modules/Visualization/R/plot_metrics.R
index ed94dfb9..fcbbd168 100644
--- a/modules/Visualization/R/plot_metrics.R
+++ b/modules/Visualization/R/plot_metrics.R
@@ -2,7 +2,8 @@ library(stringr)
 library(lubridate)
 
 plot_metrics <- function(recipe, data_cube, metrics,
-                         outdir, significance = F, output_conf) {
+                         outdir, significance = F, output_conf,
+                         logo = NULL) {
   # recipe: Auto-S2S recipe
   # archive: Auto-S2S archive
   # data_cube: s2dv_cube object with the corresponding hindcast data
@@ -474,6 +475,23 @@ plot_metrics <- function(recipe, data_cube, metrics,
                     args = c(base_args,
                              list(toptitle = toptitle,
                                   fileout = fileout)))
+            # Convert plots to user-chosen format
+            if (!is.null(recipe$Analysis$Workflow$Visualization$file_format) &&
+                tolower(recipe$Analysis$Workflow$Visualization$file_format) != "pdf") {
+              extension <- tolower(recipe$Analysis$Workflow$Visualization$file_format)
+              system_command <- paste("convert -density 300", fileout,
+                                      "-resize 40% -alpha remove",
+                                      paste0(tools::file_path_sans_ext(fileout), ".",
+                                      extension))
+              system(system_command)
+              unlink(fileout)
+              fileout <- paste0(tools::file_path_sans_ext(fileout), ".",
+                                extension)
+            }
+            # Add logo
+            if (!is.null(logo)) {
+              add_logo(file = fileout, logo = logo)
+            }
           }
         }
       }
diff --git a/modules/Visualization/R/plot_most_likely_terciles_map.R b/modules/Visualization/R/plot_most_likely_terciles_map.R
index 599e104b..35eefebe 100644
--- a/modules/Visualization/R/plot_most_likely_terciles_map.R
+++ b/modules/Visualization/R/plot_most_likely_terciles_map.R
@@ -16,7 +16,8 @@ plot_most_likely_terciles <- function(recipe,
                                       mask,
                                       dots, 
                                       outdir,
-                                      output_conf) {
+                                      output_conf,
+                                      logo = NULL) {
   
   ## TODO: Add 'anomaly' to plot title
   # Abort if frequency is daily
@@ -318,6 +319,21 @@ plot_most_likely_terciles <- function(recipe,
             }
           }
           dev.off()
+          if (!is.null(recipe$Analysis$Workflow$Visualization$file_format) &&
+              tolower(recipe$Analysis$Workflow$Visualization$file_format) != "pdf") {
+            extension <- tolower(recipe$Analysis$Workflow$Visualization$file_format)
+            system_command <- paste("convert -density 300", fileout,
+                                    "-resize 40% -alpha remove",
+                                     paste0(tools::file_path_sans_ext(fileout), ".",
+                                            extension))
+            system(system_command)
+            unlink(fileout)
+            fileout <- paste0(tools::file_path_sans_ext(fileout), ".",
+                              extension)
+          }
+          if (!is.null(logo)) {
+            add_logo(file = fileout, logo = logo)
+          }
         }
       }
     }
diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R
index 3b0d5f22..8d6557c6 100644
--- a/modules/Visualization/Visualization.R
+++ b/modules/Visualization/Visualization.R
@@ -24,7 +24,8 @@ Visualization <- function(recipe,
                           statistics = NULL,
                           probabilities = NULL,
                           significance = F,
-                          output_conf = NULL) {
+                          output_conf = NULL,
+                          logo = NULL) {
   # Try to produce and save several basic plots.
   # recipe: the auto-s2s recipe as read by read_yaml()
   # data: list containing the hcst, obs and (optional) fcst s2dv_cube objects
@@ -58,6 +59,8 @@ Visualization <- function(recipe,
       output_conf <- NULL
     }
   }
+  # Check logo
+
 
   # Get plot types and create output directories
   plots <- strsplit(recipe$Analysis$Workflow$Visualization$plots, ", | |,")[[1]]
@@ -102,21 +105,22 @@ Visualization <- function(recipe,
       if (is.logical(significance)) {
         plot_metrics(recipe = recipe, data_cube = data$hcst, 
                      metrics = skill_metrics, outdir = outdir,
-                     significance = significance, output_conf = output_conf)
+                     significance = significance, output_conf = output_conf,
+                     logo = logo)
         info(recipe$Run$logger,
              paste("Skill metrics significance set as", significance))
-        } else {
-          if (significance %in% c('both', 'dots')) {
-            plot_metrics(recipe, data$hcst, skill_metrics, outdir,
-                         significance = 'dots', output_conf = output_conf)
-            info(recipe$Run$logger, "Skill metrics significance as dots")
-            }
-          if (significance %in% c('both', 'mask')) {
-            plot_metrics(recipe, data$hcst, skill_metrics, outdir,
-                         significance = 'mask', output_conf = output_conf)
-            info(recipe$Run$logger, "Skill metrics significance as mask")
+      } else {
+        if (significance %in% c('both', 'dots')) {
+          plot_metrics(recipe, data$hcst, skill_metrics, outdir,
+                       significance = 'dots', output_conf = output_conf)
+          info(recipe$Run$logger, "Skill metrics significance as dots")
           }
+        if (significance %in% c('both', 'mask')) {
+          plot_metrics(recipe, data$hcst, skill_metrics, outdir,
+                       significance = 'mask', output_conf = output_conf)
+          info(recipe$Run$logger, "Skill metrics significance as mask")
         }
+      }
     } else {
       error(recipe$Run$logger,
             paste0("The skill metric plots have been requested, but the ",
@@ -128,7 +132,8 @@ Visualization <- function(recipe,
   if ("statistics" %in% plots) {
     if (!is.null(statistics)) {
       plot_metrics(recipe, data$hcst, statistics, outdir,
-                   significance, output_conf = output_conf)
+                   significance, output_conf = output_conf,
+                   logo = logo)
     } else {
       error(recipe$Run$logger,
             paste0("The statistics plots have been requested, but the ",
@@ -158,7 +163,8 @@ Visualization <- function(recipe,
           plot_forecast_map(recipe, data$fcst,
                             mask = NULL, dots = NULL,
                             outdir = outdir, method = method,
-                            output_conf = output_conf)
+                            output_conf = output_conf,
+                            logo = logo)
         }
         # Plots with masked
         if (recipe$Analysis$Workflow$Visualization$mask_ens %in%
@@ -176,7 +182,8 @@ Visualization <- function(recipe,
                               mask = skill_metrics$enscorr,
                               dots = NULL,
                               outdir = outdir, method = method,
-                              output_conf = output_conf)
+                              output_conf = output_conf,
+                              logo = logo)
           }
         }
         # Plots with dotted negative correlated in ens-mean-fcst
@@ -194,7 +201,8 @@ Visualization <- function(recipe,
                               mask = NULL,
                               dots = skill_metrics$enscorr,
                               outdir = outdir, method = method,
-                              output_conf = output_conf)
+                              output_conf = output_conf,
+                              logo = logo)
           }
         }
       } ## End loop over methods
@@ -223,7 +231,8 @@ Visualization <- function(recipe,
                                   probabilities,
                                   mask = NULL,
                                   dots = NULL,
-                                  outdir, output_conf = output_conf)
+                                  outdir, output_conf = output_conf,
+                                  logo = logo)
       }
       # Plots with masked terciles
       if (recipe$Analysis$Workflow$Visualization$mask_terciles %in% 
@@ -241,7 +250,8 @@ Visualization <- function(recipe,
                                     probabilities,
                                     mask = skill_metrics$rpss,
                                     dots = NULL,
-                                    outdir, output_conf = output_conf)
+                                    outdir, output_conf = output_conf,
+                                    logo = logo)
         }
       }
       # Plots with dotted terciles
@@ -259,7 +269,8 @@ Visualization <- function(recipe,
                                     probabilities,
                                     mask = NULL,
                                     dots = skill_metrics$rpss,
-                                    outdir, output_conf = output_conf)
+                                    outdir, output_conf = output_conf,
+                                    logo = logo)
         }
       }
     } else {
@@ -268,22 +279,4 @@ Visualization <- function(recipe,
                    "probabilities must be provided."))
     }
   }
-
-  # Convert plots to user-chosen format
-  if (!is.null(recipe$Analysis$Workflow$Visualization$file_format) &&
-      tolower(recipe$Analysis$Workflow$Visualization$file_format) != "pdf") {
-    extension <- tolower(recipe$Analysis$Workflow$Visualization$file_format)
-    plot_files <- list.files(path = recipe$Run$output_dir, pattern = "\\.pdf$",
-                             recursive = TRUE, full.names = TRUE)
-    for (file in plot_files) {
-      system_command <- paste("convert -density 300", file,
-                              "-resize 40% -alpha remove",
-                              paste0(tools::file_path_sans_ext(file), ".",
-                                     extension))
-      system(system_command)
-    }
-    unlink(plot_files)
-    info(recipe$Run$logger,
-         paste0("##### PLOT FILES CONVERTED TO ", toupper(extension), " #####"))
-  }
 }
diff --git a/tools/add_logo.R b/tools/add_logo.R
index 42fb87c5..03603cd1 100644
--- a/tools/add_logo.R
+++ b/tools/add_logo.R
@@ -1,15 +1,26 @@
-add_logo <- function(recipe, logo) {
-  # recipe: SUNSET recipe
+#' Function to add a logo onto an image
+#'
+#'@description This function adds a logo to the bottom right corner of the 
+# image provided in 'file'.
+#'@param file The path to the image 
+#'@param logo The path to the logo image
+#'@param logo_resize_percentage Scale factor for the logo, 1 being the original
+#' size. The default value is 0.25.
+
+add_logo <- function(file, logo, logo_resize_percentage = 0.25) {
+  # file
   # logo: URL to the logo
-  system <- list.files(paste0(recipe$Run$output_dir, "/plots/"))
-  variable <- recipe$Analysis$Variable$name
-  files <- lapply(variable, function(x) {
-                  f <- list.files(paste0(recipe$Run$output_dir, "/plots/",
-                                         system, "/", x))
-                  full_path <- paste0(recipe$Run$output_dir, "/plots/",
-                                      system, "/", x,"/", f)})[[1]]
-  dim(files) <- c(file = length(files))
-  Apply(list(files), target_dims = NULL, function(x) {
-        system(paste("composite -gravity southeast -geometry +10+10",
-                     logo, x, x))}, ncores = recipe$Analysis$ncores)
+  fig_width <- as.numeric(system(paste("identify -format '%w'", file), 
+					intern = TRUE))
+  logo_height <- as.numeric(system(paste("convert", logo, 
+		"-resize", fig_width * 0.1, "-format '%h' info:", sep=" "),
+					 intern = TRUE))
+  system(paste0("convert ", file, 
+		" -gravity south -background white -splice 0x",
+		logo_height, " extended_fig.png"))
+  system(paste0("convert extended_fig.png \\( ", 
+		logo, " -resize ", fig_width * logo_resize_percentage,
+		" \\) -gravity southeast -composite ",
+		file))     
+  file.remove("extended_fig.png")
 }
diff --git a/tools/libs.R b/tools/libs.R
index 481895fb..9727e520 100644
--- a/tools/libs.R
+++ b/tools/libs.R
@@ -38,6 +38,7 @@ source("tools/write_autosubmit_conf.R")
 source("tools/get_archive.R")
 source("tools/Utils.R")
 source("tools/restructure_recipe.R")
+source("tools/add_logo.R")
 # source("tools/add_dims.R") # Not sure if necessary yet
 
 # Settings
-- 
GitLab


From d119f3a8118b535ef298514a44cfca5d9af0ea17 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Thu, 23 Jan 2025 15:24:20 +0100
Subject: [PATCH 21/24] Bugfix: typo ('forecat_time' instead of 'forecast_time'

---
 modules/Visualization/R/plot_metrics.R     |   2 +-
 tools/{BSC_logo_95.jpg => bsc_logo_95.jpg} | Bin
 2 files changed, 1 insertion(+), 1 deletion(-)
 rename tools/{BSC_logo_95.jpg => bsc_logo_95.jpg} (100%)

diff --git a/modules/Visualization/R/plot_metrics.R b/modules/Visualization/R/plot_metrics.R
index fcbbd168..0f39f112 100644
--- a/modules/Visualization/R/plot_metrics.R
+++ b/modules/Visualization/R/plot_metrics.R
@@ -352,7 +352,7 @@ plot_metrics <- function(recipe, data_cube, metrics,
                   ## TODO: There's probably a better way to do this, like using
                   ## time_bounds$start and time_bounds$end directly.
                   forecast_time_ini <- init_month + forecast_time_ini - 1
-                  forecat_time_ini <- ifelse(forecast_time_ini > 12, forecast_time_ini - 12, forecast_time_ini) 
+                  forecast_time_ini <- ifelse(forecast_time_ini > 12, forecast_time_ini - 12, forecast_time_ini) 
                   forecast_time_ini <- month.name[forecast_time_ini]
                   forecast_time_end <- init_month + forecast_time_end - 1
                   forecast_time_end <- ifelse(forecast_time_end > 12, forecast_time_end - 12, forecast_time_end)
diff --git a/tools/BSC_logo_95.jpg b/tools/bsc_logo_95.jpg
similarity index 100%
rename from tools/BSC_logo_95.jpg
rename to tools/bsc_logo_95.jpg
-- 
GitLab


From d7a4915e8166329ddc2aa7f1a1aff9480fbc1af4 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Thu, 23 Jan 2025 15:47:05 +0100
Subject: [PATCH 22/24] Refine handling of logo; add condition not to apply
 tolower() to existing file paths

---
 modules/Visualization/Visualization.R           |   9 +++++----
 recipe_template.yml                             |   1 +
 recipes/atomic_recipes/recipe_test_multivar.yml |   1 +
 tests/recipes/recipe-seasonal_daily_1.yml       |   2 +-
 tools/{bsc_logo_95.jpg => BSC_logo_95.jpg}      | Bin
 tools/prepare_outputs.R                         |   8 +++++++-
 6 files changed, 15 insertions(+), 6 deletions(-)
 rename tools/{bsc_logo_95.jpg => BSC_logo_95.jpg} (100%)

diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R
index 8d6557c6..e5550344 100644
--- a/modules/Visualization/Visualization.R
+++ b/modules/Visualization/Visualization.R
@@ -59,9 +59,6 @@ Visualization <- function(recipe,
       output_conf <- NULL
     }
   }
-  # Check logo
-
-
   # Get plot types and create output directories
   plots <- strsplit(recipe$Analysis$Workflow$Visualization$plots, ", | |,")[[1]]
   ## TODO: Do not modify output dir here
@@ -93,11 +90,15 @@ Visualization <- function(recipe,
               "the value in the recipe. The function call value will be used.")
     }
   }
-  
   # If not specified in function call, set significance to recipe value
   if (missing(significance) && !is.null(recipe$Analysis$Workflow$Visualization$significance)) {
     significance <- recipe$Analysis$Workflow$Visualization$significance
   }
+
+  # If not specified in function call, set logo to recipe value
+  if (missing(logo) && !is.null(recipe$Analysis$Workflow$Visualization$logo)) {
+    logo <- recipe$Analysis$Workflow$Visualization$logo
+  }
   
   # Plot skill metrics
   if ("skill_metrics" %in% plots) {
diff --git a/recipe_template.yml b/recipe_template.yml
index 8adc02dd..8e35edab 100644
--- a/recipe_template.yml
+++ b/recipe_template.yml
@@ -147,6 +147,7 @@ Analysis:
       dots_terciles: yes # Whether to dot the non-significant by rpss in the most likely tercile plot. yes/true, no/false or 'both'. Default is no/false. (Optional, str)
       mask_ens: no # Whether to mask the non-significant points by rpss in the forecast ensemble mean plot. yes/true, no/false or 'both'. Default is no/false. (Optional, str)
       shapefile: # path to a shapefile (*.shp) to include in the plots. Only available for the cylindrical equidistant projection. (Optional, str)
+      logo: tools/BSC_logo_95.jpg # path to a logo (*.png or *.jpg/jpeg) to include in the plots. (Optional, str (Optional, str))
       file_format: 'PNG' # Final file format of the plots. Formats available: PNG, JPG, JPEG, EPS. Defaults to PDF.
     Scorecards:
       execute: yes # yes/no
diff --git a/recipes/atomic_recipes/recipe_test_multivar.yml b/recipes/atomic_recipes/recipe_test_multivar.yml
index 3e26ff0b..e331f02a 100644
--- a/recipes/atomic_recipes/recipe_test_multivar.yml
+++ b/recipes/atomic_recipes/recipe_test_multivar.yml
@@ -55,6 +55,7 @@ Analysis:
       dots: both
       mask_terciles: both
       shapefile: /esarchive/scratch/cdelgado/aspect_outputs/casestudy-wine/shp_spain/recintos_provinciales_inspire_peninbal_etrs89/recintos_provinciales_inspire_peninbal_etrs89.shp
+      logo: tools/BSC_logo_95.jpg
       file_format: PNG
     Indicators:
       index: no
diff --git a/tests/recipes/recipe-seasonal_daily_1.yml b/tests/recipes/recipe-seasonal_daily_1.yml
index 4378730f..926b4458 100644
--- a/tests/recipes/recipe-seasonal_daily_1.yml
+++ b/tests/recipes/recipe-seasonal_daily_1.yml
@@ -2,7 +2,7 @@ Description:
   Author: V. Agudetse
 
 Analysis:
-  Horizon: Seasonal
+  Horizon: seasonal
   Variables:
     - {name: tas, freq: daily}
   Datasets:
diff --git a/tools/bsc_logo_95.jpg b/tools/BSC_logo_95.jpg
similarity index 100%
rename from tools/bsc_logo_95.jpg
rename to tools/BSC_logo_95.jpg
diff --git a/tools/prepare_outputs.R b/tools/prepare_outputs.R
index 5bdff99e..84670d47 100644
--- a/tools/prepare_outputs.R
+++ b/tools/prepare_outputs.R
@@ -114,7 +114,13 @@ prepare_outputs <- function(recipe_file,
   reference <- recipe$Analysis$Datasets$Reference
   region <- recipe$Analysis$Region
   variable <- recipe$Analysis$Variables
-  recipe$Analysis <- rapply(recipe$Analysis, tolower,
+  recipe$Analysis <- rapply(recipe$Analysis,
+                            function(x) {
+                              if (all(!file.exists(x))) {
+                                x <- tolower(x)
+                              }
+                              return(x)
+                            },
                             how = "replace",
                             classes = "character")
   recipe$Analysis$Datasets$System <- system
-- 
GitLab


From bf893d6e1c5d35eb865a4713d546db2da0ae0ff8 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Fri, 24 Jan 2025 10:39:00 +0100
Subject: [PATCH 23/24] Finish implementation of custom dictionary for variable
 longnames, add unit test

---
 .../Visualization/R/get_variable_longname.R   |  8 +++-
 modules/Visualization/R/plot_forecast_map.R   |  3 +-
 modules/Visualization/R/plot_metrics.R        |  3 +-
 .../R/plot_most_likely_terciles_map.R         |  3 +-
 modules/Visualization/var_long_names.yml      |  8 ++++
 recipe_template.yml                           |  5 +-
 .../test-utils-get_variable_longname.R        | 46 +++++++++++++++++++
 7 files changed, 69 insertions(+), 7 deletions(-)
 create mode 100644 modules/Visualization/var_long_names.yml
 create mode 100644 tests/testthat/test-utils-get_variable_longname.R

diff --git a/modules/Visualization/R/get_variable_longname.R b/modules/Visualization/R/get_variable_longname.R
index e655f10e..1c03471c 100644
--- a/modules/Visualization/R/get_variable_longname.R
+++ b/modules/Visualization/R/get_variable_longname.R
@@ -15,11 +15,15 @@
   var_name <- data_cube$attrs$Variable$varName[[index]]
   var_long_name <- tryCatch(
     expr = {
-      var_long_name <- read_yaml("modules/Visualization/var_long_names.yml")[[format]][[var_name]]
+      dictionary <- read_yaml("modules/Visualization/var_long_names.yml")
+      dictionary_index <- grep(format, names(dictionary),
+                               ignore.case = TRUE,
+                               value = TRUE)
+      var_long_name <- dictionary[[dictionary_index]][[var_name]]
       return(var_long_name)
     },
     error = function(e) {
-      var_long_name <- str_to_title(data_cube$attrs$Variable$metadata[[var_name]]$long_name)
+      var_long_name <- str_to_title(data_cube$attrs$Variable$metadata[[var_name]]$long)
       warning("No configuration file found for this variable long name. Using s2dv_cube metadata.")
       return(var_long_name)
     }
diff --git a/modules/Visualization/R/plot_forecast_map.R b/modules/Visualization/R/plot_forecast_map.R
index 53685586..fac78092 100644
--- a/modules/Visualization/R/plot_forecast_map.R
+++ b/modules/Visualization/R/plot_forecast_map.R
@@ -78,9 +78,10 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
   # Loop over variable dimension
   for (var in 1:fcst$dims[['var']]) {
     variable <- fcst$attrs$Variable$varName[[var]]
+    custom_dictionary <- recipe$Analysis$Visualization$custom_dictionary
     var_long_name <- .get_variable_longname(data_cube = fcst,
                                             index = var,
-                                            format = NULL) ## TODO: change
+                                            format = custom_dictionary)
     units <- fcst$attrs$Variable$metadata[[variable]]$units
     # Subset ensemble mean by variable
     var_fcst <- ClimProjDiags::Subset(forecast_product,
diff --git a/modules/Visualization/R/plot_metrics.R b/modules/Visualization/R/plot_metrics.R
index 0f39f112..1ad06734 100644
--- a/modules/Visualization/R/plot_metrics.R
+++ b/modules/Visualization/R/plot_metrics.R
@@ -233,9 +233,10 @@ plot_metrics <- function(recipe, data_cube, metrics,
         
         # Get variable name and long name
         var_name <- data_cube$attrs$Variable$varName[[var]]
+        custom_dictionary <- recipe$Analysis$Visualization$custom_dictionary
         var_long_name <- .get_variable_longname(data_cube = data_cube,
                                                 index = var,
-                                                format = NULL) ## TODO: change
+                                                format = custom_dictionary)
         # Multi-panel or single-panel plots
         if (recipe$Analysis$Workflow$Visualization$multi_panel) {
           # Define titles
diff --git a/modules/Visualization/R/plot_most_likely_terciles_map.R b/modules/Visualization/R/plot_most_likely_terciles_map.R
index 35eefebe..886db225 100644
--- a/modules/Visualization/R/plot_most_likely_terciles_map.R
+++ b/modules/Visualization/R/plot_most_likely_terciles_map.R
@@ -83,9 +83,10 @@ plot_most_likely_terciles <- function(recipe,
   for (var in 1:fcst$dims[['var']]) {
     ## NOTE: Variables starting with var_* represent the variable counter.
     variable <- fcst$attrs$Variable$varName[[var]]
+    custom_dictionary <- recipe$Analysis$Visualization$custom_dictionary
     var_long_name <- .get_variable_longname(data_cube = fcst,
                                             index = var,
-                                            format = NULL) ## TODO: change
+                                            format = custom_dictionary)
     # Choose colors depending on the variable
     if (variable %in% c('prlr')) { ## add others
       cols <- list(c("#FFC473", "#FFAB38"), 
diff --git a/modules/Visualization/var_long_names.yml b/modules/Visualization/var_long_names.yml
new file mode 100644
index 00000000..a6293406
--- /dev/null
+++ b/modules/Visualization/var_long_names.yml
@@ -0,0 +1,8 @@
+Mycase:
+  tas: "2-metre air temperature"
+  tasmax: "Maximum 2-metre air temperature"
+  tasmin: "Minimum 2-metre air temperature"
+  prlr: "Total precipitation"
+  sfcWind: "Wind speed"
+  rsds: "Solar radiation"
+
diff --git a/recipe_template.yml b/recipe_template.yml
index 8e35edab..0a7b9d64 100644
--- a/recipe_template.yml
+++ b/recipe_template.yml
@@ -147,8 +147,9 @@ Analysis:
       dots_terciles: yes # Whether to dot the non-significant by rpss in the most likely tercile plot. yes/true, no/false or 'both'. Default is no/false. (Optional, str)
       mask_ens: no # Whether to mask the non-significant points by rpss in the forecast ensemble mean plot. yes/true, no/false or 'both'. Default is no/false. (Optional, str)
       shapefile: # path to a shapefile (*.shp) to include in the plots. Only available for the cylindrical equidistant projection. (Optional, str)
-      logo: tools/BSC_logo_95.jpg # path to a logo (*.png or *.jpg/jpeg) to include in the plots. (Optional, str (Optional, str))
-      file_format: 'PNG' # Final file format of the plots. Formats available: PNG, JPG, JPEG, EPS. Defaults to PDF.
+      logo: tools/BSC_logo_95.jpg # path to a logo (*.png or *.jpg/jpeg) to include in the plots. (Optional, str)
+      file_format: 'PNG' # Final file format of the plots. Formats available: PNG, JPG, JPEG, EPS. Defaults to PDF. (Optional, str)
+      custom_dictionary: 'MyCase' # Name of the desired custom dictionary entry in modules/Visualization/var_long_names.yml. This allows the user to define custom long names for each variable to be displayed in the plots. (Optional, str)
     Scorecards:
       execute: yes # yes/no
       regions: 
diff --git a/tests/testthat/test-utils-get_variable_longname.R b/tests/testthat/test-utils-get_variable_longname.R
new file mode 100644
index 00000000..2ebfa717
--- /dev/null
+++ b/tests/testthat/test-utils-get_variable_longname.R
@@ -0,0 +1,46 @@
+context("get_plot_time_labels.R")
+
+source("tools/libs.R")
+source("modules/Visualization/R/get_variable_longname.R")
+###########################################
+
+data_cube <- CSTools::lonlat_temp$exp
+index <- 1
+test_that("1. Get variable longname from dictionary", {
+  var_longname <- .get_variable_longname(data_cube = data_cube,
+                                         index = index,
+                                         format = "MyCase")
+  # Results
+  expect_equal(
+    var_longname,
+    "2-metre air temperature"
+  )
+})
+
+test_that("2. Format not present in dictionary", {
+  expect_warning(
+    var_longname <- .get_variable_longname(data_cube = data_cube,
+                                           index = index,
+                                           format = "Null"),
+    "No configuration file found for this variable long name. Using s2dv_cube metadata."
+  )
+  # Results
+  expect_equal(
+    var_longname,
+    "2 Metre Temperature"
+  )
+})
+
+test_that("3. Format is NULL", {
+  expect_warning(
+    var_longname <- .get_variable_longname(data_cube = data_cube,
+                                           index = index,
+                                           format = NULL),
+    "No configuration file found for this variable long name. Using s2dv_cube metadata."
+  )
+  # Results
+  expect_equal(
+    var_longname,
+    "2 Metre Temperature"
+  )
+})
-- 
GitLab


From 45fab9bf26ff4c27e3ee35d719a6cbe6684ce28e Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Wed, 5 Feb 2025 12:39:05 +0100
Subject: [PATCH 24/24] Adjust settings to accommodate caption

---
 modules/Visualization/R/plot_forecast_map.R | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/modules/Visualization/R/plot_forecast_map.R b/modules/Visualization/R/plot_forecast_map.R
index fac78092..090056f7 100644
--- a/modules/Visualization/R/plot_forecast_map.R
+++ b/modules/Visualization/R/plot_forecast_map.R
@@ -201,6 +201,7 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
                           subtitle_margin_scale = 2, titles = titles,
                           units = units, cols = cols, brks = brks,
                           fileout = paste0(outfile, ".pdf"),
+                          bar_extra_margin = rep(0.9, 4),
                           bar_label_digits = 4, extra_margin = rep(1, 4),
                           bar_label_scale = 1.5, axes_label_scale = 1.1,
                           shapefile = shapefile)
@@ -212,8 +213,9 @@ plot_forecast_map <- function(recipe, fcst, mask = NULL, dots = NULL,
           fun <- VizEquiMap
           base_args <- list(var = NULL, dots = NULL, mask = NULL,
                             lon = longitude, lat = latitude,
-                            dot_symbol = 20, title_scale = 0.6,
-                            font.main = 2, margin_scale = c(1, 5, 5, 5),
+                            dot_symbol = 20, title_scale = 0.7,
+                            font.main = 2, margin_scale = c(3, 5, 5, 5),
+                            bar_extra_margin = rep(0.9, 4),
                             filled.continents = F, brks = brks, cols = cols,
                             bar_label_digits = 4, bar_label_scale = 1.5,
                             axes_label_scale = 1, units = units,
-- 
GitLab