diff --git a/conf/archive.yml b/conf/archive.yml index a0caf86aa4f19b0ec1f1cf9938a6b2657d758c76..a347dc4042def2932b037fd9cfaf8d1b96613d09 100644 --- a/conf/archive.yml +++ b/conf/archive.yml @@ -187,7 +187,7 @@ esarchive: src: "recon/ecmwf/cerra/" daily_mean: {"hur":"_f3h-r2631x1113/", "ps":"_f3h-r2631x1113/", "sfcWind":"_f3h-r2631x1113/", "tas":"_f3h-r2631x1113/", "winddir":"_f3h-r2631x1113/"} - monthly_mean: {"hur":"_f3h-r2631x1113/", "ps":"_f3h-r2631x1113/", "sfcWind":"_f3h-r2631x1113/", + monthly_mean: {"hurs":"_f3h-r2631x1113/", "ps":"_f3h-r2631x1113/", "sfcWind":"_f3h-r2631x1113/", "tas":"_f3h-r2631x1113/", "winddir":"_f3h-r2631x1113/","tasmin":"_f24h-r2631x1113/","tasmax":"_f24h-r2631x1113/"} calendar: "proleptic_gregorian" reference_grid: "/esarchive/recon/ecmwf/cerra/monthly_mean/tas_f3h-r2631x1113/tas_200506.nc" diff --git a/exec_ecvs_seasonal_oper.R b/exec_ecvs_seasonal_oper.R new file mode 100644 index 0000000000000000000000000000000000000000..18f2e4938cd8d48734782717b028c4775c891316 --- /dev/null +++ b/exec_ecvs_seasonal_oper.R @@ -0,0 +1,35 @@ +rm(list=ls()) +gc() +setwd("/esarchive/scratch/nperez/git/auto-s2s") + +source("modules/Loading/Loading.R") +source("modules/Calibration/Calibration.R") +source("modules/Anomalies/Anomalies.R") +source("modules/Skill/Skill.R") +source("modules/Saving/Saving.R") +source("modules/Visualization/Visualization.R") +source("tools/prepare_outputs.R") + +# Read recipe +args = commandArgs(trailingOnly = TRUE) +recipe_file <- args[1] +recipe <- read_atomic_recipe(recipe_file) +## to test a single recipe: +#recipe_file <- "recipe_ecvs_seasonal_oper.yml" +#recipe <- prepare_outputs(recipe_file) + +# Load datasets +data <- load_datasets(recipe) +# Calibrate datasets +data <- calibrate_datasets(recipe, data) +# Compute skill metrics +skill_metrics <- compute_skill_metrics(recipe, data) +# Compute percentiles and probability bins +probabilities <- compute_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) + + diff --git a/modules/Visualization/R/plot_most_likely_terciles_map.R b/modules/Visualization/R/plot_most_likely_terciles_map.R index 1866770161f4c2a6b84583e554d509000ac5c10f..4952a90eef66477c23826b7badcce7837a3a3972 100644 --- a/modules/Visualization/R/plot_most_likely_terciles_map.R +++ b/modules/Visualization/R/plot_most_likely_terciles_map.R @@ -1,7 +1,14 @@ +source("modules/Visualization/R/tmp/PlotMostLikelyQuantileMap.R") +source("modules/Visualization/R/tmp/PlotCombinedMap.R") +source("modules/Visualization/R/tmp/GradientCatsColorBar.R") + + ## TODO: Change name plot_most_likely_terciles <- function(recipe, fcst, - probabilities, + probabilities, + mask, + dots, outdir) { ## TODO: Add 'anomaly' to plot title @@ -39,17 +46,52 @@ plot_most_likely_terciles <- function(recipe, for (var in 1:fcst$dims[['var']]) { variable <- fcst$attrs$Variable$varName[[var]] var_long_name <- fcst$attrs$Variable$metadata[[variable]]$long_name + # Choose colors depending on the variable + if (variable %in% c('prlr')) { ## add others + cols <- list(c("#FFC473", #FFAB38", + "darkorange1"), + c("#b5b5b5", "black"), + c("#A0E5E4",#41CBC9", + "deepskyblue3")) + } else { + cols <- list(c("#A0E5E4",#33BFD1", + "deepskyblue3"), + c("#b5b5b5", "black"), + c("#FFB19A", #FF764D", + "indianred3")) + } var_probs <- ClimProjDiags::Subset(probs_fcst, along = c("var"), indices = var, drop = 'selected') - var_probs <- Reorder(var_probs, c("syear", "time", "bin", "longitude", "latitude")) for (i_syear in start_date) { # Define name of output file and titles i_var_probs <- var_probs[which(start_date == i_syear), , , , ] outfile <- paste0(outdir[[var]], "forecast_most_likely_tercile-", i_syear) + # Mask + if (!is.null(mask)) { + outfile <- paste0(outfile, "_rpssmask") + var_mask <- ClimProjDiags::Subset(mask, + along = c("var"), + indices = var, + drop = 'selected') + dim_mask <- dim(var_mask) + var_mask <- as.numeric(var_mask <= 0) + dim(var_mask) <- dim_mask + } + # Dots + if (!is.null(dots)) { + outfile <- paste0(outfile, "_rpssdots") + var_dots <- ClimProjDiags::Subset(dots, + along = c("var"), + indices = var, + drop = 'selected') + dim_dots <- dim(var_dots) + var_dots <- as.numeric(var_dots <= 0) + dim(var_dots) <- dim_dots + } toptitle <- paste0(system_name, " / ", str_to_title(var_long_name), "\n", "Most Likely Tercile / Initialization: ", i_syear) @@ -66,6 +108,8 @@ plot_most_likely_terciles <- function(recipe, PlotLayout(PlotMostLikelyQuantileMap, c('bin', 'longitude', 'latitude'), cat_dim = 'bin', i_var_probs, longitude, latitude, + mask = mask, + dots = dots, coast_width = 1.5, title_scale = 0.6, title_margin_scale = 0.7, @@ -89,18 +133,36 @@ plot_most_likely_terciles <- function(recipe, forecast_time <- forecast_time + 12 } forecast_time <- sprintf("%02d", forecast_time) + # Get mask subset + if (!is.null(mask)) { + mask_i <- Subset(var_mask, along = 'time', indices = i, drop = TRUE) + } else { + mask_i <- NULL + } + # Get dots subset + if (!is.null(dots)) { + dots_i <- Subset(var_dots, along = 'time', indices = i, drop = TRUE) + } else { + dots_i <- NULL + } # Define plot title toptitle <- paste0(system_name, " / ", str_to_title(var_long_name), "\n", "Most Likely Tercile / ", months[i], " ", years[i], - " (Init.: ", i_syear, ")") + " / Start date: ", + format(as.Date(i_syear, format="%Y%m%d"), + "%d-%m-%Y")) # Plot fileout <- paste0(outfile, "_ft", forecast_time, ".png") PlotMostLikelyQuantileMap(cat_dim = 'bin', probs = i_var_probs[i, , , ], lon = longitude, lat = latitude, coast_width = 1.5, + mask = mask_i, + dots = dots_i, + col_mask = 'antiquewhite', + cols = cols, title_scale = 1, legend_scale = 0.8, cex_bar_titles = 0.9, diff --git a/modules/Visualization/R/tmp/GradientCatsColorBar.R b/modules/Visualization/R/tmp/GradientCatsColorBar.R new file mode 100644 index 0000000000000000000000000000000000000000..4cdeb0afe85ce47601228fa804ffc75011a2f249 --- /dev/null +++ b/modules/Visualization/R/tmp/GradientCatsColorBar.R @@ -0,0 +1,91 @@ +#Draws Color Bars for Categories +#A wrapper of s2dv::ColorBar to generate multiple color bars for different +#categories, and each category has different color set. +GradientCatsColorBar <- function(nmap, brks = NULL, cols = NULL, vertical = TRUE, subsampleg = NULL, + bar_limits, var_limits = NULL, + triangle_ends = NULL, plot = TRUE, + draw_separators = FALSE, + bar_titles = NULL, title_scale = 1, label_scale = 1, extra_margin = rep(0, 4), + ...) { + # bar_limits + if (!is.numeric(bar_limits) || length(bar_limits) != 2) { + stop("Parameter 'bar_limits' must be a numeric vector of length 2.") + } + + # Check brks + if (is.null(brks) || (is.numeric(brks) && length(brks) == 1)) { + num_brks <- 5 + if (is.numeric(brks)) { + num_brks <- brks + } + brks <- seq(from = bar_limits[1], to = bar_limits[2], length.out = num_brks) + } + if (!is.numeric(brks)) { + stop("Parameter 'brks' must be a numeric vector.") + } + # Check cols + col_sets <- list(c("#A1D99B", "#74C476", "#41AB5D", "#238B45"), + c("#6BAED6FF", "#4292C6FF", "#2171B5FF", "#08519CFF"), + c("#FFEDA0FF", "#FED976FF", "#FEB24CFF", "#FD8D3CFF"), + c("#FC4E2AFF", "#E31A1CFF", "#BD0026FF", "#800026FF"), + c("#FCC5C0", "#FA9FB5", "#F768A1", "#DD3497")) + if (is.null(cols)) { + if (length(col_sets) >= nmap) { + chosen_sets <- 1:nmap + chosen_sets <- chosen_sets + floor((length(col_sets) - length(chosen_sets)) / 2) + } else { + chosen_sets <- array(1:length(col_sets), nmap) + } + cols <- col_sets[chosen_sets] + } else { + if (!is.list(cols)) { + stop("Parameter 'cols' must be a list of character vectors.") + } + if (!all(sapply(cols, is.character))) { + stop("Parameter 'cols' must be a list of character vectors.") + } + if (length(cols) != nmap) { + stop("Parameter 'cols' must be a list of the same length as the number of ", + "maps in 'maps'.") + } + } + for (i in 1:length(cols)) { + if (length(cols[[i]]) != (length(brks) - 1)) { + cols[[i]] <- grDevices::colorRampPalette(cols[[i]])(length(brks) - 1) + } + } + + # Check bar_titles + if (is.null(bar_titles)) { + if (nmap == 3) { + bar_titles <- c("Below normal (%)", "Normal (%)", "Above normal (%)") + } else if (nmap == 5) { + bar_titles <- c("Low (%)", "Below normal (%)", + "Normal (%)", "Above normal (%)", "High (%)") + } else { + bar_titles <- paste0("Cat. ", 1:nmap, " (%)") + } + } + + if (plot) { +brks <- brks[-length(brks)] + for (k in 1:nmap) { +triangle_ends <- c(F,T) +col_sup <- cols[[k]][length(cols[[k]])] +cols[[k]] <- cols[[k]][-length(cols[[k]])] + s2dv::ColorBar(brks = brks, cols = cols[[k]], vertical = FALSE, + subsampleg = subsampleg, +# bar_limits = bar_limits, var_limits = var_limits, + triangle_ends = triangle_ends, plot = TRUE, + col_sup = col_sup, + draw_separators = draw_separators, + title = bar_titles[[k]], title_scale = title_scale, + label_scale = label_scale, extra_margin = extra_margin) + } + } else { + #TODO: col_inf and col_sup + return(list(brks = brks, cols = cols)) + } + +} + diff --git a/modules/Visualization/R/tmp/PlotCombinedMap.R b/modules/Visualization/R/tmp/PlotCombinedMap.R new file mode 100644 index 0000000000000000000000000000000000000000..4575f8973b638e0abeec5da0d4def5404a5e2452 --- /dev/null +++ b/modules/Visualization/R/tmp/PlotCombinedMap.R @@ -0,0 +1,445 @@ +#'Plot Multiple Lon-Lat Variables In a Single Map According to a Decision Function +#' +#'@description Plot a number a two dimensional matrices with (longitude, +#'latitude) dimensions on a single map with the cylindrical equidistant +#'latitude and longitude projection. +#' +#'@author Nicolau Manubens, \email{nicolau.manubens@bsc.es} +#'@author Veronica Torralba, \email{veronica.torralba@bsc.es} +#' +#'@param maps List of matrices to plot, each with (longitude, latitude) +#' dimensions, or 3-dimensional array with the dimensions (longitude, latitude, +#' map). Dimension names are required. +#'@param lon Vector of longitudes. Must match the length of the corresponding +#' dimension in 'maps'. +#'@param lat Vector of latitudes. Must match the length of the corresponding +#' dimension in 'maps'. +#'@param map_select_fun Function that selects, for each grid point, which value +#' to take among all the provided maps. This function receives as input a +#' vector of values for a same grid point for all the provided maps, and must +#' return a single selected value (not its index!) or NA. For example, the +#' \code{min} and \code{max} functions are accepted. +#'@param display_range Range of values to be displayed for all the maps. This +#' must be a numeric vector c(range min, range max). The values in the +#' parameter 'maps' can go beyond the limits specified in this range. If the +#' selected value for a given grid point (according to 'map_select_fun') falls +#' outside the range, it will be coloured with 'col_unknown_map'. +#'@param map_dim Optional name for the dimension of 'maps' along which the +#' multiple maps are arranged. Only applies when 'maps' is provided as a +#' 3-dimensional array. Takes the value 'map' by default. +#'@param brks Colour levels to be sent to PlotEquiMap. This parameter is +#' optional and adjusted automatically by the function. +#'@param cols List of vectors of colours to be sent to PlotEquiMap for the +#' colour bar of each map. This parameter is optional and adjusted +#' automatically by the function (up to 5 maps). The colours provided for each +#' colour bar will be automatically interpolated to match the number of breaks. +#' Each item in this list can be named, and the name will be used as title for +#' the corresponding colour bar (equivalent to the parameter 'bar_titles'). +#'@param col_unknown_map Colour to use to paint the grid cells for which a map +#' is not possible to be chosen according to 'map_select_fun' or for those +#' values that go beyond 'display_range'. Takes the value 'white' by default. +#'@param mask Optional numeric array with dimensions (latitude, longitude), with +#' values in the range [0, 1], indicating the opacity of the mask over each +#' grid point. Cells with a 0 will result in no mask, whereas cells with a 1 +#' will result in a totally opaque superimposed pixel coloured in 'col_mask'. +#'@param col_mask Colour to be used for the superimposed mask (if specified in +#' 'mask'). Takes the value 'grey' by default. +#'@param dots Array of same dimensions as 'var' or with dimensions +#' c(n, dim(var)), where n is the number of dot/symbol layers to add to the +#' plot. A value of TRUE at a grid cell will draw a dot/symbol on the +#' corresponding square of the plot. By default all layers provided in 'dots' +#' are plotted with dots, but a symbol can be specified for each of the +#' layers via the parameter 'dot_symbol'. +#'@param bar_titles Optional vector of character strings providing the titles to +#' be shown on top of each of the colour bars. +#'@param legend_scale Scale factor for the size of the colour bar labels. Takes +#' 1 by default. +#'@param cex_bar_titles Scale factor for the sizes of the bar titles. Takes 1.5 +#' by default. +#'@param plot_margin Numeric vector of length 4 for the margin sizes in the +#' following order: bottom, left, top, and right. If not specified, use the +#' default of par("mar"), c(5.1, 4.1, 4.1, 2.1). Used as 'margin_scale' in +#' s2dv::PlotEquiMap. +#'@param fileout File where to save the plot. If not specified (default) a +#' graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, bmp +#' and tiff +#'@param width File width, in the units specified in the parameter size_units +#' (inches by default). Takes 8 by default. +#'@param height File height, in the units specified in the parameter size_units +#' (inches by default). Takes 5 by default. +#'@param size_units Units of the size of the device (file or window) to plot in. +#' Inches ('in') by default. See ?Devices and the creator function of the +#' corresponding device. +#'@param res Resolution of the device (file or window) to plot in. See ?Devices +#' and the creator function of the corresponding device. +#'@param drawleg Where to draw the common colour bar. Can take values TRUE, +#' FALSE or:\cr +#' 'up', 'u', 'U', 'top', 't', 'T', 'north', 'n', 'N'\cr +#' 'down', 'd', 'D', 'bottom', 'b', 'B', 'south', 's', 'S' (default)\cr +#' 'right', 'r', 'R', 'east', 'e', 'E'\cr +#' 'left', 'l', 'L', 'west', 'w', 'W' +#'@param ... Additional parameters to be passed on to \code{PlotEquiMap}. +#' +#'@examples +#'# Simple example +#'x <- array(1:(20 * 10), dim = c(lat = 10, lon = 20)) / 200 +#'a <- x * 0.6 +#'b <- (1 - x) * 0.6 +#'c <- 1 - (a + b) +#'lons <- seq(0, 359.5, length = 20) +#'lats <- seq(-89.5, 89.5, length = 10) +#'\dontrun{ +#'PlotCombinedMap(list(a, b, c), lons, lats, +#' toptitle = 'Maximum map', +#' map_select_fun = max, +#' display_range = c(0, 1), +#' bar_titles = paste('% of belonging to', c('a', 'b', 'c')), +#' brks = 20, width = 12, height = 10) +#'} +#' +#'Lon <- c(0:40, 350:359) +#'Lat <- 51:26 +#'data <- rnorm(51 * 26 * 3) +#'dim(data) <- c(map = 3, lon = 51, lat = 26) +#'mask <- sample(c(0,1), replace = TRUE, size = 51 * 26) +#'dim(mask) <- c(lat = 26, lon = 51) +#'\dontrun{ +#'PlotCombinedMap(data, lon = Lon, lat = Lat, map_select_fun = max, +#' display_range = range(data), mask = mask, +#' width = 14, height = 10) +#'} +#' +#'@seealso \code{PlotCombinedMap} and \code{PlotEquiMap} +#' +#'@importFrom s2dv PlotEquiMap ColorBar +#'@importFrom maps map +#'@importFrom graphics box image layout mtext par plot.new +#'@importFrom grDevices adjustcolor bmp colorRampPalette dev.cur dev.new dev.off +#' hcl jpeg pdf png postscript svg tiff +#'@export +PlotCombinedMap <- function(maps, lon, lat, + map_select_fun, display_range, + map_dim = 'map', + brks = NULL, cols = NULL, + col_unknown_map = 'white', + mask = NULL, col_mask = 'grey', + dots = NULL, + bar_titles = NULL, legend_scale = 1, + cex_bar_titles = 1.5, + plot_margin = NULL, + fileout = NULL, width = 8, height = 5, + size_units = 'in', res = 100, drawleg = T, + ...) { + args <- list(...) + # If there is any filenames to store the graphics, process them + # to select the right device + if (!is.null(fileout)) { + .SelectDevice <- utils::getFromNamespace(".SelectDevice", "s2dv") + deviceInfo <- .SelectDevice(fileout = fileout, width = width, height = height, + units = size_units, res = res) + saveToFile <- deviceInfo$fun + fileout <- deviceInfo$files + } + + # Check probs + error <- FALSE + # Change list into an array + if (is.list(maps)) { + if (length(maps) < 1) { + stop("Parameter 'maps' must be of length >= 1 if provided as a list.") + } + check_fun <- function(x) { + is.numeric(x) && (length(dim(x)) == 2) + } + if (!all(sapply(maps, check_fun))) { + error <- TRUE + } + ref_dims <- dim(maps[[1]]) + equal_dims <- all(sapply(maps, function(x) identical(dim(x), ref_dims))) + if (!equal_dims) { + stop("All arrays in parameter 'maps' must have the same dimension ", + "sizes and names when 'maps' is provided as a list of arrays.") + } + num_maps <- length(maps) + maps <- unlist(maps) + dim(maps) <- c(ref_dims, map = num_maps) + map_dim <- 'map' + } + if (!is.numeric(maps)) { + error <- TRUE + } + if (is.null(dim(maps))) { + error <- TRUE + } + if (length(dim(maps)) != 3) { + error <- TRUE + } + if (error) { + stop("Parameter 'maps' must be either a numeric array with 3 dimensions ", + " or a list of numeric arrays of the same size with the 'lon' and ", + "'lat' dimensions.") + } + dimnames <- names(dim(maps)) + + # Check map_dim + if (is.character(map_dim)) { + if (is.null(dimnames)) { + stop("Specified a dimension name in 'map_dim' but no dimension names provided ", + "in 'maps'.") + } + map_dim <- which(dimnames == map_dim) + if (length(map_dim) < 1) { + stop("Dimension 'map_dim' not found in 'maps'.") + } else { + map_dim <- map_dim[1] + } + } else if (!is.numeric(map_dim)) { + stop("Parameter 'map_dim' must be either a numeric value or a ", + "dimension name.") + } + if (length(map_dim) != 1) { + stop("Parameter 'map_dim' must be of length 1.") + } + map_dim <- round(map_dim) + + # Work out lon_dim and lat_dim + lon_dim <- NULL + if (!is.null(dimnames)) { + lon_dim <- which(dimnames %in% c('lon', 'longitude'))[1] + } + if (length(lon_dim) < 1) { + lon_dim <- (1:3)[-map_dim][1] + } + lon_dim <- round(lon_dim) + + lat_dim <- NULL + if (!is.null(dimnames)) { + lat_dim <- which(dimnames %in% c('lat', 'latitude'))[1] + } + if (length(lat_dim) < 1) { + lat_dim <- (1:3)[-map_dim][2] + } + lat_dim <- round(lat_dim) + + # Check lon + if (!is.numeric(lon)) { + stop("Parameter 'lon' must be a numeric vector.") + } + if (length(lon) != dim(maps)[lon_dim]) { + stop("Parameter 'lon' does not match the longitude dimension in 'maps'.") + } + + # Check lat + if (!is.numeric(lat)) { + stop("Parameter 'lat' must be a numeric vector.") + } + if (length(lat) != dim(maps)[lat_dim]) { + stop("Parameter 'lat' does not match the longitude dimension in 'maps'.") + } + + # Check map_select_fun + if (is.numeric(map_select_fun)) { + if (length(dim(map_select_fun)) != 2) { + stop("Parameter 'map_select_fun' must be an array with dimensions ", + "'lon' and 'lat' if provided as an array.") + } + if (!identical(dim(map_select_fun), dim(maps)[-map_dim])) { + stop("The dimensions 'lon' and 'lat' in the 'map_select_fun' array must ", + "have the same size, name and order as in the 'maps' parameter.") + } + } + if (!is.function(map_select_fun)) { + stop("The parameter 'map_select_fun' must be a function or a numeric array.") + } + + # Generate the desired brks and cols. Only nmap, brks, cols, bar_limits, and + # bar_titles matter here because plot = F. + colorbar <- GradientCatsColorBar(nmap = dim(maps)[map_dim], + brks = brks, cols = cols, vertical = FALSE, + subsampleg = NULL, bar_limits = display_range, var_limits = NULL, + triangle_ends = NULL, plot = FALSE, draw_separators = TRUE, + bar_titles = bar_titles, title_scale = cex_bar_titles, label_scale = legend_scale * 1.5, + extra_margin = c(2, 0, 2, 0)) + + # Check legend_scale + if (!is.numeric(legend_scale)) { + stop("Parameter 'legend_scale' must be numeric.") + } + + # Check col_unknown_map + if (!is.character(col_unknown_map)) { + stop("Parameter 'col_unknown_map' must be a character string.") + } + + # Check col_mask + if (!is.character(col_mask)) { + stop("Parameter 'col_mask' must be a character string.") + } + + # Check mask + if (!is.null(mask)) { + if (!is.numeric(mask)) { + stop("Parameter 'mask' must be numeric.") + } + if (length(dim(mask)) != 2) { + stop("Parameter 'mask' must have two dimensions.") + } + if ((dim(mask)[1] != dim(maps)[lat_dim]) || + (dim(mask)[2] != dim(maps)[lon_dim])) { + stop("Parameter 'mask' must have dimensions c(lat, lon).") + } + } + # Check dots + if (!is.null(dots)) { + if (length(dim(dots)) != 2) { + stop("Parameter 'mask' must have two dimensions.") + } + if ((dim(dots)[1] != dim(maps)[lat_dim]) || + (dim(dots)[2] != dim(maps)[lon_dim])) { + stop("Parameter 'mask' must have dimensions c(lat, lon).") + } + } + + #---------------------- + # Identify the most likely map + #---------------------- + brks_norm <- seq(0, 1, length.out = length(colorbar$brks)) + if (is.function(map_select_fun)) { + range_width <- display_range[2] - display_range[1] + ml_map <- apply(maps, c(lat_dim, lon_dim), function(x) { + if (any(is.na(x))) { + res <- NA + } else { + res <- which(x == map_select_fun(x)) + if (length(res) > 0) { + res <- res[1] + if (map_select_fun(x) < display_range[1] || + map_select_fun(x) > display_range[2]) { + res <- -0.5 + } else { + res <- res + (map_select_fun(x) - display_range[1]) / range_width + if (map_select_fun(x) == display_range[1]) { + res <- res + brks_norm[2] / (length(brks_norm) * 2) + } + } + } else { + res <- -0.5 + } + } + res + }) + } else { + stop("Providing 'map_select_fun' as array not implemented yet.") + ml_map <- map_select_fun + } + nmap <- dim(maps)[map_dim] + nlat <- length(lat) + nlon <- length(lon) + + #---------------------- + # Set latitudes from minimum to maximum + #---------------------- + if (lat[1] > lat[nlat]){ + lat <- lat[nlat:1] + indices <- list(nlat:1, TRUE) + ml_map <- do.call("[", c(list(x = ml_map), indices)) + if (!is.null(mask)){ + mask <- mask[nlat:1, ] + } + if (!is.null(dots)){ + dots <- dots[nlat:1,] + } + } + + #---------------------- + # Set layout and parameters + #---------------------- + # Open connection to graphical device + if (!is.null(fileout)) { + saveToFile(fileout) + } else if (names(dev.cur()) == 'null device') { + dev.new(units = size_units, res = res, width = width, height = height) + } + #NOTE: I think plot.new() is not necessary in any case. +# plot.new() + #TODO: Don't hardcoded. Let users decide. + par(font.main = 1) + # If colorbars need to be plotted, re-define layout. + if (drawleg) { + layout(matrix(c(rep(1, nmap),2:(nmap + 1)), 2, nmap, byrow = TRUE), heights = c(6, 1.5)) + } + + #---------------------- + # Set colors and breaks and then PlotEquiMap + #---------------------- + tcols <- c(col_unknown_map, colorbar$cols[[1]]) + for (k in 2:nmap) { + tcols <- append(tcols, c(col_unknown_map, colorbar$cols[[k]])) + } + tbrks <- c(-1, brks_norm + rep(1:nmap, each = length(brks_norm))) + + if (is.null(plot_margin)) { + plot_margin <- c(5, 4, 4, 2) + 0.1 # default of par()$mar + } + + PlotEquiMap(var = ml_map, lon = lon, lat = lat, + brks = tbrks, cols = tcols, drawleg = FALSE, + filled.continents = FALSE, dots = dots, margin_scale = plot_margin, ...) + + #---------------------- + # Add overplot on top + #---------------------- + if (!is.null(mask)) { + dims_mask <- dim(mask) + latb <- sort(lat, index.return = TRUE) + dlon <- lon[2:dims_mask[2]] - lon[1:(dims_mask[2] - 1)] + wher <- which(dlon > (mean(dlon) + 1)) + if (length(wher) > 0) { + lon[(wher + 1):dims_mask[2]] <- lon[(wher + 1):dims_mask[2]] - 360 + } + lonb <- sort(lon, index.return = TRUE) + + cols_mask <- sapply(seq(from = 0, to = 1, length.out = 10), + function(x) adjustcolor(col_mask, alpha.f = x)) + image(lonb$x, latb$x, t(mask)[lonb$ix, latb$ix], + axes = FALSE, col = cols_mask, + breaks = seq(from = 0, to = 1, by = 0.1), + xlab='', ylab='', add = TRUE, xpd = TRUE) + if (!exists('coast_color')) { + coast_color <- 'black' + } + if (min(lon) < 0) { + map('world', interior = FALSE, add = TRUE, lwd = 1, col = coast_color) # Low resolution world map (lon -180 to 180). + } else { + map('world2', interior = FALSE, add = TRUE, lwd = 1, col = coast_color) # Low resolution world map (lon 0 to 360). + } + box() + } + + #---------------------- + # Add colorbars + #---------------------- + if ('toptitle' %in% names(args)) { + size_title <- 1 + if ('title_scale' %in% names(args)) { + size_title <- args[['title_scale']] + } + old_mar <- par('mar') + old_mar[3] <- old_mar[3] - (2 * size_title + 1) + par(mar = old_mar) + } + + if (drawleg) { + GradientCatsColorBar(nmap = dim(maps)[map_dim], + brks = colorbar$brks, cols = colorbar$cols, vertical = FALSE, + subsampleg = NULL, bar_limits = display_range, var_limits = NULL, + triangle_ends = NULL, plot = TRUE, draw_separators = TRUE, + bar_titles = bar_titles, title_scale = cex_bar_titles, label_scale = legend_scale * 1.5, + extra_margin = c(2, 0, 2, 0)) + } + + # If the graphic was saved to file, close the connection with the device + if (!is.null(fileout)) dev.off() +} + + diff --git a/modules/Visualization/R/tmp/PlotMostLikelyQuantileMap.R b/modules/Visualization/R/tmp/PlotMostLikelyQuantileMap.R new file mode 100644 index 0000000000000000000000000000000000000000..a81515e986cb78a0b52339f580ab9ff918053036 --- /dev/null +++ b/modules/Visualization/R/tmp/PlotMostLikelyQuantileMap.R @@ -0,0 +1,223 @@ +#'Plot Maps of Most Likely Quantiles +#' +#'@author Veronica Torralba, \email{veronica.torralba@bsc.es}, Nicolau Manubens, +#'\email{nicolau.manubens@bsc.es} +#'@description This function receives as main input (via the parameter +#'\code{probs}) a collection of longitude-latitude maps, each containing the +#'probabilities (from 0 to 1) of the different grid cells of belonging to a +#'category. As many categories as maps provided as inputs are understood to +#'exist. The maps of probabilities must be provided on a common rectangular +#'regular grid, and a vector with the longitudes and a vector with the latitudes +#'of the grid must be provided. The input maps can be provided in two forms, +#'either as a list of multiple two-dimensional arrays (one for each category) or +#'as a three-dimensional array, where one of the dimensions corresponds to the +#'different categories. +#' +#'@param probs A list of bi-dimensional arrays with the named dimensions +#' 'latitude' (or 'lat') and 'longitude' (or 'lon'), with equal size and in the +#' same order, or a single tri-dimensional array with an additional dimension +#' (e.g. 'bin') for the different categories. The arrays must contain +#' probability values between 0 and 1, and the probabilities for all categories +#' of a grid cell should not exceed 1 when added. +#'@param lon A numeric vector with the longitudes of the map grid, in the same +#' order as the values along the corresponding dimension in \code{probs}. +#'@param lat A numeric vector with the latitudes of the map grid, in the same +#' order as the values along the corresponding dimension in \code{probs}. +#'@param cat_dim The name of the dimension along which the different categories +#' are stored in \code{probs}. This only applies if \code{probs} is provided in +#' the form of 3-dimensional array. The default expected name is 'bin'. +#'@param bar_titles Vector of character strings with the names to be drawn on +#' top of the color bar for each of the categories. As many titles as +#' categories provided in \code{probs} must be provided. +#'@param col_unknown_cat Character string with a colour representation of the +#' colour to be used to paint the cells for which no category can be clearly +#' assigned. Takes the value 'white' by default. +#'@param drawleg Where to draw the common colour bar. Can take values TRUE, +#' FALSE or:\cr +#' 'up', 'u', 'U', 'top', 't', 'T', 'north', 'n', 'N'\cr +#' 'down', 'd', 'D', 'bottom', 'b', 'B', 'south', 's', 'S' (default)\cr +#' 'right', 'r', 'R', 'east', 'e', 'E'\cr +#' 'left', 'l', 'L', 'west', 'w', 'W' +#'@param ... Additional parameters to be sent to \code{PlotCombinedMap} and +#' \code{PlotEquiMap}. +#'@seealso \code{PlotCombinedMap} and \code{PlotEquiMap} +#'@examples +#'# Simple example +#'x <- array(1:(20 * 10), dim = c(lat = 10, lon = 20)) / 200 +#'a <- x * 0.6 +#'b <- (1 - x) * 0.6 +#'c <- 1 - (a + b) +#'lons <- seq(0, 359.5, length = 20) +#'lats <- seq(-89.5, 89.5, length = 10) +#'\dontrun{ +#'PlotMostLikelyQuantileMap(list(a, b, c), lons, lats, +#' toptitle = 'Most likely tercile map', +#' bar_titles = paste('% of belonging to', c('a', 'b', 'c')), +#' brks = 20, width = 10, height = 8) +#'} +#' +#'# More complex example +#'n_lons <- 40 +#'n_lats <- 20 +#'n_timesteps <- 100 +#'n_bins <- 4 +#' +#'# 1. Generation of sample data +#'lons <- seq(0, 359.5, length = n_lons) +#'lats <- seq(-89.5, 89.5, length = n_lats) +#' +#'# This function builds a 3-D gaussian at a specified point in the map. +#'make_gaussian <- function(lon, sd_lon, lat, sd_lat) { +#' w <- outer(lons, lats, function(x, y) dnorm(x, lon, sd_lon) * dnorm(y, lat, sd_lat)) +#' min_w <- min(w) +#' w <- w - min_w +#' w <- w / max(w) +#' w <- t(w) +#' names(dim(w)) <- c('lat', 'lon') +#' w +#'} +#' +#'# This function generates random time series (with values ranging 1 to 5) +#'# according to 2 input weights. +#'gen_data <- function(w1, w2, n) { +#' r <- sample(1:5, n, +#' prob = c(.05, .9 * w1, .05, .05, .9 * w2), +#' replace = TRUE) +#' r <- r + runif(n, -0.5, 0.5) +#' dim(r) <- c(time = n) +#' r +#'} +#' +#'# We build two 3-D gaussians. +#'w1 <- make_gaussian(120, 80, 20, 30) +#'w2 <- make_gaussian(260, 60, -10, 40) +#' +#'# We generate sample data (with dimensions time, lat, lon) according +#'# to the generated gaussians +#'sample_data <- multiApply::Apply(list(w1, w2), NULL, +#' gen_data, n = n_timesteps)$output1 +#' +#'# 2. Binning sample data +#'prob_thresholds <- 1:n_bins / n_bins +#'prob_thresholds <- prob_thresholds[1:(n_bins - 1)] +#'thresholds <- quantile(sample_data, prob_thresholds) +#' +#'binning <- function(x, thresholds) { +#' n_samples <- length(x) +#' n_bins <- length(thresholds) + 1 +#' +#' thresholds <- c(thresholds, max(x)) +#' result <- 1:n_bins +#' lower_threshold <- min(x) - 1 +#' for (i in 1:n_bins) { +#' result[i] <- sum(x > lower_threshold & x <= thresholds[i]) / n_samples +#' lower_threshold <- thresholds[i] +#' } +#' +#' dim(result) <- c(bin = n_bins) +#' result +#'} +#' +#'bins <- multiApply::Apply(sample_data, 'time', binning, thresholds)$output1 +#' +#'# 3. Plotting most likely quantile/bin +#'\dontrun{ +#'PlotMostLikelyQuantileMap(bins, lons, lats, +#' toptitle = 'Most likely quantile map', +#' bar_titles = paste('% of belonging to', letters[1:n_bins]), +#' mask = 1 - (w1 + w2 / max(c(w1, w2))), +#' brks = 20, width = 10, height = 8) +#'} +#'@importFrom maps map +#'@importFrom graphics box image layout mtext par plot.new +#'@importFrom grDevices adjustcolor bmp colorRampPalette dev.cur dev.new dev.off hcl jpeg pdf png postscript svg tiff +#'@export +PlotMostLikelyQuantileMap <- function(probs, lon, lat, cat_dim = 'bin', + bar_titles = NULL, + col_unknown_cat = 'white', drawleg = T, + ...) { + # Check probs + error <- FALSE + if (is.list(probs)) { + if (length(probs) < 1) { + stop("Parameter 'probs' must be of length >= 1 if provided as a list.") + } + check_fun <- function(x) { + is.numeric(x) && (length(dim(x)) == 2) + } + if (!all(sapply(probs, check_fun))) { + error <- TRUE + } + ref_dims <- dim(probs[[1]]) + equal_dims <- all(sapply(probs, function(x) identical(dim(x), ref_dims))) + if (!equal_dims) { + stop("All arrays in parameter 'probs' must have the same dimension ", + "sizes and names when 'probs' is provided as a list of arrays.") + } + num_probs <- length(probs) + probs <- unlist(probs) + dim(probs) <- c(ref_dims, map = num_probs) + cat_dim <- 'map' + } + if (!is.numeric(probs)) { + error <- TRUE + } + if (is.null(dim(probs))) { + error <- TRUE + } + if (length(dim(probs)) != 3) { + error <- TRUE + } + if (error) { + stop("Parameter 'probs' must be either a numeric array with 3 dimensions ", + " or a list of numeric arrays of the same size with the 'lon' and ", + "'lat' dimensions.") + } + dimnames <- names(dim(probs)) + + # Check cat_dim + if (is.character(cat_dim)) { + if (is.null(dimnames)) { + stop("Specified a dimension name in 'cat_dim' but no dimension names provided ", + "in 'probs'.") + } + cat_dim <- which(dimnames == cat_dim) + if (length(cat_dim) < 1) { + stop("Dimension 'cat_dim' not found in 'probs'.") + } + cat_dim <- cat_dim[1] + } else if (!is.numeric(cat_dim)) { + stop("Parameter 'cat_dim' must be either a numeric value or a ", + "dimension name.") + } + if (length(cat_dim) != 1) { + stop("Parameter 'cat_dim' must be of length 1.") + } + cat_dim <- round(cat_dim) + nprobs <- dim(probs)[cat_dim] + + # Check bar_titles + if (is.null(bar_titles)) { + if (nprobs == 3) { + bar_titles <- c("Below normal (%)", "Normal (%)", "Above normal (%)") + } else if (nprobs == 5) { + bar_titles <- c("Low (%)", "Below normal (%)", + "Normal (%)", "Above normal (%)", "High (%)") + } else { + bar_titles <- paste0("Cat. ", 1:nprobs, " (%)") + } + } + + minimum_value <- ceiling(1 / nprobs * 10 * 1.1) * 10 + + # By now, the PlotCombinedMap function is included below in this file. + # In the future, PlotCombinedMap will be part of s2dverification and will + # be properly imported. + PlotCombinedMap(probs * 100, lon, lat, map_select_fun = max, + display_range = c(minimum_value, 100), + map_dim = cat_dim, + bar_titles = bar_titles, + col_unknown_map = col_unknown_cat, + drawleg = drawleg, ...) +} + diff --git a/modules/Visualization/R/tmp/PlotRobinson.R b/modules/Visualization/R/tmp/PlotRobinson.R index 6a3412f53d8261748f6daa31171a0096dd99f614..d8210258275d481c9e2724e0e860281fc2e6e048 100644 --- a/modules/Visualization/R/tmp/PlotRobinson.R +++ b/modules/Visualization/R/tmp/PlotRobinson.R @@ -314,8 +314,19 @@ PlotRobinson <- function(data, lon, lat, lon_dim = NULL, lat_dim = NULL, 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 brks to triangles - brks_ggplot <- c(min(data, na.rm = T), brks, max(data, na.rm = T)) + + # Add triangles to brks + brks_ggplot <- brks + if (max(data, na.rm = T) > 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 (min(data, na.rm = T) < 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 @@ -328,8 +339,15 @@ PlotRobinson <- function(data, lon, lat, lon_dim = NULL, lat_dim = NULL, 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)) } @@ -436,12 +454,18 @@ PlotRobinson <- function(data, lon, lat, lon_dim = NULL, lat_dim = NULL, #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 = range(data_df$long), ylim = range(data_df$lat), + coord_sf(xlim = coord_sf_lim[1:2], ylim = coord_sf_lim[3:4], expand = TRUE, datum = target_proj) if (!is.null(dots)) { diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R index 03a4d5a25926a21bdb1adc3b815eeb7fc539c226..5fb3b946b3251d2c13c68d671dd8937dbc4de7fd 100644 --- a/modules/Visualization/Visualization.R +++ b/modules/Visualization/Visualization.R @@ -39,7 +39,10 @@ plot_data <- function(recipe, so there is no data that can be plotted.") stop() } - + # Set default single-panel plots if not specified + if (is.null(recipe$Analysis$Workflow$Visualization$multi_panel)) { + recipe$Analysis$Workflow$Visualization$multi_panel <- FALSE + } # Plot skill metrics if ("skill_metrics" %in% plots) { if (!is.null(skill_metrics)) { @@ -66,8 +69,61 @@ plot_data <- function(recipe, # Plot Most Likely Terciles if ("most_likely_terciles" %in% plots) { if ((!is.null(probabilities)) && (!is.null(data$fcst))) { - plot_most_likely_terciles(recipe, data$fcst, - probabilities, outdir) + # Assign default parameters + if (is.null(recipe$Analysis$Workflow$Visualization$mask_terciles)) { + recipe$Analysis$Workflow$Visualization$mask_terciles <- FALSE + } + if (is.null(recipe$Analysis$Workflow$Visualization$dots)) { + recipe$Analysis$Workflow$Visualization$dots <- FALSE + } + # Plots without masked terciles/dots + if ((recipe$Analysis$Workflow$Visualization$mask_terciles + %in% c('both', FALSE)) || + (recipe$Analysis$Workflow$Visualization$dots + %in% c('both', FALSE))) { + plot_most_likely_terciles(recipe, data$fcst, + probabilities, + mask = NULL, + dots = NULL, + outdir) + } + # Plots with masked terciles + if (recipe$Analysis$Workflow$Visualization$mask_terciles %in% + c('both', TRUE)) { + if (is.null(skill_metrics)) { + error(recipe$Run$logger, + paste0("For the most likely terciles plot, skill_metrics ", + "need to be provided to be masked.")) + } else if (!('rpss' %in% names(skill_metrics))) { + error(recipe$Run$logger, + paste0("For the most likely terciles plot, rpss metric ", + "need to be provided to be masked")) + } else { + plot_most_likely_terciles(recipe, data$fcst, + probabilities, + mask = skill_metrics$rpss, + dots = NULL, + outdir) + } + } + # Plots with dotted terciles + if (recipe$Analysis$Workflow$Visualization$dots %in% c('both', TRUE)) { + if (is.null(skill_metrics)) { + error(recipe$Run$logger, + paste0("For the most likely terciles plot, skill_metrics ", + "need to be provided for the dots")) + } else if (!('rpss' %in% names(skill_metrics))) { + error(recipe$Run$logger, + paste0("For the most likely terciles plot, rpss metric ", + "needs to be provided for the dots")) + } else { + plot_most_likely_terciles(recipe, data$fcst, + probabilities, + mask = NULL, + dots = skill_metrics$rpss, + outdir) + } + } } else { error(recipe$Run$logger, paste0("For the most likely terciles plot, both the fcst and the ", @@ -75,4 +131,3 @@ plot_data <- function(recipe, } } } - diff --git a/recipe_ecvs_seasonal_oper.yml b/recipe_ecvs_seasonal_oper.yml new file mode 100644 index 0000000000000000000000000000000000000000..d47fd1593749cb07c04dd8166a73075ac2058d76 --- /dev/null +++ b/recipe_ecvs_seasonal_oper.yml @@ -0,0 +1,73 @@ +Description: + Author: nperez + Info: ECVs Oper ESS ECMWF SEAS5 Seasonal Forecast recipe (monthly mean, tas) + +Analysis: + Horizon: seasonal # Mandatory, str: either subseasonal, seasonal, or decadal + Variables: + - {name: tas, freq: monthly_mean} + - {name: prlr, freq: monthly_mean} + Datasets: + System: + - {name: ECMWF-SEAS5.1} # system21_m1 system35c3s + Multimodel: no # Mandatory, bool: Either yes/true or no/false + Reference: + - {name: ERA5} # Mandatory, str: Reference codename. See docu. + Time: + sdate: '0701' ## MMDD + fcst_year: '2023' # Optional, int: Forecast year 'YYYY' + hcst_start: '1993' # Mandatory, int: Hindcast start year 'YYYY' + hcst_end: '2016' # Mandatory, int: Hindcast end year 'YYYY' + ftime_min: 1 # Mandatory, int: First leadtime time step in months + ftime_max: 6 # Mandatory, int: Last leadtime time step in months + Region: + - {name: "UE", latmin: 20, latmax: 80, lonmin: -20, lonmax: 40} + Regrid: + method: bilinear # Mandatory, str: Interpolation method. See docu. + type: "to_system" + #type: /esarchive/scratch/nmilders/gitlab/git_clones/auto-s2s/conf/grid_description.txt #'r360x180' # Mandatory, str: to_system, to_reference, or CDO-accepted grid. + Workflow: + Anomalies: + compute: no + cross_validation: no + save: none + Calibration: + method: evmos # Mandatory, str: Calibration method. See docu. + cross_validation: yes + save: none + Skill: + metric: mean_bias EnsCorr rpss crpss bss10 bss90 + save: 'all' + cross_validation: yes + Probabilities: + percentiles: [[1/3, 2/3], [1/10, 9/10]] # frac: Quantile thresholds. + save: 'all' + Indicators: + index: no + Visualization: + plots: skill_metrics forecast_ensemble_mean most_likely_terciles + multi_panel: no + projection: lambert_europe + ncores: 4 # Optional, int: number of cores, defaults to 1 + remove_NAs: # Optional, bool: Whether NAs are removed, defaults to FALSE + Output_format: scorecards + logo: yes +Run: + Loglevel: INFO + Terminal: yes + filesystem: esarchive + output_dir: /esarchive/scratch/nperez/cs_oper/seasonal/ # replace with the directory where you want to save the outputs + code_dir: /esarchive/scratch/nperez/git/auto-s2s/ # replace with the directory where your code is + autosubmit: yes + # fill only if using autosubmit + auto_conf: + script: /esarchive/scratch/nperez/git/auto-s2s/exec_ecvs_seasonal_oper.R # replace with the path to your script + expid: a68v # replace with your EXPID + hpc_user: bsc32339 # replace with your hpc username + wallclock: 02:00 # hh:mm + processors_per_job: 4 + platform: nord3v2 + email_notifications: yes # enable/disable email notifications. Change it if you want to. + email_address: nuria.perez@bsc.es # replace with your email address + notify_completed: yes # notify me by email when a job finishes + notify_failed: yes # notify me by email when a job fails diff --git a/recipe_prlr_seasonal_oper.yml b/recipe_prlr_seasonal_oper.yml new file mode 100644 index 0000000000000000000000000000000000000000..fb6d1a1cd039df755abc45da9aae0d4486abdc5c --- /dev/null +++ b/recipe_prlr_seasonal_oper.yml @@ -0,0 +1,62 @@ +Description: + Author: nperez + Info: ECVs Oper ESS ECMWF SEAS5 Seasonal Forecast recipe (monthly mean, tas) + +Analysis: + Horizon: seasonal # Mandatory, str: either subseasonal, seasonal, or decadal + Variables: + name: prlr + freq: monthly_mean + Datasets: + System: + name: ECMWF-SEAS5.1 # Mandatory, str: system5c3s system21_m1 system35c3s + Multimodel: no # Mandatory, bool: Either yes/true or no/false + Reference: + name: ERA5 # Mandatory, str: Reference codename. See docu. + Time: + sdate: '0701' ## MMDD + fcst_year: '2023' # Optional, int: Forecast year 'YYYY' + hcst_start: '1993' # Mandatory, int: Hindcast start year 'YYYY' + hcst_end: '2016' # Mandatory, int: Hindcast end year 'YYYY' + ftime_min: 1 # Mandatory, int: First leadtime time step in months + ftime_max: 6 # Mandatory, int: Last leadtime time step in months + Region: + latmin: 20 # Mandatory, int: minimum latitude + latmax: 80 # Mandatory, int: maximum latitude + lonmin: -20 # Mandatory, int: minimum longitude + lonmax: 40 # Mandatory, int: maximum longitude + Regrid: + method: bilinear # Mandatory, str: Interpolation method. See docu. + type: "to_system" + #type: /esarchive/scratch/nmilders/gitlab/git_clones/auto-s2s/conf/grid_description.txt #'r360x180' # Mandatory, str: to_system, to_reference, or CDO-accepted grid. + Workflow: + Anomalies: + compute: no + cross_validation: no + save: none + Calibration: + method: evmos # Mandatory, str: Calibration method. See docu. + cross_validation: yes + save: none + Skill: + metric: mean_bias EnsCorr rpss crpss bss10 bss90 + save: 'all' + cross_validation: yes + Probabilities: + percentiles: [[1/3, 2/3], [1/10, 9/10]] # frac: Quantile thresholds. + save: 'all' + Indicators: + index: no + Visualization: + plots: skill_metrics forecast_ensemble_mean most_likely_terciles + multi_panel: no + dots: both + ncores: 4 # Optional, int: number of cores, defaults to 1 + remove_NAs: # Optional, bool: Whether NAs are removed, defaults to FALSE + Output_format: scorecards + logo: yes +Run: + Loglevel: INFO + Terminal: TRUE + output_dir: /esarchive/scratch/nperez/cs_oper/ + code_dir: /esarchive/scratch/nperez/git/s2s-suite/ diff --git a/recipe_tas_seasonal_oper.yml b/recipe_tas_seasonal_oper.yml new file mode 100644 index 0000000000000000000000000000000000000000..75918265ba4d12a20226bd69511ef53eb1de1d4e --- /dev/null +++ b/recipe_tas_seasonal_oper.yml @@ -0,0 +1,62 @@ +Description: + Author: nperez + Info: ECVs Oper ESS ECMWF SEAS5 Seasonal Forecast recipe (monthly mean, tas) + +Analysis: + Horizon: seasonal # Mandatory, str: either subseasonal, seasonal, or decadal + Variables: + name: tas + freq: monthly_mean + Datasets: + System: + name: ECMWF-SEAS5.1 # Mandatory, str: system5c3s system21_m1 system35c3s + Multimodel: no # Mandatory, bool: Either yes/true or no/false + Reference: + name: ERA5 # Mandatory, str: Reference codename. See docu. + Time: + sdate: '0601' ## MMDD + fcst_year: '2023' # Optional, int: Forecast year 'YYYY' + hcst_start: '1993' # Mandatory, int: Hindcast start year 'YYYY' + hcst_end: '2016' # Mandatory, int: Hindcast end year 'YYYY' + ftime_min: 1 # Mandatory, int: First leadtime time step in months + ftime_max: 6 # Mandatory, int: Last leadtime time step in months + Region: + latmin: 20 # Mandatory, int: minimum latitude + latmax: 80 # Mandatory, int: maximum latitude + lonmin: -20 # Mandatory, int: minimum longitude + lonmax: 40 # Mandatory, int: maximum longitude + Regrid: + method: bilinear # Mandatory, str: Interpolation method. See docu. + type: "to_system" + #type: /esarchive/scratch/nmilders/gitlab/git_clones/auto-s2s/conf/grid_description.txt #'r360x180' # Mandatory, str: to_system, to_reference, or CDO-accepted grid. + Workflow: + Anomalies: + compute: no + cross_validation: no + save: none + Calibration: + method: evmos # Mandatory, str: Calibration method. See docu. + cross_validation: yes + save: none + Skill: + metric: mean_bias EnsCorr rpss crpss bss10 bss90 + save: 'all' + cross_validation: yes + Probabilities: + percentiles: [[1/3, 2/3], [1/10, 9/10]] # frac: Quantile thresholds. + save: 'all' + Indicators: + index: no + Visualization: + plots: skill_metrics forecast_ensemble_mean most_likely_terciles + multi_panel: no + mask_terciles: both + ncores: 4 # Optional, int: number of cores, defaults to 1 + remove_NAs: # Optional, bool: Whether NAs are removed, defaults to FALSE + Output_format: scorecards + logo: yes +Run: + Loglevel: INFO + Terminal: TRUE + output_dir: /esarchive/scratch/nperez/cs_oper/ + code_dir: /esarchive/scratch/nperez/git/s2s-suite/ diff --git a/recipes/atomic_recipes/recipe_test_multivar.yml b/recipes/atomic_recipes/recipe_test_multivar.yml index 94e41223aed39ffd5757d8f636e8a260ce3fe7b0..8eb3e962e0ae5f05ead34ffa8dd23aa88d9b3293 100644 --- a/recipes/atomic_recipes/recipe_test_multivar.yml +++ b/recipes/atomic_recipes/recipe_test_multivar.yml @@ -41,6 +41,12 @@ Analysis: Probabilities: percentiles: [[1/3, 2/3], [1/10, 9/10]] save: 'all' + Visualization: + plots: most_likely_terciles + multi_panel: no + projection: cylindrical_equidistant + dots: both + mask_terciles: both Indicators: index: no ncores: 7 diff --git a/tools/check_recipe.R b/tools/check_recipe.R index c105353455c984dd9d3e412b0dd0275c979f1624..cfac94e397c18849102e2458fde0a55819bc49ea 100644 --- a/tools/check_recipe.R +++ b/tools/check_recipe.R @@ -400,22 +400,47 @@ check_recipe <- function(recipe) { } } # Check multi_panel option - if ((is.null(recipe$Analysis$Workflow$Visualization$multi_panel)) || - (!is.logical(recipe$Analysis$Workflow$Visualization$multi_panel))) { + ## TODO: Multi-panel + if (is.null(recipe$Analysis$Workflow$Visualization$multi_panel)) { + warn(recipe$Run$logger, + paste0("Visualization:multi_panel not specified for the plots, the", + " default is 'no/False'.")) + } else if (!is.logical(recipe$Analysis$Workflow$Visualization$multi_panel)) { error(recipe$Run$logger, - paste0("The 'multi_panel' element must be defined under ", - "Visualization, the options are 'yes' or 'no'.")) + paste0("Parameter 'Visualization:multi_panel' must be a logical ", + "value: either 'yes/True' or 'no/False'")) error_status <- T } # Check projection if (is.null(recipe$Analysis$Workflow$Visualization$projection)) { warn(recipe$Run$logger, - paste0("No projection has been specified for the plots, the ", + paste0("Visualization:projection not specified for the plots, the ", "default projection is cylindrical equidistant.")) - ## TODO: This will not work - recipe$Analysis$Workflow$Visualization$projection <- "cylindrical_equidistant" } ## TODO: Add significance? + if ("most_likely_terciles" %in% plots) { + if (is.null(recipe$Analysis$Workflow$Visualization$mask_terciles)) { + warn(recipe$Run$logger, + paste0("Visualization:mask_terciles not set for tercile plots,", + " the default setting is: 'no/False'.")) + } else if (!(recipe$Analysis$Workflow$Visualization$mask_terciles %in% + c(TRUE, FALSE, "both"))) { + error(recipe$Run$logger, + paste0("Parameter Visualization:mask_terciles must be one of: ", + "yes/True, no/False, 'both'")) + error_status <- T + } + if (is.null(recipe$Analysis$Workflow$Visualization$dots)) { + warn(recipe$Run$logger, + paste0("Visualization:dots not set for tercile plots, the default", + " setting is: 'no/False'.")) + } else if (!(recipe$Analysis$Workflow$Visualization$dots %in% + c(TRUE, FALSE, "both"))) { + error(recipe$Run$logger, + paste0("Parameter Visualization:plots must be one of: ", + "yes/True, no/False, 'both'")) + } + } } # --------------------------------------------------------------------- @@ -550,6 +575,5 @@ check_recipe <- function(recipe) { " startup.log file.") } else { info(recipe$Run$logger, "##### RECIPE CHECK SUCCESSFULL #####") - # return(append(nverifications, fcst.sdate)) } }