From 52a478d2f452d1527fa5193a6631a348e4876b4e Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 28 May 2024 11:30:47 +0200 Subject: [PATCH] Mask and/or dots for metrics visualisation --- example_scripts/exec_units.R | 65 +++++++++++++++++-- modules/Visualization/R/plot_metrics.R | 55 ++++++++++++---- modules/Visualization/R/tmp/PlotRobinson.R | 8 ++- modules/Visualization/Visualization.R | 23 ++++++- .../examples/recipe_tas_seasonal_units.yml | 4 +- 5 files changed, 128 insertions(+), 27 deletions(-) diff --git a/example_scripts/exec_units.R b/example_scripts/exec_units.R index 819121c9..f4085c11 100644 --- a/example_scripts/exec_units.R +++ b/example_scripts/exec_units.R @@ -1,6 +1,5 @@ rm(list=ls()) gc() -setwd("/esarchive/scratch/nperez/git/auto-s2s") source("modules/Loading/Loading.R") source("modules/Calibration/Calibration.R") @@ -21,18 +20,72 @@ source("tools/prepare_outputs.R") recipe <- prepare_outputs(recipe_file) # Load datasets -data <- load_datasets(recipe) +data <- Loading(recipe) # Units transformation source("modules/Units/Units.R") test <- Units(recipe, data) # Calibrate datasets -data <- calibrate_datasets(recipe, test) +data <- Calibration(recipe, test) # Compute skill metrics -skill_metrics <- compute_skill_metrics(recipe, data) +skill_metrics <- Skill(recipe, data) # Compute percentiles and probability bins -probabilities <- compute_probabilities(recipe, data) +probabilities <- Probabilities(recipe, data) # Export all data to netCDF ## TODO: Fix plotting # save_data(recipe, data, skill_metrics, probabilities) # Plot data -plot_data(recipe, data, skill_metrics, probabilities, significance = T) +Visualization(recipe = recipe, data = data, skill_metrics = skill_metrics, + probabilities = probabilities, significance = F) +list.files(paste0(recipe$Run$output_dir, "/plots/ECMWF-SEAS51/ERA5/evmos/tas/")) +# signif TRUE: +Visualization(recipe = recipe, data = data, skill_metrics = skill_metrics, + probabilities = probabilities, significance = T) +list.files(paste0(recipe$Run$output_dir, "/plots/ECMWF-SEAS51/ERA5/evmos/tas/")) +# signif mask: this option is not available in PlotEquiMap, it return dots +Visualization(recipe = recipe, data = data, skill_metrics = skill_metrics, + probabilities = probabilities, significance = 'mask') +list.files(paste0(recipe$Run$output_dir, "/plots/ECMWF-SEAS51/ERA5/evmos/tas/")) +## The difference between TRUE and mask for PlotEquiMap is the file name + +# signif dots: +visualization(recipe = recipe, data = data, skill_metrics = skill_metrics, + probabilities = probabilities, significance = 'dots') +list.files(paste0(recipe$run$output_dir, "/plots/ecmwf-seas51/era5/evmos/tas/")) + +# In this case duplicates files but simplifying the code makes it very complex +Visualization(recipe = recipe, data = data, skill_metrics = skill_metrics, + probabilities = probabilities, significance = 'both') + +recipe$Analysis$Workflow$Visualization$multi_panel <- TRUE +recipe$Analysis$Workflow$Visualization$mask_terciles <- NULL +Visualization(recipe = recipe, data = data, skill_metrics = skill_metrics, + probabilities = probabilities, significance = F) + +Error in PlotCombinedMap(probs * 100, lon, lat, map_select_fun = max, : + Parameter 'mask' must have two dimensions. + +Visualization(recipe = recipe, data = data, skill_metrics = skill_metrics, + probabilities = probabilities, significance = T) + +Visualization(recipe = recipe, data = data, skill_metrics = skill_metrics, + probabilities = probabilities, significance = 'mask') + +recipe$Analysis$Workflow$Visualization$multi_panel <- FALSE +recipe$Analysis$Workflow$Visualization$projection <- 'Robinson' +Visualization(recipe = recipe, data = data, skill_metrics = skill_metrics, + probabilities = probabilities, significance = F) +Visualization(recipe = recipe, data = data, skill_metrics = skill_metrics, + probabilities = probabilities, significance = T) +Visualization(recipe = recipe, data = data, skill_metrics = skill_metrics, + probabilities = probabilities, significance = 'dots') +Visualization(recipe = recipe, data = data, skill_metrics = skill_metrics, + probabilities = probabilities, significance = 'mask') +Visualization(recipe = recipe, data = data, skill_metrics = skill_metrics, + probabilities = probabilities, significance = 'both') + + + + + + + diff --git a/modules/Visualization/R/plot_metrics.R b/modules/Visualization/R/plot_metrics.R index 2df06e8b..6976e6ac 100644 --- a/modules/Visualization/R/plot_metrics.R +++ b/modules/Visualization/R/plot_metrics.R @@ -175,19 +175,23 @@ plot_metrics <- function(recipe, data_cube, metrics, # Reorder dimensions metric <- Reorder(metric, c("time", "longitude", "latitude")) # If the significance has been requested and the variable has it, - # retrieve it and reorder its dimensions. - significance_name <- paste0(name, "_significance") - if ((significance) && (significance_name %in% names(metrics))) { - metric_significance <- var_metric[[significance_name]] - metric_significance <- Reorder(metric_significance, c("time", - "longitude", - "latitude")) - # Split significance into list of lists, along the time dimension - # This allows for plotting the significance dots correctly. - metric_significance <- ClimProjDiags::ArrayToList(metric_significance, + # retrieve it and reorder its dimensions. + 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]] + metric_significance <- Reorder(metric_significance, c("time", + "longitude", + "latitude")) + # Split significance into list of lists, along the time dimension + # This allows for plotting the significance dots correctly. + metric_significance <- ClimProjDiags::ArrayToList(metric_significance, dim = "time", level = "sublist", names = "dots") + } else { + metric_significance <- NULL + } } else { metric_significance <- NULL } @@ -274,9 +278,10 @@ plot_metrics <- function(recipe, data_cube, metrics, lon = longitude, lat = latitude, lon_dim = 'longitude', lat_dim = 'latitude', target_proj = target_proj, legend = 's2dv', - style = 'point', brks = brks, cols = cols, + style = 'point', dots = NULL, brks = brks, + cols = cols, col_inf = col_inf, col_sup = col_sup, - units = units) + units = units) } # Loop over forecast times for (i in 1:dim(metric)[['time']]) { @@ -324,10 +329,31 @@ plot_metrics <- function(recipe, data_cube, metrics, # Modify base arguments base_args[[1]] <- metric[i, , ] if (!is.null(metric_significance)) { - base_args[[2]] <- metric_significance[[i]][[1]] + sign_file_label <- NULL + if (is.logical(significance)) { + if (significance) { + base_args[[2]] <- metric_significance[[i]][[1]] + sign_file_lable <- '_mask' + } + } else { + if (significance == 'dots') { + if (projection != 'cylindrical_equidistant') { + base_args[[10]] <- metric_significance[[i]][[1]] + } else { + # The position of arguments in base_args requires this cond + # so PlotEquiMap plots dots when requested + base_args[[2]] <- metric_significance[[i]][[1]] + } + sign_file_label <- '_dots' + } else if (significance == 'mask') { + base_args[[2]] <- metric_significance[[i]][[1]] + sign_file_label <- '_mask' + } + } significance_caption <- "alpha = 0.05" } else { significance_caption <- NULL + sign_file_label <- NULL } if (identical(fun, PlotRobinson)) { ## TODO: Customize alpha and other technical details depending on the metric @@ -338,7 +364,8 @@ plot_metrics <- function(recipe, data_cube, metrics, "Reference: ", recipe$Analysis$Datasets$Reference, "\n", significance_caption) } - fileout <- paste0(outfile, "_ft", forecast_time, ".pdf") + fileout <- paste0(outfile, "_ft", forecast_time, + sign_file_label, ".pdf") # Plot info(recipe$Run$logger, paste("Plotting", display_name)) diff --git a/modules/Visualization/R/tmp/PlotRobinson.R b/modules/Visualization/R/tmp/PlotRobinson.R index bd427448..67ca034d 100644 --- a/modules/Visualization/R/tmp/PlotRobinson.R +++ b/modules/Visualization/R/tmp/PlotRobinson.R @@ -118,11 +118,11 @@ PlotRobinson <- function(data, lon, lat, lon_dim = NULL, lat_dim = NULL, target_proj = 54030, legend = 's2dv', 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 = clim.palette(), bar_extra_margin = rep(0, 4), vertical = TRUE, + color_fun = clim.palette(), bar_extra_margin = c(3.5, 0, 3.5, 0), vertical = TRUE, toptitle = NULL, caption = NULL, units = NULL, crop_coastlines = NULL, - point_size = "auto", title_size = 16, dots_size = 0.5, + point_size = "auto", title_size = 10, dots_size = 0.2, dots_shape = 47, coastlines_width = 0.3, - fileout = NULL, width = 8, height = 4, size_units = "in", + fileout = NULL, width = 8, height = 5, size_units = "in", res = 300) { # Sanity check @@ -234,12 +234,14 @@ PlotRobinson <- function(data, lon, lat, lon_dim = NULL, lat_dim = NULL, 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.") } diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R index 3750baea..b124c311 100644 --- a/modules/Visualization/Visualization.R +++ b/modules/Visualization/Visualization.R @@ -77,8 +77,27 @@ Visualization <- function(recipe, # Plot skill metrics if ("skill_metrics" %in% plots) { if (!is.null(skill_metrics)) { - plot_metrics(recipe, data$hcst, skill_metrics, outdir, - significance, output_conf = output_conf) + if (is.logical(significance)) { + plot_metrics(recipe, data$hcst, skill_metrics, outdir, + significance, output_conf = output_conf) + 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 { error(recipe$Run$logger, paste0("The skill metric plots have been requested, but the ", diff --git a/recipes/examples/recipe_tas_seasonal_units.yml b/recipes/examples/recipe_tas_seasonal_units.yml index d2b25321..17df1efd 100644 --- a/recipes/examples/recipe_tas_seasonal_units.yml +++ b/recipes/examples/recipe_tas_seasonal_units.yml @@ -59,5 +59,5 @@ Analysis: Run: Loglevel: INFO Terminal: TRUE - output_dir: /esarchive/scratch/nperez/cs_oper/ - code_dir: /esarchive/scratch/nperez/git/s2s-suite/ + output_dir: /esarchive/scratch/nperez/ + code_dir: /esarchive/scratch/nperez/git6/sunset/ -- GitLab