diff --git a/modules/Loading/testing_recipes/recipe_seasonal-tests.yml b/modules/Loading/testing_recipes/recipe_seasonal-tests.yml index 61177b7162e55b07e8471e59a2f4a493e7cfb42e..4e8128e74342358c40e0bac5dd44f3e27c05a7e8 100644 --- a/modules/Loading/testing_recipes/recipe_seasonal-tests.yml +++ b/modules/Loading/testing_recipes/recipe_seasonal-tests.yml @@ -8,40 +8,40 @@ Analysis: freq: monthly_mean Datasets: System: - name: system7c3s + name: system21_m1 Multimodel: False Reference: name: era5 Time: sdate: '1101' - fcst_year: '2020' + fcst_year: hcst_start: '2000' hcst_end: '2015' ftime_min: 1 - ftime_max: 2 + ftime_max: 6 Region: - latmin: -10 - latmax: 10 - lonmin: 0 - lonmax: 20 + latmin: 30 + latmax: 50 + lonmin: -10 + lonmax: 30 Regrid: method: bilinear type: to_system Workflow: Calibration: - method: mse_min + method: raw Anomalies: compute: yes cross_validation: yes Skill: - metric: RPS RPSS CRPS CRPSS FRPSS BSS10 BSS90 EnsCorr Corr mean_bias mean_bias_SS + metric: mean_bias EnsCorr RPSS CRPSS EnsSprErr # RPS RPSS CRPS CRPSS FRPSS BSS10 BSS90 EnsCorr Corr mean_bias mean_bias_SS Probabilities: percentiles: [[1/3, 2/3], [1/10, 9/10], [1/4, 2/4, 3/4]] Indicators: index: no - ncores: 7 + ncores: 14 remove_NAs: yes - Output_format: S2S4E + Output_format: Scorecards Run: Loglevel: INFO Terminal: yes diff --git a/modules/Loading/testing_recipes/recipe_testing_nadia.yml b/modules/Loading/testing_recipes/recipe_testing_nadia.yml new file mode 100644 index 0000000000000000000000000000000000000000..6071198148b887a99598a265b9ddd62b2f007565 --- /dev/null +++ b/modules/Loading/testing_recipes/recipe_testing_nadia.yml @@ -0,0 +1,49 @@ +Description: + Author: V. Agudetse + +Analysis: + Horizon: Seasonal + Variables: + name: tas + freq: monthly_mean + Datasets: + System: + name: system5c3s + Multimodel: False + Reference: + name: era5 + Time: + sdate: '1101' + fcst_year: + hcst_start: '2010' + hcst_end: '2015' + ftime_min: 1 + ftime_max: 6 + Region: + latmin: 30 + latmax: 50 + lonmin: -10 + lonmax: 30 + Regrid: + method: bilinear + type: to_system + Workflow: + Calibration: + method: raw + Anomalies: + compute: yes + cross_validation: yes + Skill: + metric: mean_bias EnsCorr RPSS CRPSS EnsSprErr + Probabilities: + percentiles: [[1/3, 2/3], [1/10, 9/10], [1/4, 2/4, 3/4]] + Indicators: + index: no + ncores: 7 + remove_NAs: yes + Output_format: scorecards +Run: + Loglevel: INFO + Terminal: yes + output_dir: /esarchive/scratch/vagudets/repos/auto-s2s/out-logs/ + code_dir: /esarchive/scratch/vagudets/repos/auto-s2s/ diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R index 7de55fb334fc1defe2de6cfc4a13a775ab016534..a8664569047bb60185733e27a1250ec985cd481d 100644 --- a/modules/Visualization/Visualization.R +++ b/modules/Visualization/Visualization.R @@ -1,12 +1,11 @@ #G# TODO: Remove once released in s2dv/CSTools source("modules/Visualization/tmp/PlotMostLikelyQuantileMap.R") source("modules/Visualization/tmp/PlotCombinedMap.R") +source("modules/Visualization/tmp/clim.palette.R") ## TODO: Add the possibility to read the data directly from netCDF ## TODO: Adapt to multi-model case ## TODO: Add param 'raw'? -## TODO: Reduce colorbar size and increase colorbar label size -## Param: bar_label_scale and ???? plot_data <- function(recipe, data, @@ -63,7 +62,14 @@ plot_data <- function(recipe, plot_skill_metrics <- function(recipe, archive, data_cube, skill_metrics, outdir, significance = F) { - + # recipe: Auto-S2S recipe + # archive: Auto-S2S archive + # data_cube: s2dv_cube object with the corresponding hindcast data + # skill_metrics: list of named skill metrics arrays + # outdir: output directory + # significance: T/F, whether to display the significance dots in the plots + + ## TODO: OPTION for CERISE: Using PuOr # Abort if frequency is daily if (recipe$Analysis$Variables$freq == "daily_mean") { error(recipe$Run$logger, "Visualization functions not yet implemented @@ -83,15 +89,24 @@ plot_skill_metrics <- function(recipe, archive, data_cube, skill_metrics, init_month <- lubridate::month(as.numeric(substr(recipe$Analysis$Time$sdate, start = 1, stop = 2)), label = T, abb = T) + # Define color palette and number of breaks according to output format + ## TODO: Make separate function + if (tolower(recipe$Analysis$Output_format) %in% c("scorecards", "cerise")) { + diverging_palette <- "purpleorange" + sequential_palette <- "Oranges" + } else { + diverging_palette <- "bluered" + sequential_palette <- "Reds" + } # Group different metrics by type skill_scores <- c("rpss", "bss90", "bss10", "frpss", "crpss", "mean_bias_ss", "enscorr", "rpss_specs", "bss90_specs", "bss10_specs", "enscorr_specs", "rmsss") scores <- c("rps", "frps", "crps", "frps_specs") - + # Assign colorbar to each metric type + ## TODO: Triangle ends for (name in c(skill_scores, scores, "mean_bias", "enssprerr")) { - if (name %in% names(skill_metrics)) { # Define plot characteristics and metric name to display in plot if (name %in% c("rpss", "bss90", "bss10", "frpss", "crpss", @@ -99,96 +114,104 @@ plot_skill_metrics <- function(recipe, archive, data_cube, skill_metrics, "rmsss")) { display_name <- toupper(strsplit(name, "_")[[1]][1]) skill <- skill_metrics[[name]] - brks <- seq(-1, 1, by = 0.1) - col2 <- grDevices::hcl.colors(length(brks) - 1, "RdBu", rev = TRUE) - + brks <- seq(-1, 1, by = 0.2) + colorbar <- clim.colors(length(brks) + 1, diverging_palette) + cols <- colorbar[2:(length(colorbar) - 1)] + col_inf <- colorbar[1] + col_sup <- NULL } else if (name == "mean_bias_ss") { display_name <- "Mean Bias Skill Score" skill <- skill_metrics[[name]] - brks <- seq(-1, 1, by = 0.1) - col2 <- grDevices::hcl.colors(length(brks) - 1, "RdBu", rev = TRUE) - + brks <- seq(-1, 1, by = 0.2) + colorbar <- clim.colors(length(brks) + 1, diverging_palette) + cols <- colorbar[2:(length(colorbar) - 1)] + col_inf <- colorbar[1] + col_sup <- NULL } else if (name %in% c("enscorr", "enscorr_specs")) { display_name <- "Ensemble Mean Correlation" skill <- skill_metrics[[name]] - brks <- seq(-1, 1, by = 0.1) - col2 <- grDevices::hcl.colors(length(brks) - 1, "RdBu", rev = TRUE) - + brks <- seq(-1, 1, by = 0.2) + cols <- clim.colors(length(brks) - 1, diverging_palette) + col_inf <- NULL + col_sup <- NULL } else if (name %in% scores) { skill <- skill_metrics[[name]] display_name <- toupper(strsplit(name, "_")[[1]][1]) brks <- seq(0, 1, by = 0.1) - col2 <- grDevices::hcl.colors(length(brks) - 1, "Reds") - + colorbar <- grDevices::hcl.colors(length(brks), sequential_palette) + cols <- colorbar[1:(length(colorbar) - 1)] + col_inf <- NULL + col_sup <- colorbar[length(colorbar)] } else if (name == "enssprerr") { - ## TODO: Adjust colorbar parameters skill <- skill_metrics[[name]] display_name <- "Spread-to-Error Ratio" - brks <- pretty(0:max(skill, na.rm = T), n = 20, min.n = 10) - col2 <- grDevices::hcl.colors(length(brks) - 1, "RdBu", rev = TRUE) - + brks <- c(0, 0.6, 0.7, 0.8, 0.9, 1, 1.2, 1.4, 1.6, 1.8, 2) + colorbar <- clim.colors(length(brks), diverging_palette) + cols <- colorbar[1:length(colorbar) - 1] + col_inf <- NULL + col_sup <- colorbar[length(colorbar)] } else if (name == "mean_bias") { skill <- skill_metrics[[name]] display_name <- "Mean Bias" - max_value <- max(abs(skill)) - ugly_intervals <- seq(-max_value, max_value, (max_value*2)/10) - brks <- pretty(ugly_intervals, n = 20, min.n = 10) - col2 <- grDevices::hcl.colors(length(brks) - 1, "RdBu", rev = TRUE) + max_value <- max(abs(quantile(skill, 0.02, na.rm = T)), + abs(quantile(skill, 0.98, na.rm = T))) + brks <- max_value * seq(-1, 1, by = 0.2) + colorbar <- clim.colors(length(brks) + 1, diverging_palette) + cols <- colorbar[2:(length(colorbar) - 1)] + col_inf <- colorbar[1] + col_sup <- colorbar[length(colorbar)] } - - options(bitmapType = "cairo") - - # Reorder dimensions - skill <- Reorder(skill, 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(skill_metrics))) { + options(bitmapType = "cairo") + # Reorder dimensions + skill <- Reorder(skill, 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") - skill_significance <- skill_metrics[[significance_name]] - skill_significance <- Reorder(skill_significance, c("time", - "longitude", - "latitude")) - # Split skill significance into list of lists, along the time dimension - # This allows for plotting the significance dots correctly. - skill_significance <- ClimProjDiags::ArrayToList(skill_significance, - dim = 'time', - level = "sublist", - names = "dots") - } else { - skill_significance <- NULL - } - # Define output file name and titles - outfile <- paste0(outdir, name, ".png") - toptitle <- paste(display_name, "-", data_cube$Variable$varName, - "-", system_name, "-", init_month, hcst_period) - months <- unique(lubridate::month(data_cube$Dates$start, - label = T, abb = F)) - titles <- as.vector(months) - # Plot - suppressWarnings( - PlotLayout(PlotEquiMap, c('longitude', 'latitude'), - asplit(skill, MARGIN=1), # Splitting array into a list - longitude, latitude, - special_args = skill_significance, - dot_symbol = 20, - toptitle = toptitle, - title_scale = 0.6, - titles = titles, - filled.continents=F, - brks = brks, - cols = col2, - col_inf = col2[1], - col_sup = col2[length(col2)], - fileout = outfile, - bar_label_digits = 3, - bar_extra_margin = rep(0.9, 4), - bar_label_scale = 1.5, - axes_label_scale = 1.3) - ) - } + if ((significance) && (significance_name %in% names(skill_metrics))) { + significance_name <- paste0(name, "_significance") + skill_significance <- skill_metrics[[significance_name]] + skill_significance <- Reorder(skill_significance, c("time", + "longitude", + "latitude")) + # Split skill significance into list of lists, along the time dimension + # to avoid overlapping of significance dots. + skill_significance <- ClimProjDiags::ArrayToList(skill_significance, + dim = 'time', + level = "sublist", + names = "dots") + } else { + skill_significance <- NULL + } + # Define output file name and titles + outfile <- paste0(outdir, name, ".png") + toptitle <- paste(display_name, "-", data_cube$Variable$varName, + "-", system_name, "-", init_month, hcst_period) + months <- unique(lubridate::month(data_cube$Dates$start, + label = T, abb = F)) + titles <- as.vector(months) + # Plot + suppressWarnings( + PlotLayout(PlotEquiMap, c('longitude', 'latitude'), + asplit(skill, MARGIN=1), # Splitting array into a list + longitude, latitude, + special_args = skill_significance, + dot_symbol = 20, + toptitle = toptitle, + title_scale = 0.6, + titles = titles, + filled.continents=F, + brks = brks, + cols = cols, + col_inf = col_inf, + col_sup = col_sup, + fileout = outfile, + bar_label_digits = 3, + bar_extra_margin = rep(0.9, 4), + bar_label_scale = 1.5, + axes_label_scale = 1.3) + ) + } } - info(recipe$Run$logger, "##### SKILL METRIC PLOTS SAVED TO OUTPUT DIRECTORY #####") } @@ -208,7 +231,6 @@ plot_ensemble_mean <- function(recipe, archive, fcst, outdir) { units <- attr(fcst$Variable, "variable")$units start_date <- paste0(recipe$Analysis$Time$fcst_year, recipe$Analysis$Time$sdate) - # Compute ensemble mean ensemble_mean <- s2dv::MeanDims(fcst$data, 'ensemble') # Drop extra dims, add time dim if missing: @@ -218,9 +240,14 @@ plot_ensemble_mean <- function(recipe, archive, fcst, outdir) { dim(ensemble_mean) <- c("time" = 1, dim(ensemble_mean)) } if (!'syear' %in% names(dim(ensemble_mean))) { - ensemble_mean <- Reorder(ensemble_mean, c("time", "longitude", "latitude")) + ensemble_mean <- Reorder(ensemble_mean, c("time", + "longitude", + "latitude")) } else { - ensemble_mean <- Reorder(ensemble_mean, c("syear", "time", "longitude", "latitude")) + ensemble_mean <- Reorder(ensemble_mean, c("syear", + "time", + "longitude", + "latitude")) } ## TODO: Redefine column colors, possibly depending on variable if (variable == 'prlr') { @@ -230,10 +257,18 @@ plot_ensemble_mean <- function(recipe, archive, fcst, outdir) { palette = "RdBu" rev = T } - - brks <- pretty(range(ensemble_mean, na.rm = T), n = 15, min.n = 8) - col2 <- grDevices::hcl.colors(length(brks) - 1, palette, rev = rev) - # color <- colorRampPalette(col2)(length(brks) - 1) + # Define brks, centered on in the case of anomalies + ## + if (grepl("anomaly", + attr(fcst$Variable, "variable")$long_name)) { + variable <- paste(variable, "anomaly") + max_value <- max(abs(ensemble_mean)) + ugly_intervals <- seq(-max_value, max_value, max_value/20) + brks <- pretty(ugly_intervals, n = 12, min.n = 8) + } else { + brks <- pretty(range(ensemble_mean, na.rm = T), n = 15, min.n = 8) + } + cols <- grDevices::hcl.colors(length(brks) - 1, palette, rev = rev) options(bitmapType = "cairo") for (i_syear in start_date) { @@ -258,7 +293,7 @@ plot_ensemble_mean <- function(recipe, archive, fcst, outdir) { title_scale = 0.6, titles = titles, units = units, - cols = col2, + cols = cols, brks = brks, fileout = outfile, bar_label_digits = 4, @@ -266,7 +301,6 @@ plot_ensemble_mean <- function(recipe, archive, fcst, outdir) { bar_label_scale = 1.5, axes_label_scale = 1.3) } - info(recipe$Run$logger, "##### FCST ENSEMBLE MEAN PLOT SAVED TO OUTPUT DIRECTORY #####") } diff --git a/modules/Visualization/tmp/clim.palette.R b/modules/Visualization/tmp/clim.palette.R new file mode 100644 index 0000000000000000000000000000000000000000..b23ff8428f201dbe4cb129e5f202ab2f1924532f --- /dev/null +++ b/modules/Visualization/tmp/clim.palette.R @@ -0,0 +1,69 @@ +#'Generate Climate Color Palettes +#' +#'Generates a colorblind friendly color palette with color ranges useful in +#'climate temperature variable plotting. +#' +#'@param palette Which type of palette to generate: from blue through white +#' to red ('bluered'), from red through white to blue ('redblue'), from +#' yellow through orange to red ('yellowred'), from red through orange to +#' red ('redyellow'), from purple through white to orange ('purpleorange'), +#' and from orange through white to purple ('orangepurple'). +#'@param n Number of colors to generate. +#' +#'@examples +#'lims <- seq(-1, 1, length.out = 21) +#' +#'ColorBar(lims, color_fun = clim.palette('redyellow')) +#' +#'cols <- clim.colors(20) +#'ColorBar(lims, cols) +#' +#'@rdname clim.palette +#'@importFrom grDevices colorRampPalette +#'@export +clim.palette <- 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 clim.palette +#'@export +clim.colors <- function(n, palette = "bluered") { + clim.palette(palette)(n) +}