From 2157419a88fd11831468a7e382e43e92953da4a6 Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 28 Jul 2023 17:40:43 +0200 Subject: [PATCH 01/21] First approach --- modules/Units/R/transform_units_preasure.R | 2 + .../Units/R/transform_units_precipitation.R | 18 ++++++++ modules/Units/R/transform_units_temperature.R | 10 +++++ modules/Units/R/units_transform.R | 25 +++++++++++ modules/Units/Units.R | 45 +++++++++++++++++++ 5 files changed, 100 insertions(+) create mode 100644 modules/Units/R/transform_units_preasure.R create mode 100644 modules/Units/R/transform_units_precipitation.R create mode 100644 modules/Units/R/transform_units_temperature.R create mode 100644 modules/Units/R/units_transform.R create mode 100644 modules/Units/Units.R diff --git a/modules/Units/R/transform_units_preasure.R b/modules/Units/R/transform_units_preasure.R new file mode 100644 index 00000000..30694114 --- /dev/null +++ b/modules/Units/R/transform_units_preasure.R @@ -0,0 +1,2 @@ +transform_units_preasure <- function(data, data, original_units, new_units){ +} diff --git a/modules/Units/R/transform_units_precipitation.R b/modules/Units/R/transform_units_precipitation.R new file mode 100644 index 00000000..81e900b0 --- /dev/null +++ b/modules/Units/R/transform_units_precipitation.R @@ -0,0 +1,18 @@ +transform_units_precipitation <- function(data, original_units, new_units){ + if (!substr(original_units, 1, 2) == substr(new_units, 1, 2)){ + if (substr(original_units, 1, 2) == 'mm' & (substr(new_units, 1, 2) == 'm ' | substr(new_units, 1, 2) == 'm/')){ + data <- data / 1000 + } else if ((substr(original_units, 1, 2) == 'm ' | substr(original_units, 1, 2) == 'm/') & substr(new_units, 1, 2) == 'mm'){ + data <- data * 1000 + } else {stop(paste0("not prepared to transform from ", original_units, " to ", new_units))} + } + if (('s' %in% strsplit(original_units, " |/|per") | 'second' %in% strsplit(original_units, " |/|per")) & 'month' %in% strsplit(new_units, " |/|per")){ + data <- data * 3600 * 24 * 365.25/12 + } else if ('day' %in% strsplit(original_units, " |/|per")[[1]] & 'month' %in% strsplit(new_units, " |/|per")[[1]]){ + data <- data * 365.25/12 + } else {stop(paste0("not prepared to transform from ", original_units, " to ", new_units))} + + return(data) +} + + diff --git a/modules/Units/R/transform_units_temperature.R b/modules/Units/R/transform_units_temperature.R new file mode 100644 index 00000000..859ae479 --- /dev/null +++ b/modules/Units/R/transform_units_temperature.R @@ -0,0 +1,10 @@ +transform_units_temperature <- function(data, original_units, new_units){ + if (original_units == 'C' & new_units == 'K'){ + data <- data + 273.15 + } + if (original_units == 'K' & new_units == 'C'){ + data <- data - 273.15 + } + return(data) +} + diff --git a/modules/Units/R/units_transform.R b/modules/Units/R/units_transform.R new file mode 100644 index 00000000..088d3d0b --- /dev/null +++ b/modules/Units/R/units_transform.R @@ -0,0 +1,25 @@ +# data is a s2dv_cube object +# units as character +units_transform <- function(data, orig_units, user_units) { + if (orig_units %in% c("C", "k")) { + if (user_units %in% c("C", "k")) { + trans <- transform_units_temperature(data, orig_units, user_units) + } else { + stop("Transformation temperature units not available.") + } if else (orig_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { + if (user_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { + trans <- transform_units_precipitation(data, orig_units, user_units) + } else { + stop("Transformation precipitation units not available.") + } + } if else (orig_units %in% c("pa", "hpa", "mb")) { + if (user_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { + trans <- transform_units_pressure(data, orig_units, user_units) + } else { + stop("Transformation precipitation units not available.") + } + } else { + stop("Transformation unknown.") + } +} + diff --git a/modules/Units/Units.R b/modules/Units/Units.R new file mode 100644 index 00000000..15a7e628 --- /dev/null +++ b/modules/Units/Units.R @@ -0,0 +1,45 @@ +# This function aims to convert units +## +source("modules/R/units_transform.R") +source("modules/R/transform_units_temperature.R") +source("modules/R/transform_units_precipitation.R") +source("modules/R/transform_units_preasure.R") +Units <- function(recipe, data) { + # from recipe read the user defined units + # from data the original units + # deaccumulate option for CDS accumulated variables? + ## Do we need to convert other than ECVs? + user_units <- recipe$Analysis$Variables$units + var_names <- lapply(data, function(x) {x$attrs$Variable$varName}) + orig_units <- lapply(1:length(data), function(x) { + data[[x]]$attrs$Variable$metadata[[var_names[[x]]]]$units}) + # remove spaces, "**", "*" and "per" from units + user_units <- tolower(user_units) + user_units <- gsub(" ", "", user_units) + user_units <- gsub("**", "", user_units) + user_units <- gsub("*", "", user_units) + orig_units <- lapply(orig_units, function(x) { + x <- tolower(x) + x <- gsub(" ", "", x) + x <- gsub("**", "", x) + x <- gsub("*", "", x)}) + # if "/" appears substitute by -1 in at the end of next unit. How to know? + + if (all(orig_units == user_units)) { + # no transformation needed + res <- data + } else { + obj2trans <- which(orig_units != user_units) + res <- lapply(1:length(data), function(x) { + if (x %in% obj2trans) { + result <- units_transform(data[x], + orig_units = orig_units[x], + user_units = user_units) + } else { + result <- data[x] + } + return(result) + }) + } + return(res) +} -- GitLab From 581782a9e34a783f9b0e3a13c866bc5e7ae04227 Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 28 Jul 2023 17:41:25 +0200 Subject: [PATCH 02/21] example --- exec_units.R | 40 ++++++++++++++++++++++ recipe_tas_seasonal_units.yml | 63 +++++++++++++++++++++++++++++++++++ 2 files changed, 103 insertions(+) create mode 100644 exec_units.R create mode 100644 recipe_tas_seasonal_units.yml diff --git a/exec_units.R b/exec_units.R new file mode 100644 index 00000000..bc242e78 --- /dev/null +++ b/exec_units.R @@ -0,0 +1,40 @@ +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_tas_seasonal_units.yml" +recipe <- prepare_outputs(recipe_file) + +# Load datasets +data <- load_datasets(recipe) +# Units transformation +source("modules/Units/Units.R") +test <- Units(recipe, data) +# 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/recipe_tas_seasonal_units.yml b/recipe_tas_seasonal_units.yml new file mode 100644 index 00000000..4c1eccc1 --- /dev/null +++ b/recipe_tas_seasonal_units.yml @@ -0,0 +1,63 @@ +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 + units: C + 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: '20133' # 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/ -- GitLab From 759a0ea7c29f89850cc3c9f20547d498b5e66d70 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 31 Jul 2023 13:25:06 +0200 Subject: [PATCH 03/21] ongoing precipitation units --- exec_units.R | 4 +- modules/Loading/Loading.R | 40 ++++++------ modules/Units/R/transform_units_preasure.R | 3 +- .../Units/R/transform_units_precipitation.R | 62 ++++++++++++++---- modules/Units/R/transform_units_temperature.R | 15 +++-- modules/Units/R/units_transform.R | 22 ++++--- modules/Units/Units.R | 26 ++++---- recipe_prlr_seasonal_units.yml | 63 +++++++++++++++++++ recipe_tas_seasonal_units.yml | 10 +-- 9 files changed, 179 insertions(+), 66 deletions(-) create mode 100644 recipe_prlr_seasonal_units.yml diff --git a/exec_units.R b/exec_units.R index bc242e78..09d6e711 100644 --- a/exec_units.R +++ b/exec_units.R @@ -15,7 +15,9 @@ source("tools/prepare_outputs.R") #recipe_file <- args[1] #recipe <- read_atomic_recipe(recipe_file) ## to test a single recipe: -recipe_file <- "recipe_tas_seasonal_units.yml" + # recipe_file <- "recipe_tas_seasonal_units.yml" + # recipe_file <- "recipe_prlr_seasonal_units.yml" + recipe <- prepare_outputs(recipe_file) # Load datasets diff --git a/modules/Loading/Loading.R b/modules/Loading/Loading.R index d945d5af..8cc46dc7 100644 --- a/modules/Loading/Loading.R +++ b/modules/Loading/Loading.R @@ -365,26 +365,26 @@ load_datasets <- function(recipe) { # Convert prlr from m/s to mm/day ## TODO: Make a unit conversion function - if (variable[[var_idx]] == "prlr") { - # Verify that the units are m/s and the same in obs and hcst - if (((obs$attrs$Variable$metadata[[var_name]]$units == "m s-1") || - (obs$attrs$Variable$metadata[[var_name]]$units == "m s**-1")) && - ((hcst$attrs$Variable$metadata[[var_name]]$units == "m s-1") || - (hcst$attrs$Variable$metadata[[var_name]]$units == "m s**-1"))) { - info(recipe$Run$logger, "Converting precipitation from m/s to mm/day.") - obs$data[, var_idx, , , , , , , ] <- - obs$data[, var_idx, , , , , , , ]*86400*1000 - obs$attrs$Variable$metadata[[var_name]]$units <- "mm/day" - hcst$data[, var_idx, , , , , , , ] <- - hcst$data[, var_idx, , , , , , , ]*86400*1000 - hcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day" - if (!is.null(fcst)) { - fcst$data[, var_idx, , , , , , , ] <- - fcst$data[, var_idx, , , , , , , ]*86400*1000 - fcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day" - } - } - } +# if (variable[[var_idx]] == "prlr") { +# # Verify that the units are m/s and the same in obs and hcst +# if (((obs$attrs$Variable$metadata[[var_name]]$units == "m s-1") || +# (obs$attrs$Variable$metadata[[var_name]]$units == "m s**-1")) && +# ((hcst$attrs$Variable$metadata[[var_name]]$units == "m s-1") || +# (hcst$attrs$Variable$metadata[[var_name]]$units == "m s**-1"))) { +# info(recipe$Run$logger, "Converting precipitation from m/s to mm/day.") +# obs$data[, var_idx, , , , , , , ] <- +# obs$data[, var_idx, , , , , , , ]*86400*1000 +# obs$attrs$Variable$metadata[[var_name]]$units <- "mm/day" +# hcst$data[, var_idx, , , , , , , ] <- +# hcst$data[, var_idx, , , , , , , ]*86400*1000 +# hcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day" +# if (!is.null(fcst)) { +# fcst$data[, var_idx, , , , , , , ] <- +# fcst$data[, var_idx, , , , , , , ]*86400*1000 +# fcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day" +# } +# } +# } } # Compute anomalies if requested # Print a summary of the loaded data for the user, for each object diff --git a/modules/Units/R/transform_units_preasure.R b/modules/Units/R/transform_units_preasure.R index 30694114..01ca271f 100644 --- a/modules/Units/R/transform_units_preasure.R +++ b/modules/Units/R/transform_units_preasure.R @@ -1,2 +1,3 @@ -transform_units_preasure <- function(data, data, original_units, new_units){ +transform_units_preasure <- function(data, original_units, new_units) { + } diff --git a/modules/Units/R/transform_units_precipitation.R b/modules/Units/R/transform_units_precipitation.R index 81e900b0..56df5266 100644 --- a/modules/Units/R/transform_units_precipitation.R +++ b/modules/Units/R/transform_units_precipitation.R @@ -1,18 +1,54 @@ -transform_units_precipitation <- function(data, original_units, new_units){ - if (!substr(original_units, 1, 2) == substr(new_units, 1, 2)){ - if (substr(original_units, 1, 2) == 'mm' & (substr(new_units, 1, 2) == 'm ' | substr(new_units, 1, 2) == 'm/')){ - data <- data / 1000 - } else if ((substr(original_units, 1, 2) == 'm ' | substr(original_units, 1, 2) == 'm/') & substr(new_units, 1, 2) == 'mm'){ - data <- data * 1000 - } else {stop(paste0("not prepared to transform from ", original_units, " to ", new_units))} +transform_units_precipitation <- function(data, original_units, new_units, + var_name) { + if (original_units == "ms-1") { + if (new_units == "mm") { + data[[1]]$data <- data[[1]]$data * 3600 * 24 * 1000 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'mm' + } else if (new_units == "m") { + data[[1]]$data <- data[[1]]$data * 3600 * 24 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm' + } else if (new_units == "kgm-2s-1") { + data[[1]]$data <- data[[1]]$data * 1000 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'kg m-2 s-1' + } else { + stop(paste("Unknown transformation from", original_units, "to", + new_units)) } - if (('s' %in% strsplit(original_units, " |/|per") | 'second' %in% strsplit(original_units, " |/|per")) & 'month' %in% strsplit(new_units, " |/|per")){ - data <- data * 3600 * 24 * 365.25/12 - } else if ('day' %in% strsplit(original_units, " |/|per")[[1]] & 'month' %in% strsplit(new_units, " |/|per")[[1]]){ - data <- data * 365.25/12 - } else {stop(paste0("not prepared to transform from ", original_units, " to ", new_units))} + } else if (original_units == "mm") { + if (new_units == "ms-1") { + data[[1]]$data <- data[[1]]$data / (3600 * 24 * 1000) + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm s-1' + } else if (new_units == "m") { + data[[1]]$data <- data[[1]]$data / 1000 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'mm' + } else if (new_units == "kgm-2s-1") { +# data[[1]]$data <- data[[1]]$data / (3600 * 24 * 1000) + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'kg m-2 s-1' + } else { + stop(paste("Unknown transformation from", original_units, "to", + new_units)) + } + } else if (original_units == "m") { + if (new_units == "ms-1") { + } else if (new_units == "mm") { + } else if (new_units == "kgm-2s-1") { + } else { + stop(paste("Unknown transformation from", original_units, "to", + new_units)) + } + } else if (original_units == "kgm-2s-1") { + if (new_units == "ms-1") { + } else if (new_units == "mm") { + } else if (new_units == "kgm-2s-1") { + } else { + stop(paste("Unknown transformation from", original_units, "to", + new_units)) + } + } else { + stop("Unknown precipitation units transformation") + } - return(data) + return(data) } diff --git a/modules/Units/R/transform_units_temperature.R b/modules/Units/R/transform_units_temperature.R index 859ae479..09496b0f 100644 --- a/modules/Units/R/transform_units_temperature.R +++ b/modules/Units/R/transform_units_temperature.R @@ -1,10 +1,13 @@ -transform_units_temperature <- function(data, original_units, new_units){ - if (original_units == 'C' & new_units == 'K'){ - data <- data + 273.15 +transform_units_temperature <- function(data, original_units, new_units, + var_name) { + if (original_units == 'c' & new_units == 'k') { + data[[1]]$data <- data[[1]]$data + 273.15 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'K' } - if (original_units == 'K' & new_units == 'C'){ - data <- data - 273.15 + if (original_units == 'k' & new_units == 'c') { + data[[1]]$data <- data[[1]]$data - 273.15 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "C" + } return(data) } - diff --git a/modules/Units/R/units_transform.R b/modules/Units/R/units_transform.R index 088d3d0b..8cecdb6d 100644 --- a/modules/Units/R/units_transform.R +++ b/modules/Units/R/units_transform.R @@ -1,25 +1,31 @@ # data is a s2dv_cube object # units as character -units_transform <- function(data, orig_units, user_units) { - if (orig_units %in% c("C", "k")) { - if (user_units %in% c("C", "k")) { - trans <- transform_units_temperature(data, orig_units, user_units) +units_transform <- function(data, orig_units, user_units, var_name) { + var_name <- unlist(var_name) + if (orig_units %in% c("c", "k")) { + if (user_units %in% c("c", "k")) { + trans <- transform_units_temperature(data, orig_units, user_units, + var_name) } else { stop("Transformation temperature units not available.") - } if else (orig_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { + } + } else if (orig_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { if (user_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { - trans <- transform_units_precipitation(data, orig_units, user_units) + trans <- transform_units_precipitation(data, orig_units, user_units, + var_name) } else { stop("Transformation precipitation units not available.") } - } if else (orig_units %in% c("pa", "hpa", "mb")) { + } else if (orig_units %in% c("pa", "hpa", "mb")) { if (user_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { - trans <- transform_units_pressure(data, orig_units, user_units) + trans <- transform_units_pressure(data, orig_units, user_units, + var_name) } else { stop("Transformation precipitation units not available.") } } else { stop("Transformation unknown.") } + return(trans) } diff --git a/modules/Units/Units.R b/modules/Units/Units.R index 15a7e628..d4e057cf 100644 --- a/modules/Units/Units.R +++ b/modules/Units/Units.R @@ -1,9 +1,9 @@ # This function aims to convert units ## -source("modules/R/units_transform.R") -source("modules/R/transform_units_temperature.R") -source("modules/R/transform_units_precipitation.R") -source("modules/R/transform_units_preasure.R") +source("modules/Units/R/units_transform.R") +source("modules/Units/R/transform_units_temperature.R") +source("modules/Units/R/transform_units_precipitation.R") +source("modules/Units/R/transform_units_preasure.R") Units <- function(recipe, data) { # from recipe read the user defined units # from data the original units @@ -16,30 +16,32 @@ Units <- function(recipe, data) { # remove spaces, "**", "*" and "per" from units user_units <- tolower(user_units) user_units <- gsub(" ", "", user_units) - user_units <- gsub("**", "", user_units) - user_units <- gsub("*", "", user_units) + user_units <- gsub("\\**", "", user_units) + user_units <- gsub("\\*", "", user_units) orig_units <- lapply(orig_units, function(x) { x <- tolower(x) x <- gsub(" ", "", x) - x <- gsub("**", "", x) - x <- gsub("*", "", x)}) + x <- gsub("\\**", "", x) + x <- gsub("\\*", "", x)}) # if "/" appears substitute by -1 in at the end of next unit. How to know? - +browser() if (all(orig_units == user_units)) { # no transformation needed res <- data - } else { + } else { obj2trans <- which(orig_units != user_units) res <- lapply(1:length(data), function(x) { if (x %in% obj2trans) { result <- units_transform(data[x], - orig_units = orig_units[x], - user_units = user_units) + orig_units = orig_units[[x]], + user_units = user_units, + var_names[x]) } else { result <- data[x] } return(result) }) + } return(res) } diff --git a/recipe_prlr_seasonal_units.yml b/recipe_prlr_seasonal_units.yml new file mode 100644 index 00000000..ee68ad39 --- /dev/null +++ b/recipe_prlr_seasonal_units.yml @@ -0,0 +1,63 @@ +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 + units: mm + 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: '2014' # 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: 30 # Mandatory, int: minimum latitude + latmax: 50 # Mandatory, int: maximum latitude + lonmin: -10 # Mandatory, int: minimum longitude + lonmax: 10 # 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/recipe_tas_seasonal_units.yml b/recipe_tas_seasonal_units.yml index 4c1eccc1..c217b9fd 100644 --- a/recipe_tas_seasonal_units.yml +++ b/recipe_tas_seasonal_units.yml @@ -17,15 +17,15 @@ Analysis: Time: sdate: '0601' ## MMDD fcst_year: '2023' # Optional, int: Forecast year 'YYYY' - hcst_start: '20133' # Mandatory, int: Hindcast start year 'YYYY' + hcst_start: '2014' # 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 + latmin: 30 # Mandatory, int: minimum latitude + latmax: 50 # Mandatory, int: maximum latitude + lonmin: -10 # Mandatory, int: minimum longitude + lonmax: 10 # Mandatory, int: maximum longitude Regrid: method: bilinear # Mandatory, str: Interpolation method. See docu. type: "to_system" -- GitLab From 7df27415f8884b518e63ec6e1f834c6a76361cba Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 3 Aug 2023 13:06:08 +0200 Subject: [PATCH 04/21] Monthly precipitation --- modules/Units/R/transform_units_precipitation.R | 11 ++++++++++- modules/Units/Units.R | 16 +++++++++++----- recipe_prlr_seasonal_units.yml | 2 +- 3 files changed, 22 insertions(+), 7 deletions(-) diff --git a/modules/Units/R/transform_units_precipitation.R b/modules/Units/R/transform_units_precipitation.R index 56df5266..f57deca9 100644 --- a/modules/Units/R/transform_units_precipitation.R +++ b/modules/Units/R/transform_units_precipitation.R @@ -1,5 +1,6 @@ transform_units_precipitation <- function(data, original_units, new_units, var_name) { +browser() if (original_units == "ms-1") { if (new_units == "mm") { data[[1]]$data <- data[[1]]$data * 3600 * 24 * 1000 @@ -47,7 +48,15 @@ transform_units_precipitation <- function(data, original_units, new_units, } else { stop("Unknown precipitation units transformation") } - + # If monthly data compute the total monthly accumlated + if (recipe$Analysis$Variable$freq == "monthly_mean") { # could it not be mean? + data[[1]] <- Apply(list(data[[1]]$data, data[[1]]$attrs$Dates), + target_dim = 'time', fun = function(x, y) { + date <- as.Date(y, "%Y-%m-%d") + num_days <- lubridate::days_in_month(date) + res <- x * num_days + })$output1 + } return(data) } diff --git a/modules/Units/Units.R b/modules/Units/Units.R index d4e057cf..212c66f0 100644 --- a/modules/Units/Units.R +++ b/modules/Units/Units.R @@ -9,10 +9,14 @@ Units <- function(recipe, data) { # from data the original units # deaccumulate option for CDS accumulated variables? ## Do we need to convert other than ECVs? - user_units <- recipe$Analysis$Variables$units var_names <- lapply(data, function(x) {x$attrs$Variable$varName}) orig_units <- lapply(1:length(data), function(x) { data[[x]]$attrs$Variable$metadata[[var_names[[x]]]]$units}) + if (is.null(recipe$Analysis$Variables$units)) { + user_units <- orig_units[[which(!is.null(orig_units))[1]]] + } else { + user_units <- recipe$Analysis$Variables$units + } # remove spaces, "**", "*" and "per" from units user_units <- tolower(user_units) user_units <- gsub(" ", "", user_units) @@ -22,15 +26,17 @@ Units <- function(recipe, data) { x <- tolower(x) x <- gsub(" ", "", x) x <- gsub("\\**", "", x) - x <- gsub("\\*", "", x)}) + x <- gsub("\\*", "", x) + x <- unlist(ifelse(length(x)==0, user_units, x)) #when fcst is NULL + }) # if "/" appears substitute by -1 in at the end of next unit. How to know? -browser() + if (all(orig_units == user_units)) { # no transformation needed res <- data } else { obj2trans <- which(orig_units != user_units) - res <- lapply(1:length(data), function(x) { + res <- sapply(1:length(data), function(x) { if (x %in% obj2trans) { result <- units_transform(data[x], orig_units = orig_units[[x]], @@ -40,7 +46,7 @@ browser() result <- data[x] } return(result) - }) + }, simplify = TRUE) # instead of lapply to get the named list directly } return(res) diff --git a/recipe_prlr_seasonal_units.yml b/recipe_prlr_seasonal_units.yml index ee68ad39..e03428ac 100644 --- a/recipe_prlr_seasonal_units.yml +++ b/recipe_prlr_seasonal_units.yml @@ -16,7 +16,7 @@ Analysis: name: ERA5 # Mandatory, str: Reference codename. See docu. Time: sdate: '0601' ## MMDD - fcst_year: '2023' # Optional, int: Forecast year 'YYYY' + fcst_year: # Optional, int: Forecast year 'YYYY' hcst_start: '2014' # 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 -- GitLab From 054af4f4c5aee0f20ff2ba2184bfa6f4c1c59d96 Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 3 Aug 2023 13:53:39 +0200 Subject: [PATCH 05/21] review transform ppt --- .../Units/R/transform_units_precipitation.R | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/modules/Units/R/transform_units_precipitation.R b/modules/Units/R/transform_units_precipitation.R index f57deca9..00938ffd 100644 --- a/modules/Units/R/transform_units_precipitation.R +++ b/modules/Units/R/transform_units_precipitation.R @@ -1,6 +1,7 @@ transform_units_precipitation <- function(data, original_units, new_units, var_name) { browser() + if (original_units == "ms-1") { if (new_units == "mm") { data[[1]]$data <- data[[1]]$data * 3600 * 24 * 1000 @@ -21,9 +22,9 @@ browser() data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm s-1' } else if (new_units == "m") { data[[1]]$data <- data[[1]]$data / 1000 - data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'mm' + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm' } else if (new_units == "kgm-2s-1") { -# data[[1]]$data <- data[[1]]$data / (3600 * 24 * 1000) + data[[1]]$data <- data[[1]]$data / (3600 * 24 ) data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'kg m-2 s-1' } else { stop(paste("Unknown transformation from", original_units, "to", @@ -31,16 +32,28 @@ browser() } } else if (original_units == "m") { if (new_units == "ms-1") { + data[[1]]$data <- data[[1]]$data / (3600 * 24) + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm s-1' } else if (new_units == "mm") { + data[[1]]$data <- data[[1]]$data * 1000 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'mm' } else if (new_units == "kgm-2s-1") { + data[[1]]$data <- data[[1]]$data * 1000 / (3600 * 24) + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'kg m-2 s-1' } else { stop(paste("Unknown transformation from", original_units, "to", new_units)) } } else if (original_units == "kgm-2s-1") { if (new_units == "ms-1") { + data[[1]]$data <- data[[1]]$data / 1000 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm s-1' } else if (new_units == "mm") { - } else if (new_units == "kgm-2s-1") { + data[[1]]$data <- data[[1]]$data * 3600 * 24 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'mm' + } else if (new_units == "m") { + data[[1]]$data <- data[[1]]$data * 3600 * 24 / 1000 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm' } else { stop(paste("Unknown transformation from", original_units, "to", new_units)) -- GitLab From f0f2eff7d6fb35be7a7d19ba3fe3d302a7af080e Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 3 Aug 2023 14:08:41 +0200 Subject: [PATCH 06/21] pressure transform --- modules/Units/R/transform_units_preasure.R | 3 --- modules/Units/R/transform_units_pressure.R | 27 ++++++++++++++++++++++ modules/Units/R/units_transform.R | 2 +- modules/Units/Units.R | 3 ++- 4 files changed, 30 insertions(+), 5 deletions(-) delete mode 100644 modules/Units/R/transform_units_preasure.R create mode 100644 modules/Units/R/transform_units_pressure.R diff --git a/modules/Units/R/transform_units_preasure.R b/modules/Units/R/transform_units_preasure.R deleted file mode 100644 index 01ca271f..00000000 --- a/modules/Units/R/transform_units_preasure.R +++ /dev/null @@ -1,3 +0,0 @@ -transform_units_preasure <- function(data, original_units, new_units) { - -} diff --git a/modules/Units/R/transform_units_pressure.R b/modules/Units/R/transform_units_pressure.R new file mode 100644 index 00000000..6ebd39ae --- /dev/null +++ b/modules/Units/R/transform_units_pressure.R @@ -0,0 +1,27 @@ +transform_units_pressure <- function(data, original_units, new_units, var_name) { + + if (original_units == 'pa') { + if (new_units == 'hpa') { + data[[1]]$data <- data[[1]]$data /100 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'hPa' + } else if (new_units == 'mb') { + data[[1]]$data <- data[[1]]$data /100 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'mb' + } + } else if (original_units == 'hpa') { + if (new_units == 'pa') { + data[[1]]$data <- data[[1]]$data * 100 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "Pa" + } else if (new_units == "mb") { + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "mb" + } + } else if (original_units == "mb") { + if (new_units == 'pa') { + data[[1]]$data <- data[[1]]$data * 100 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "Pa" + } else if (new_units == "hPa") { + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "hPa" + } + } + return(data) +} diff --git a/modules/Units/R/units_transform.R b/modules/Units/R/units_transform.R index 8cecdb6d..10ab15cc 100644 --- a/modules/Units/R/units_transform.R +++ b/modules/Units/R/units_transform.R @@ -17,7 +17,7 @@ units_transform <- function(data, orig_units, user_units, var_name) { stop("Transformation precipitation units not available.") } } else if (orig_units %in% c("pa", "hpa", "mb")) { - if (user_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { + if (user_units %in% c("pa", "hpa", "mb")) { trans <- transform_units_pressure(data, orig_units, user_units, var_name) } else { diff --git a/modules/Units/Units.R b/modules/Units/Units.R index 212c66f0..b6afb44b 100644 --- a/modules/Units/Units.R +++ b/modules/Units/Units.R @@ -3,7 +3,7 @@ source("modules/Units/R/units_transform.R") source("modules/Units/R/transform_units_temperature.R") source("modules/Units/R/transform_units_precipitation.R") -source("modules/Units/R/transform_units_preasure.R") +source("modules/Units/R/transform_units_pressure.R") Units <- function(recipe, data) { # from recipe read the user defined units # from data the original units @@ -51,3 +51,4 @@ Units <- function(recipe, data) { } return(res) } + -- GitLab From 3c4ff96d709bc933523f425b2b4a603a2429b039 Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 3 Aug 2023 15:55:06 +0200 Subject: [PATCH 07/21] deaccumulate function --- modules/Units/R/deaccumulate.R | 6 ++++++ modules/Units/R/transform_units_precipitation.R | 6 ++---- modules/Units/R/units_transform.R | 4 ++-- modules/Units/Units.R | 3 ++- 4 files changed, 12 insertions(+), 7 deletions(-) create mode 100644 modules/Units/R/deaccumulate.R diff --git a/modules/Units/R/deaccumulate.R b/modules/Units/R/deaccumulate.R new file mode 100644 index 00000000..474b7382 --- /dev/null +++ b/modules/Units/R/deaccumulate.R @@ -0,0 +1,6 @@ +deaccumulate <- function(data) { + data[[1]]$data <- Apply(list(data[[1]]$data), target_dim = 'time', + fun = function(x) { + c(x[1], diff(x)) + })$output1 +} diff --git a/modules/Units/R/transform_units_precipitation.R b/modules/Units/R/transform_units_precipitation.R index 00938ffd..9de12c41 100644 --- a/modules/Units/R/transform_units_precipitation.R +++ b/modules/Units/R/transform_units_precipitation.R @@ -1,7 +1,5 @@ transform_units_precipitation <- function(data, original_units, new_units, - var_name) { -browser() - + var_name, freq) { if (original_units == "ms-1") { if (new_units == "mm") { data[[1]]$data <- data[[1]]$data * 3600 * 24 * 1000 @@ -62,7 +60,7 @@ browser() stop("Unknown precipitation units transformation") } # If monthly data compute the total monthly accumlated - if (recipe$Analysis$Variable$freq == "monthly_mean") { # could it not be mean? + if (freq == "monthly_mean") { # could it not be mean? data[[1]] <- Apply(list(data[[1]]$data, data[[1]]$attrs$Dates), target_dim = 'time', fun = function(x, y) { date <- as.Date(y, "%Y-%m-%d") diff --git a/modules/Units/R/units_transform.R b/modules/Units/R/units_transform.R index 10ab15cc..01008524 100644 --- a/modules/Units/R/units_transform.R +++ b/modules/Units/R/units_transform.R @@ -1,6 +1,6 @@ # data is a s2dv_cube object # units as character -units_transform <- function(data, orig_units, user_units, var_name) { +units_transform <- function(data, orig_units, user_units, var_name, freq) { var_name <- unlist(var_name) if (orig_units %in% c("c", "k")) { if (user_units %in% c("c", "k")) { @@ -12,7 +12,7 @@ units_transform <- function(data, orig_units, user_units, var_name) { } else if (orig_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { if (user_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { trans <- transform_units_precipitation(data, orig_units, user_units, - var_name) + var_name, freq) } else { stop("Transformation precipitation units not available.") } diff --git a/modules/Units/Units.R b/modules/Units/Units.R index b6afb44b..bd9f9045 100644 --- a/modules/Units/Units.R +++ b/modules/Units/Units.R @@ -10,6 +10,7 @@ Units <- function(recipe, data) { # deaccumulate option for CDS accumulated variables? ## Do we need to convert other than ECVs? var_names <- lapply(data, function(x) {x$attrs$Variable$varName}) + freq <- recipe$Analysis$Variable$freq orig_units <- lapply(1:length(data), function(x) { data[[x]]$attrs$Variable$metadata[[var_names[[x]]]]$units}) if (is.null(recipe$Analysis$Variables$units)) { @@ -41,7 +42,7 @@ Units <- function(recipe, data) { result <- units_transform(data[x], orig_units = orig_units[[x]], user_units = user_units, - var_names[x]) + var_names[x], freq = freq) } else { result <- data[x] } -- GitLab From 01fcad0ca0a2d60f1ab0349b475f580b16e552f0 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 21 Aug 2023 15:32:09 +0200 Subject: [PATCH 08/21] Include parameter ncores --- modules/Units/R/deaccumulate.R | 4 ++-- modules/Units/R/transform_units_precipitation.R | 14 +++++++------- modules/Units/R/units_transform.R | 5 +++-- modules/Units/Units.R | 4 +++- 4 files changed, 15 insertions(+), 12 deletions(-) diff --git a/modules/Units/R/deaccumulate.R b/modules/Units/R/deaccumulate.R index 474b7382..e1fcd6f5 100644 --- a/modules/Units/R/deaccumulate.R +++ b/modules/Units/R/deaccumulate.R @@ -1,6 +1,6 @@ -deaccumulate <- function(data) { +deaccumulate <- function(data, ncores = NULL) { data[[1]]$data <- Apply(list(data[[1]]$data), target_dim = 'time', fun = function(x) { c(x[1], diff(x)) - })$output1 + }, ncores = ncores)$output1 } diff --git a/modules/Units/R/transform_units_precipitation.R b/modules/Units/R/transform_units_precipitation.R index 9de12c41..20921147 100644 --- a/modules/Units/R/transform_units_precipitation.R +++ b/modules/Units/R/transform_units_precipitation.R @@ -1,5 +1,5 @@ transform_units_precipitation <- function(data, original_units, new_units, - var_name, freq) { + var_name, freq, ncores = NULL) { if (original_units == "ms-1") { if (new_units == "mm") { data[[1]]$data <- data[[1]]$data * 3600 * 24 * 1000 @@ -61,12 +61,12 @@ transform_units_precipitation <- function(data, original_units, new_units, } # If monthly data compute the total monthly accumlated if (freq == "monthly_mean") { # could it not be mean? - data[[1]] <- Apply(list(data[[1]]$data, data[[1]]$attrs$Dates), - target_dim = 'time', fun = function(x, y) { - date <- as.Date(y, "%Y-%m-%d") - num_days <- lubridate::days_in_month(date) - res <- x * num_days - })$output1 + data[[1]]$data <- Apply(list(data[[1]]$data, data[[1]]$attrs$Dates), + target_dim = 'time', fun = function(x, y) { + date <- as.Date(y, "%Y-%m-%d") + num_days <- lubridate::days_in_month(date) + res <- x * num_days + }, ncores = ncores)$output1 } return(data) } diff --git a/modules/Units/R/units_transform.R b/modules/Units/R/units_transform.R index 01008524..a9232647 100644 --- a/modules/Units/R/units_transform.R +++ b/modules/Units/R/units_transform.R @@ -1,6 +1,7 @@ # data is a s2dv_cube object # units as character -units_transform <- function(data, orig_units, user_units, var_name, freq) { +units_transform <- function(data, orig_units, user_units, var_name, freq, + ncores = NULL) { var_name <- unlist(var_name) if (orig_units %in% c("c", "k")) { if (user_units %in% c("c", "k")) { @@ -12,7 +13,7 @@ units_transform <- function(data, orig_units, user_units, var_name, freq) { } else if (orig_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { if (user_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { trans <- transform_units_precipitation(data, orig_units, user_units, - var_name, freq) + var_name, freq, ncores = ncores) } else { stop("Transformation precipitation units not available.") } diff --git a/modules/Units/Units.R b/modules/Units/Units.R index bd9f9045..386730f5 100644 --- a/modules/Units/Units.R +++ b/modules/Units/Units.R @@ -8,6 +8,7 @@ Units <- function(recipe, data) { # from recipe read the user defined units # from data the original units # deaccumulate option for CDS accumulated variables? + ncores <- recipe$Analysis$ncores ## Do we need to convert other than ECVs? var_names <- lapply(data, function(x) {x$attrs$Variable$varName}) freq <- recipe$Analysis$Variable$freq @@ -42,7 +43,8 @@ Units <- function(recipe, data) { result <- units_transform(data[x], orig_units = orig_units[[x]], user_units = user_units, - var_names[x], freq = freq) + var_names[x], freq = freq, + ncores = ncores) } else { result <- data[x] } -- GitLab From b9da976b9c36c280aa349eb61d136c7dd82424c0 Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 24 Aug 2023 17:28:23 +0200 Subject: [PATCH 09/21] flux parameter in recipe --- exec_units.R | 2 +- .../Units/R/transform_units_precipitation.R | 47 +++++++++++++++---- modules/Units/R/units_transform.R | 4 +- modules/Units/Units.R | 11 +++++ recipe_tas_seasonal_units.yml | 4 +- 5 files changed, 54 insertions(+), 14 deletions(-) diff --git a/exec_units.R b/exec_units.R index 09d6e711..fa95767b 100644 --- a/exec_units.R +++ b/exec_units.R @@ -26,7 +26,7 @@ data <- load_datasets(recipe) source("modules/Units/Units.R") test <- Units(recipe, data) # Calibrate datasets -data <- calibrate_datasets(recipe, data) +data <- calibrate_datasets(recipe, test) # Compute skill metrics skill_metrics <- compute_skill_metrics(recipe, data) # Compute percentiles and probability bins diff --git a/modules/Units/R/transform_units_precipitation.R b/modules/Units/R/transform_units_precipitation.R index 20921147..bc840d15 100644 --- a/modules/Units/R/transform_units_precipitation.R +++ b/modules/Units/R/transform_units_precipitation.R @@ -1,5 +1,7 @@ transform_units_precipitation <- function(data, original_units, new_units, - var_name, freq, ncores = NULL) { + var_name, freq, flux = FALSE, ncores = NULL) { +## TODO consider higher frequencies (e.g. 6hourly) +## could create a constant hours <- 24 or hours <- 6 and use the same code if (original_units == "ms-1") { if (new_units == "mm") { data[[1]]$data <- data[[1]]$data * 3600 * 24 * 1000 @@ -59,16 +61,43 @@ transform_units_precipitation <- function(data, original_units, new_units, } else { stop("Unknown precipitation units transformation") } - # If monthly data compute the total monthly accumlated - if (freq == "monthly_mean") { # could it not be mean? - data[[1]]$data <- Apply(list(data[[1]]$data, data[[1]]$attrs$Dates), - target_dim = 'time', fun = function(x, y) { - date <- as.Date(y, "%Y-%m-%d") - num_days <- lubridate::days_in_month(date) - res <- x * num_days - }, ncores = ncores)$output1 + if (flux) { + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- paste0( + data[[1]]$attrs$Variable$metadata[[var_name]]$units, "/day") + } else { + if (freq == "monthly_mean") { # could it not be mean? + time_pos <- which(lapply(data[[1]]$attrs$Variable$metadata[[var_name]]$dim, + function(x) {x$name}) == 'time') + cal <- tolower(data[[1]]$attrs$Variable$metadata[[var_name]]$dim[[time_pos]]$calendar) + data[[1]]$data <- Apply(list(data[[1]]$data, data[[1]]$attrs$Dates), + target_dim = 'time', fun = function(x, y) { + date <- as.Date(y, "%Y-%m-%d") + num_days <- .days_in_month(date, cal = cal) + res <- x * num_days + }, ncores = ncores)$output1 + } } return(data) } +.days_in_month <- function (x, cal) { + if (cal %in% c('gregorian', 'standard', 'proleptic_gregorian')) { + N_DAYS_IN_MONTHS <- lubridate:::N_DAYS_IN_MONTHS + if (leap_year(year(x))) { + N_DAYS_IN_MONTHS[2] <- N_DAYS_IN_MONTHS[2] + 1 + } + } else if (cal %in% c('360', '360_day')) { + N_DAYS_IN_MONTHS <- rep(30, 12) + names(N_DAYS_IN_MONTHS) <- month.abb + } else if (cal %in% c('365', '365_day')) { + N_DAYS_IN_MONTHS <- lubridate:::N_DAYS_IN_MONTHS + } else { + stop("Unknown calendar") + } + + month_x <- month(x, label = TRUE, locale = "C") + n_days <- N_DAYS_IN_MONTHS[month_x] + n_days[month_x == "Feb" & leap_year(x)] <- 29L + n_days +} diff --git a/modules/Units/R/units_transform.R b/modules/Units/R/units_transform.R index a9232647..fd435fa8 100644 --- a/modules/Units/R/units_transform.R +++ b/modules/Units/R/units_transform.R @@ -1,6 +1,6 @@ # data is a s2dv_cube object # units as character -units_transform <- function(data, orig_units, user_units, var_name, freq, +units_transform <- function(data, orig_units, user_units, var_name, freq, flux = FALSE, ncores = NULL) { var_name <- unlist(var_name) if (orig_units %in% c("c", "k")) { @@ -13,7 +13,7 @@ units_transform <- function(data, orig_units, user_units, var_name, freq, } else if (orig_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { if (user_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { trans <- transform_units_precipitation(data, orig_units, user_units, - var_name, freq, ncores = ncores) + var_name, freq, flux, ncores = ncores) } else { stop("Transformation precipitation units not available.") } diff --git a/modules/Units/Units.R b/modules/Units/Units.R index 386730f5..5e7ea88c 100644 --- a/modules/Units/Units.R +++ b/modules/Units/Units.R @@ -4,6 +4,11 @@ source("modules/Units/R/units_transform.R") source("modules/Units/R/transform_units_temperature.R") source("modules/Units/R/transform_units_precipitation.R") source("modules/Units/R/transform_units_pressure.R") +# Requires recipe to include: +# recipe$Analysis$Variables: +# $freq +# $flux (if precipitation) +# $units (optional) Units <- function(recipe, data) { # from recipe read the user defined units # from data the original units @@ -12,6 +17,11 @@ Units <- function(recipe, data) { ## Do we need to convert other than ECVs? var_names <- lapply(data, function(x) {x$attrs$Variable$varName}) freq <- recipe$Analysis$Variable$freq + if (is.null(recipe$Analysis$Variable$flux)) { + flux <- FALSE + } else { + flux <- recipe$Analysis$Variable$flux + } orig_units <- lapply(1:length(data), function(x) { data[[x]]$attrs$Variable$metadata[[var_names[[x]]]]$units}) if (is.null(recipe$Analysis$Variables$units)) { @@ -44,6 +54,7 @@ Units <- function(recipe, data) { orig_units = orig_units[[x]], user_units = user_units, var_names[x], freq = freq, + flux = flux, ncores = ncores) } else { result <- data[x] diff --git a/recipe_tas_seasonal_units.yml b/recipe_tas_seasonal_units.yml index c217b9fd..d2b25321 100644 --- a/recipe_tas_seasonal_units.yml +++ b/recipe_tas_seasonal_units.yml @@ -7,7 +7,7 @@ Analysis: Variables: name: tas freq: monthly_mean - units: C + units: K Datasets: System: name: ECMWF-SEAS5.1 # Mandatory, str: system5c3s system21_m1 system35c3s @@ -17,7 +17,7 @@ Analysis: Time: sdate: '0601' ## MMDD fcst_year: '2023' # Optional, int: Forecast year 'YYYY' - hcst_start: '2014' # Mandatory, int: Hindcast start year 'YYYY' + hcst_start: '2009' # 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 -- GitLab From 15a25f013bc2a0ede30d1da7db4a1b220b015565 Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 5 Sep 2023 11:51:01 +0200 Subject: [PATCH 10/21] Improve error message --- modules/Units/R/units_transform.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/Units/R/units_transform.R b/modules/Units/R/units_transform.R index fd435fa8..694bd3a7 100644 --- a/modules/Units/R/units_transform.R +++ b/modules/Units/R/units_transform.R @@ -25,7 +25,7 @@ units_transform <- function(data, orig_units, user_units, var_name, freq, flux = stop("Transformation precipitation units not available.") } } else { - stop("Transformation unknown.") + stop(paste("Transformation unknown from", orig_units, "to", user_units)) } return(trans) } -- GitLab From aba972a89fa9402b5edbd23f3ce4b4229597cdb8 Mon Sep 17 00:00:00 2001 From: Victoria Agudetse Roures Date: Tue, 5 Sep 2023 14:45:44 +0200 Subject: [PATCH 11/21] Adapt Units module to multivar case --- .../Units/R/transform_units_precipitation.R | 66 ++++++++++++------- modules/Units/R/transform_units_pressure.R | 15 +++-- modules/Units/R/transform_units_temperature.R | 23 ++++--- modules/Units/R/units_transform.R | 54 ++++++++------- modules/Units/Units.R | 28 +++++--- .../atomic_recipes/recipe_test_multivar.yml | 1 + 6 files changed, 113 insertions(+), 74 deletions(-) diff --git a/modules/Units/R/transform_units_precipitation.R b/modules/Units/R/transform_units_precipitation.R index bc840d15..afe29d4f 100644 --- a/modules/Units/R/transform_units_precipitation.R +++ b/modules/Units/R/transform_units_precipitation.R @@ -1,16 +1,20 @@ transform_units_precipitation <- function(data, original_units, new_units, - var_name, freq, flux = FALSE, ncores = NULL) { + var_name, freq, flux = FALSE, ncores = NULL, + var_index = 1) { ## TODO consider higher frequencies (e.g. 6hourly) ## could create a constant hours <- 24 or hours <- 6 and use the same code if (original_units == "ms-1") { if (new_units == "mm") { - data[[1]]$data <- data[[1]]$data * 3600 * 24 * 1000 + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] * 3600 * 24 * 1000 data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'mm' } else if (new_units == "m") { - data[[1]]$data <- data[[1]]$data * 3600 * 24 + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] * 3600 * 24 data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm' } else if (new_units == "kgm-2s-1") { - data[[1]]$data <- data[[1]]$data * 1000 + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] * 1000 data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'kg m-2 s-1' } else { stop(paste("Unknown transformation from", original_units, "to", @@ -18,13 +22,16 @@ transform_units_precipitation <- function(data, original_units, new_units, } } else if (original_units == "mm") { if (new_units == "ms-1") { - data[[1]]$data <- data[[1]]$data / (3600 * 24 * 1000) + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] / (3600 * 24 * 1000) data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm s-1' } else if (new_units == "m") { - data[[1]]$data <- data[[1]]$data / 1000 + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] / 1000 data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm' } else if (new_units == "kgm-2s-1") { - data[[1]]$data <- data[[1]]$data / (3600 * 24 ) + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] / (3600 * 24 ) data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'kg m-2 s-1' } else { stop(paste("Unknown transformation from", original_units, "to", @@ -32,13 +39,16 @@ transform_units_precipitation <- function(data, original_units, new_units, } } else if (original_units == "m") { if (new_units == "ms-1") { - data[[1]]$data <- data[[1]]$data / (3600 * 24) + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] / (3600 * 24) data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm s-1' } else if (new_units == "mm") { - data[[1]]$data <- data[[1]]$data * 1000 + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] * 1000 data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'mm' } else if (new_units == "kgm-2s-1") { - data[[1]]$data <- data[[1]]$data * 1000 / (3600 * 24) + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] * 1000 / (3600 * 24) data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'kg m-2 s-1' } else { stop(paste("Unknown transformation from", original_units, "to", @@ -46,13 +56,16 @@ transform_units_precipitation <- function(data, original_units, new_units, } } else if (original_units == "kgm-2s-1") { if (new_units == "ms-1") { - data[[1]]$data <- data[[1]]$data / 1000 + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] / 1000 data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm s-1' } else if (new_units == "mm") { - data[[1]]$data <- data[[1]]$data * 3600 * 24 + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] * 3600 * 24 data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'mm' } else if (new_units == "m") { - data[[1]]$data <- data[[1]]$data * 3600 * 24 / 1000 + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] * 3600 * 24 / 1000 data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'm' } else { stop(paste("Unknown transformation from", original_units, "to", @@ -65,21 +78,27 @@ transform_units_precipitation <- function(data, original_units, new_units, data[[1]]$attrs$Variable$metadata[[var_name]]$units <- paste0( data[[1]]$attrs$Variable$metadata[[var_name]]$units, "/day") } else { + ## TODO: Shouldn't use time dimension, need it for Compute() if (freq == "monthly_mean") { # could it not be mean? time_pos <- which(lapply(data[[1]]$attrs$Variable$metadata[[var_name]]$dim, function(x) {x$name}) == 'time') cal <- tolower(data[[1]]$attrs$Variable$metadata[[var_name]]$dim[[time_pos]]$calendar) - data[[1]]$data <- Apply(list(data[[1]]$data, data[[1]]$attrs$Dates), - target_dim = 'time', fun = function(x, y) { - date <- as.Date(y, "%Y-%m-%d") - num_days <- .days_in_month(date, cal = cal) - res <- x * num_days - }, ncores = ncores)$output1 + data_subset <- Subset(data[[1]]$data, along = "var", indices = var_index, drop = 'selected') + data[[1]]$data[ , var_index, , , , , , , ] <- + Apply(list(data_subset, data[[1]]$attrs$Dates), + target_dim = list(c('time'), c('time')), + extra_info = list(cal = cal, days_in_month = .days_in_month), + fun = function(x, y) { + date <- as.Date(y, "%Y-%m-%d") + num_days <- .days_in_month(date, cal = .cal) + res <- x * num_days + }, ncores = ncores)$output1 } } return(data) } -.days_in_month <- function (x, cal) { + +.days_in_month <- function(x, cal) { if (cal %in% c('gregorian', 'standard', 'proleptic_gregorian')) { N_DAYS_IN_MONTHS <- lubridate:::N_DAYS_IN_MONTHS if (leap_year(year(x))) { @@ -92,12 +111,9 @@ transform_units_precipitation <- function(data, original_units, new_units, N_DAYS_IN_MONTHS <- lubridate:::N_DAYS_IN_MONTHS } else { stop("Unknown calendar") - } - + } month_x <- month(x, label = TRUE, locale = "C") n_days <- N_DAYS_IN_MONTHS[month_x] n_days[month_x == "Feb" & leap_year(x)] <- 29L - n_days + return(n_days) } - - diff --git a/modules/Units/R/transform_units_pressure.R b/modules/Units/R/transform_units_pressure.R index 6ebd39ae..db3f88f1 100644 --- a/modules/Units/R/transform_units_pressure.R +++ b/modules/Units/R/transform_units_pressure.R @@ -1,23 +1,28 @@ -transform_units_pressure <- function(data, original_units, new_units, var_name) { +transform_units_pressure <- function(data, original_units, new_units, var_name, + var_index = 1) { if (original_units == 'pa') { if (new_units == 'hpa') { - data[[1]]$data <- data[[1]]$data /100 + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] /100 data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'hPa' } else if (new_units == 'mb') { - data[[1]]$data <- data[[1]]$data /100 + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] /100 data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'mb' } } else if (original_units == 'hpa') { if (new_units == 'pa') { - data[[1]]$data <- data[[1]]$data * 100 + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] * 100 data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "Pa" } else if (new_units == "mb") { data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "mb" } } else if (original_units == "mb") { if (new_units == 'pa') { - data[[1]]$data <- data[[1]]$data * 100 + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] * 100 data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "Pa" } else if (new_units == "hPa") { data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "hPa" diff --git a/modules/Units/R/transform_units_temperature.R b/modules/Units/R/transform_units_temperature.R index 09496b0f..bcd27524 100644 --- a/modules/Units/R/transform_units_temperature.R +++ b/modules/Units/R/transform_units_temperature.R @@ -1,13 +1,16 @@ transform_units_temperature <- function(data, original_units, new_units, - var_name) { - if (original_units == 'c' & new_units == 'k') { - data[[1]]$data <- data[[1]]$data + 273.15 - data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'K' - } - if (original_units == 'k' & new_units == 'c') { - data[[1]]$data <- data[[1]]$data - 273.15 - data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "C" + var_name, var_index = 1, + var_dim = "var") { + if (original_units == 'c' & new_units == 'k') { + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] + 273.15 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- 'K' + } + if (original_units == 'k' & new_units == 'c') { + data[[1]]$data[ , var_index, , , , , , , ] <- + data[[1]]$data[ , var_index, , , , , , , ] - 273.15 + data[[1]]$attrs$Variable$metadata[[var_name]]$units <- "C" - } - return(data) + } + return(data) } diff --git a/modules/Units/R/units_transform.R b/modules/Units/R/units_transform.R index 694bd3a7..3da48638 100644 --- a/modules/Units/R/units_transform.R +++ b/modules/Units/R/units_transform.R @@ -2,31 +2,37 @@ # units as character units_transform <- function(data, orig_units, user_units, var_name, freq, flux = FALSE, ncores = NULL) { - var_name <- unlist(var_name) - if (orig_units %in% c("c", "k")) { - if (user_units %in% c("c", "k")) { - trans <- transform_units_temperature(data, orig_units, user_units, - var_name) - } else { - stop("Transformation temperature units not available.") + ## TODO: Change how argument 'flux' works + for (i in 1:length(var_name)) { + if (!(orig_units[i] == user_units[i])) { + if (orig_units[i] %in% c("c", "k")) { + if (user_units[i] %in% c("c", "k")) { + data <- transform_units_temperature(data, orig_units[i], user_units[i], + var_name[i], var_index = i) + } else { + stop("Transformation temperature units not available.") + } + } else if (orig_units[i] %in% c("ms-1", "kgm-2s-1", "mm", "m")) { + if (user_units[i] %in% c("ms-1", "kgm-2s-1", "mm", "m")) { + print("Converting precip units") + data <- transform_units_precipitation(data, orig_units[i], user_units[i], + var_name[i], freq, flux, + ncores = ncores, var_index = i) + } else { + stop("Transformation precipitation units not available.") + } + } else if (orig_units[i] %in% c("pa", "hpa", "mb")) { + if (user_units[i] %in% c("pa", "hpa", "mb")) { + data <- transform_units_pressure(data, orig_units[i], user_units[i], + var_name[i]) + } else { + stop("Transformation precipitation units not available.") + } + } else { + stop(paste("Transformation unknown from", orig_units[i], "to", user_units[i])) + } } - } else if (orig_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { - if (user_units %in% c("ms-1", "kgm-2s-1", "mm", "m")) { - trans <- transform_units_precipitation(data, orig_units, user_units, - var_name, freq, flux, ncores = ncores) - } else { - stop("Transformation precipitation units not available.") - } - } else if (orig_units %in% c("pa", "hpa", "mb")) { - if (user_units %in% c("pa", "hpa", "mb")) { - trans <- transform_units_pressure(data, orig_units, user_units, - var_name) - } else { - stop("Transformation precipitation units not available.") - } - } else { - stop(paste("Transformation unknown from", orig_units, "to", user_units)) } - return(trans) + return(data) } diff --git a/modules/Units/Units.R b/modules/Units/Units.R index 5e7ea88c..07b830d5 100644 --- a/modules/Units/Units.R +++ b/modules/Units/Units.R @@ -22,13 +22,21 @@ Units <- function(recipe, data) { } else { flux <- recipe$Analysis$Variable$flux } - orig_units <- lapply(1:length(data), function(x) { - data[[x]]$attrs$Variable$metadata[[var_names[[x]]]]$units}) + orig_units <- list() + for (element in names(var_names)) { + orig_units[[element]] <- c() + for (x in var_names[[element]]) { + orig_units[[element]] <- c(orig_units[[element]], + data[[element]]$attrs$Variable$metadata[[x]]$units) + } + } if (is.null(recipe$Analysis$Variables$units)) { user_units <- orig_units[[which(!is.null(orig_units))[1]]] } else { user_units <- recipe$Analysis$Variables$units } + ## TODO: How to handle spaces in multi-var case? + user_units <- strsplit(recipe$Analysis$Variables$units, ", | |,")[[1]] # remove spaces, "**", "*" and "per" from units user_units <- tolower(user_units) user_units <- gsub(" ", "", user_units) @@ -39,21 +47,21 @@ Units <- function(recipe, data) { x <- gsub(" ", "", x) x <- gsub("\\**", "", x) x <- gsub("\\*", "", x) - x <- unlist(ifelse(length(x)==0, user_units, x)) #when fcst is NULL + x <- unlist(ifelse(rep(length(x)==0, length(user_units)), user_units, x)) #when fcst is NULL }) - # if "/" appears substitute by -1 in at the end of next unit. How to know? - - if (all(orig_units == user_units)) { + ## TODO: + ## if "/" appears substitute by -1 in at the end of next unit. How to know? + identicalValue <- function(x,y) if (identical(x,y)) TRUE else FALSE + if (Reduce(identicalValue, c(orig_units, user_units))) { # no transformation needed res <- data } else { - obj2trans <- which(orig_units != user_units) res <- sapply(1:length(data), function(x) { - if (x %in% obj2trans) { + if (!all(orig_units[x] == user_units)) { result <- units_transform(data[x], orig_units = orig_units[[x]], user_units = user_units, - var_names[x], freq = freq, + var_names[[x]], freq = freq, flux = flux, ncores = ncores) } else { @@ -61,7 +69,7 @@ Units <- function(recipe, data) { } return(result) }, simplify = TRUE) # instead of lapply to get the named list directly - + info(recipe$Run$logger, "##### UNIT CONVERSION COMPLETE #####") } return(res) } diff --git a/recipes/atomic_recipes/recipe_test_multivar.yml b/recipes/atomic_recipes/recipe_test_multivar.yml index 8eb3e962..9c626a51 100644 --- a/recipes/atomic_recipes/recipe_test_multivar.yml +++ b/recipes/atomic_recipes/recipe_test_multivar.yml @@ -6,6 +6,7 @@ Analysis: Variables: name: tas prlr freq: monthly_mean + units: C mm Datasets: System: name: ECMWF-SEAS5 -- GitLab From abf6675dc963482990e6669ffb87de6aa9dee5f4 Mon Sep 17 00:00:00 2001 From: Victoria Agudetse Roures Date: Tue, 5 Sep 2023 14:49:47 +0200 Subject: [PATCH 12/21] Add a TODO --- modules/Units/R/transform_units_precipitation.R | 1 + modules/Units/R/transform_units_pressure.R | 1 + modules/Units/R/transform_units_temperature.R | 1 + 3 files changed, 3 insertions(+) diff --git a/modules/Units/R/transform_units_precipitation.R b/modules/Units/R/transform_units_precipitation.R index afe29d4f..43d00a6a 100644 --- a/modules/Units/R/transform_units_precipitation.R +++ b/modules/Units/R/transform_units_precipitation.R @@ -1,6 +1,7 @@ transform_units_precipitation <- function(data, original_units, new_units, var_name, freq, flux = FALSE, ncores = NULL, var_index = 1) { + ## TODO: Hard-coded subsetting ## TODO consider higher frequencies (e.g. 6hourly) ## could create a constant hours <- 24 or hours <- 6 and use the same code if (original_units == "ms-1") { diff --git a/modules/Units/R/transform_units_pressure.R b/modules/Units/R/transform_units_pressure.R index db3f88f1..9712e9fe 100644 --- a/modules/Units/R/transform_units_pressure.R +++ b/modules/Units/R/transform_units_pressure.R @@ -1,6 +1,7 @@ transform_units_pressure <- function(data, original_units, new_units, var_name, var_index = 1) { + ## TODO: Hard-coded subsetting if (original_units == 'pa') { if (new_units == 'hpa') { data[[1]]$data[ , var_index, , , , , , , ] <- diff --git a/modules/Units/R/transform_units_temperature.R b/modules/Units/R/transform_units_temperature.R index bcd27524..366f0d34 100644 --- a/modules/Units/R/transform_units_temperature.R +++ b/modules/Units/R/transform_units_temperature.R @@ -1,6 +1,7 @@ transform_units_temperature <- function(data, original_units, new_units, var_name, var_index = 1, var_dim = "var") { + ## TODO: Hard-coded subsetting if (original_units == 'c' & new_units == 'k') { data[[1]]$data[ , var_index, , , , , , , ] <- data[[1]]$data[ , var_index, , , , , , , ] + 273.15 -- GitLab From bddc54dfe86cfb67f56fc77e364dc79111e70754 Mon Sep 17 00:00:00 2001 From: Victoria Agudetse Roures Date: Tue, 5 Sep 2023 14:58:43 +0200 Subject: [PATCH 13/21] Display units in data summary --- tools/data_summary.R | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tools/data_summary.R b/tools/data_summary.R index 714c5a0b..11f365cd 100644 --- a/tools/data_summary.R +++ b/tools/data_summary.R @@ -34,8 +34,11 @@ data_summary <- function(data_cube, recipe) { info(recipe$Run$logger, paste0("Statistical summary of the data in ", object_name, ":")) for (var_index in 1:data_cube$dims[['var']]) { + variable_name <- data_cube$attrs$Variable$varName[var_index] + variable_units <- data_cube$attrs$Variable$metadata[[variable_name]]$units info(recipe$Run$logger, - paste("Variable:", data_cube$attrs$Variable$varName[var_index])) + paste0("Variable: ", variable_name, + " (units: ", variable_units, ")")) output_string <- capture.output(summary(Subset(data_cube$data, along = "var", indices = var_index))) -- GitLab From ad5e156cb066a6231c9a05fb04b7a40a22ef9de4 Mon Sep 17 00:00:00 2001 From: Victoria Agudetse Roures Date: Tue, 5 Sep 2023 16:26:47 +0200 Subject: [PATCH 14/21] Remove printed message --- modules/Units/R/units_transform.R | 1 - 1 file changed, 1 deletion(-) diff --git a/modules/Units/R/units_transform.R b/modules/Units/R/units_transform.R index 3da48638..83f5f700 100644 --- a/modules/Units/R/units_transform.R +++ b/modules/Units/R/units_transform.R @@ -14,7 +14,6 @@ units_transform <- function(data, orig_units, user_units, var_name, freq, flux = } } else if (orig_units[i] %in% c("ms-1", "kgm-2s-1", "mm", "m")) { if (user_units[i] %in% c("ms-1", "kgm-2s-1", "mm", "m")) { - print("Converting precip units") data <- transform_units_precipitation(data, orig_units[i], user_units[i], var_name[i], freq, flux, ncores = ncores, var_index = i) -- GitLab From 36b7be80fd5498b76d2683b88ef445d7dce559fb Mon Sep 17 00:00:00 2001 From: Victoria Agudetse Roures Date: Tue, 5 Sep 2023 16:37:17 +0200 Subject: [PATCH 15/21] Change margin dimension in transform_units_precipitation() from 'time' to 'syear' --- modules/Units/R/transform_units_precipitation.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/Units/R/transform_units_precipitation.R b/modules/Units/R/transform_units_precipitation.R index 43d00a6a..a24ef851 100644 --- a/modules/Units/R/transform_units_precipitation.R +++ b/modules/Units/R/transform_units_precipitation.R @@ -87,7 +87,7 @@ transform_units_precipitation <- function(data, original_units, new_units, data_subset <- Subset(data[[1]]$data, along = "var", indices = var_index, drop = 'selected') data[[1]]$data[ , var_index, , , , , , , ] <- Apply(list(data_subset, data[[1]]$attrs$Dates), - target_dim = list(c('time'), c('time')), + target_dim = list(c('syear'), c('syear')), extra_info = list(cal = cal, days_in_month = .days_in_month), fun = function(x, y) { date <- as.Date(y, "%Y-%m-%d") -- GitLab From 10047b33e3b67270623f3f080971dda21a134916 Mon Sep 17 00:00:00 2001 From: Victoria Agudetse Roures Date: Fri, 8 Sep 2023 15:29:46 +0200 Subject: [PATCH 16/21] Handle multivar units in recipe, fix bug in unit equality assessment logic --- modules/Units/Units.R | 14 +++++++++----- recipes/atomic_recipes/recipe_test_multivar.yml | 6 +++--- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/modules/Units/Units.R b/modules/Units/Units.R index 07b830d5..f2e8413f 100644 --- a/modules/Units/Units.R +++ b/modules/Units/Units.R @@ -35,8 +35,12 @@ Units <- function(recipe, data) { } else { user_units <- recipe$Analysis$Variables$units } + if (!is.null(names(user_units))) { + # Ensure that the units are in the correct order + user_units <- unlist(user_units[var_names[[1]]]) + } ## TODO: How to handle spaces in multi-var case? - user_units <- strsplit(recipe$Analysis$Variables$units, ", | |,")[[1]] + # user_units <- strsplit(recipe$Analysis$Variables$units, ", | |,")[[1]] # remove spaces, "**", "*" and "per" from units user_units <- tolower(user_units) user_units <- gsub(" ", "", user_units) @@ -51,9 +55,10 @@ Units <- function(recipe, data) { }) ## TODO: ## if "/" appears substitute by -1 in at the end of next unit. How to know? - identicalValue <- function(x,y) if (identical(x,y)) TRUE else FALSE - if (Reduce(identicalValue, c(orig_units, user_units))) { - # no transformation needed + # Check if all units are equal (if so, no conversion needed; else convert) + browser() + if (length(unique(c(unlist(orig_units), user_units))) == length(user_units)) { + info(recipe$Run$logger, "##### NO UNIT CONVERSION NEEDED #####") res <- data } else { res <- sapply(1:length(data), function(x) { @@ -73,4 +78,3 @@ Units <- function(recipe, data) { } return(res) } - diff --git a/recipes/atomic_recipes/recipe_test_multivar.yml b/recipes/atomic_recipes/recipe_test_multivar.yml index 9c626a51..a8c2e71b 100644 --- a/recipes/atomic_recipes/recipe_test_multivar.yml +++ b/recipes/atomic_recipes/recipe_test_multivar.yml @@ -4,9 +4,9 @@ Description: Analysis: Horizon: Seasonal Variables: - name: tas prlr + name: tas, prlr freq: monthly_mean - units: C mm + units: {tas: C, prlr: mm} Datasets: System: name: ECMWF-SEAS5 @@ -19,7 +19,7 @@ Analysis: hcst_start: '1999' hcst_end: '2010' ftime_min: 1 - ftime_max: 2 + ftime_max: 1 Region: latmin: -10 latmax: 10 -- GitLab From d357047e82392771e7ea810234efb9f433ed95c7 Mon Sep 17 00:00:00 2001 From: Victoria Agudetse Roures Date: Fri, 8 Sep 2023 15:43:55 +0200 Subject: [PATCH 17/21] Add warning in the case that filesystem: 'esarchive' and the units between s2dv_cube objects are not consistent --- modules/Units/Units.R | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/modules/Units/Units.R b/modules/Units/Units.R index f2e8413f..137ad56b 100644 --- a/modules/Units/Units.R +++ b/modules/Units/Units.R @@ -56,11 +56,16 @@ Units <- function(recipe, data) { ## TODO: ## if "/" appears substitute by -1 in at the end of next unit. How to know? # Check if all units are equal (if so, no conversion needed; else convert) - browser() if (length(unique(c(unlist(orig_units), user_units))) == length(user_units)) { info(recipe$Run$logger, "##### NO UNIT CONVERSION NEEDED #####") res <- data } else { + if (recipe$Run$filesystem == 'esarchive' && + !(length(unique(orig_units)) == 1)) { + warn(recipe$Run$logger, + paste("The units in", names(orig_units), "are not all equal.", + "If this is not expected, please contact the ES data team.")) + } res <- sapply(1:length(data), function(x) { if (!all(orig_units[x] == user_units)) { result <- units_transform(data[x], -- GitLab From 537866bab24a0914129ba56c26d5314a9085c416 Mon Sep 17 00:00:00 2001 From: Victoria Agudetse Roures Date: Fri, 8 Sep 2023 15:58:36 +0200 Subject: [PATCH 18/21] Fix indentation --- modules/Units/R/units_transform.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/Units/R/units_transform.R b/modules/Units/R/units_transform.R index 83f5f700..11e05f31 100644 --- a/modules/Units/R/units_transform.R +++ b/modules/Units/R/units_transform.R @@ -28,7 +28,7 @@ units_transform <- function(data, orig_units, user_units, var_name, freq, flux = stop("Transformation precipitation units not available.") } } else { - stop(paste("Transformation unknown from", orig_units[i], "to", user_units[i])) + stop(paste("Transformation unknown from", orig_units[i], "to", user_units[i])) } } } -- GitLab From c7c70ec0021d165436ca6390bc2b5fd6c39bb546 Mon Sep 17 00:00:00 2001 From: Victoria Agudetse Roures Date: Fri, 8 Sep 2023 16:32:53 +0200 Subject: [PATCH 19/21] Add Units() --- example_scripts/test_seasonal.R | 3 +++ 1 file changed, 3 insertions(+) diff --git a/example_scripts/test_seasonal.R b/example_scripts/test_seasonal.R index 0c1c1d2b..9734aeaf 100644 --- a/example_scripts/test_seasonal.R +++ b/example_scripts/test_seasonal.R @@ -6,6 +6,7 @@ # Load modules source("modules/Loading/Loading.R") +source("modules/Units/Units.R") source("modules/Calibration/Calibration.R") source("modules/Anomalies/Anomalies.R") source("modules/Skill/Skill.R") @@ -18,6 +19,8 @@ recipe <- prepare_outputs(recipe_file) # Load datasets data <- load_datasets(recipe) +# Change units +data <- Units(recipe, data) # Calibrate datasets data <- calibrate_datasets(recipe, data) # Compute anomalies -- GitLab From 1b6fd9e2f693b400a630c9c65feb57b426f2108b Mon Sep 17 00:00:00 2001 From: Victoria Agudetse Roures Date: Tue, 12 Sep 2023 12:04:35 +0200 Subject: [PATCH 20/21] Improve logic when reading Units from recipe --- modules/Units/Units.R | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/modules/Units/Units.R b/modules/Units/Units.R index 137ad56b..07e17543 100644 --- a/modules/Units/Units.R +++ b/modules/Units/Units.R @@ -34,13 +34,14 @@ Units <- function(recipe, data) { user_units <- orig_units[[which(!is.null(orig_units))[1]]] } else { user_units <- recipe$Analysis$Variables$units + if (!is.null(names(user_units))) { + # Change to vector and ensure that the units are in the correct order + user_units <- unlist(user_units[var_names[[1]]]) + } else if (length(var_names[[1]] > 1)) { + ## TODO: How to handle spaces in multi-var case? + user_units <- strsplit(recipe$Analysis$Variables$units, ", |,")[[1]] + } } - if (!is.null(names(user_units))) { - # Ensure that the units are in the correct order - user_units <- unlist(user_units[var_names[[1]]]) - } - ## TODO: How to handle spaces in multi-var case? - # user_units <- strsplit(recipe$Analysis$Variables$units, ", | |,")[[1]] # remove spaces, "**", "*" and "per" from units user_units <- tolower(user_units) user_units <- gsub(" ", "", user_units) -- GitLab From e4cbd88a4b8ef1632e88b3cdeba14b043ea4963b Mon Sep 17 00:00:00 2001 From: Victoria Agudetse Roures Date: Tue, 12 Sep 2023 12:37:22 +0200 Subject: [PATCH 21/21] Move testing recipes and scripts --- example_scripts/exec_ecvs_seasonal_oper.R | 35 +++++++++ example_scripts/exec_units.R | 42 +++++++++++ .../examples/recipe_ecvs_seasonal_oper.yml | 73 +++++++++++++++++++ .../examples/recipe_prlr_seasonal_oper.yml | 62 ++++++++++++++++ .../examples/recipe_prlr_seasonal_units.yml | 63 ++++++++++++++++ recipes/examples/recipe_tas_seasonal_oper.yml | 62 ++++++++++++++++ .../examples/recipe_tas_seasonal_units.yml | 63 ++++++++++++++++ 7 files changed, 400 insertions(+) create mode 100644 example_scripts/exec_ecvs_seasonal_oper.R create mode 100644 example_scripts/exec_units.R create mode 100644 recipes/examples/recipe_ecvs_seasonal_oper.yml create mode 100644 recipes/examples/recipe_prlr_seasonal_oper.yml create mode 100644 recipes/examples/recipe_prlr_seasonal_units.yml create mode 100644 recipes/examples/recipe_tas_seasonal_oper.yml create mode 100644 recipes/examples/recipe_tas_seasonal_units.yml diff --git a/example_scripts/exec_ecvs_seasonal_oper.R b/example_scripts/exec_ecvs_seasonal_oper.R new file mode 100644 index 00000000..18f2e493 --- /dev/null +++ b/example_scripts/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/example_scripts/exec_units.R b/example_scripts/exec_units.R new file mode 100644 index 00000000..bcd443c5 --- /dev/null +++ b/example_scripts/exec_units.R @@ -0,0 +1,42 @@ +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 <- "recipes/examples/recipe_tas_seasonal_units.yml" + # recipe_file <- "recipes/examples/recipe_prlr_seasonal_units.yml" + +recipe <- prepare_outputs(recipe_file) + +# Load datasets +data <- load_datasets(recipe) +# Units transformation +source("modules/Units/Units.R") +test <- Units(recipe, data) +# Calibrate datasets +data <- calibrate_datasets(recipe, test) +# 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/recipes/examples/recipe_ecvs_seasonal_oper.yml b/recipes/examples/recipe_ecvs_seasonal_oper.yml new file mode 100644 index 00000000..d47fd159 --- /dev/null +++ b/recipes/examples/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/recipes/examples/recipe_prlr_seasonal_oper.yml b/recipes/examples/recipe_prlr_seasonal_oper.yml new file mode 100644 index 00000000..fb6d1a1c --- /dev/null +++ b/recipes/examples/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/recipes/examples/recipe_prlr_seasonal_units.yml b/recipes/examples/recipe_prlr_seasonal_units.yml new file mode 100644 index 00000000..e03428ac --- /dev/null +++ b/recipes/examples/recipe_prlr_seasonal_units.yml @@ -0,0 +1,63 @@ +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 + units: mm + 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: # Optional, int: Forecast year 'YYYY' + hcst_start: '2014' # 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: 30 # Mandatory, int: minimum latitude + latmax: 50 # Mandatory, int: maximum latitude + lonmin: -10 # Mandatory, int: minimum longitude + lonmax: 10 # 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/examples/recipe_tas_seasonal_oper.yml b/recipes/examples/recipe_tas_seasonal_oper.yml new file mode 100644 index 00000000..75918265 --- /dev/null +++ b/recipes/examples/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/examples/recipe_tas_seasonal_units.yml b/recipes/examples/recipe_tas_seasonal_units.yml new file mode 100644 index 00000000..d2b25321 --- /dev/null +++ b/recipes/examples/recipe_tas_seasonal_units.yml @@ -0,0 +1,63 @@ +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 + units: K + 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: '2009' # 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: 30 # Mandatory, int: minimum latitude + latmax: 50 # Mandatory, int: maximum latitude + lonmin: -10 # Mandatory, int: minimum longitude + lonmax: 10 # 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/ -- GitLab