diff --git a/modules/Calibration/test_calib.R b/modules/Calibration/test_calib.R index fa3f1359ae8b39fb2485d0097ae039e3ede092f4..8c6600f3e60bf4986c56597265531030cd4e8b16 100644 --- a/modules/Calibration/test_calib.R +++ b/modules/Calibration/test_calib.R @@ -4,6 +4,11 @@ source("https://earth.bsc.es/gitlab/external/cstools/-/raw/develop-CalibrationFo ## Entry params data and recipe? calibrate_datasets <- function(data, recipe) { + # Function that calibrates the hindcast using the method stated in the + # recipe. If the forecast is not null, it calibrates it as well. + # + # data: list of s2dv_cube objects containing the hcst, obs and fcst. + # recipe: object obtained when passing the .yml recipe file to read_yaml() hcst <- data$hcst obs <- data$obs fcst <- data$fcst @@ -32,10 +37,11 @@ calibrate_datasets <- function(data, recipe) { ## TODO: implement other calibration methods ## TODO: Restructure the code? if (!(method %in% CST_CALIB_METHODS)) { - stop("Calibration method is not implemented in CSTools") + 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 and forecast + # Calibrate the hindcast hcst_calibrated <- CST_Calibration(hcst, obs, cal.method = method, eval.method = "leave-one-out", @@ -68,11 +74,13 @@ calibrate_datasets <- function(data, recipe) { } else if (recipe$Analysis$Variables$freq == "daily_mean") { # Daily data calibration using Quantile Mapping - warning("Data is at daily frequency. Quantile Mapping will be used", - " as the calibration method.") + 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(hcst, obs, exp_cor = NULL, - ## TODO: Tidy up sample_dims = c("syear", "time", "ensemble"), @@ -82,6 +90,7 @@ calibrate_datasets <- function(data, recipe) { na.rm = na.rm) if (!is.null(fcst)) { + # Calibrate the forecast fcst_calibrated <- CST_QuantileMapping(hcst, obs, exp_cor = fcst, sample_dims = c("syear", @@ -96,6 +105,8 @@ calibrate_datasets <- function(data, recipe) { } } - print("###### CALIBRATION COMPLETE ######") - return(list(hcst_calib = hcst_calibrated, fcst_calib = fcst_calibrated)) + print("##### CALIBRATION COMPLETE #####") + ## TODO: Return observations too? + ## TODO: Change naming convention? + return(list(hcst = hcst_calibrated, fcst = fcst_calibrated)) } diff --git a/modules/data_load/load.R b/modules/data_load/load.R index 685881562eee6295e2181ecc7156fbc4e24baf47..3941f5dc0902a0ef2a9441735f0814dbc90153b3 100755 --- a/modules/data_load/load.R +++ b/modules/data_load/load.R @@ -5,60 +5,6 @@ source("/esarchive/scratch/vagudets/repos/cstools/R/s2dv_cube.R") source("modules/data_load/dates2load.R") source("tools/libs.R") -path <- "/esarchive/scratch/lpalma/git/startR-test-attrDims/R/" -ff <- lapply(list.files(path), function(x) paste0(path, x)) -invisible(lapply(ff, source)) - -## TODO: Move to CSOperational -# Conversion from startR_array to s2dv_array -#------------------------------------------------------------------- -attr2array <- function(attr){ - return(array(as.vector(attr), dim(attr))) -} - -startR_to_s2dv <- function(startR_array){ - - dates_dims <- dim(Subset(startR_array, - along=c('dat', 'var', - 'latitude', 'longitude', 'ensemble'), - list(1,1,1,1,1), drop="selected")) - - #sdates_dims <- dim(Subset(startR_array, - # along=c('dat', 'var', 'time', 'sweek', 'sday', - # 'latitude', 'longitude', 'ensemble'), - # list(1,1,1,1,1,1,1,1), drop="selected")) - - dates_start <- attr(startR_array, 'Variables')$common$time - dates_end <- attr(startR_array, 'Variables')$common$time - #sdates <- unlist(attr(startR_array, 'FileSelectors')$dat1$file_date) - - dim(dates_start) <- dates_dims - dim(dates_end) <- dates_dims - #dim(sdates) <- sdates_dims - - ## TODO: change this line when time attributes work correctly? - time_dims <- c("time", "sday", "sweek", "syear") - - s2dv_object <- s2dv_cube(data = attr2array(startR_array), - lon = attr2array(attr(startR_array, - 'Variables')$dat1$longitude), - lat = attr2array(attr(startR_array, - 'Variables')$dat1$latitude), - Variable = list(varName=names(attr(startR_array, - 'Variables')$common)[2], - level = NULL), - Dates = list(start = dates_start, - end = dates_end), - #sdate = sdates), - time_dims = time_dims, - when = Sys.time(), - source_files = attr2array(attr(startR_array, "Files")) - # Datasets=list(exp1=list(InitializationsDates=list( Member_1="01011990", - # Members="Member_1"))) - ) - return(s2dv_object) -} - # RECIPE FOR TESTING # -------------------------------------------------------------------------------- # recipe_file <- "modules/data_load/testing_recipes/recipe_3.yml" @@ -199,7 +145,7 @@ load_datasets <- function(recipe_file){ default_dims[names(dim(hcst))] <- dim(hcst) dim(hcst) <- default_dims } - hcst <- startR_to_s2dv(hcst) + hcst <- as.s2dv_cube(hcst) # Load forecast #------------------------------------------------------------------- @@ -230,7 +176,6 @@ load_datasets <- function(recipe_file){ split_multiselected_dims = split_multiselected_dims, retrieve = TRUE) - ## TODO: Add 'sday' and 'sweek' dims to fcst if (recipe$Analysis$Variables$freq == "daily_mean") { # Adjusts dims for daily case, could be removed if startR allows # multidim split @@ -241,8 +186,8 @@ load_datasets <- function(recipe_file){ default_dims[names(dim(fcst))] <- dim(fcst) dim(fcst) <- default_dims } - # names(dim(fcst))[which(names(dim(fcst)) == 'file_date')] <- "syear" - fcst <- startR_to_s2dv(fcst) + + fcst <- as.s2dv_cube(fcst) } else { fcst <- NULL @@ -324,7 +269,7 @@ load_datasets <- function(recipe_file){ default_dims[names(dim(obs))] <- dim(obs) dim(obs) <- default_dims - obs <- startR_to_s2dv(obs) + obs <- as.s2dv_cube(obs) ############################################################################ # diff --git a/modules/data_load/testing_recipes/recipe_1.yml b/modules/data_load/testing_recipes/recipe_1.yml index 68715462779ad4f4a8b2a67dc80c75945df8ba11..d1aa13e6ba9cd5700003a42975664bd3b43fced5 100644 --- a/modules/data_load/testing_recipes/recipe_1.yml +++ b/modules/data_load/testing_recipes/recipe_1.yml @@ -2,37 +2,37 @@ Description: Author: V. Agudetse Analysis: - Horizon: Seasonal + Horizon: Seasonal # Mandatory, str: either subseasonal, seasonal, or decadal Variables: - name: tas - freq: daily_mean + name: tas # Mandatory, str: tas, + freq: daily_mean # Mandatory, str: either monthly_mean or daily_mean Datasets: System: - name: system5c3s - Multimodel: no + name: system5c3s # Mandatory, str: System codename + Multimodel: no # Mandatory, bool: Either yes/true or no/false Reference: - name: era5 + name: era5 # Mandatory, str: Reference codename Time: - sdate: - fcst_syear: '2020' - fcst_sday: '1101' - hcst_start: '1993' - hcst_end: '2016' - leadtimemin: 1 - leadtimemax: 4 + sdate: + fcst_syear: '2020' # Optional, int: Forecast year 'YYYY' + fcst_sday: '1101' # Mandatory, int: Forecast startdate, 'MMDD' + hcst_start: '1993' # Mandatory, int: Hindcast start year 'YYYY' + hcst_end: '2016' # Mandatory, int: Hindcast end year 'YYYY' + leadtimemin: 1 # Mandatory, int: First leadtime time step in months + leadtimemax: 4 # Mandatory, int: Last leadtime time step in months Region: - latmin: -10 - latmax: 10 - lonmin: 0 - lonmax: 20 + latmin: -10 # Mandatory, int: minimum latitude + latmax: 10 # Mandatory, int: maximum latitude + lonmin: 0 # Mandatory, int: minimum longitude + lonmax: 20 # Mandatory, int: maximum longitude Regrid: - method: bilinear - type: to_system + method: bilinear # Mandatory, str: Interpolation method. See docu. + type: to_system # Mandatory, str: to_system, to_reference, or CDO-compatible grid. Workflow: Calibration: - method: SBC + method: qmap # Mandatory, str: Calibration method. See docu. Skill: - metric: RPSS + metric: RPSS #Mandatory, str: Skill metric or list of skill metrics. See docu. Indicators: index: no Output_format: S2S4E diff --git a/modules/data_load/testing_recipes/recipe_3.yml b/modules/data_load/testing_recipes/recipe_3.yml index da1bd78a1a3298320439124016f6f773bb2c0637..4488d43d9282c7bc3190e48524c06e6f7f8e3933 100644 --- a/modules/data_load/testing_recipes/recipe_3.yml +++ b/modules/data_load/testing_recipes/recipe_3.yml @@ -30,7 +30,7 @@ Analysis: type: to_system Workflow: Calibration: - method: mse_min + method: qmap Skill: metric: RPSS Indicators: diff --git a/modules/test_victoria.R b/modules/test_victoria.R index 96a0e6bb6972a93a8dc0b192b103ff41888bc0c4..6b90d2e93b2cbe2265a8b47783d232bb47c04f14 100644 --- a/modules/test_victoria.R +++ b/modules/test_victoria.R @@ -12,6 +12,6 @@ recipe <- read_yaml(recipe_file) # Calibrate datasets calibrated_data <- calibrate_datasets(data, recipe) # Compute skill metric -skill_metrics <- compute_skill_metrics(calibrated_data$hcst_calib, data$obs, +skill_metrics <- compute_skill_metrics(calibrated_data$hcst, data$obs, recipe, na.rm = T, ncores = 4)