diff --git a/modules/Indices/R/compute_nao.R b/modules/Indices/R/compute_nao.R
index 5bb76e4edd2a6fa10434e3872ca1bf7da2bb63c1..bdf2bb858fea35ac7f3b0b349345d3697f404523 100644
--- a/modules/Indices/R/compute_nao.R
+++ b/modules/Indices/R/compute_nao.R
@@ -182,7 +182,7 @@ compute_nao <- function(data, recipe, obsproj, plot_ts, plot_sp,
         plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/NAO_",
                            system, "_", recipe$Analysis$Datasets$Reference$name,
                            "_s", recipe$Analysis$Time$sdate, "_ftime",
-                           sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min), ".png") 
+                           sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min), ".pdf") 
         caption <- paste0("NAO method: ",
                           ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
                                  "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
@@ -199,7 +199,7 @@ compute_nao <- function(data, recipe, obsproj, plot_ts, plot_sp,
         plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/NAO_",
                            system, "_", recipe$Analysis$Datasets$Reference$name,
                            "_ftime",
-                           sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min), ".png") 
+                           sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min), ".pdf") 
         caption <- paste0("NAO method: ",
                           ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
                                  "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
@@ -277,7 +277,7 @@ compute_nao <- function(data, recipe, obsproj, plot_ts, plot_sp,
                            recipe$Analysis$Variable$name, "_",
                            recipe$Analysis$Datasets$Reference$name,
                            "_s", recipe$Analysis$Time$sdate,
-                           "_ftime", fmonth, ".png")
+                           "_ftime", fmonth, ".pdf")
         caption <-  paste0("NAO method: ",
                            ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
                                   "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
@@ -295,7 +295,7 @@ compute_nao <- function(data, recipe, obsproj, plot_ts, plot_sp,
         plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/NAO_correlation_",
                            recipe$Analysis$Variable$name, "_",
                            recipe$Analysis$Datasets$Reference$name,
-                           "_ftime", fmonth, ".png")
+                           "_ftime", fmonth, ".pdf")
         caption <-  paste0("NAO method: ",
                            ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
                                   "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
@@ -335,7 +335,7 @@ compute_nao <- function(data, recipe, obsproj, plot_ts, plot_sp,
                            recipe$Analysis$Variable$name, "_ensmean_",
                            system,
                            "_s", recipe$Analysis$Time$sdate,
-                           "_ftime", fmonth, ".png")
+                           "_ftime", fmonth, ".pdf")
         caption <-  paste0("NAO method: ",
                            ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
                                   "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
@@ -355,7 +355,7 @@ compute_nao <- function(data, recipe, obsproj, plot_ts, plot_sp,
                            recipe$Analysis$Variable$name, "_ensmean_",
                            system, 
                            recipe$Analysis$Datasets$Reference$name,
-                           "_ftime", fmonth, ".png")
+                           "_ftime", fmonth, ".pdf")
         caption <-  paste0("NAO method: ",
                            ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
                                   "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
@@ -392,7 +392,7 @@ compute_nao <- function(data, recipe, obsproj, plot_ts, plot_sp,
                            recipe$Analysis$Variable$name, "_member_",
                            system,
                            "_s", recipe$Analysis$Time$sdate,
-                           "_ftime", fmonth, ".png")
+                           "_ftime", fmonth, ".pdf")
         caption <-  paste0("NAO method: ",
                            ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
                                   "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
@@ -412,7 +412,7 @@ compute_nao <- function(data, recipe, obsproj, plot_ts, plot_sp,
                            recipe$Analysis$Variable$name, "_member_",
                            system, 
                            recipe$Analysis$Datasets$Reference$name,
-                           "_ftime", fmonth, ".png")
+                           "_ftime", fmonth, ".pdf")
         caption <-  paste0("NAO method: ",
                            ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
                                   "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
diff --git a/modules/Indices/R/compute_nino.R b/modules/Indices/R/compute_nino.R
index 915dc9cedd826e9733f0b8495dbb7f72ee8edcbb..b47d9ee238a582d894ae1872b9a94bbf4464dc36 100644
--- a/modules/Indices/R/compute_nino.R
+++ b/modules/Indices/R/compute_nino.R
@@ -170,7 +170,7 @@ compute_nino <- function(data, recipe, region, standardised = TRUE,
                            system, "_", recipe$Analysis$Datasets$Reference$name,
                            "_s", recipe$Analysis$Time$sdate, "_ftime",
                            sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min), 
-                           ".png")
+                           ".pdf")
         caption <-  paste0("Nominal start date: 1st of ",
                            month.name[as.numeric(substr(recipe$Analysis$Time$sdate, 1,2))],
                            "\n",
@@ -184,7 +184,7 @@ compute_nino <- function(data, recipe, region, standardised = TRUE,
         plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/", nino_name, "_",
                            system, "_", recipe$Analysis$Datasets$Reference$name,
                            "_ftime",
-                           sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min), ".png")
+                           sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min), ".pdf")
        caption <- paste0("Start date month: ",
                           month.name[get_archive(recipe)$System[[recipe$Analysis$Datasets$System$name]]$initial_month],
                           "\n",
@@ -272,7 +272,7 @@ compute_nino <- function(data, recipe, region, standardised = TRUE,
                            recipe$Analysis$Variable$name, "_",
                            recipe$Analysis$Datasets$Reference$name,
                            "_s", recipe$Analysis$Time$sdate,
-                           "_ftime", fmonth, ".png")
+                           "_ftime", fmonth, ".pdf")
         caption <-  paste0("Nominal start date: 1st of ",
                 month.name[as.numeric(substr(recipe$Analysis$Time$sdate, 1,2))],
                          "\n",
@@ -289,7 +289,7 @@ compute_nino <- function(data, recipe, region, standardised = TRUE,
                            nino_name, "_correlation_",
                            recipe$Analysis$Variable$name, "_",
                            recipe$Analysis$Datasets$Reference$name,
-                           "_ftime", fmonth, ".png")
+                           "_ftime", fmonth, ".pdf")
         caption <-  paste0("Start date: month ",
                            month.name[get_archive(recipe)$System[[recipe$Analysis$Datasets$System$name]]$initial_month],
                            "\n",
@@ -330,7 +330,7 @@ compute_nino <- function(data, recipe, region, standardised = TRUE,
                            recipe$Analysis$Variable$name, "_ensmean_",
                            system,
                            "_s", recipe$Analysis$Time$sdate,
-                           "_ftime", fmonth, ".png")
+                           "_ftime", fmonth, ".pdf")
         caption <-  paste0("Ensemble mean\n",
                            "Nominal start date: 1st of ",
                 month.name[as.numeric(substr(recipe$Analysis$Time$sdate, 1,2))],
@@ -348,7 +348,7 @@ compute_nino <- function(data, recipe, region, standardised = TRUE,
                            recipe$Analysis$Variable$name, "_ensmean_",
                            system,
                            recipe$Analysis$Datasets$Reference$name,
-                           "_ftime", fmonth, ".png")
+                           "_ftime", fmonth, ".pdf")
         caption <-  paste0("Correlation ensemble mean\n",
                            "Start date month: ",
                            month.name[get_archive(recipe)$System[[recipe$Analysis$Datasets$System$name]]$initial_month],
@@ -384,7 +384,7 @@ compute_nino <- function(data, recipe, region, standardised = TRUE,
                            recipe$Analysis$Variable$name, "_member_",
                            system,
                            "_s", recipe$Analysis$Time$sdate,
-                           "_ftime", fmonth, ".png")
+                           "_ftime", fmonth, ".pdf")
         caption <-  paste0("Individual members\n",
                            "Nominal start date: 1st of ",
                 month.name[as.numeric(substr(recipe$Analysis$Time$sdate, 1,2))],
@@ -402,7 +402,7 @@ compute_nino <- function(data, recipe, region, standardised = TRUE,
                          recipe$Analysis$Variable$name, "_ensmean_",
                          system,
                          recipe$Analysis$Datasets$Reference$name,
-                         "_ftime", fmonth, ".png")
+                         "_ftime", fmonth, ".pdf")
       caption <-  paste0("Correlation ensemble mean\n",
                          "Start date month: ",
                          month.name[get_archive(recipe)$System[[recipe$Analysis$Datasets$System$name]]$initial_month],
diff --git a/modules/Visualization/R/plot_ensemble_mean.R b/modules/Visualization/R/plot_ensemble_mean.R
index 3d00742dd4de5443c9cafddd13a7990ffafc1dd4..b39a8f438f6382057424b1e5e6464cab27d89e8a 100644
--- a/modules/Visualization/R/plot_ensemble_mean.R
+++ b/modules/Visualization/R/plot_ensemble_mean.R
@@ -134,7 +134,7 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
                    units = units,
                    cols = cols,
                    brks = brks,
-                   fileout = paste0(outfile, ".png"),
+                   fileout = paste0(outfile, ".pdf"),
                    bar_label_digits = 4,
                    extra_margin = rep(1, 4),
                    bar_label_scale = 1.5,
@@ -204,7 +204,7 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
           }
           # Modify base arguments
           base_args[[1]] <- i_var_ens_mean[i, , ]
-          fileout <- paste0(outfile, "_ft", sprintf("%02d", i), ".png")
+          fileout <- paste0(outfile, "_ft", sprintf("%02d", i), ".pdf")
           base_args$mask <- mask_i
           base_args$dots <- dots_i
           # Plot
diff --git a/modules/Visualization/R/plot_most_likely_terciles_map.R b/modules/Visualization/R/plot_most_likely_terciles_map.R
index 5d60f8c1747c8ec8b5435902a3c44107d9a4873f..b4ea493ed8fe6f5d092877dcc40727c083dc2b67 100644
--- a/modules/Visualization/R/plot_most_likely_terciles_map.R
+++ b/modules/Visualization/R/plot_most_likely_terciles_map.R
@@ -130,7 +130,7 @@ plot_most_likely_terciles <- function(recipe,
                    legend_scale = 0.8, #cex_bar_titles = 0.6,
                    toptitle = toptitle,
                    titles = titles,
-                   fileout = paste0(outfile, ".png"),
+                   fileout = paste0(outfile, ".pdf"),
                    bar_label_digits = 2,
                    bar_scale = rep(0.7, 4),
                    extra_margin = rep(1, 4),
@@ -191,7 +191,7 @@ plot_most_likely_terciles <- function(recipe,
                              format(as.Date(i_syear, format="%Y%m%d"), 
                                             "%d-%m-%Y"))
           # Plot
-          fileout <- paste0(outfile, "_ft", forecast_time, ".png")
+          fileout <- paste0(outfile, "_ft", forecast_time, ".pdf")
           base_args$probs <- i_var_probs[i, , , ]
           base_args$mask <- mask_i
           base_args$dots <- dots_i
diff --git a/modules/Visualization/R/plot_skill_metrics.R b/modules/Visualization/R/plot_skill_metrics.R
index b4c2b273dffb70f3d413b33b9575b0c7e1662bc6..9ef90bf947a2355b2d202c32d02eb56304616912 100644
--- a/modules/Visualization/R/plot_skill_metrics.R
+++ b/modules/Visualization/R/plot_skill_metrics.R
@@ -170,7 +170,7 @@ plot_skill_metrics <- function(recipe, data_cube, skill_metrics,
                      cols = cols,
                      col_inf = col_inf,
                      col_sup = col_sup,
-                     fileout = paste0(outfile, ".png"),
+                     fileout = paste0(outfile, ".pdf"),
                      bar_label_digits = 3,
                      bar_extra_margin = rep(0.9, 4),
                      extra_margin = rep(1, 4),
@@ -241,7 +241,7 @@ plot_skill_metrics <- function(recipe, data_cube, skill_metrics,
                        "Reference: ", recipe$Analysis$Datasets$Reference, "\n",
                        significance_caption)
             }
-            fileout <- paste0(outfile, "_ft", forecast_time, ".png")
+            fileout <- paste0(outfile, "_ft", forecast_time, ".pdf")
             # Plot
 
             do.call(fun,
diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R
index ae94ed3d01526f8cf98932a2adf6f6bce6d7359d..fb12c66a6c3aa4b496bea15bdf704d24538fac5d 100644
--- a/modules/Visualization/Visualization.R
+++ b/modules/Visualization/Visualization.R
@@ -55,6 +55,7 @@ Visualization <- function(recipe,
 
   # Get plot types and create output directories
   plots <- strsplit(recipe$Analysis$Workflow$Visualization$plots, ", | |,")[[1]]
+  ## TODO: Do not modify output dir here
   recipe$Run$output_dir <- paste0(recipe$Run$output_dir, "/plots/")
   outdir <- get_dir(recipe = recipe,
                     variable = data$hcst$attrs$Variable$varName)
@@ -209,5 +210,22 @@ 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/recipes/atomic_recipes/recipe_test_multivar.yml b/recipes/atomic_recipes/recipe_test_multivar.yml
index 01625d40c13a0a3c9cdc6879ff7736c015c36408..86bc59f92bf51da6405e79ed9a9c36e3d418849b 100644
--- a/recipes/atomic_recipes/recipe_test_multivar.yml
+++ b/recipes/atomic_recipes/recipe_test_multivar.yml
@@ -43,11 +43,12 @@ Analysis:
       percentiles: [[1/3, 2/3], [1/10, 9/10]]
       save: 'all'
     Visualization:
-      plots: most_likely_terciles
+      plots: most_likely_terciles, skill_metrics, forecast_ensemble_mean
       multi_panel: no
       projection: cylindrical_equidistant
       dots: both
       mask_terciles: both
+      file_format: PNG
     Indicators:
       index: no
   ncores: 7
diff --git a/tests/testthat/test-decadal_monthly_1.R b/tests/testthat/test-decadal_monthly_1.R
index 3346e529ae76f9e58c08371f95cf0ded95fb6533..c33bd51160e98d921a2e281547396655bed35ef0 100644
--- a/tests/testthat/test-decadal_monthly_1.R
+++ b/tests/testthat/test-decadal_monthly_1.R
@@ -273,8 +273,8 @@ 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.png", "forecast_most_likely_tercile-2021.png",
-  "rpss.png")),
+c("forecast_ensemble_mean-2021.pdf", "forecast_most_likely_tercile-2021.pdf",
+  "rpss.pdf")),
 TRUE
 )
 expect_equal(
diff --git a/tests/testthat/test-decadal_monthly_2.R b/tests/testthat/test-decadal_monthly_2.R
index 9adc16b65deb3b84f9a825d93c175504f6304f2b..02f42c0fe0c12fbe976de0c9cb9d4ce96a6857f4 100644
--- a/tests/testthat/test-decadal_monthly_2.R
+++ b/tests/testthat/test-decadal_monthly_2.R
@@ -267,7 +267,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.png", "enscorr_specs.png", "forecast_ensemble_mean-2020.png", "forecast_ensemble_mean-2021.png", "forecast_most_likely_tercile-2020.png", "forecast_most_likely_tercile-2021.png", "frps_specs.png", "frps.png", "rpss_specs.png")
+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")
 ),
 TRUE
 )
diff --git a/tests/testthat/test-seasonal_NAO.R b/tests/testthat/test-seasonal_NAO.R
index 9582bdfcfdb89ac11b4fa8f3c7dbb79f983fac04..90b446c35ab0a9f8d14290741472f8d80971a517 100644
--- a/tests/testthat/test-seasonal_NAO.R
+++ b/tests/testthat/test-seasonal_NAO.R
@@ -190,10 +190,10 @@ outputs <- paste0(recipe$Run$output_dir, "/plots/")
 expect_equal(
 all(list.files(outputs, recursive = T) %in%
 paste0("Indices/",
-  c("NAO_ECMWF-SEAS5_ERA5_s0301_ftime02.png",
-    "NAO_correlation_psl_ERA5_s0301_ftime02.png",
-    "NAO_correlation_psl_ensmean_ECMWF-SEAS5_s0301_ftime02.png",
-    "NAO_correlation_psl_member_ECMWF-SEAS5_s0301_ftime02.png"))),
+  c("NAO_ECMWF-SEAS5_ERA5_s0301_ftime02.pdf",
+    "NAO_correlation_psl_ERA5_s0301_ftime02.pdf",
+    "NAO_correlation_psl_ensmean_ECMWF-SEAS5_s0301_ftime02.pdf",
+    "NAO_correlation_psl_member_ECMWF-SEAS5_s0301_ftime02.pdf"))),
 TRUE)
 
 expect_equal(
diff --git a/tests/testthat/test-seasonal_monthly.R b/tests/testthat/test-seasonal_monthly.R
index d9b7b3c3857f4d79df127f4be135c25c0b77818c..18bde4c01ec55b2fb29e618e588af74dc167c5f7 100644
--- a/tests/testthat/test-seasonal_monthly.R
+++ b/tests/testthat/test-seasonal_monthly.R
@@ -235,9 +235,9 @@ test_that("5. Visualization", {
 plots <- paste0(recipe$Run$output_dir, "/plots/")
 expect_equal(
 all(basename(list.files(plots, recursive = T)) %in%
-c("crpss-november.png", "enscorr_specs-november.png", "enscorr-november.png",
-  "forecast_ensemble_mean-20201101.png", "forecast_most_likely_tercile-20201101.png",
-  "rpss-november.png")),
+c("crpss-november.pdf", "enscorr_specs-november.pdf", "enscorr-november.pdf",
+  "forecast_ensemble_mean-20201101.pdf", "forecast_most_likely_tercile-20201101.pdf",
+  "rpss-november.pdf")),
 TRUE
 )
 expect_equal(
diff --git a/tools/check_recipe.R b/tools/check_recipe.R
index 01b790b8255daacce1482357cd174ec26479c039..cf469975505390864616a8a5a822035fe35339d3 100644
--- a/tools/check_recipe.R
+++ b/tools/check_recipe.R
@@ -577,6 +577,16 @@ check_recipe <- function(recipe) {
                      "yes/True, no/False, 'both'"))
       }
     }
+    # Check user-defined file format for plots
+    PLOT_FILE_FORMATS <- c("pdf", "png", "jpg", "jpeg", "eps")
+    if (!is.null(recipe$Analysis$Workflow$Visualization$file_format) &&
+        !(tolower(recipe$Analysis$Workflow$Visualization$file_format)) %in%
+        PLOT_FILE_FORMATS) {
+      error(recipe$Run$logger,
+            paste0("Visualization:file_format must be one of the following: ",
+                   paste(PLOT_FILE_FORMATS, collapse = ", "), "."))
+      error_status <- TRUE
+    }
   }
   # Scorecards
   if ("Scorecards" %in% names(recipe$Analysis$Workflow)) {