From ba62ffc23c2ea862bc4ef8e70c82b6b27b6982e3 Mon Sep 17 00:00:00 2001 From: Jaume Ramon Date: Wed, 26 Oct 2022 14:45:01 +0200 Subject: [PATCH 1/8] Downscaling recipie --- modules/Downscaling/Downscaling.R | 147 ++++++++++++++++++ .../testing_recipes/recipe_system5c3s-tas.yml | 60 ++++--- modules/test_seasonal.R | 3 +- modules/test_seasonal_downscaling.R | 27 ++++ 4 files changed, 212 insertions(+), 25 deletions(-) create mode 100644 modules/Downscaling/Downscaling.R create mode 100644 modules/test_seasonal_downscaling.R diff --git a/modules/Downscaling/Downscaling.R b/modules/Downscaling/Downscaling.R new file mode 100644 index 00000000..11c57a2a --- /dev/null +++ b/modules/Downscaling/Downscaling.R @@ -0,0 +1,147 @@ + +## TODO: Remove once CSDownscale is on CRAN +source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Interpolation.R') +source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Intbc.R') +source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Intlr.R') +source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Analogs.R') +source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Utils.R') +source("https://earth.bsc.es/gitlab/external/cstools/-/raw/master/R/CST_BiasCorrection.R") + + +## Entry params data and recipe? +downscale_datasets <- function(data, recipe) { + # Function that downscale the hindcast using the method stated in the + # recipe. For the moment, forecast must be null. + # + # data: list of s2dv_cube objects containing the hcst, obs and fcst. + # recipe: object obtained when passing the .yml recipe file to read_yaml() + + type <- tolower(recipe$Analysis$Workflow$Downscaling[["type"]]) + + if (!is.null(data$fcst)) { + warning("The downscaling will be only performed to the hindcast data") + data$fcst <- NULL + } + + if (type == "none") { + + hcst_calibrated <- data$hcst + DOWNSCAL_MSG <- "##### NO DOWNSCALING PERFORMED #####" + + } else { + # Downscaling function params + int_method <- tolower(recipe$Analysis$Workflow$Downscaling[["int_method"]]) + bc_method <- tolower(recipe$Analysis$Workflow$Downscaling[["bc_method"]]) + lr_method <- tolower(recipe$Analysis$Workflow$Downscaling[["lr_method"]]) + target_grid <- tolower(recipe$Analysis$Workflow$Downscaling[["target_grid"]]) + nanalogs <- tolower(recipe$Analysis$Workflow$Downscaling[["nanalogs"]]) + + if (is.null(recipe$Analysis$ncores)) { + ncores <- 1 + } else { + ncores <- recipe$Analysis$ncores + } + + if (is.null(recipe$Analysis$loocv)) { + loocv <- TRUE + } else { + loocv <- recipe$Analysis$loocv + } + + if (is.null(recipe$Analysis$remove_NAs)) { + na.rm = F + } else { + na.rm = recipe$Analysis$remove_NAs + } + + CALIB_MSG <- "##### CALIBRATION COMPLETE #####" + # Replicate observation array for the multi-model case + if (mm) { + obs.mm <- obs$data + for(dat in 1:(dim(data$hcst$data)['dat'][[1]]-1)) { + obs.mm <- abind(obs.mm, data$obs$data, + along=which(names(dim(data$obs$data)) == 'dat')) + } + names(dim(obs.mm)) <- names(dim(obs$data)) + data$obs$data <- obs.mm + remove(obs.mm) + } + + if (recipe$Analysis$Variables$freq == "monthly_mean") { + + CST_CALIB_METHODS <- c("bias", "evmos", "mse_min", "crps_min", "rpc-based") + ## TODO: implement other calibration methods + ## TODO: Restructure the code? + if (!(method %in% CST_CALIB_METHODS)) { + stop("Calibration method in the recipe is not available for monthly", + " data.") + } else { + ## Alba's version of CST_Calibration (pending merge) is being used + # Calibrate the hindcast + hcst_calibrated <- CST_Calibration(data$hcst, data$obs, + cal.method = method, + eval.method = "leave-one-out", + multi.model = mm, + na.fill = TRUE, + na.rm = na.rm, + apply_to = NULL, + alpha = NULL, + memb_dim = "ensemble", + sdate_dim = "syear", + ncores = ncores) + if (!is.null(data$fcst)) { + # Calibrate the forecast + fcst_calibrated <- CST_Calibration(data$hcst, data$obs, data$fcst, + cal.method = method, + eval.method = "leave-one-out", + multi.model = mm, + na.fill = TRUE, + na.rm = na.rm, + apply_to = NULL, + alpha = NULL, + memb_dim = "ensemble", + sdate_dim = "syear", + ncores = ncores) + } else { + fcst_calibrated <- NULL + } + } + } else if (recipe$Analysis$Variables$freq == "daily_mean") { + # Daily data calibration using Quantile Mapping + if (!(method %in% c("qmap"))) { + stop("Calibration method in the recipe is not available at daily ", + "frequency. Only quantile mapping 'qmap' is implemented.") + } + # Calibrate the hindcast + hcst_calibrated <- CST_QuantileMapping(data$hcst, data$obs, + exp_cor = NULL, + sample_dims = c("syear", + "time", + "ensemble"), + sample_length = NULL, + method = "QUANT", + wet.day = FALSE, + ncores = ncores, + na.rm = na.rm) + + if (!is.null(data$fcst)) { + # Calibrate the forecast + fcst_calibrated <- CST_QuantileMapping(data$hcst, data$obs, + exp_cor = data$fcst, + sample_dims = c("syear", + "time", + "ensemble"), + sample_length = NULL, + method = "QUANT", + wet.day = FALSE, + ncores = ncores, + na.rm = na.rm) + } else { + fcst_calibrated <- NULL + } + } + } +print(CALIB_MSG) + ## TODO: Return observations too? + return(list(hcst = hcst_calibrated, fcst = fcst_calibrated)) +} diff --git a/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml b/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml index 27cdccdc..15159553 100644 --- a/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml +++ b/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml @@ -1,11 +1,22 @@ +# ▄████▄ ██████ ▓█████▄ ▒█████ █ █░███▄ █ ██████ ▄████▄ ▄▄▄ ██▓ ▓█████ +#▒██▀ ▀█ ▒██ ▒ ▒██▀ ██▌▒██▒ ██▒▓█░ █ ░█░██ ▀█ █ ▒██ ▒ ▒██▀ ▀█ ▒████▄ ▓██▒ ▓█ ▀ +#▒▓█ ▄ ░ ▓██▄ ░██ █▌▒██░ ██▒▒█░ █ ░█▓██ ▀█ ██▒░ ▓██▄ ▒▓█ ▄ ▒██ ▀█▄ ▒██░ ▒███ +#▒▓▓▄ ▄██▒ ▒ ██▒░▓█▄ ▌▒██ ██░░█░ █ ░█▓██▒ ▐▌██▒ ▒ ██▒▒▓▓▄ ▄██▒░██▄▄▄▄██ ▒██░ ▒▓█ ▄ +#▒ ▓███▀ ░▒██████▒▒░▒████▓ ░ ████▓▒░░░██▒██▓▒██░ ▓██░▒██████▒▒▒ ▓███▀ ░ ▓█ ▓██▒░██████▒░▒████▒ +#░ ░▒ ▒ ░▒ ▒▓▒ ▒ ░ ▒▒▓ ▒ ░ ▒░▒░▒░ ░ ▓░▒ ▒ ░ ▒░ ▒ ▒ ▒ ▒▓▒ ▒ ░░ ░▒ ▒ ░ ▒▒ ▓▒█░░ ▒░▓ ░░░ ▒░ ░ +# ░ ▒ ░ ░▒ ░ ░ ░ ▒ ▒ ░ ▒ ▒░ ▒ ░ ░ ░ ░░ ░ ▒░░ ░▒ ░ ░ ░ ▒ ▒ ▒▒ ░░ ░ ▒ ░ ░ ░ ░ +#░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ▒ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ▒ ░ ░ ░ +#░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ +#░ ░ ░ Description: - Author: V. Agudetse + Author: J. Ramon +Info: downscaling of seasonal predictions from coarse to fine grids Analysis: Horizon: Seasonal Variables: name: tas - freq: monthly_mean + freq: daily_mean Datasets: System: name: system5c3s @@ -13,32 +24,33 @@ Analysis: Reference: name: era5 Time: - sdate: '0601' - fcst_year: '2020' - hcst_start: '1993' - hcst_end: '2006' + sdate: '0501' + hcst_start: '2000' + hcst_end: '2004' ftime_min: 1 ftime_max: 3 +#sometimes we want the region for the obs to be bigger than for hcst Region: - latmin: -10 - latmax: 10 - lonmin: 0 - lonmax: 20 - Regrid: - method: bilinear - type: to_system - Workflow: - Calibration: - method: raw - Skill: - metric: RPSS_specs BSS90_specs EnsCorr_specs FRPS_specs FRPSS_specs BSS10_specs FRPS - Probabilities: - percentiles: [[1/3, 2/3]] - Indicators: - index: no + latmin: 34.1 + latmax: 45.1 + lonmin: -12.5 + lonmax: 6.35 + Regrid: + type: none + Workflow: + Downscaling: + # Assumption 1: leave-one-out cross-validation is always applied + # TO DO: add downscaling to point locations + # TO DO: how can we pass the parameter 'predictors'? + type: # mandatory, 'none', 'int', 'intbc', 'intlr', 'analogs', 'logres' + int_method: # optional, regridding method accepted by CDO + bc_method: # optional, 'simple_bias', 'calibration', 'quantile_mapping' + lr_method: # optional, 'basic', 'large_scale', '4nn' + target_grid: # optional, nc file or grid accepted by CDO + nanalogs: # optional, number of analogs to be searched Output_format: S2S4E Run: Loglevel: INFO Terminal: yes - output_dir: /esarchive/scratch/vagudets/repos/auto-s2s/out-logs/ - code_dir: /esarchive/scratch/vagudets/repos/auto-s2s/ + output_dir: /esarchive/scratch/jramon/GitLab_jramon/auto-s2s/out-logs/ + code_dir: /esarchive/scratch/jramon/GitLab_jramon/auto-s2s/ diff --git a/modules/test_seasonal.R b/modules/test_seasonal.R index 5f59794f..b9f94b13 100644 --- a/modules/test_seasonal.R +++ b/modules/test_seasonal.R @@ -1,11 +1,12 @@ +setwd("/esarchive/scratch/jramon/GitLab_jramon/auto-s2s/") source("modules/Loading/Loading.R") source("modules/Calibration/Calibration.R") source("modules/Skill/Skill.R") source("modules/Saving/Saving.R") source("modules/Visualization/Visualization.R") -recipe_file <- "modules/Loading/testing_recipes/recipe_system7c3s-tas.yml" +recipe_file <- "modules/Loading/testing_recipes/recipe_system5c3s-tas.yml" recipe <- read_yaml(recipe_file) archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive.yml"))$archive diff --git a/modules/test_seasonal_downscaling.R b/modules/test_seasonal_downscaling.R new file mode 100644 index 00000000..ec1d2f2e --- /dev/null +++ b/modules/test_seasonal_downscaling.R @@ -0,0 +1,27 @@ + +# DO source ../MODULES before opening the R session + +setwd("/esarchive/scratch/jramon/GitLab_jramon/auto-s2s/") +source("modules/Loading/Loading.R") +source("modules/Calibration/Calibration.R") +source("modules/Skill/Skill.R") +source("modules/Saving/Saving.R") +source("modules/Visualization/Visualization.R") + +recipe_file <- "modules/Loading/testing_recipes/recipe_system5c3s-tas.yml" +recipe <- read_yaml(recipe_file) +archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive.yml"))$archive + +# Load datasets +data <- load_datasets(recipe_file) +# Calibrate datasets +calibrated_data <- calibrate_datasets(data, recipe) +# Compute skill metrics +skill_metrics <- compute_skill_metrics(calibrated_data$hcst, data$obs, recipe) +# Compute percentiles and probability bins +probabilities <- compute_probabilities(calibrated_data$hcst, recipe) +# Export all data to netCDF +save_data(recipe, archive, data, calibrated_data, skill_metrics, probabilities) +# Plot data +plot_data(recipe, archive, data, calibrated_data, skill_metrics, + probabilities, significance = T) -- GitLab From 0052980a4b978b76823d93afe0b716419dbed043 Mon Sep 17 00:00:00 2001 From: Jaume Ramon Date: Fri, 28 Oct 2022 15:09:28 +0200 Subject: [PATCH 2/8] Downscaling module v0 --- modules/Downscaling/Downscaling.R | 299 ++++++++++++++++++++---------- 1 file changed, 198 insertions(+), 101 deletions(-) diff --git a/modules/Downscaling/Downscaling.R b/modules/Downscaling/Downscaling.R index 11c57a2a..8ff8dee9 100644 --- a/modules/Downscaling/Downscaling.R +++ b/modules/Downscaling/Downscaling.R @@ -1,4 +1,13 @@ - +# ▄████▄ ██████ ▓█████▄ ▒█████ █ █░███▄ █ ██████ ▄████▄ ▄▄▄ ██▓ ▓█████ +#▒██▀ ▀█ ▒██ ▒ ▒██▀ ██▌▒██▒ ██▒▓█░ █ ░█░██ ▀█ █ ▒██ ▒ ▒██▀ ▀█ ▒████▄ ▓██▒ ▓█ ▀ +#▒▓█ ▄ ░ ▓██▄ ░██ █▌▒██░ ██▒▒█░ █ ░█▓██ ▀█ ██▒░ ▓██▄ ▒▓█ ▄ ▒██ ▀█▄ ▒██░ ▒███ +#▒▓▓▄ ▄██▒ ▒ ██▒░▓█▄ ▌▒██ ██░░█░ █ ░█▓██▒ ▐▌██▒ ▒ ██▒▒▓▓▄ ▄██▒░██▄▄▄▄██ ▒██░ ▒▓█ ▄ +#▒ ▓███▀ ░▒██████▒▒░▒████▓ ░ ████▓▒░░░██▒██▓▒██░ ▓██░▒██████▒▒▒ ▓███▀ ░ ▓█ ▓██▒░██████▒░▒████▒ +#░ ░▒ ▒ ░▒ ▒▓▒ ▒ ░ ▒▒▓ ▒ ░ ▒░▒░▒░ ░ ▓░▒ ▒ ░ ▒░ ▒ ▒ ▒ ▒▓▒ ▒ ░░ ░▒ ▒ ░ ▒▒ ▓▒█░░ ▒░▓ ░░░ ▒░ ░ +# ░ ▒ ░ ░▒ ░ ░ ░ ▒ ▒ ░ ▒ ▒░ ▒ ░ ░ ░ ░░ ░ ▒░░ ░▒ ░ ░ ░ ▒ ▒ ▒▒ ░░ ░ ▒ ░ ░ ░ ░ +#░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ▒ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ▒ ░ ░ ░ +#░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ +#░ ░ ░ ## TODO: Remove once CSDownscale is on CRAN source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Interpolation.R') source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Intbc.R') @@ -16,7 +25,7 @@ downscale_datasets <- function(data, recipe) { # data: list of s2dv_cube objects containing the hcst, obs and fcst. # recipe: object obtained when passing the .yml recipe file to read_yaml() - type <- tolower(recipe$Analysis$Workflow$Downscaling[["type"]]) + type <- tolower(recipe$Analysis$Workflow$Downscaling$type) if (!is.null(data$fcst)) { warning("The downscaling will be only performed to the hindcast data") @@ -30,118 +39,206 @@ downscale_datasets <- function(data, recipe) { } else { # Downscaling function params - int_method <- tolower(recipe$Analysis$Workflow$Downscaling[["int_method"]]) - bc_method <- tolower(recipe$Analysis$Workflow$Downscaling[["bc_method"]]) - lr_method <- tolower(recipe$Analysis$Workflow$Downscaling[["lr_method"]]) - target_grid <- tolower(recipe$Analysis$Workflow$Downscaling[["target_grid"]]) - nanalogs <- tolower(recipe$Analysis$Workflow$Downscaling[["nanalogs"]]) - + int_method <- tolower(recipe$Analysis$Workflow$Downscaling$int_method) + bc_method <- tolower(recipe$Analysis$Workflow$Downscaling$bc_method) + lr_method <- tolower(recipe$Analysis$Workflow$Downscaling$lr_method) + log_reg_method <- tolower(recipe$Analysis$Workflow$Downscaling$log_reg_method) + target_grid <- tolower(recipe$Analysis$Workflow$Downscaling$target_grid) + nanalogs <- tolower(recipe$Analysis$Workflow$Downscaling$nanalog) + if (is.null(recipe$Analysis$ncores)) { ncores <- 1 } else { ncores <- recipe$Analysis$ncores } - + + #TO DO: add the parametre loocv where it corresponds if (is.null(recipe$Analysis$loocv)) { loocv <- TRUE } else { loocv <- recipe$Analysis$loocv } - if (is.null(recipe$Analysis$remove_NAs)) { - na.rm = F - } else { - na.rm = recipe$Analysis$remove_NAs - } - - CALIB_MSG <- "##### CALIBRATION COMPLETE #####" - # Replicate observation array for the multi-model case - if (mm) { - obs.mm <- obs$data - for(dat in 1:(dim(data$hcst$data)['dat'][[1]]-1)) { - obs.mm <- abind(obs.mm, data$obs$data, - along=which(names(dim(data$obs$data)) == 'dat')) - } - names(dim(obs.mm)) <- names(dim(obs$data)) - data$obs$data <- obs.mm - remove(obs.mm) - } - - if (recipe$Analysis$Variables$freq == "monthly_mean") { + DOWNSCAL_TYPES <- c("none", "int", "intbc", "intlr", "analogs", "logres") + BC_METHODS <- c("simple_bias", "calibration", "quantile_mapping") + LR_METHODS <- c("basic", "large_scale", "4nn") + LOG_REG_METHODS <- c("ens_mean", "ens_mean_sd", "sorted_members") + + if (!(type %in% DOWNSCAL_TYPES)) { + stop("Downscaling type in the recipe is not available. Accepted types ", + "are 'none', 'int', 'intbc', 'intlr', 'analogs', 'logres'.") + } + + if (type == "int") { + if (is.null(int_method)) { + stop("Please provide one interpolation method in the recipe.") + } + + if (is.null(target_grid)) { + stop("Please provide the target grid in the recipe.") + } + + hcst_downscal <- CST_Interpolation(data$hcst$data, + points = NULL, + method_remap = int_method, + target_grid = target_grid, + lat_dim = "latitude", + lon_dim = "longitude", + region = NULL, + method_point_interp = NULL) + + + DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" + } else if (type == "intbc") { + if (is.null(int_method)) { + stop("Please provide one interpolation method in the recipe.") + } + + + if (is.null(bc_method)) { + stop("Please provide one bias-correction method in the recipe. Accepted ", + "methods are 'simple_bias', 'calibration', 'quantile_mapping'.") + } + + if (is.null(target_grid)) { + stop("Please provide the target grid in the recipe.") + } + + if (!(bc_method %in% BC_METHODS)) { + stop("Bias-correction method in the recipe is not available. Accepted methods ", + "are 'simple_bias', 'calibration', 'quantile_mapping'.") + } + + hcst_downscal <- CST_Intbc(data$hcst$data, data$obs$data, + target_grid = target_grid, + bc_method = bc_method, + int_method = int_method, + points = NULL, + method_point_interp = NULL, + lat_dim = "latitude", + lon_dim = "longitude", + sdate_dim = "syear", + member_dim = "ensemble", + region = NULL, + ncores = ncores) - CST_CALIB_METHODS <- c("bias", "evmos", "mse_min", "crps_min", "rpc-based") - ## TODO: implement other calibration methods - ## TODO: Restructure the code? - if (!(method %in% CST_CALIB_METHODS)) { - stop("Calibration method in the recipe is not available for monthly", - " data.") - } else { - ## Alba's version of CST_Calibration (pending merge) is being used - # Calibrate the hindcast - hcst_calibrated <- CST_Calibration(data$hcst, data$obs, - cal.method = method, - eval.method = "leave-one-out", - multi.model = mm, - na.fill = TRUE, - na.rm = na.rm, - apply_to = NULL, - alpha = NULL, - memb_dim = "ensemble", - sdate_dim = "syear", - ncores = ncores) - if (!is.null(data$fcst)) { - # Calibrate the forecast - fcst_calibrated <- CST_Calibration(data$hcst, data$obs, data$fcst, - cal.method = method, - eval.method = "leave-one-out", - multi.model = mm, - na.fill = TRUE, - na.rm = na.rm, - apply_to = NULL, - alpha = NULL, - memb_dim = "ensemble", - sdate_dim = "syear", - ncores = ncores) - } else { - fcst_calibrated <- NULL + DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" + } else if (type == "intlr") { + + if (is.null(int_method)) { + stop("Please provide one interpolation method in the recipe.") + } + + if (is.null(lr_method)) { + stop("Please provide one linear regression method in the recipe. Accepted ", + "methods are 'basic', 'large_scale', '4nn'.") + } + + if (is.null(target_grid)) { + stop("Please provide the target grid in the recipe.") + } + + if (!(lr_method %in% LR_METHODS)) { + stop("Linear regression method in the recipe is not available. Accepted methods ", + "are 'basic', 'large_scale', '4nn'.") + } + + # TO DO: add the possibility to have the element 'pred' in 'data' + if (lr_method == "large_scale") { + if (is.null(data$pred$data)) { + stop("Please provide the large scale predictors in the element 'data$pred$data'.") } - } - } else if (recipe$Analysis$Variables$freq == "daily_mean") { - # Daily data calibration using Quantile Mapping - if (!(method %in% c("qmap"))) { - stop("Calibration method in the recipe is not available at daily ", - "frequency. Only quantile mapping 'qmap' is implemented.") - } - # Calibrate the hindcast - hcst_calibrated <- CST_QuantileMapping(data$hcst, data$obs, - exp_cor = NULL, - sample_dims = c("syear", - "time", - "ensemble"), - sample_length = NULL, - method = "QUANT", - wet.day = FALSE, - ncores = ncores, - na.rm = na.rm) - - if (!is.null(data$fcst)) { - # Calibrate the forecast - fcst_calibrated <- CST_QuantileMapping(data$hcst, data$obs, - exp_cor = data$fcst, - sample_dims = c("syear", - "time", - "ensemble"), - sample_length = NULL, - method = "QUANT", - wet.day = FALSE, - ncores = ncores, - na.rm = na.rm) } else { - fcst_calibrated <- NULL + data$pred$data <- NULL + } + + if (is.null(target_grid)) { + stop("Please provide the target grid in the recipe.") + } + if (is.null(target_grid)) { + stop("Please provide the target grid in the recipe.") } + + + hcst_downscal <- CST_Intlr(data$hcst$data, data$obs$data, + lr_method = lr_method, + target_grid = target_grid, + points = NULL, + int_method = int_method, + method_point_interp = NULL, + predictors = data$pred$data, + lat_dim = "latitude", + lon_dim = "longitude", + sdate_dim = "syear", + time_dim = "time", + large_scale_predictor_dimname = 'vars', + loocv = loocv, + region = NULL, + ncores = ncores) + + DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" + } else if (type == "analogs") { + + if (is.null(nanalogs)) { + warning("The number of analogs for searching has not been provided in the ", + "recipe. Setting it to 3.") + nanalogs <- 3 + } + + hcst_downscal <- CST_Analogs(data$hcst$data, data$obs$data, + grid_exp = data$hcst$source_files[1], + nanalogs = nanalogs, + fun_analog = NULL, + lat_dim = "latitude", + lon_dim = "longitude", + sdate_dim = "syear", + time_dim = "time", + region = NULL, + return_indices = FALSE, + loocv_window = loocv, + ncores = ncores) + + DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" + } else if (type == "logres") { + + if (is.null(int_method)) { + stop("Please provide one interpolation method in the recipe.") + } + + if (is.null(log_reg_method)) { + stop("Please provide one logistic regression method in the recipe. Accepted ", + "methods are 'ens_mean', 'ens_mean_sd', 'sorted_members'.") + } + + if (is.null(target_grid)) { + stop("Please provide the target grid in the recipe.") + } + + if (!(log_reg_method %in% LOG_REG_METHODS)) { + stop("Logistic regression method in the recipe is not available. Accepted methods ", + "are 'ens_mean', 'ens_mean_sd', 'sorted_members'.") + } + + hcst_downscal <- CST_LogisticReg(data$hcst$data, data$obs$data, + target_grid = target_grid, + int_method = int_method, + log_reg_method = log_reg_method, + probs_cat = c(1/3,2/3), + return_most_likely_cat = FALSE, + points = NULL, + method_point_interp = NULL, + lat_dim = "latitude", + lon_dim = "longitude", + sdate_dim = "syear", + member_dim = "ensemble", + region = NULL, + loocv = loocv, + ncores = ncores) + + DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" } - } -print(CALIB_MSG) - ## TODO: Return observations too? - return(list(hcst = hcst_calibrated, fcst = fcst_calibrated)) -} + + print(DOWNSCAL_MSG) + return(list(hcst = hcst_downscal, fcst = NULL)) +} + -- GitLab From 4a59cf73766e0d92544c069cd89fe6914330ecb2 Mon Sep 17 00:00:00 2001 From: Jaume Ramon Date: Wed, 2 Nov 2022 11:01:29 +0100 Subject: [PATCH 3/8] . --- .../testing_recipes/recipe_system5c3s-tas.yml | 61 ++++++++----------- .../recipe_system5c3s-tas_downscaling.yml | 57 +++++++++++++++++ modules/test_seasonal.R | 2 +- 3 files changed, 83 insertions(+), 37 deletions(-) create mode 100644 modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml diff --git a/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml b/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml index 15159553..38e23c8b 100644 --- a/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml +++ b/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml @@ -1,22 +1,11 @@ -# ▄████▄ ██████ ▓█████▄ ▒█████ █ █░███▄ █ ██████ ▄████▄ ▄▄▄ ██▓ ▓█████ -#▒██▀ ▀█ ▒██ ▒ ▒██▀ ██▌▒██▒ ██▒▓█░ █ ░█░██ ▀█ █ ▒██ ▒ ▒██▀ ▀█ ▒████▄ ▓██▒ ▓█ ▀ -#▒▓█ ▄ ░ ▓██▄ ░██ █▌▒██░ ██▒▒█░ █ ░█▓██ ▀█ ██▒░ ▓██▄ ▒▓█ ▄ ▒██ ▀█▄ ▒██░ ▒███ -#▒▓▓▄ ▄██▒ ▒ ██▒░▓█▄ ▌▒██ ██░░█░ █ ░█▓██▒ ▐▌██▒ ▒ ██▒▒▓▓▄ ▄██▒░██▄▄▄▄██ ▒██░ ▒▓█ ▄ -#▒ ▓███▀ ░▒██████▒▒░▒████▓ ░ ████▓▒░░░██▒██▓▒██░ ▓██░▒██████▒▒▒ ▓███▀ ░ ▓█ ▓██▒░██████▒░▒████▒ -#░ ░▒ ▒ ░▒ ▒▓▒ ▒ ░ ▒▒▓ ▒ ░ ▒░▒░▒░ ░ ▓░▒ ▒ ░ ▒░ ▒ ▒ ▒ ▒▓▒ ▒ ░░ ░▒ ▒ ░ ▒▒ ▓▒█░░ ▒░▓ ░░░ ▒░ ░ -# ░ ▒ ░ ░▒ ░ ░ ░ ▒ ▒ ░ ▒ ▒░ ▒ ░ ░ ░ ░░ ░ ▒░░ ░▒ ░ ░ ░ ▒ ▒ ▒▒ ░░ ░ ▒ ░ ░ ░ ░ -#░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ▒ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ▒ ░ ░ ░ -#░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ -#░ ░ ░ Description: - Author: J. Ramon -Info: downscaling of seasonal predictions from coarse to fine grids + Author: V. Agudetse Analysis: Horizon: Seasonal Variables: name: tas - freq: daily_mean + freq: monthly_mean Datasets: System: name: system5c3s @@ -24,33 +13,33 @@ Analysis: Reference: name: era5 Time: - sdate: '0501' - hcst_start: '2000' - hcst_end: '2004' + sdate: '0601' + fcst_year: '2020' + hcst_start: '1993' + hcst_end: '2006' ftime_min: 1 ftime_max: 3 -#sometimes we want the region for the obs to be bigger than for hcst Region: - latmin: 34.1 - latmax: 45.1 - lonmin: -12.5 - lonmax: 6.35 - Regrid: - type: none - Workflow: - Downscaling: - # Assumption 1: leave-one-out cross-validation is always applied - # TO DO: add downscaling to point locations - # TO DO: how can we pass the parameter 'predictors'? - type: # mandatory, 'none', 'int', 'intbc', 'intlr', 'analogs', 'logres' - int_method: # optional, regridding method accepted by CDO - bc_method: # optional, 'simple_bias', 'calibration', 'quantile_mapping' - lr_method: # optional, 'basic', 'large_scale', '4nn' - target_grid: # optional, nc file or grid accepted by CDO - nanalogs: # optional, number of analogs to be searched + latmin: -10 + latmax: 10 + lonmin: 0 + lonmax: 20 + Regrid: + method: bilinear + type: to_system + Workflow: + Calibration: + method: raw + Skill: + metric: RPSS_specs BSS90_specs EnsCorr_specs FRPS_specs FRPSS_specs BSS10_specs FRPS + Probabilities: + percentiles: [[1/3, 2/3]] + Indicators: + index: no Output_format: S2S4E Run: Loglevel: INFO Terminal: yes - output_dir: /esarchive/scratch/jramon/GitLab_jramon/auto-s2s/out-logs/ - code_dir: /esarchive/scratch/jramon/GitLab_jramon/auto-s2s/ + output_dir: /esarchive/scratch/vagudets/repos/auto-s2s/out-logs/ + code_dir: /esarchive/scratch/vagudets/repos/auto-s2s/ + diff --git a/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml b/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml new file mode 100644 index 00000000..13f4e156 --- /dev/null +++ b/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml @@ -0,0 +1,57 @@ +# ▄████▄ ██████ ▓█████▄ ▒█████ █ █░███▄ █ ██████ ▄████▄ ▄▄▄ ██▓ ▓█████ +#▒██▀ ▀█ ▒██ ▒ ▒██▀ ██▌▒██▒ ██▒▓█░ █ ░█░██ ▀█ █ ▒██ ▒ ▒██▀ ▀█ ▒████▄ ▓██▒ ▓█ ▀ +#▒▓█ ▄ ░ ▓██▄ ░██ █▌▒██░ ██▒▒█░ █ ░█▓██ ▀█ ██▒░ ▓██▄ ▒▓█ ▄ ▒██ ▀█▄ ▒██░ ▒███ +#▒▓▓▄ ▄██▒ ▒ ██▒░▓█▄ ▌▒██ ██░░█░ █ ░█▓██▒ ▐▌██▒ ▒ ██▒▒▓▓▄ ▄██▒░██▄▄▄▄██ ▒██░ ▒▓█ ▄ +#▒ ▓███▀ ░▒██████▒▒░▒████▓ ░ ████▓▒░░░██▒██▓▒██░ ▓██░▒██████▒▒▒ ▓███▀ ░ ▓█ ▓██▒░██████▒░▒████▒ +#░ ░▒ ▒ ░▒ ▒▓▒ ▒ ░ ▒▒▓ ▒ ░ ▒░▒░▒░ ░ ▓░▒ ▒ ░ ▒░ ▒ ▒ ▒ ▒▓▒ ▒ ░░ ░▒ ▒ ░ ▒▒ ▓▒█░░ ▒░▓ ░░░ ▒░ ░ +# ░ ▒ ░ ░▒ ░ ░ ░ ▒ ▒ ░ ▒ ▒░ ▒ ░ ░ ░ ░░ ░ ▒░░ ░▒ ░ ░ ░ ▒ ▒ ▒▒ ░░ ░ ▒ ░ ░ ░ ░ +#░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ▒ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ▒ ░ ░ ░ +#░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ ░ +#░ ░ ░ +Description: + Author: J. Ramon +Info: downscaling of seasonal predictions from coarse to fine grids + +Analysis: + Horizon: Seasonal + Variables: + name: tas + freq: daily_mean + Datasets: + System: + name: system5c3s + Multimodel: no + Reference: + name: era5 + Time: + sdate: '0501' + hcst_start: '2000' + hcst_end: '2004' + ftime_min: 1 + ftime_max: 3 +#sometimes we want the region for the obs to be bigger than for hcst + Region: + latmin: 34.1 + latmax: 45.1 + lonmin: -12.5 + lonmax: 6.35 + Regrid: + type: none + Workflow: + Downscaling: + # Assumption 1: leave-one-out cross-validation is always applied + # TO DO: add downscaling to point locations + # TO DO: how can we pass the parameter 'predictors'? + type: # mandatory, 'none', 'int', 'intbc', 'intlr', 'analogs', 'logreg' + int_method: # optional, regridding method accepted by CDO + bc_method: # optional, 'simple_bias', 'calibration', 'quantile_mapping' + lr_method: # optional, 'basic', 'large_scale', '4nn' + log_reg_method: # optional, 'ens_mean', 'ens_mean_sd', 'sorted_members' + target_grid: # optional, nc file or grid accepted by CDO + nanalogs: # optional, number of analogs to be searched + Output_format: S2S4E +Run: + Loglevel: INFO + Terminal: yes + output_dir: /esarchive/scratch/jramon/GitLab_jramon/auto-s2s/out-logs/ + code_dir: /esarchive/scratch/jramon/GitLab_jramon/auto-s2s/ diff --git a/modules/test_seasonal.R b/modules/test_seasonal.R index b9f94b13..58adece6 100644 --- a/modules/test_seasonal.R +++ b/modules/test_seasonal.R @@ -6,7 +6,7 @@ source("modules/Skill/Skill.R") source("modules/Saving/Saving.R") source("modules/Visualization/Visualization.R") -recipe_file <- "modules/Loading/testing_recipes/recipe_system5c3s-tas.yml" +recipe_file <- "modules/Loading/testing_recipes/recipe_system7c3s-tas.yml" recipe <- read_yaml(recipe_file) archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive.yml"))$archive -- GitLab From fac2975020bc4efc5504a8fd0ed70bb1f5f56cc3 Mon Sep 17 00:00:00 2001 From: Jaume Ramon Date: Wed, 2 Nov 2022 17:30:21 +0100 Subject: [PATCH 4/8] Downscaling recipe with the downscaling module --- .../testing_recipes/recipe_system5c3s-tas_downscaling.yml | 8 +++++++- modules/test_seasonal_downscaling.R | 7 +++++-- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml b/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml index 13f4e156..8b56f3aa 100644 --- a/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml +++ b/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml @@ -16,7 +16,7 @@ Analysis: Horizon: Seasonal Variables: name: tas - freq: daily_mean + freq: monthly_mean Datasets: System: name: system5c3s @@ -25,6 +25,7 @@ Analysis: name: era5 Time: sdate: '0501' + fcst_year: '2020' hcst_start: '2000' hcst_end: '2004' ftime_min: 1 @@ -36,8 +37,13 @@ Analysis: lonmin: -12.5 lonmax: 6.35 Regrid: + method: bilinear type: none Workflow: + Calibration: + method: raw + Skill: + metric: RPSS Downscaling: # Assumption 1: leave-one-out cross-validation is always applied # TO DO: add downscaling to point locations diff --git a/modules/test_seasonal_downscaling.R b/modules/test_seasonal_downscaling.R index ec1d2f2e..9043c2a9 100644 --- a/modules/test_seasonal_downscaling.R +++ b/modules/test_seasonal_downscaling.R @@ -3,17 +3,20 @@ setwd("/esarchive/scratch/jramon/GitLab_jramon/auto-s2s/") source("modules/Loading/Loading.R") +source("modules/Loading/Downscaling.R") source("modules/Calibration/Calibration.R") source("modules/Skill/Skill.R") source("modules/Saving/Saving.R") source("modules/Visualization/Visualization.R") -recipe_file <- "modules/Loading/testing_recipes/recipe_system5c3s-tas.yml" +recipe_file <- "modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml" recipe <- read_yaml(recipe_file) archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive.yml"))$archive # Load datasets -data <- load_datasets(recipe_file) +data <- load_datasets(recipe) +# Downscale datasets +downscaled_data <- downscale_datasets(data, recipe) # Calibrate datasets calibrated_data <- calibrate_datasets(data, recipe) # Compute skill metrics -- GitLab From 521489d25a4a1c4b72b73fd112b254b08b2c2b65 Mon Sep 17 00:00:00 2001 From: Jaume Ramon Date: Thu, 10 Nov 2022 12:53:13 +0100 Subject: [PATCH 5/8] module working for Interpolation and Intbc --- modules/Downscaling/Downscaling.R | 223 ++++++++++-------- .../recipe_system5c3s-tas_downscaling.yml | 14 +- modules/test_seasonal_downscaling.R | 24 +- 3 files changed, 134 insertions(+), 127 deletions(-) diff --git a/modules/Downscaling/Downscaling.R b/modules/Downscaling/Downscaling.R index 8ff8dee9..4d3edd2e 100644 --- a/modules/Downscaling/Downscaling.R +++ b/modules/Downscaling/Downscaling.R @@ -18,7 +18,7 @@ source("https://earth.bsc.es/gitlab/external/cstools/-/raw/master/R/CST_BiasCorr ## Entry params data and recipe? -downscale_datasets <- function(data, recipe) { +downscale_datasets <- function(recipe, data) { # Function that downscale the hindcast using the method stated in the # recipe. For the moment, forecast must be null. # @@ -34,17 +34,17 @@ downscale_datasets <- function(data, recipe) { if (type == "none") { - hcst_calibrated <- data$hcst + hcst_downscal <- data$hcst DOWNSCAL_MSG <- "##### NO DOWNSCALING PERFORMED #####" } else { # Downscaling function params - int_method <- tolower(recipe$Analysis$Workflow$Downscaling$int_method) - bc_method <- tolower(recipe$Analysis$Workflow$Downscaling$bc_method) - lr_method <- tolower(recipe$Analysis$Workflow$Downscaling$lr_method) - log_reg_method <- tolower(recipe$Analysis$Workflow$Downscaling$log_reg_method) - target_grid <- tolower(recipe$Analysis$Workflow$Downscaling$target_grid) - nanalogs <- tolower(recipe$Analysis$Workflow$Downscaling$nanalog) + int_methods <- tolower(recipe$Analysis$Workflow$Downscaling$int_method) + bc_methods <- tolower(recipe$Analysis$Workflow$Downscaling$bc_method) + lr_methods <- tolower(recipe$Analysis$Workflow$Downscaling$lr_method) + log_reg_methods <- tolower(recipe$Analysis$Workflow$Downscaling$log_reg_method) + target_grid <- tolower(recipe$Analysis$Workflow$Downscaling$target_grid) + nanalogs <- as.numeric(recipe$Analysis$Workflow$Downscaling$nanalogs) if (is.null(recipe$Analysis$ncores)) { ncores <- 1 @@ -70,32 +70,36 @@ downscale_datasets <- function(data, recipe) { } if (type == "int") { - if (is.null(int_method)) { - stop("Please provide one interpolation method in the recipe.") + if (is.null(int_methods)) { + stop("Please provide one or more interpolation methods in the recipe.") } if (is.null(target_grid)) { stop("Please provide the target grid in the recipe.") } - hcst_downscal <- CST_Interpolation(data$hcst$data, - points = NULL, - method_remap = int_method, - target_grid = target_grid, - lat_dim = "latitude", - lon_dim = "longitude", - region = NULL, - method_point_interp = NULL) - + hcst_downscal <- list() + for (int_method in strsplit(int_methods, ", | |,")[[1]]) { + + hcst_d <- CST_Interpolation(data$hcst, + points = NULL, + method_remap = int_method, + target_grid = target_grid, + lat_dim = "latitude", + lon_dim = "longitude", + region = NULL, + method_point_interp = NULL) + + hcst_downscal[[ int_method ]] <- hcst_d + } DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" } else if (type == "intbc") { - if (is.null(int_method)) { - stop("Please provide one interpolation method in the recipe.") + if (is.null(int_methods) | (length(strsplit(int_methods, ", | |,")[[1]]) != 1)) { + stop("Please provide one (and only one) interpolation method in the recipe.") } - - if (is.null(bc_method)) { + if (is.null(bc_methods)) { stop("Please provide one bias-correction method in the recipe. Accepted ", "methods are 'simple_bias', 'calibration', 'quantile_mapping'.") } @@ -104,32 +108,38 @@ downscale_datasets <- function(data, recipe) { stop("Please provide the target grid in the recipe.") } - if (!(bc_method %in% BC_METHODS)) { - stop("Bias-correction method in the recipe is not available. Accepted methods ", - "are 'simple_bias', 'calibration', 'quantile_mapping'.") - } - - hcst_downscal <- CST_Intbc(data$hcst$data, data$obs$data, - target_grid = target_grid, - bc_method = bc_method, - int_method = int_method, - points = NULL, - method_point_interp = NULL, - lat_dim = "latitude", - lon_dim = "longitude", - sdate_dim = "syear", - member_dim = "ensemble", - region = NULL, - ncores = ncores) - + hcst_downscal <- list() + for (bc_method in strsplit(bc_methods, ", | |,")[[1]]) { + + if (!(bc_method %in% BC_METHODS)) { + stop(paste0(bc_method, " method in the recipe is not available. Accepted methods ", + "are 'simple_bias', 'calibration', 'quantile_mapping'.")) + } + + hcst_d <- CST_Intbc(data$hcst, data$obs, + target_grid = target_grid, + bc_method = bc_method, + int_method = int_methods, + points = NULL, + method_point_interp = NULL, + lat_dim = "latitude", + lon_dim = "longitude", + sdate_dim = "syear", + member_dim = "ensemble", + region = NULL, + ncores = ncores) + + hcst_downscal[[ bc_method ]] <- hcst_d + } + DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" } else if (type == "intlr") { - if (is.null(int_method)) { - stop("Please provide one interpolation method in the recipe.") + if (is.null(int_methods) | (length(strsplit(int_methods, ", | |,")[[1]]) != 1)) { + stop("Please provide one (and only one) interpolation method in the recipe.") } - if (is.null(lr_method)) { + if (is.null(lr_methods)) { stop("Please provide one linear regression method in the recipe. Accepted ", "methods are 'basic', 'large_scale', '4nn'.") } @@ -138,44 +148,42 @@ downscale_datasets <- function(data, recipe) { stop("Please provide the target grid in the recipe.") } - if (!(lr_method %in% LR_METHODS)) { - stop("Linear regression method in the recipe is not available. Accepted methods ", - "are 'basic', 'large_scale', '4nn'.") - } + hcst_downscal <- list() + for (lr_method in strsplit(lr_methods, ", | |,")[[1]]) { - # TO DO: add the possibility to have the element 'pred' in 'data' - if (lr_method == "large_scale") { - if (is.null(data$pred$data)) { - stop("Please provide the large scale predictors in the element 'data$pred$data'.") + if (!(lr_method %in% LR_METHODS)) { + stop(paste0(lr_method, " method in the recipe is not available. Accepted methods ", + "are 'basic', 'large_scale', '4nn'.")) } - } else { - data$pred$data <- NULL - } - if (is.null(target_grid)) { - stop("Please provide the target grid in the recipe.") - } - if (is.null(target_grid)) { - stop("Please provide the target grid in the recipe.") - } + # TO DO: add the possibility to have the element 'pred' in 'data' + if (lr_method == "large_scale") { + if (is.null(data$pred$data)) { + stop("Please provide the large scale predictors in the element 'data$pred$data'.") + } + } else { + data$pred$data <- NULL + } + hcst_d <- CST_Intlr(data$hcst, data$obs, + lr_method = lr_method, + target_grid = target_grid, + points = NULL, + int_method = int_methods, + method_point_interp = NULL, + predictors = data$pred$data, + lat_dim = "latitude", + lon_dim = "longitude", + sdate_dim = "syear", + time_dim = "time", + large_scale_predictor_dimname = 'vars', + loocv = loocv, + region = NULL, + ncores = ncores) - hcst_downscal <- CST_Intlr(data$hcst$data, data$obs$data, - lr_method = lr_method, - target_grid = target_grid, - points = NULL, - int_method = int_method, - method_point_interp = NULL, - predictors = data$pred$data, - lat_dim = "latitude", - lon_dim = "longitude", - sdate_dim = "syear", - time_dim = "time", - large_scale_predictor_dimname = 'vars', - loocv = loocv, - region = NULL, - ncores = ncores) - + hcst_downscal[[ lr_method ]] <- hcst_d + } + DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" } else if (type == "analogs") { @@ -184,8 +192,9 @@ downscale_datasets <- function(data, recipe) { "recipe. Setting it to 3.") nanalogs <- 3 } - - hcst_downscal <- CST_Analogs(data$hcst$data, data$obs$data, + + + hcst_downscal <- CST_Analogs(data$hcst, data$obs, grid_exp = data$hcst$source_files[1], nanalogs = nanalogs, fun_analog = NULL, @@ -201,11 +210,11 @@ downscale_datasets <- function(data, recipe) { DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" } else if (type == "logres") { - if (is.null(int_method)) { - stop("Please provide one interpolation method in the recipe.") + if (is.null(int_methods) | (length(strsplit(int_methods, ", | |,")[[1]]) != 1)) { + stop("Please provide one (and only one) interpolation method in the recipe.") } - if (is.null(log_reg_method)) { + if (is.null(log_reg_methods)) { stop("Please provide one logistic regression method in the recipe. Accepted ", "methods are 'ens_mean', 'ens_mean_sd', 'sorted_members'.") } @@ -214,31 +223,37 @@ downscale_datasets <- function(data, recipe) { stop("Please provide the target grid in the recipe.") } - if (!(log_reg_method %in% LOG_REG_METHODS)) { - stop("Logistic regression method in the recipe is not available. Accepted methods ", - "are 'ens_mean', 'ens_mean_sd', 'sorted_members'.") - } + hcst_downscal <- list() + for (log_reg_method in strsplit(log_reg_methods, ", | |,")[[1]]) { + + if (!(log_reg_method %in% LOG_REG_METHODS)) { + stop(paste0(log_reg_method, " method in the recipe is not available. Accepted methods ", + "are 'ens_mean', 'ens_mean_sd', 'sorted_members'.")) + } - hcst_downscal <- CST_LogisticReg(data$hcst$data, data$obs$data, - target_grid = target_grid, - int_method = int_method, - log_reg_method = log_reg_method, - probs_cat = c(1/3,2/3), - return_most_likely_cat = FALSE, - points = NULL, - method_point_interp = NULL, - lat_dim = "latitude", - lon_dim = "longitude", - sdate_dim = "syear", - member_dim = "ensemble", - region = NULL, - loocv = loocv, - ncores = ncores) + hcst_d <- CST_LogisticReg(data$hcst, data$obs, + target_grid = target_grid, + int_method = int_methods, + log_reg_method = log_reg_method, + probs_cat = c(1/3,2/3), + return_most_likely_cat = FALSE, + points = NULL, + method_point_interp = NULL, + lat_dim = "latitude", + lon_dim = "longitude", + sdate_dim = "syear", + member_dim = "ensemble", + region = NULL, + loocv = loocv, + ncores = ncores) + + hcst_downscal[[ log_reg_method ]] <- hcst_d + } DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" } - - print(DOWNSCAL_MSG) - return(list(hcst = hcst_downscal, fcst = NULL)) + } + print(DOWNSCAL_MSG) + return(list(hcst = hcst_downscal, fcst = NULL)) } diff --git a/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml b/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml index 8b56f3aa..a6549863 100644 --- a/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml +++ b/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml @@ -25,11 +25,10 @@ Analysis: name: era5 Time: sdate: '0501' - fcst_year: '2020' hcst_start: '2000' hcst_end: '2004' ftime_min: 1 - ftime_max: 3 + ftime_max: 1 #sometimes we want the region for the obs to be bigger than for hcst Region: latmin: 34.1 @@ -43,17 +42,16 @@ Analysis: Calibration: method: raw Skill: - metric: RPSS + metric: RPSS EnsCorr BSS90 Mean_Bias_SS Downscaling: # Assumption 1: leave-one-out cross-validation is always applied # TO DO: add downscaling to point locations - # TO DO: how can we pass the parameter 'predictors'? - type: # mandatory, 'none', 'int', 'intbc', 'intlr', 'analogs', 'logreg' - int_method: # optional, regridding method accepted by CDO - bc_method: # optional, 'simple_bias', 'calibration', 'quantile_mapping' + type: intbc # mandatory, 'none', 'int', 'intbc', 'intlr', 'analogs', 'logreg' + int_method: bilinear # optional, regridding method accepted by CDO + bc_method: simple_bias # optional, 'simple_bias', 'calibration', 'quantile_mapping' lr_method: # optional, 'basic', 'large_scale', '4nn' log_reg_method: # optional, 'ens_mean', 'ens_mean_sd', 'sorted_members' - target_grid: # optional, nc file or grid accepted by CDO + target_grid: /esarchive/recon/ecmwf/era5/monthly_mean/tas_f1h/tas_200002.nc # optional, nc file or grid accepted by CDO nanalogs: # optional, number of analogs to be searched Output_format: S2S4E Run: diff --git a/modules/test_seasonal_downscaling.R b/modules/test_seasonal_downscaling.R index 9043c2a9..2f9ffbd0 100644 --- a/modules/test_seasonal_downscaling.R +++ b/modules/test_seasonal_downscaling.R @@ -1,30 +1,24 @@ -# DO source ../MODULES before opening the R session +# DO source /esarchive/scratch/jramon/GitLab_jramon/auto-s2s/MODULES before opening the R session setwd("/esarchive/scratch/jramon/GitLab_jramon/auto-s2s/") source("modules/Loading/Loading.R") -source("modules/Loading/Downscaling.R") +source("modules/Downscaling/Downscaling.R") source("modules/Calibration/Calibration.R") source("modules/Skill/Skill.R") source("modules/Saving/Saving.R") source("modules/Visualization/Visualization.R") recipe_file <- "modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml" -recipe <- read_yaml(recipe_file) -archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive.yml"))$archive +recipe <- prepare_outputs(recipe_file) +#archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive.yml"))$archive # Load datasets data <- load_datasets(recipe) # Downscale datasets -downscaled_data <- downscale_datasets(data, recipe) -# Calibrate datasets -calibrated_data <- calibrate_datasets(data, recipe) -# Compute skill metrics -skill_metrics <- compute_skill_metrics(calibrated_data$hcst, data$obs, recipe) -# Compute percentiles and probability bins -probabilities <- compute_probabilities(calibrated_data$hcst, recipe) -# Export all data to netCDF -save_data(recipe, archive, data, calibrated_data, skill_metrics, probabilities) +downscaled_data <- downscale_datasets(recipe, data) +# Compute skill +skill_metrics <- compute_skill_metrics(recipe, downscaled_data$hcst$simple_bias, data$obs) # Plot data -plot_data(recipe, archive, data, calibrated_data, skill_metrics, - probabilities, significance = T) +plot_data(recipe, calibrated_data = NULL, skill_metrics = NULL, + probabilities = NULL, significance = F) -- GitLab From 99bf98713140707e445a4c29182d17cd70497671 Mon Sep 17 00:00:00 2001 From: Jaume Ramon Date: Tue, 29 Nov 2022 11:57:10 +0100 Subject: [PATCH 6/8] Workflow work arounded --- modules/test_seasonal_downscaling.R | 45 ++++++++++++++++++++++++++--- 1 file changed, 41 insertions(+), 4 deletions(-) diff --git a/modules/test_seasonal_downscaling.R b/modules/test_seasonal_downscaling.R index 2f9ffbd0..898af7eb 100644 --- a/modules/test_seasonal_downscaling.R +++ b/modules/test_seasonal_downscaling.R @@ -9,16 +9,53 @@ source("modules/Skill/Skill.R") source("modules/Saving/Saving.R") source("modules/Visualization/Visualization.R") +# use this development to run the analogs function. Also change the Interpolation.R script to +# not load CDORemap from s2dv. +library(easyNCDF) +library(ncdf4) +library(ClimProjDiags) +source('https://earth.bsc.es/gitlab/es/s2dv/-/raw/master/R/Utils.R') +source('https://earth.bsc.es/gitlab/es/s2dv/-/raw/master/R/CDORemap.R') + recipe_file <- "modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml" recipe <- prepare_outputs(recipe_file) -#archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive.yml"))$archive # Load datasets data <- load_datasets(recipe) + # Downscale datasets downscaled_data <- downscale_datasets(recipe, data) + # Compute skill -skill_metrics <- compute_skill_metrics(recipe, downscaled_data$hcst$simple_bias, data$obs) +# NOTE: LogReg method does not have member dimension because it provides directly the +# predicted categories +# The following line should work only for Interpolation, Intbc, Intlr, Analogs +stopifnot(recipe$Analysis$Workflow$Downscaling$type %in% c('int', 'intbc', 'intlr', 'analogs')) +skill_metrics <- lapply(downscaled_data$hcst, function(x) compute_skill_metrics(recipe, x$exp, x$obs)) + +# Plot downscaled data +lapply(downscaled_data$hcst, function(x) plot_data(recipe, downscaled_data = x, data = NULL)) + +# Plot skill data +my_data <- list() +my_data$hcst <- downscaled_data$hcst$bilinear$exp +plot_data(recipe, skill_metrics = skill_metrics$bilinear, data = my_data) #Plot skill metrics + +################################## +# EXAMPLE FOR SEASONAL CALIBRATION +################################## +#data <- load_datasets(recipe) +# Calibrate datasets +#calibrated_data <- calibrate_datasets(recipe, data) +# Compute skill metrics +#skill_metrics <- compute_skill_metrics(recipe, calibrated_data$hcst, data$obs) +# Compute percentiles and probability bins +#probabilities <- compute_probabilities(recipe, calibrated_data$hcst) +# Export all data to netCDF +#save_data(recipe, data, calibrated_data, skill_metrics, probabilities) # Plot data -plot_data(recipe, calibrated_data = NULL, skill_metrics = NULL, - probabilities = NULL, significance = F) +#plot_data(recipe, data, calibrated_data, skill_metrics, probabilities, +# significance = T) + + + -- GitLab From 35476e7d5dd22e3f0f03136bf183ac5d938a18f9 Mon Sep 17 00:00:00 2001 From: Jaume Ramon Date: Tue, 29 Nov 2022 17:18:49 +0100 Subject: [PATCH 7/8] Scripts adapted to work with several downscaling methods --- modules/Downscaling/Downscaling.R | 80 ++++++++++++++----- .../recipe_system5c3s-tas_downscaling.yml | 11 +-- modules/Saving/paths2save.R | 24 ++++-- modules/Visualization/Visualization.R | 29 +++++-- modules/test_seasonal_downscaling.R | 7 +- 5 files changed, 115 insertions(+), 36 deletions(-) diff --git a/modules/Downscaling/Downscaling.R b/modules/Downscaling/Downscaling.R index 4d3edd2e..c3d063b4 100644 --- a/modules/Downscaling/Downscaling.R +++ b/modules/Downscaling/Downscaling.R @@ -13,6 +13,7 @@ source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Interp source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Intbc.R') source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Intlr.R') source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Analogs.R') +source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/LogisticReg.R') source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Utils.R') source("https://earth.bsc.es/gitlab/external/cstools/-/raw/master/R/CST_BiasCorrection.R") @@ -59,14 +60,14 @@ downscale_datasets <- function(recipe, data) { loocv <- recipe$Analysis$loocv } - DOWNSCAL_TYPES <- c("none", "int", "intbc", "intlr", "analogs", "logres") + DOWNSCAL_TYPES <- c("none", "int", "intbc", "intlr", "analogs", "logreg") BC_METHODS <- c("simple_bias", "calibration", "quantile_mapping") LR_METHODS <- c("basic", "large_scale", "4nn") LOG_REG_METHODS <- c("ens_mean", "ens_mean_sd", "sorted_members") if (!(type %in% DOWNSCAL_TYPES)) { stop("Downscaling type in the recipe is not available. Accepted types ", - "are 'none', 'int', 'intbc', 'intlr', 'analogs', 'logres'.") + "are 'none', 'int', 'intbc', 'intlr', 'analogs', 'logreg'.") } if (type == "int") { @@ -81,21 +82,43 @@ downscale_datasets <- function(recipe, data) { hcst_downscal <- list() for (int_method in strsplit(int_methods, ", | |,")[[1]]) { + # Ensure that observations are in the same grid as experiments + # Only needed for this method because the others already return the + # observations + latmin <- data$hcst$lat[1] + lonmin <- data$hcst$lon[1] + latmax <- data$hcst$lat[length(data$hcst$lat)] + lonmax <- data$hcst$lon[length(data$hcst$lon)] hcst_d <- CST_Interpolation(data$hcst, points = NULL, method_remap = int_method, target_grid = target_grid, lat_dim = "latitude", lon_dim = "longitude", - region = NULL, + region = c(lonmin, lonmax, latmin, latmax), method_point_interp = NULL) + obs_d <- CST_Interpolation(data$obs, + points = NULL, + method_remap = int_method, + target_grid = target_grid, + lat_dim = "latitude", + lon_dim = "longitude", + region = c(lonmin, lonmax, latmin, latmax), + method_point_interp = NULL) + + + hcst_d$obs <- obs_d$exp hcst_downscal[[ int_method ]] <- hcst_d } DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" } else if (type == "intbc") { - if (is.null(int_methods) | (length(strsplit(int_methods, ", | |,")[[1]]) != 1)) { + if (length(int_methods) == 0) { + stop("Please provide one (and only one) interpolation method in the recipe.") + } + + if (length(strsplit(int_methods, ", | |,")[[1]]) != 1) { stop("Please provide one (and only one) interpolation method in the recipe.") } @@ -134,8 +157,11 @@ downscale_datasets <- function(recipe, data) { DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" } else if (type == "intlr") { + if (length(int_methods) == 0) { + stop("Please provide one (and only one) interpolation method in the recipe.") + } - if (is.null(int_methods) | (length(strsplit(int_methods, ", | |,")[[1]]) != 1)) { + if (length(strsplit(int_methods, ", | |,")[[1]]) != 1) { stop("Please provide one (and only one) interpolation method in the recipe.") } @@ -176,6 +202,7 @@ downscale_datasets <- function(recipe, data) { lon_dim = "longitude", sdate_dim = "syear", time_dim = "time", + member_dim = "ensemble", large_scale_predictor_dimname = 'vars', loocv = loocv, region = NULL, @@ -193,24 +220,29 @@ downscale_datasets <- function(recipe, data) { nanalogs <- 3 } - - hcst_downscal <- CST_Analogs(data$hcst, data$obs, - grid_exp = data$hcst$source_files[1], - nanalogs = nanalogs, - fun_analog = NULL, - lat_dim = "latitude", - lon_dim = "longitude", - sdate_dim = "syear", - time_dim = "time", - region = NULL, - return_indices = FALSE, - loocv_window = loocv, - ncores = ncores) + hcst_downscal <- list() + hcst_downscal[["analogs"]] <- CST_Analogs(data$hcst, data$obs, + grid_exp = data$hcst$source_files[1], + nanalogs = nanalogs, + fun_analog = "min", + lat_dim = "latitude", + lon_dim = "longitude", + sdate_dim = "syear", + time_dim = "time", + member_dim = "ensemble", + region = NULL, + return_indices = FALSE, + loocv_window = loocv, + ncores = ncores) DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" - } else if (type == "logres") { + } else if (type == "logreg") { - if (is.null(int_methods) | (length(strsplit(int_methods, ", | |,")[[1]]) != 1)) { + if (length(int_methods) == 0) { + stop("Please provide one (and only one) interpolation method in the recipe.") + } + + if (length(strsplit(int_methods, ", | |,")[[1]]) != 1) { stop("Please provide one (and only one) interpolation method in the recipe.") } @@ -222,6 +254,14 @@ downscale_datasets <- function(recipe, data) { if (is.null(target_grid)) { stop("Please provide the target grid in the recipe.") } + + # Since we are forcing to create three categories, and applying cross-validation, + # we need at least six years of data for the logistic regression function to not + # crash + if (dim(data$hcst$data)[names(dim(data$hcst$data)) == "syear"] <= 5) { + stop("The number of start dates is insufficient for the logisitic regression method. ", + "Please provide six or more.") + } hcst_downscal <- list() for (log_reg_method in strsplit(log_reg_methods, ", | |,")[[1]]) { diff --git a/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml b/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml index a6549863..400448b7 100644 --- a/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml +++ b/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml @@ -26,7 +26,7 @@ Analysis: Time: sdate: '0501' hcst_start: '2000' - hcst_end: '2004' + hcst_end: '2005' ftime_min: 1 ftime_max: 1 #sometimes we want the region for the obs to be bigger than for hcst @@ -45,14 +45,15 @@ Analysis: metric: RPSS EnsCorr BSS90 Mean_Bias_SS Downscaling: # Assumption 1: leave-one-out cross-validation is always applied + # Assumption 2: for analogs, we select the best analog (minimum distance) # TO DO: add downscaling to point locations - type: intbc # mandatory, 'none', 'int', 'intbc', 'intlr', 'analogs', 'logreg' - int_method: bilinear # optional, regridding method accepted by CDO - bc_method: simple_bias # optional, 'simple_bias', 'calibration', 'quantile_mapping' + type: int # mandatory, 'none', 'int', 'intbc', 'intlr', 'analogs', 'logreg' + int_method: bilinear nn conservative # optional, regridding method accepted by CDO + bc_method: # optional, 'simple_bias', 'calibration', 'quantile_mapping' lr_method: # optional, 'basic', 'large_scale', '4nn' log_reg_method: # optional, 'ens_mean', 'ens_mean_sd', 'sorted_members' target_grid: /esarchive/recon/ecmwf/era5/monthly_mean/tas_f1h/tas_200002.nc # optional, nc file or grid accepted by CDO - nanalogs: # optional, number of analogs to be searched + nanalogs: 1 # optional, number of analogs to be searched Output_format: S2S4E Run: Loglevel: INFO diff --git a/modules/Saving/paths2save.R b/modules/Saving/paths2save.R index f48ebe7b..e34a6cf8 100644 --- a/modules/Saving/paths2save.R +++ b/modules/Saving/paths2save.R @@ -60,13 +60,27 @@ get_dir <- function(recipe, agg = "global") { calib.method <- tolower(recipe$Analysis$Workflow$Calibration$method) store.freq <- recipe$Analysis$Variables$freq - + downsc.type <- tolower(recipe$Analysis$Workflow$Downscaling$type) + if (downsc.type == "int") { + downsc.method <- strsplit(tolower(recipe$Analysis$Workflow$Downscaling$int_method), ", | |,")[[1]] + } else if (downsc.type == "intbc") { + downsc.method <- strsplit(tolower(recipe$Analysis$Workflow$Downscaling$bc_method), ", | |,")[[1]] + } else if (downsc.type == "intlr") { + downsc.method <- strsplit(tolower(recipe$Analysis$Workflow$Downscaling$lr_method), ", | |,")[[1]] + } else if (downsc.type == "analogs") { + downsc.method <- "analogs" + } else if (downsc.type == "logreg") { + downsc.method <- strsplit(tolower(recipe$Analysis$Workflow$Downscaling$log_reg_method), ", | |,")[[1]] + } + switch(tolower(agg), - "country" = {dir <- paste0(outdir, "/", calib.method, "-", - store.freq, "/", variable, + "country" = {dir <- paste0(outdir, calib.method, "-", + store.freq, "/", downsc.type, "/", + downsc.method, "/", variable, "_country/", fcst.sdate, "/")}, - "global" = {dir <- paste0(outdir, "/", calib.method, "-", - store.freq, "/", variable, "/", + "global" = {dir <- paste0(outdir, calib.method, "-", + store.freq, "/", downsc.type, "/", + downsc.method, "/", variable, "/", fcst.sdate, "/")}) return(dir) diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R index ff0e9fd4..17f58196 100644 --- a/modules/Visualization/Visualization.R +++ b/modules/Visualization/Visualization.R @@ -8,7 +8,8 @@ source("modules/Visualization/tmp/PlotCombinedMap.R") plot_data <- function(recipe, data, - calibrated_data = NULL, + downscaled_data = NULL, + calibrated_data = NULL, skill_metrics = NULL, probabilities = NULL, archive = NULL, @@ -23,11 +24,12 @@ plot_data <- function(recipe, # skill_metrics: list of arrays containing the computed skill metrics # significance: Bool. Whether to include significance dots where applicable - outdir <- paste0(get_dir(recipe), "/plots/") - dir.create(outdir, showWarnings = FALSE, recursive = TRUE) + outdir <- paste0(get_dir(recipe), "plots/") + sapply(outdir, dir.create, showWarnings = FALSE, recursive = TRUE) + #dir.create(outdir, showWarnings = FALSE, recursive = TRUE) if ((is.null(skill_metrics)) && (is.null(calibrated_data)) && - (is.null(data$fcst))) { + (is.null(data$fcst)) && (is.null(downscaled_data))) { stop("The Visualization module has been called, but there is no data ", "that can be plotted.") } @@ -47,6 +49,12 @@ plot_data <- function(recipe, plot_skill_metrics(recipe, archive, data$hcst, skill_metrics, outdir, significance) } + + # Plot forecast ensemble mean + if (!is.null(downscaled_data)) { + plot_downscaled_data(recipe, archive, downscaled_data, outdir) + #plot_ensemble_mean(recipe, archive, downscaled_data$exp, outdir) + } # Plot forecast ensemble mean if (!is.null(calibrated_data$fcst)) { @@ -69,6 +77,16 @@ plot_data <- function(recipe, } } +plot_downscaled_data <- function(recipe, archive, data_ls, outdir) { + + methods <- names(data_ls) + for (method in methods) { + outdir_method <- grep(pattern = method, x = outdir, value = TRUE) + print(outdir_method) + plot_ensemble_mean(recipe, archive, fcst = data_ls[[method]]$exp, outdir = outdir_method) + } +} + plot_skill_metrics <- function(recipe, archive, data_cube, skill_metrics, outdir, significance = F) { @@ -80,9 +98,10 @@ plot_skill_metrics <- function(recipe, archive, data_cube, skill_metrics, if (!is.list(skill_metrics) || is.null(names(skill_metrics))) { stop("The element 'skill_metrics' must be a list of named arrays.") } - + latitude <- data_cube$lat longitude <- data_cube$lon + system_name <- archive$System[[recipe$Analysis$Datasets$System$name]]$name hcst_period <- paste0(recipe$Analysis$Time$hcst_start, "-", recipe$Analysis$Time$hcst_end) diff --git a/modules/test_seasonal_downscaling.R b/modules/test_seasonal_downscaling.R index 898af7eb..12026bba 100644 --- a/modules/test_seasonal_downscaling.R +++ b/modules/test_seasonal_downscaling.R @@ -34,7 +34,12 @@ stopifnot(recipe$Analysis$Workflow$Downscaling$type %in% c('int', 'intbc', 'intl skill_metrics <- lapply(downscaled_data$hcst, function(x) compute_skill_metrics(recipe, x$exp, x$obs)) # Plot downscaled data -lapply(downscaled_data$hcst, function(x) plot_data(recipe, downscaled_data = x, data = NULL)) +plot_data(recipe, downscaled_data = downscaled_data$hcst, data = NULL) + +# Plot raw data +my_data <- list() +my_data$fcst <- data$hcst +plot_data(recipe, data = my_data) # Plot skill data my_data <- list() -- GitLab From 07d3f53501fd29cd11364ae6a3dcd3bd03ea0eb3 Mon Sep 17 00:00:00 2001 From: Jaume Ramon Date: Thu, 16 Feb 2023 15:06:04 +0100 Subject: [PATCH 8/8] last changes before git pull --- modules/Downscaling/Downscaling.R | 269 ++++++++---------- .../testing_recipes/recipe_system5c3s-tas.yml | 2 +- .../recipe_system5c3s-tas_downscaling.yml | 8 +- modules/Visualization/Visualization.R | 24 +- modules/test_seasonal_downscaling.R | 65 ++++- 5 files changed, 182 insertions(+), 186 deletions(-) diff --git a/modules/Downscaling/Downscaling.R b/modules/Downscaling/Downscaling.R index c3d063b4..1f814b74 100644 --- a/modules/Downscaling/Downscaling.R +++ b/modules/Downscaling/Downscaling.R @@ -15,7 +15,7 @@ source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Intlr. source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Analogs.R') source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/LogisticReg.R') source('/esarchive/scratch/jramon/GitLab_jramon/downscaling/csdownscale/R/Utils.R') -source("https://earth.bsc.es/gitlab/external/cstools/-/raw/master/R/CST_BiasCorrection.R") +#source("https://earth.bsc.es/gitlab/external/cstools/-/raw/master/R/CST_BiasCorrection.R") ## Entry params data and recipe? @@ -40,12 +40,12 @@ downscale_datasets <- function(recipe, data) { } else { # Downscaling function params - int_methods <- tolower(recipe$Analysis$Workflow$Downscaling$int_method) - bc_methods <- tolower(recipe$Analysis$Workflow$Downscaling$bc_method) - lr_methods <- tolower(recipe$Analysis$Workflow$Downscaling$lr_method) - log_reg_methods <- tolower(recipe$Analysis$Workflow$Downscaling$log_reg_method) - target_grid <- tolower(recipe$Analysis$Workflow$Downscaling$target_grid) - nanalogs <- as.numeric(recipe$Analysis$Workflow$Downscaling$nanalogs) + int_method <- tolower(recipe$Analysis$Workflow$Downscaling$int_method) + bc_method <- tolower(recipe$Analysis$Workflow$Downscaling$bc_method) + lr_method <- tolower(recipe$Analysis$Workflow$Downscaling$lr_method) + log_reg_method <- tolower(recipe$Analysis$Workflow$Downscaling$log_reg_method) + target_grid <- tolower(recipe$Analysis$Workflow$Downscaling$target_grid) + nanalogs <- as.numeric(recipe$Analysis$Workflow$Downscaling$nanalogs) if (is.null(recipe$Analysis$ncores)) { ncores <- 1 @@ -71,58 +71,48 @@ downscale_datasets <- function(recipe, data) { } if (type == "int") { - if (is.null(int_methods)) { - stop("Please provide one or more interpolation methods in the recipe.") + if (is.null(int_method)) { + stop("Please provide one interpolation method in the recipe.") } if (is.null(target_grid)) { stop("Please provide the target grid in the recipe.") } - hcst_downscal <- list() - for (int_method in strsplit(int_methods, ", | |,")[[1]]) { - - # Ensure that observations are in the same grid as experiments - # Only needed for this method because the others already return the - # observations - latmin <- data$hcst$lat[1] - lonmin <- data$hcst$lon[1] - latmax <- data$hcst$lat[length(data$hcst$lat)] - lonmax <- data$hcst$lon[length(data$hcst$lon)] - hcst_d <- CST_Interpolation(data$hcst, - points = NULL, - method_remap = int_method, - target_grid = target_grid, - lat_dim = "latitude", - lon_dim = "longitude", - region = c(lonmin, lonmax, latmin, latmax), - method_point_interp = NULL) - - obs_d <- CST_Interpolation(data$obs, - points = NULL, - method_remap = int_method, - target_grid = target_grid, - lat_dim = "latitude", - lon_dim = "longitude", - region = c(lonmin, lonmax, latmin, latmax), - method_point_interp = NULL) - - - hcst_d$obs <- obs_d$exp - hcst_downscal[[ int_method ]] <- hcst_d - } + # Ensure that observations are in the same grid as experiments + # Only needed for this method because the others already return the + # observations + latmin <- data$hcst$lat[1] + lonmin <- data$hcst$lon[1] + latmax <- data$hcst$lat[length(data$hcst$lat)] + lonmax <- data$hcst$lon[length(data$hcst$lon)] + hcst_downscal <- CST_Interpolation(data$hcst, + points = NULL, + method_remap = int_method, + target_grid = target_grid, + lat_dim = "latitude", + lon_dim = "longitude", + region = c(lonmin, lonmax, latmin, latmax), + method_point_interp = NULL) + + obs_downscal <- CST_Interpolation(data$obs, + points = NULL, + method_remap = int_method, + target_grid = target_grid, + lat_dim = "latitude", + lon_dim = "longitude", + region = c(lonmin, lonmax, latmin, latmax), + method_point_interp = NULL) + + hcst_downscal$obs <- obs_downscal$exp DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" } else if (type == "intbc") { - if (length(int_methods) == 0) { - stop("Please provide one (and only one) interpolation method in the recipe.") - } - - if (length(strsplit(int_methods, ", | |,")[[1]]) != 1) { + if (length(int_method) == 0) { stop("Please provide one (and only one) interpolation method in the recipe.") } - if (is.null(bc_methods)) { + if (is.null(bc_method)) { stop("Please provide one bias-correction method in the recipe. Accepted ", "methods are 'simple_bias', 'calibration', 'quantile_mapping'.") } @@ -131,41 +121,31 @@ downscale_datasets <- function(recipe, data) { stop("Please provide the target grid in the recipe.") } - hcst_downscal <- list() - for (bc_method in strsplit(bc_methods, ", | |,")[[1]]) { - - if (!(bc_method %in% BC_METHODS)) { - stop(paste0(bc_method, " method in the recipe is not available. Accepted methods ", - "are 'simple_bias', 'calibration', 'quantile_mapping'.")) - } + if (!(bc_method %in% BC_METHODS)) { + stop(paste0(bc_method, " method in the recipe is not available. Accepted methods ", + "are 'simple_bias', 'calibration', 'quantile_mapping'.")) + } - hcst_d <- CST_Intbc(data$hcst, data$obs, - target_grid = target_grid, - bc_method = bc_method, - int_method = int_methods, - points = NULL, - method_point_interp = NULL, - lat_dim = "latitude", - lon_dim = "longitude", - sdate_dim = "syear", - member_dim = "ensemble", - region = NULL, - ncores = ncores) + hcst_downscal <- CST_Intbc(data$hcst, data$obs, + target_grid = target_grid, + bc_method = bc_method, + int_method = int_method, + points = NULL, + method_point_interp = NULL, + lat_dim = "latitude", + lon_dim = "longitude", + sdate_dim = "syear", + member_dim = "ensemble", + region = NULL, + ncores = ncores) - hcst_downscal[[ bc_method ]] <- hcst_d - } - DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" } else if (type == "intlr") { - if (length(int_methods) == 0) { - stop("Please provide one (and only one) interpolation method in the recipe.") - } - - if (length(strsplit(int_methods, ", | |,")[[1]]) != 1) { + if (length(int_method) == 0) { stop("Please provide one (and only one) interpolation method in the recipe.") } - if (is.null(lr_methods)) { + if (is.null(lr_method)) { stop("Please provide one linear regression method in the recipe. Accepted ", "methods are 'basic', 'large_scale', '4nn'.") } @@ -174,43 +154,37 @@ downscale_datasets <- function(recipe, data) { stop("Please provide the target grid in the recipe.") } - hcst_downscal <- list() - for (lr_method in strsplit(lr_methods, ", | |,")[[1]]) { - - if (!(lr_method %in% LR_METHODS)) { - stop(paste0(lr_method, " method in the recipe is not available. Accepted methods ", - "are 'basic', 'large_scale', '4nn'.")) - } + if (!(lr_method %in% LR_METHODS)) { + stop(paste0(lr_method, " method in the recipe is not available. Accepted methods ", + "are 'basic', 'large_scale', '4nn'.")) + } - # TO DO: add the possibility to have the element 'pred' in 'data' - if (lr_method == "large_scale") { - if (is.null(data$pred$data)) { - stop("Please provide the large scale predictors in the element 'data$pred$data'.") - } - } else { - data$pred$data <- NULL + # TO DO: add the possibility to have the element 'pred' in 'data' + if (lr_method == "large_scale") { + if (is.null(data$pred$data)) { + stop("Please provide the large scale predictors in the element 'data$pred$data'.") } + } else { + data$pred$data <- NULL + } - hcst_d <- CST_Intlr(data$hcst, data$obs, - lr_method = lr_method, - target_grid = target_grid, - points = NULL, - int_method = int_methods, - method_point_interp = NULL, - predictors = data$pred$data, - lat_dim = "latitude", - lon_dim = "longitude", - sdate_dim = "syear", - time_dim = "time", - member_dim = "ensemble", - large_scale_predictor_dimname = 'vars', - loocv = loocv, - region = NULL, - ncores = ncores) + hcst_downscal <- CST_Intlr(data$hcst, data$obs, + lr_method = lr_method, + target_grid = target_grid, + points = NULL, + int_method = int_method, + method_point_interp = NULL, + predictors = data$pred$data, + lat_dim = "latitude", + lon_dim = "longitude", + sdate_dim = "syear", + time_dim = "time", + member_dim = "ensemble", + large_scale_predictor_dimname = 'vars', + loocv = loocv, + region = NULL, + ncores = ncores) - hcst_downscal[[ lr_method ]] <- hcst_d - } - DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" } else if (type == "analogs") { @@ -220,33 +194,28 @@ downscale_datasets <- function(recipe, data) { nanalogs <- 3 } - hcst_downscal <- list() - hcst_downscal[["analogs"]] <- CST_Analogs(data$hcst, data$obs, - grid_exp = data$hcst$source_files[1], - nanalogs = nanalogs, - fun_analog = "min", - lat_dim = "latitude", - lon_dim = "longitude", - sdate_dim = "syear", - time_dim = "time", - member_dim = "ensemble", - region = NULL, - return_indices = FALSE, - loocv_window = loocv, - ncores = ncores) + hcst_downscal <- CST_Analogs(data$hcst, data$obs, + grid_exp = data$hcst$source_files[1], + nanalogs = nanalogs, + fun_analog = "min", + lat_dim = "latitude", + lon_dim = "longitude", + sdate_dim = "syear", + time_dim = "time", + member_dim = "ensemble", + region = NULL, + return_indices = FALSE, + loocv_window = loocv, + ncores = ncores) DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" } else if (type == "logreg") { - if (length(int_methods) == 0) { + if (length(int_method) == 0) { stop("Please provide one (and only one) interpolation method in the recipe.") } - if (length(strsplit(int_methods, ", | |,")[[1]]) != 1) { - stop("Please provide one (and only one) interpolation method in the recipe.") - } - - if (is.null(log_reg_methods)) { + if (is.null(log_reg_method)) { stop("Please provide one logistic regression method in the recipe. Accepted ", "methods are 'ens_mean', 'ens_mean_sd', 'sorted_members'.") } @@ -263,33 +232,27 @@ downscale_datasets <- function(recipe, data) { "Please provide six or more.") } - hcst_downscal <- list() - for (log_reg_method in strsplit(log_reg_methods, ", | |,")[[1]]) { - - if (!(log_reg_method %in% LOG_REG_METHODS)) { - stop(paste0(log_reg_method, " method in the recipe is not available. Accepted methods ", - "are 'ens_mean', 'ens_mean_sd', 'sorted_members'.")) - } - - hcst_d <- CST_LogisticReg(data$hcst, data$obs, - target_grid = target_grid, - int_method = int_methods, - log_reg_method = log_reg_method, - probs_cat = c(1/3,2/3), - return_most_likely_cat = FALSE, - points = NULL, - method_point_interp = NULL, - lat_dim = "latitude", - lon_dim = "longitude", - sdate_dim = "syear", - member_dim = "ensemble", - region = NULL, - loocv = loocv, - ncores = ncores) - - hcst_downscal[[ log_reg_method ]] <- hcst_d + if (!(log_reg_method %in% LOG_REG_METHODS)) { + stop(paste0(log_reg_method, " method in the recipe is not available. Accepted methods ", + "are 'ens_mean', 'ens_mean_sd', 'sorted_members'.")) } + hcst_downscal <- CST_LogisticReg(data$hcst, data$obs, + target_grid = target_grid, + int_method = int_method, + log_reg_method = log_reg_method, + probs_cat = c(1/3,2/3), + return_most_likely_cat = FALSE, + points = NULL, + method_point_interp = NULL, + lat_dim = "latitude", + lon_dim = "longitude", + sdate_dim = "syear", + member_dim = "ensemble", + region = NULL, + loocv = loocv, + ncores = ncores) + DOWNSCAL_MSG <- "##### DOWNSCALING COMPLETE #####" } } diff --git a/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml b/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml index 38e23c8b..9ed0bf80 100644 --- a/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml +++ b/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml @@ -40,6 +40,6 @@ Analysis: Run: Loglevel: INFO Terminal: yes - output_dir: /esarchive/scratch/vagudets/repos/auto-s2s/out-logs/ + output_dir: /esarchive/scratch/jramon/GitLab_jramon/auto-s2s/out-logs/ code_dir: /esarchive/scratch/vagudets/repos/auto-s2s/ diff --git a/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml b/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml index 400448b7..0cad69ec 100644 --- a/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml +++ b/modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml @@ -42,18 +42,18 @@ Analysis: Calibration: method: raw Skill: - metric: RPSS EnsCorr BSS90 Mean_Bias_SS + metric: BSS10 BSS90 Mean_Bias_SS RPSS ROCSS Downscaling: # Assumption 1: leave-one-out cross-validation is always applied # Assumption 2: for analogs, we select the best analog (minimum distance) # TO DO: add downscaling to point locations - type: int # mandatory, 'none', 'int', 'intbc', 'intlr', 'analogs', 'logreg' - int_method: bilinear nn conservative # optional, regridding method accepted by CDO + type: # mandatory, 'none', 'int', 'intbc', 'intlr', 'analogs', 'logreg' + int_method: # optional, regridding method accepted by CDO bc_method: # optional, 'simple_bias', 'calibration', 'quantile_mapping' lr_method: # optional, 'basic', 'large_scale', '4nn' log_reg_method: # optional, 'ens_mean', 'ens_mean_sd', 'sorted_members' target_grid: /esarchive/recon/ecmwf/era5/monthly_mean/tas_f1h/tas_200002.nc # optional, nc file or grid accepted by CDO - nanalogs: 1 # optional, number of analogs to be searched + nanalogs: 4 # optional, number of analogs to be searched Output_format: S2S4E Run: Loglevel: INFO diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R index 17f58196..61f5de9a 100644 --- a/modules/Visualization/Visualization.R +++ b/modules/Visualization/Visualization.R @@ -8,8 +8,8 @@ source("modules/Visualization/tmp/PlotCombinedMap.R") plot_data <- function(recipe, data, - downscaled_data = NULL, calibrated_data = NULL, + downscaled_data = NULL, skill_metrics = NULL, probabilities = NULL, archive = NULL, @@ -25,8 +25,7 @@ plot_data <- function(recipe, # significance: Bool. Whether to include significance dots where applicable outdir <- paste0(get_dir(recipe), "plots/") - sapply(outdir, dir.create, showWarnings = FALSE, recursive = TRUE) - #dir.create(outdir, showWarnings = FALSE, recursive = TRUE) + dir.create(outdir, showWarnings = FALSE, recursive = TRUE) if ((is.null(skill_metrics)) && (is.null(calibrated_data)) && (is.null(data$fcst)) && (is.null(downscaled_data))) { @@ -49,12 +48,11 @@ plot_data <- function(recipe, plot_skill_metrics(recipe, archive, data$hcst, skill_metrics, outdir, significance) } - - # Plot forecast ensemble mean + + # Plot downscaled data if (!is.null(downscaled_data)) { - plot_downscaled_data(recipe, archive, downscaled_data, outdir) - #plot_ensemble_mean(recipe, archive, downscaled_data$exp, outdir) - } + plot_ensemble_mean(recipe, archive, downscaled_data$hcst$exp, outdir) + } # Plot forecast ensemble mean if (!is.null(calibrated_data$fcst)) { @@ -77,16 +75,6 @@ plot_data <- function(recipe, } } -plot_downscaled_data <- function(recipe, archive, data_ls, outdir) { - - methods <- names(data_ls) - for (method in methods) { - outdir_method <- grep(pattern = method, x = outdir, value = TRUE) - print(outdir_method) - plot_ensemble_mean(recipe, archive, fcst = data_ls[[method]]$exp, outdir = outdir_method) - } -} - plot_skill_metrics <- function(recipe, archive, data_cube, skill_metrics, outdir, significance = F) { diff --git a/modules/test_seasonal_downscaling.R b/modules/test_seasonal_downscaling.R index 12026bba..5f375abe 100644 --- a/modules/test_seasonal_downscaling.R +++ b/modules/test_seasonal_downscaling.R @@ -17,12 +17,27 @@ library(ClimProjDiags) source('https://earth.bsc.es/gitlab/es/s2dv/-/raw/master/R/Utils.R') source('https://earth.bsc.es/gitlab/es/s2dv/-/raw/master/R/CDORemap.R') +# Loop all downscaling methods and types +#types <- c('int', 'intbc', 'intlr', 'analogs', 'logreg') +#int_methods <- c('con', 'bil', 'bic', 'nn', 'con2') +#bc_methods <- c('simple_bias', 'calibration', 'quantile_mapping') +#lr_methods <- c('basic', 'large_scale', '4nn') +#logreg_methods <- c('ens_mean', 'ens_mean_sd', 'sorted_members') + +# Read recipe recipe_file <- "modules/Loading/testing_recipes/recipe_system5c3s-tas_downscaling.yml" recipe <- prepare_outputs(recipe_file) # Load datasets data <- load_datasets(recipe) +# Modify the recipe manually +recipe$Analysis$Workflow$Downscaling$type <- 'analogs' +recipe$Analysis$Workflow$Downscaling$int_method <- NULL # con, bil, bic, nn, con2 +recipe$Analysis$Workflow$Downscaling$bc_method <- NULL +recipe$Analysis$Workflow$Downscaling$lr_method <- NULL +recipe$Analysis$Workflow$Downscaling$log_reg_method <- NULL + # Downscale datasets downscaled_data <- downscale_datasets(recipe, data) @@ -30,21 +45,51 @@ downscaled_data <- downscale_datasets(recipe, data) # NOTE: LogReg method does not have member dimension because it provides directly the # predicted categories # The following line should work only for Interpolation, Intbc, Intlr, Analogs -stopifnot(recipe$Analysis$Workflow$Downscaling$type %in% c('int', 'intbc', 'intlr', 'analogs')) -skill_metrics <- lapply(downscaled_data$hcst, function(x) compute_skill_metrics(recipe, x$exp, x$obs)) +#stopifnot(recipe$Analysis$Workflow$Downscaling$type %in% c('int', 'intbc', 'intlr', 'analogs')) +#skill_metrics <- compute_skill_metrics(recipe, downscaled_data$hcst$exp, downscaled_data$hcst$obs) + +# Save downscaled data to later intercompare them +saveRDS(downscaled_data, file = paste0('/esarchive/scratch/jramon/downscaling/data/FOCUS/', + recipe$Analysis$Workflow$Downscaling$type, '-', + recipe$Analysis$Workflow$Downscaling$int_method, '-', + recipe$Analysis$Workflow$Downscaling$bc_method, '-', + recipe$Analysis$Workflow$Downscaling$lr_method, '-', + recipe$Analysis$Workflow$Downscaling$log_reg_method, '.RDS')) # Plot downscaled data -plot_data(recipe, downscaled_data = downscaled_data$hcst, data = NULL) +plot_data(recipe, data = NULL, downscaled_data = downscaled_data) -# Plot raw data +# Plot skill metrics my_data <- list() -my_data$fcst <- data$hcst -plot_data(recipe, data = my_data) +my_data$hcst <- downscaled_data$hcst$exp +plot_data(recipe, data = my_data, skill_metrics = skill_metrics) -# Plot skill data -my_data <- list() -my_data$hcst <- downscaled_data$hcst$bilinear$exp -plot_data(recipe, skill_metrics = skill_metrics$bilinear, data = my_data) #Plot skill metrics +# Plot Taylor diagram intercomparison +source('/esarchive/scratch/jramon/GitLab_jramon/nc2ts/taylor_diagram.R') +rds_files <- list.files(path = '/esarchive/scratch/jramon/downscaling/data/FOCUS', full.names = T) +down_data <- list() +methods <- c(list.files(path = '/esarchive/scratch/jramon/downscaling/data/FOCUS')) +for (nfile in seq(rds_files)) { + + # to plot all points in the same Taylor diagram + flag <- ifelse(nfile == 1, F, T) + + down_data[[nfile]] <- readRDS(rds_files[[nfile]]) + + #compute model ensemble mean + ens_mean <- Apply(down_data[[nfile]]$hcst$exp$data, 'ensemble', mean)$output1 + + # remove ensemble dimension in obs + obs <- Subset(down_data[[nfile]]$hcst$obs$data, along = 'ensemble', indices = 1, drop = 'selected') + + ens_mean <- .reorder_dims(arr_ref = obs, arr_to_reorder = ens_mean) + + plot <- taylor.diagram(ref = as.vector(obs), model = as.vector(ens_mean), normalize = TRUE, + add = flag, pch = nfile, sd.arcs = TRUE, maxsd = 1.5) + +} + +legend(1.3, 1.5, cex = 0.8, pt.cex = 1, legend = methods, pch = 1:length(rds_files), ncol = 1) ################################## # EXAMPLE FOR SEASONAL CALIBRATION -- GitLab