From c2fa1b6f2769f0482442a700a2f9646646f4794b Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 11 Jun 2025 16:17:18 +0200 Subject: [PATCH 1/6] editing Crossval_calibration to use Apply --- modules/Crossval/Crossval_calibration.R | 504 ++++++++++++++++++++++++ modules/Crossval/R/tmp/MergeDims.R | 163 ++++++++ 2 files changed, 667 insertions(+) create mode 100644 modules/Crossval/Crossval_calibration.R create mode 100644 modules/Crossval/R/tmp/MergeDims.R diff --git a/modules/Crossval/Crossval_calibration.R b/modules/Crossval/Crossval_calibration.R new file mode 100644 index 00000000..26e7e31c --- /dev/null +++ b/modules/Crossval/Crossval_calibration.R @@ -0,0 +1,504 @@ +# take the output of Flor/s2s/subseasonal_loading.R + +source("modules/Crossval/R/tmp/GetProbs.R") +source("modules/Crossval/R/tmp/MergeDims.R") +Crossval_calibration <- function(recipe, data, correct_negative = FALSE) { + + cal_method <- recipe$Analysis$Workflow$Calibration$method + cross.method <- recipe$Analysis$cross.method + # TODO move check + if (is.null(cross.method)) { + cross.method <- 'leave-one-out' + } + categories <- recipe$Analysis$Workflow$Probabilities$percentiles + categories <- lapply(categories, function (x) { + sapply(x, function(y) { + eval(parse(text = y))})}) + ncores <- recipe$Analysis$ncores + na.rm <- recipe$Analysis$remove_NAs + ## data dimensions + sdate_dim <- dim(data$hcst$data)['syear'] + orig_dims <- names(dim(data$hcst$data)) + # spatial dims + if ('latitude' %in% names(dim(data$hcst$data))) { + nlats <- dim(data$hcst$data)['latitude'] + nlons <- dim(data$hcst$data)['longitude'] + agg = 'global' + } else if ('region' %in% names(dim(data$hcst$data))) { + agg = 'region' + nregions <- dim(data$hcst$data)['region'] + } + # TODO fix it to use new version https://earth.bsc.es/gitlab/external/cstools/-/blob/dev-cross-indices/R/CST_Calibration.R#L570 + cross <- CSTools:::.make.eval.train.dexes(eval.method = cross.method, + amt.points = sdate_dim, + amt.points_cor = NULL) # k = ? + .cross_calibration <- function(indices, exp, obs, cross, categories, recipe) { + ## Objects to return: + # lims_cal_hcst_tr: hcst percentiles (per category) + # lims_cal_obs_tr: obs percentiles (per category) + # cal_hcst_ev: hcst cross-validated calibrated + # cal_hcst_tr: hc + # cal_obs_tr: obs training +source("modules/Crossval/R/tmp/MergeDims.R") + ## training indices + cross <- cross[[indices]] + print(paste("Crossval Index:", cross$eval.dexes)) + obs_tr <- Subset(obs, along = 'syear', + indices = cross$train.dexes) + hcst_tr <- Subset(exp, along = 'syear', + indices = cross$train.dexes) + ## evaluation indices + hcst_ev <- Subset(exp, along = 'syear', + indices = cross$eval.dexes) + obs_ev <- Subset(obs, along = 'syear', + indices = cross$eval.dexes) + + if (tolower(recipe$Analysis$Horizon) == 'subseasonal') { + central_day <- (dim(exp)['sday'] + 1)/2 + hcst_tr <- MergeDims(hcst_tr, merge_dims = c('sday', 'syear'), + rename_dim = 'syear', na.rm = na.rm) + + obs_tr <- MergeDims(obs_tr, merge_dims = c('sday', 'syear'), + rename_dim = 'syear', na.rm = na.rm) + # 'sday' dim to select the central day + hcst_ev <- Subset(hcst_ev, along = 'sday', indices = central_day, + drop = 'selected') + + } + if (cal_method %in% c('PTF', 'DIST', 'RQUANT', 'QUANT', 'SSPLIN')) { + cal_hcst_tr <- QuantileMapping(exp = hcst_tr, obs = obs_tr, + method = cal_method, + memb_dim = 'ensemble', + sdate_dim = 'syear', + window_dim = NULL, + na.rm = na.rm, qstep = 0.1, + wet.day = FALSE, + ncores = ncores) + cal_hcst_ev <- QuantileMapping(exp = hcst_tr, obs = obs_tr, + exp_cor = hcst_ev, + method = cal_method, + memb_dim = 'ensemble', + sdate_dim = 'syear', + window_dim = NULL, + wet.day = FALSE, + na.rm = na.rm, qstep = 0.1, + ncores = ncores) + } else { + cal_hcst_tr <- Calibration(exp = hcst_tr, obs = obs_tr, + cal.method = cal_method, + memb_dim = 'ensemble', sdate_dim = 'syear', + eval.method = 'in-sample', + na.fill = TRUE, na.rm = na.rm, + apply_to = NULL, + alpha = NULL, ncores = ncores) + cal_hcst_ev <- Calibration(exp = hcst_tr, obs = obs_tr, exp_cor = hcst_ev, + cal.method = cal_method, + memb_dim = 'ensemble', sdate_dim = 'syear', + eval.method = 'in-sample', + na.fill = TRUE, na.rm = na.rm, + apply_to = NULL, + alpha = NULL, ncores = ncores) + } + # If precipitation is negative: + if (correct_negative) { + cal_hcst_tr[cal_hcst_tr < 0] <- 0 + cal_hcst_ev[cal_hcst_ev < 0] <- 0 + } + lims_cal_hcst_tr <- Apply(cal_hcst_tr, target_dims = c('syear', 'ensemble'), + fun = function(x, prob_lims) { + as.array(quantile(as.vector(x), prob_lims, + na.rm = TRUE))}, + prob_lims = unlist(categories), + ncores = ncores)$output1 + names(dim(lims_cal_hcst_tr))[1] <- 'bin' + lims_obs_tr <- Apply(obs_tr, target_dims = c('syear'),#, 'ensemble'), + fun = function(x, prob_lims) { + as.array(quantile(as.vector(x), prob_lims, + na.rm = TRUE))}, + prob_lims = unlist(categories), + ncores = ncores)$output1 + names(dim(lims_obs_tr))[1] <- 'bin' + if (tolower(recipe$Analysis$Horizon) == 'subseasonal') { + cal_hcst_tr <- SplitDim(cal_hcst_tr, split_dim = 'syear', + new_dim_name = 'sday', + indices = rep(1:dim(data$hcst$data)['sday'], + length(cross$train.dexes))) + obs_tr <- SplitDim(obs_tr, split_dim = 'syear', new_dim_name = 'sday', + indices = rep(1:dim(data$hcst$data)['sday'], + length(cross$train.dexes))) + } + # check order before storing: + # Add back 'sday' dim if needed or select central 'sday' + if ('sday' %in% names(dim(cal_hcst_ev))) { + if (dim(cal_hcst_ev)['sday'] != 1) { + # To fail because it should not appear + } + } else { + cal_hcst_ev <- InsertDim(cal_hcst_ev, 'sday', len = 1, pos = 1) + } + # remove syear dim from ev objects + if ('syear' %in% names(dim(cal_hcst_ev))) { + if (dim(cal_hcst_ev)['syear'] == 1) { + cal_hcst_ev <- Subset(cal_hcst_ev, along = 'syear', indices = 1, + drop = 'selected') + } + } + # ---- + if ('sday' %in% names(dim(cal_hcst_tr))) { + if (dim(cal_hcst_tr)['sday'] > 1) { + # 'sday' dim to select the central day + cal_hcst_tr <- Subset(cal_hcst_tr, along = 'sday', + indices = central_day) + } + } else { + cal_hcst_tr <- InsertDim(cal_hcst_tr, 'sday', len = 1, pos = 1) + } + cal_hcst_tr <- Reorder(cal_hcst_tr, + names(dim(data$hcst$data))) + # ---- + if ('sday' %in% names(dim(obs_tr))) { + if (dim(obs_tr)['sday'] > 1) { + # 'sday' dim to select the central day + obs_tr <- Subset(obs_tr, along = 'sday', + indices = central_day) + } + } else { + obs_tr <- InsertDim(obs_tr, 'sday', len = 1, pos = 1) + } + obs_tr <- Reorder(obs_tr, names(dim(data$obs$data))) + # ---- + if ('sday' %in% names(dim(lims_obs_tr))) { + if (dim(lims_obs_tr)['sday'] > 1) { + # 'sday' dim to select the central day + lims_obs_tr <- Subset(lims_obs_tr, along = 'sday', + indices = central_day) + } + } else { + lims_obs_tr <- InsertDim(lims_obs_tr, 'sday', len = 1, pos = 1) + } + # lims cannot have ensemble dim: + if ('ensemble' %in% names(dim(lims_obs_tr))) { + if (dim(lims_obs_tr)['ensemble'] == 1) { + lims_obs_tr <- Subset(lims_obs_tr, along = 'ensemble', + indices = 1, drop = 'selected') + } else { + # error: + stop("Why it has memb_dim") + } + } + # ---- + if ('sday' %in% names(dim(lims_cal_hcst_tr))) { + if (dim(lims_cal_hcst_tr)['sday'] > 1) { + # 'sday' dim to select the central day + lims_cal_hcst_tr <- Subset(lims_cal_hcst_tr, along = 'sday', + indices = central_day) + } + } else { + lims_cal_hcst_tr <- InsertDim(lims_cal_hcst_tr, 'sday', len = 1, pos = 1) + } + # lims cannot have ensemble dim: + if ('ensemble' %in% names(dim(lims_cal_hcst_tr))) { + if (dim(lims_cal_hcst_tr)['ensemble'] == 1) { + lims_cal_hcst_tr <- Subset(lims_cal_hcst_tr, along = 'ensemble', + indices = 1, drop = 'selected') + } else { + # error: + stop("Why it has memb_dim") + } + } + return(list(cal_hcst_ev = cal_hcst_ev, + cal_hcst_tr = cal_hcst_tr, + obs_tr = obs_tr, + lims_cal_hcst_tr = lims_cal_hcst_tr, + lims_obs_tr = lims_obs_tr)) + } + indices <- array(1:length(cross), dim = c(syear = length(cross))) + res <- Apply(data = list(indices = indices, + exp = data$hcst$data, + obs = data$obs$data), + target_dims = list(indices = NULL, + obs = names(dim(data$obs$data)), + exp = names(dim(data$hcst$data))), + fun = .cross_calibration, + ncores = 1, ## TODO: Explore options for core distribution + cross = cross, + categories = categories, + recipe = recipe) + + info(recipe$Run$logger, + "#### Calibration Cross-validation loop ended #####") + # Rearrange limits into a list: + lims_tmp_hcst <- list() + lims_tmp_obs <- list() + last_index <- 0 + for (ps in 1:length(categories)) { + indices <- last_index + 1:length(categories[[ps]]) + lims_category_hcst <- Subset(res$lims_cal_hcst_tr, along = "bin", + indices = list(indices), + drop = FALSE) + lims_tmp_hcst[[ps]] <- lims_category_hcst + lims_category_obs <- Subset(res$lims_obs_tr, along = "bin", + indices = list(indices), + drop = FALSE) + lims_tmp_obs[[ps]] <- lims_category_obs + last_index <- indices[length(indices)] + } + res$lims_cal_hcst_tr <- lims_tmp_hcst + res$lims_obs_tr <- lims_tmp_obs + rm(lims_tmp_hcst, lims_tmp_obs) + gc() + + + recipe$Run$output_dir <- paste0(recipe$Run$output_dir, + "/", "outputs/Calibration/") + + # Forecast calibration: + if (!is.null(data$fcst)) { + if (tolower(recipe$Analysis$Horizon) %in% c('subseasonal')) { + # merge sample dimensions and select central week + hcst <- MergeDims(data$hcst$data, merge_dims = c('sday', 'syear'), + rename_dim = 'syear', na.rm = FALSE) + hcst <- Subset(hcst, along = 'sweek', + indices = (dim(data$hcst$data)['sweek'] + 1) / 2) + obs <- MergeDims(data$obs$data, merge_dims = c('sday', 'syear'), + rename_dim = 'syear', na.rm = FALSE) + obs <- Subset(obs, along = 'sweek', + indices = (dim(data$obs$data)['sweek'] + 1) / 2) + fcst <- Subset(data$fcst$data, along = 'sday', indices = 1, + drop = 'selected') + } + if (cal_method %in% c('PTF', 'DIST', 'RQUANT', 'QUANT', 'SSPLIN')) { + data$fcst$data <- QuantileMapping(exp = data$hcst$data, + obs = data$obs$data, + exp_cor = data$fcst$data, + method = cal_method, + memb_dim = 'ensemble', + sdate_dim = 'syear', + window_dim = NULL, + na.rm = na.rm, qstep = 0.1, + wet.day = FALSE, + ncores = ncores) + hcst_cal <- QuantileMapping(exp = data$hcst$data, + obs = data$obs$data, + method = cal_method, + memb_dim = 'ensemble', + sdate_dim = 'syear', + window_dim = NULL, + wet.day = FALSE, + na.rm = na.rm, qstep = 0.1, + ncores = ncores) + } else { + data$fcst$data <- Calibration(exp = data$hcst$data, + obs = data$obs$data, + exp_cor = data$fcst$data, + cal.method = cal_method, + multi.model = FALSE, + na.fill = TRUE, na.rm = na.rm, + apply_to = NULL, + alpha = NULL, memb_dim = 'ensemble', + sdate_dim = 'syear', + dat_dim = NULL, ncores = ncores) + hcst_cal <- Calibration(exp = data$hcst$data, + obs = data$obs$data, + cal.method = cal_method, + eval.method = 'in-sample', + multi.model = FALSE, + na.fill = TRUE, na.rm = na.rm, + apply_to = NULL, + alpha = NULL, memb_dim = 'ensemble', + sdate_dim = 'syear', + dat_dim = NULL, ncores = ncores) + } + if (!('sday' %in% names(dim(data$fcst$data)))) { + data$fcst$data <- InsertDim(data$fcst$data, + len = 1, pos = 3, name = 'sday') + } + if (correct_negative) { + hcst_cal[hcst_cal < 0] <- 0 + data$fcst$data[data$fcst$data < 0] <- 0 + } + # Terciles limits using the whole hindcast period: + lims_fcst <- Apply(hcst_cal, target_dims = c('syear', 'ensemble'), + fun = function(x, prob_lims) { + lapply(prob_lims, function(ps) { + as.array(quantile(as.vector(x), ps, + na.rm = TRUE))})}, + output_dims = lapply(categories, function(x) {'bin'}), + prob_lims = categories, + ncores = ncores) + lims <- Apply(obs, target_dims = c('syear', 'ensemble'), + fun = function(x, prob_lims) { + lapply(prob_lims, function(ps) { + as.array(quantile(as.vector(x), ps, + na.rm = TRUE))})}, + output_dims = lapply(categories, function(x) {'bin'}), + prob_lims = categories, + ncores = ncores) + } + tmp_lims2 <- list() + + # SAVING 'lims' which are the observed category limits: + # TODO saving: + recipe$Run$output_dir <- paste0(recipe$Run$output_dir, + "/outputs/Calibration/") + + # Make categories rounded number to use as names: + categories <- recipe$Analysis$Workflow$Probabilities$percentiles + categories <- lapply(categories, function (x) { + sapply(x, function(y) { + round(eval(parse(text = y)),2)})}) + + # Compute Probabilities + hcst_probs_ev <- list() + obs_probs_ev <- list() + fcst_probs <- list() + for (ps in 1:length(categories)) { + hcst_probs_ev[[ps]] <- GetProbs(res$cal_hcst_ev, time_dim = 'syear', + prob_thresholds = NULL, + bin_dim_abs = 'bin', + indices_for_quantiles = NULL, + memb_dim = 'ensemble', + abs_thresholds = res$lims_cal_hcst_tr[[ps]], + ncores = ncores) + central_day <- (dim(data$obs$data)['sday'] + 1)/2 + tmp <- Subset(data$obs$data, along = 'sday', indices = central_day) + obs_probs_ev[[ps]] <- GetProbs(tmp, time_dim = 'syear', + prob_thresholds = NULL, + bin_dim_abs = 'bin', + indices_for_quantiles = NULL, + memb_dim = 'ensemble', + abs_thresholds = res$lims_obs_tr[[ps]], + ncores = ncores) + if (!is.null(data$fcst)) { + fcst_probs[[ps]] <- GetProbs(data$fcst$data, time_dim = 'syear', + prob_thresholds = NULL, + bin_dim_abs = 'bin', + indices_for_quantiles = NULL, + memb_dim = 'ensemble', + abs_thresholds = lims_fcst[[ps]], + ncores = ncores) + } + } + if (recipe$Analysis$Workflow$Probabilities$save == 'all') { + for (ps in 1:length(categories)) { + ## TODO: use .drop_dims() + tmp_lims3 <- .drop_dims(lims[[ps]]) + for (l in 1:dim(lims[[ps]])['bin']) { + tmp_lims <- ClimProjDiags::Subset(x = tmp_lims3, + along = "bin", + indices = l, + drop = "selected") + tmp_lims2 <- append(tmp_lims2, list(tmp_lims)) + names(tmp_lims2)[length(tmp_lims2)] <- + paste0("p_", as.character(categories[[ps]][l]*100)) + } + } + info(recipe$Run$logger, + "SAVING OBSERVED CATEGORY LIMITS") + save_percentiles(recipe = recipe, percentiles = tmp_lims2, + data_cube = data$obs, + agg = "global", outdir = NULL) + #} + } + + # Convert to s2dv_cubes the resulting calibrated + hcst <- data$hcst + hcst$data <- res$cal_hcst_ev + + if (tolower(recipe$Analysis$Horizon) == 'subseasonal') { + # Keep only de central sday position from the dimension + if (dim(hcst$data)['sday'] > 1) { + hcst <- CST_Subset(hcst, along = 'sday', + indices = (dim(hcst$data)['sday'] + 1) / 2) + } + if (dim(data$obs$data)['sday'] > 1) { + data$obs <- CST_Subset(data$obs, along = 'sday', + indices = (dim(data$obs$data)['sday'] + 1) / 2) + } + if (dim(data$hcst$data)['sday'] > 1) { + data$hcst <- CST_Subset(data$hcst, along = 'sday', + indices = (dim(data$hcst$data)['sday'] + 1) / 2) + } + } + info(recipe$Run$logger, + "#### Calibrated and Probabilities Done #####") + if (recipe$Analysis$Workflow$Calibration$save != FALSE) { + info(recipe$Run$logger, "##### START SAVING CALIBRATED #####") + # Save forecast + if ((recipe$Analysis$Workflow$Calibration$save %in% + c('all', 'exp_only', 'fcst_only')) && !is.null(data$fcst)) { + save_forecast(recipe = recipe, data_cube = data$fcst, type = 'fcst') + } + } + # Save probability bins + probs_hcst <- list() + probs_fcst <- list() + probs_obs <- list() + all_names <- NULL + + for (ps in 1:length(categories)) { + for (perc in 1:(length(categories[[ps]]) + 1)) { + if (perc == 1) { + name_elem <- paste0("below_", categories[[ps]][perc]*100) + } else if (perc == length(categories[[ps]]) + 1) { + name_elem <- paste0("above_", categories[[ps]][perc-1]*100) + } else { + name_elem <- paste0("from_", categories[[ps]][perc-1]*100, + "_to_", categories[[ps]][perc]*100) + } + probs_hcst <- append(probs_hcst, + list(Subset(hcst_probs_ev[[ps]], + along = 'bin', + indices = perc, + drop = 'selected'))) + probs_obs <- append(probs_obs, + list(Subset(obs_probs_ev[[ps]], + along = 'bin', + indices = perc, + drop = 'selected'))) + if (!is.null(data$fcst)) { + probs_fcst <- append(probs_fcst, + list(Subset(fcst_probs[[ps]], + along = 'bin', + indices = perc, + drop = 'selected'))) + } + all_names <- c(all_names, name_elem) + } + } + ## TODO: Apply .drop_dims() directly to hcst_probs_ev etc; and move + ## reorganizing and renaming to save_probabilities() + names(probs_hcst) <- all_names + probs_hcst <- lapply(probs_hcst, .drop_dims) + names(probs_obs) <- all_names + probs_obs <- lapply(probs_obs, .drop_dims) + if (!is.null(data$fcst)) { + names(probs_fcst) <- all_names + probs_fcst <- lapply(probs_fcst, .drop_dims) + } + if (recipe$Analysis$Workflow$Probabilities$save %in% + c('all', 'bins_only')) { + save_probabilities(recipe = recipe, probs = probs_fcst, + data_cube = data$fcst, agg = agg, + type = "fcst") + } + + return(list(hcst = hcst, obs = data$obs, fcst = data$fcst, + hcst.full_val = data$hcst, obs.full_val = data$obs, + probs = list(hcst_ev = hcst_probs_ev, + obs_ev = obs_probs_ev, + fcst = fcst_probs), + ref_obs_tr = res$obs_tr)) + +} + +# Outputs: +# hcst is the calibrated hcst for eval indices +# obs are the original obs +# fcst is the calibrated forecast +# .full_val are the original data to calculate bias (although it may be 0 in calibrated data) +# hcst_ev and obs_ev probabilities are needed for the RPS +# fcst_probs for an operaitonal forecast +#ref_obs_tr is needed for the CRPSS +# observerd category limits for the operational forecast + diff --git a/modules/Crossval/R/tmp/MergeDims.R b/modules/Crossval/R/tmp/MergeDims.R new file mode 100644 index 00000000..b0952e5b --- /dev/null +++ b/modules/Crossval/R/tmp/MergeDims.R @@ -0,0 +1,163 @@ +#'Function to Merge Dimensions +#' +#'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} +#' +#'@description This function merges two dimensions of the array \code{data} in a +#''s2dv_cube' object into one. The user can select the dimensions to merge and +#'provide the final name of the dimension. The user can select to remove NA +#'values or keep them. +#' +#'@param data An 's2dv_cube' object +#'@param merge_dims A character vector indicating the names of the dimensions to +#' merge. +#'@param rename_dim a character string indicating the name of the output +#' dimension. If left at NULL, the first dimension name provided in parameter +#' \code{merge_dims} will be used. +#'@param na.rm A logical indicating if the NA values should be removed or not. +#' +#'@examples +#'data <- 1 : c(2 * 3 * 4 * 5 * 6 * 7) +#'dim(data) <- c(time = 7, lat = 2, lon = 3, monthly = 4, member = 6, +#' dataset = 5, var = 1) +#'data[2,,,,,,] <- NA +#'data[c(3,27)] <- NA +#'data <- list(data = data) +#'class(data) <- 's2dv_cube' +#'new_data <- CST_MergeDims(data, merge_dims = c('time', 'monthly')) +#'new_data <- CST_MergeDims(data, merge_dims = c('lon', 'lat'), rename_dim = 'grid') +#'new_data <- CST_MergeDims(data, merge_dims = c('time', 'monthly'), na.rm = TRUE) +#'@export +CST_MergeDims <- function(data, merge_dims = c('ftime', 'monthly'), + rename_dim = NULL, na.rm = FALSE) { + # Check 's2dv_cube' + if (!inherits(data, 's2dv_cube')) { + stop("Parameter 'data' must be of the class 's2dv_cube'.") + } + if (is.null(rename_dim)) { + rename_dim <- merge_dims[1] + } + # data + data$data <- MergeDims(data$data, merge_dims = merge_dims, + rename_dim = rename_dim, na.rm = na.rm) + # dims + data$dims <- dim(data$data) + + # rename_dim + if (length(rename_dim) > 1) { + rename_dim <- as.character(rename_dim[1]) + } + # coords + data$coords[merge_dims] <- NULL + data$coords[[rename_dim]] <- 1:dim(data$data)[rename_dim] + attr(data$coords[[rename_dim]], 'indices') <- TRUE + + # attrs + if (all(merge_dims %in% names(dim(data$attrs$Dates)))) { + original_timezone <- attr(data$attrs$Dates[1], "tzone") + data$attrs$Dates <- MergeDims(data$attrs$Dates, merge_dims = merge_dims, + rename_dim = rename_dim, na.rm = na.rm) + # Transform dates back to POSIXct + data$attrs$Dates <- as.POSIXct(data$attrs$Dates, + origin = "1970-01-01", + tz = original_timezone) + + } else if (any(merge_dims %in% names(dim(data$attrs$Dates)))) { + warning("The dimensions of 'Dates' array will be different from ", + "the temporal dimensions in 'data'. Parameter 'merge_dims' ", + "only includes one temporal dimension of 'Dates'.") + } + return(data) +} +#'Function to Split Dimension +#' +#'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} +#' +#'@description This function merges two dimensions of an array into one. The +#'user can select the dimensions to merge and provide the final name of the +#'dimension. The user can select to remove NA values or keep them. +#' +#'@param data An n-dimensional array with named dimensions +#'@param merge_dims A character vector indicating the names of the dimensions to +#' merge. +#'@param rename_dim A character string indicating the name of the output +#' dimension. If left at NULL, the first dimension name provided in parameter +#' \code{merge_dims} will be used. +#'@param na.rm A logical indicating if the NA values should be removed or not. +#' +#'@examples +#'data <- 1 : 20 +#'dim(data) <- c(time = 10, lat = 2) +#'new_data <- MergeDims(data, merge_dims = c('time', 'lat')) +#'@import abind +#'@importFrom ClimProjDiags Subset +#'@export +MergeDims <- function(data, merge_dims = c('time', 'monthly'), + rename_dim = NULL, na.rm = FALSE) { + # check data + if (is.null(data)) { + stop("Parameter 'data' cannot be NULL.") + } + if (is.null(dim(data))) { + stop("Parameter 'data' must have dimensions.") + } + if (is.null(names(dim(data)))) { + stop("Parameter 'data' must have dimension names.") + } + dims <- dim(data) + # check merge_dims + if (is.null(merge_dims)) { + stop("Parameter 'merge_dims' cannot be NULL.") + } + if (!is.character(merge_dims)) { + stop("Parameter 'merge_dims' must be a character vector ", + "indicating the names of the dimensions to be merged.") + } + if (length(merge_dims) > 2) { + warning("Only two dimensions can be merge, only the first two ", + "dimension will be used. To merge further dimensions ", + "consider to use this function multiple times.") + merge_dims <- merge_dims[1 : 2] + } else if (length(merge_dims) < 2) { + stop("Parameter 'merge_dims' must be of length two.") + } + if (is.null(rename_dim)) { + rename_dim <- merge_dims[1] + } + if (length(rename_dim) > 1) { + warning("Parameter 'rename_dim' has length greater than 1 ", + "and only the first element will be used.") + rename_dim <- as.character(rename_dim[1]) + } + if (!any(names(dims) %in% merge_dims)) { + stop("Parameter 'merge_dims' must match with dimension ", + "names in parameter 'data'.") + } + pos1 <- which(names(dims) == merge_dims[1]) + pos2 <- which(names(dims) == merge_dims[2]) + if (length(pos1) == 0 | length(pos2) == 0) { + stop("Parameter 'merge_dims' must match with dimension ", + "names in parameter 'data'.") + } + if (pos1 > pos2) { + pos1 <- pos1 - 1 + } + data <- lapply(1:dims[pos2], function(x) {Subset(data, along = pos2, + indices = x, drop = 'selected')}) + data <- abind(data, along = pos1) + names(dim(data)) <- names(dims)[-pos2] + if (!is.null(rename_dim)) { + names(dim(data))[pos1] <- rename_dim + } + if (na.rm) { + nas <- which(is.na(Subset(data, along = -pos1, indices = 1))) + if (length(nas) != 0) { + nas <- unlist(lapply(nas, function(x) { + if(all(is.na(Subset(data, along = pos1, + indices = x)))) { + return(x)}})) + data <- Subset(data, along = pos1, indices = -nas) + } + } + return(data) +} + -- GitLab From 40be5516219ecc29a1033bfaf542b1b98ce83071 Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 11 Jun 2025 19:54:52 +0200 Subject: [PATCH 2/6] fix dimensions --- modules/Crossval/Crossval_calibration.R | 79 ++++++------ modules/Crossval/Crossval_metrics.R | 1 + modules/Crossval/R/tmp/MergeDims.R | 163 ------------------------ 3 files changed, 41 insertions(+), 202 deletions(-) delete mode 100644 modules/Crossval/R/tmp/MergeDims.R diff --git a/modules/Crossval/Crossval_calibration.R b/modules/Crossval/Crossval_calibration.R index 26e7e31c..d4c86ac1 100644 --- a/modules/Crossval/Crossval_calibration.R +++ b/modules/Crossval/Crossval_calibration.R @@ -1,7 +1,7 @@ # take the output of Flor/s2s/subseasonal_loading.R source("modules/Crossval/R/tmp/GetProbs.R") -source("modules/Crossval/R/tmp/MergeDims.R") +source("modules/Crossval/R/tmp/CST_MergeDims.R") Crossval_calibration <- function(recipe, data, correct_negative = FALSE) { cal_method <- recipe$Analysis$Workflow$Calibration$method @@ -37,8 +37,7 @@ Crossval_calibration <- function(recipe, data, correct_negative = FALSE) { # lims_cal_hcst_tr: hcst percentiles (per category) # lims_cal_obs_tr: obs percentiles (per category) # cal_hcst_ev: hcst cross-validated calibrated - # cal_hcst_tr: hc - # cal_obs_tr: obs training + # obs_tr: obs training source("modules/Crossval/R/tmp/MergeDims.R") ## training indices cross <- cross[[indices]] @@ -144,18 +143,6 @@ source("modules/Crossval/R/tmp/MergeDims.R") } } # ---- - if ('sday' %in% names(dim(cal_hcst_tr))) { - if (dim(cal_hcst_tr)['sday'] > 1) { - # 'sday' dim to select the central day - cal_hcst_tr <- Subset(cal_hcst_tr, along = 'sday', - indices = central_day) - } - } else { - cal_hcst_tr <- InsertDim(cal_hcst_tr, 'sday', len = 1, pos = 1) - } - cal_hcst_tr <- Reorder(cal_hcst_tr, - names(dim(data$hcst$data))) - # ---- if ('sday' %in% names(dim(obs_tr))) { if (dim(obs_tr)['sday'] > 1) { # 'sday' dim to select the central day @@ -165,7 +152,18 @@ source("modules/Crossval/R/tmp/MergeDims.R") } else { obs_tr <- InsertDim(obs_tr, 'sday', len = 1, pos = 1) } - obs_tr <- Reorder(obs_tr, names(dim(data$obs$data))) + # rename dimension to be ensemble + if ('ensemble' %in% names(dim(obs_tr))) { + if (dim(obs_tr)['ensemble'] == 1) { + obs_tr <- Subset(obs_tr, along = 'ensemble', indices = 1, + drop = 'selected') + } else { + stop("why it is longer than 1?") + } + } + if ('syear' %in% names(dim(obs_tr))) { # should be true + names(dim(obs_tr))[which(names(dim(obs_tr)) == 'syear')] <- 'ensemble' + } # ---- if ('sday' %in% names(dim(lims_obs_tr))) { if (dim(lims_obs_tr)['sday'] > 1) { @@ -207,7 +205,6 @@ source("modules/Crossval/R/tmp/MergeDims.R") } } return(list(cal_hcst_ev = cal_hcst_ev, - cal_hcst_tr = cal_hcst_tr, obs_tr = obs_tr, lims_cal_hcst_tr = lims_cal_hcst_tr, lims_obs_tr = lims_obs_tr)) @@ -224,7 +221,6 @@ source("modules/Crossval/R/tmp/MergeDims.R") cross = cross, categories = categories, recipe = recipe) - info(recipe$Run$logger, "#### Calibration Cross-validation loop ended #####") # Rearrange limits into a list: @@ -264,29 +260,34 @@ source("modules/Crossval/R/tmp/MergeDims.R") rename_dim = 'syear', na.rm = FALSE) obs <- Subset(obs, along = 'sweek', indices = (dim(data$obs$data)['sweek'] + 1) / 2) - fcst <- Subset(data$fcst$data, along = 'sday', indices = 1, + fcst <- Subset(data$fcst$data, along = 'sday', indices = 1, drop = 'selected') - } + } else { + # we need to keep original objects (to later compute bias) + hcst <- data$hcst$data + obs <- data$obs$data + fcst <- data$fcst$data + } if (cal_method %in% c('PTF', 'DIST', 'RQUANT', 'QUANT', 'SSPLIN')) { - data$fcst$data <- QuantileMapping(exp = data$hcst$data, - obs = data$obs$data, - exp_cor = data$fcst$data, - method = cal_method, - memb_dim = 'ensemble', - sdate_dim = 'syear', - window_dim = NULL, - na.rm = na.rm, qstep = 0.1, - wet.day = FALSE, - ncores = ncores) - hcst_cal <- QuantileMapping(exp = data$hcst$data, - obs = data$obs$data, - method = cal_method, - memb_dim = 'ensemble', - sdate_dim = 'syear', - window_dim = NULL, - wet.day = FALSE, - na.rm = na.rm, qstep = 0.1, - ncores = ncores) + data$fcst$data <- QuantileMapping(exp = hcst, + obs = obs, + exp_cor = fcst, + method = cal_method, + memb_dim = 'ensemble', + sdate_dim = 'syear', + window_dim = NULL, + na.rm = na.rm, qstep = 0.1, + wet.day = FALSE, + ncores = ncores) + hcst_cal <- QuantileMapping(exp = hcst, + obs = obs, + method = cal_method, + memb_dim = 'ensemble', + sdate_dim = 'syear', + window_dim = NULL, + wet.day = FALSE, + na.rm = na.rm, qstep = 0.1, + ncores = ncores) } else { data$fcst$data <- Calibration(exp = data$hcst$data, obs = data$obs$data, diff --git a/modules/Crossval/Crossval_metrics.R b/modules/Crossval/Crossval_metrics.R index af1a4425..5a3a2239 100644 --- a/modules/Crossval/Crossval_metrics.R +++ b/modules/Crossval/Crossval_metrics.R @@ -128,6 +128,7 @@ Crossval_metrics <- function(recipe, data_crossval, data_crossval$obs <- CST_MergeDims(data_crossval$obs, merge_dims = c('sweek', 'syear'), rename_dim = 'syear', na.rm = FALSE) +browser() data_crossval$ref_obs_tr <- MergeDims(data_crossval$ref_obs_tr, merge_dims = c('sweek', 'syear'), rename_dim = 'syear', na.rm = FALSE) diff --git a/modules/Crossval/R/tmp/MergeDims.R b/modules/Crossval/R/tmp/MergeDims.R deleted file mode 100644 index b0952e5b..00000000 --- a/modules/Crossval/R/tmp/MergeDims.R +++ /dev/null @@ -1,163 +0,0 @@ -#'Function to Merge Dimensions -#' -#'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} -#' -#'@description This function merges two dimensions of the array \code{data} in a -#''s2dv_cube' object into one. The user can select the dimensions to merge and -#'provide the final name of the dimension. The user can select to remove NA -#'values or keep them. -#' -#'@param data An 's2dv_cube' object -#'@param merge_dims A character vector indicating the names of the dimensions to -#' merge. -#'@param rename_dim a character string indicating the name of the output -#' dimension. If left at NULL, the first dimension name provided in parameter -#' \code{merge_dims} will be used. -#'@param na.rm A logical indicating if the NA values should be removed or not. -#' -#'@examples -#'data <- 1 : c(2 * 3 * 4 * 5 * 6 * 7) -#'dim(data) <- c(time = 7, lat = 2, lon = 3, monthly = 4, member = 6, -#' dataset = 5, var = 1) -#'data[2,,,,,,] <- NA -#'data[c(3,27)] <- NA -#'data <- list(data = data) -#'class(data) <- 's2dv_cube' -#'new_data <- CST_MergeDims(data, merge_dims = c('time', 'monthly')) -#'new_data <- CST_MergeDims(data, merge_dims = c('lon', 'lat'), rename_dim = 'grid') -#'new_data <- CST_MergeDims(data, merge_dims = c('time', 'monthly'), na.rm = TRUE) -#'@export -CST_MergeDims <- function(data, merge_dims = c('ftime', 'monthly'), - rename_dim = NULL, na.rm = FALSE) { - # Check 's2dv_cube' - if (!inherits(data, 's2dv_cube')) { - stop("Parameter 'data' must be of the class 's2dv_cube'.") - } - if (is.null(rename_dim)) { - rename_dim <- merge_dims[1] - } - # data - data$data <- MergeDims(data$data, merge_dims = merge_dims, - rename_dim = rename_dim, na.rm = na.rm) - # dims - data$dims <- dim(data$data) - - # rename_dim - if (length(rename_dim) > 1) { - rename_dim <- as.character(rename_dim[1]) - } - # coords - data$coords[merge_dims] <- NULL - data$coords[[rename_dim]] <- 1:dim(data$data)[rename_dim] - attr(data$coords[[rename_dim]], 'indices') <- TRUE - - # attrs - if (all(merge_dims %in% names(dim(data$attrs$Dates)))) { - original_timezone <- attr(data$attrs$Dates[1], "tzone") - data$attrs$Dates <- MergeDims(data$attrs$Dates, merge_dims = merge_dims, - rename_dim = rename_dim, na.rm = na.rm) - # Transform dates back to POSIXct - data$attrs$Dates <- as.POSIXct(data$attrs$Dates, - origin = "1970-01-01", - tz = original_timezone) - - } else if (any(merge_dims %in% names(dim(data$attrs$Dates)))) { - warning("The dimensions of 'Dates' array will be different from ", - "the temporal dimensions in 'data'. Parameter 'merge_dims' ", - "only includes one temporal dimension of 'Dates'.") - } - return(data) -} -#'Function to Split Dimension -#' -#'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} -#' -#'@description This function merges two dimensions of an array into one. The -#'user can select the dimensions to merge and provide the final name of the -#'dimension. The user can select to remove NA values or keep them. -#' -#'@param data An n-dimensional array with named dimensions -#'@param merge_dims A character vector indicating the names of the dimensions to -#' merge. -#'@param rename_dim A character string indicating the name of the output -#' dimension. If left at NULL, the first dimension name provided in parameter -#' \code{merge_dims} will be used. -#'@param na.rm A logical indicating if the NA values should be removed or not. -#' -#'@examples -#'data <- 1 : 20 -#'dim(data) <- c(time = 10, lat = 2) -#'new_data <- MergeDims(data, merge_dims = c('time', 'lat')) -#'@import abind -#'@importFrom ClimProjDiags Subset -#'@export -MergeDims <- function(data, merge_dims = c('time', 'monthly'), - rename_dim = NULL, na.rm = FALSE) { - # check data - if (is.null(data)) { - stop("Parameter 'data' cannot be NULL.") - } - if (is.null(dim(data))) { - stop("Parameter 'data' must have dimensions.") - } - if (is.null(names(dim(data)))) { - stop("Parameter 'data' must have dimension names.") - } - dims <- dim(data) - # check merge_dims - if (is.null(merge_dims)) { - stop("Parameter 'merge_dims' cannot be NULL.") - } - if (!is.character(merge_dims)) { - stop("Parameter 'merge_dims' must be a character vector ", - "indicating the names of the dimensions to be merged.") - } - if (length(merge_dims) > 2) { - warning("Only two dimensions can be merge, only the first two ", - "dimension will be used. To merge further dimensions ", - "consider to use this function multiple times.") - merge_dims <- merge_dims[1 : 2] - } else if (length(merge_dims) < 2) { - stop("Parameter 'merge_dims' must be of length two.") - } - if (is.null(rename_dim)) { - rename_dim <- merge_dims[1] - } - if (length(rename_dim) > 1) { - warning("Parameter 'rename_dim' has length greater than 1 ", - "and only the first element will be used.") - rename_dim <- as.character(rename_dim[1]) - } - if (!any(names(dims) %in% merge_dims)) { - stop("Parameter 'merge_dims' must match with dimension ", - "names in parameter 'data'.") - } - pos1 <- which(names(dims) == merge_dims[1]) - pos2 <- which(names(dims) == merge_dims[2]) - if (length(pos1) == 0 | length(pos2) == 0) { - stop("Parameter 'merge_dims' must match with dimension ", - "names in parameter 'data'.") - } - if (pos1 > pos2) { - pos1 <- pos1 - 1 - } - data <- lapply(1:dims[pos2], function(x) {Subset(data, along = pos2, - indices = x, drop = 'selected')}) - data <- abind(data, along = pos1) - names(dim(data)) <- names(dims)[-pos2] - if (!is.null(rename_dim)) { - names(dim(data))[pos1] <- rename_dim - } - if (na.rm) { - nas <- which(is.na(Subset(data, along = -pos1, indices = 1))) - if (length(nas) != 0) { - nas <- unlist(lapply(nas, function(x) { - if(all(is.na(Subset(data, along = pos1, - indices = x)))) { - return(x)}})) - data <- Subset(data, along = pos1, indices = -nas) - } - } - return(data) -} - -- GitLab From cd9390c27b7c2fd9b6900687d4a08dc7d98ca676 Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 12 Jun 2025 17:33:10 +0200 Subject: [PATCH 3/6] fix seasonal and remove browser --- modules/Crossval/Crossval_calibration.R | 41 +++++++++++++------------ modules/Crossval/Crossval_metrics.R | 1 - 2 files changed, 21 insertions(+), 21 deletions(-) diff --git a/modules/Crossval/Crossval_calibration.R b/modules/Crossval/Crossval_calibration.R index d4c86ac1..5b0694c3 100644 --- a/modules/Crossval/Crossval_calibration.R +++ b/modules/Crossval/Crossval_calibration.R @@ -38,7 +38,6 @@ Crossval_calibration <- function(recipe, data, correct_negative = FALSE) { # lims_cal_obs_tr: obs percentiles (per category) # cal_hcst_ev: hcst cross-validated calibrated # obs_tr: obs training -source("modules/Crossval/R/tmp/MergeDims.R") ## training indices cross <- cross[[indices]] print(paste("Crossval Index:", cross$eval.dexes)) @@ -249,6 +248,7 @@ source("modules/Crossval/R/tmp/MergeDims.R") "/", "outputs/Calibration/") # Forecast calibration: + lims <- list() if (!is.null(data$fcst)) { if (tolower(recipe$Analysis$Horizon) %in% c('subseasonal')) { # merge sample dimensions and select central week @@ -289,9 +289,9 @@ source("modules/Crossval/R/tmp/MergeDims.R") na.rm = na.rm, qstep = 0.1, ncores = ncores) } else { - data$fcst$data <- Calibration(exp = data$hcst$data, - obs = data$obs$data, - exp_cor = data$fcst$data, + data$fcst$data <- Calibration(exp = hcst, + obs = obs, + exp_cor = fcst, cal.method = cal_method, multi.model = FALSE, na.fill = TRUE, na.rm = na.rm, @@ -299,16 +299,16 @@ source("modules/Crossval/R/tmp/MergeDims.R") alpha = NULL, memb_dim = 'ensemble', sdate_dim = 'syear', dat_dim = NULL, ncores = ncores) - hcst_cal <- Calibration(exp = data$hcst$data, - obs = data$obs$data, - cal.method = cal_method, - eval.method = 'in-sample', - multi.model = FALSE, - na.fill = TRUE, na.rm = na.rm, - apply_to = NULL, - alpha = NULL, memb_dim = 'ensemble', - sdate_dim = 'syear', - dat_dim = NULL, ncores = ncores) + hcst_cal <- Calibration(exp = hcst, + obs = obs, + cal.method = cal_method, + eval.method = 'in-sample', + multi.model = FALSE, + na.fill = TRUE, na.rm = na.rm, + apply_to = NULL, + alpha = NULL, memb_dim = 'ensemble', + sdate_dim = 'syear', + dat_dim = NULL, ncores = ncores) } if (!('sday' %in% names(dim(data$fcst$data)))) { data$fcst$data <- InsertDim(data$fcst$data, @@ -327,11 +327,11 @@ source("modules/Crossval/R/tmp/MergeDims.R") output_dims = lapply(categories, function(x) {'bin'}), prob_lims = categories, ncores = ncores) - lims <- Apply(obs, target_dims = c('syear', 'ensemble'), - fun = function(x, prob_lims) { - lapply(prob_lims, function(ps) { - as.array(quantile(as.vector(x), ps, - na.rm = TRUE))})}, + lims <- Apply(obs, target_dims = c('syear', 'ensemble'), + fun = function(x, prob_lims) { + lapply(prob_lims, function(ps) { + as.array(quantile(as.vector(x), ps, + na.rm = TRUE))})}, output_dims = lapply(categories, function(x) {'bin'}), prob_lims = categories, ncores = ncores) @@ -486,6 +486,7 @@ source("modules/Crossval/R/tmp/MergeDims.R") return(list(hcst = hcst, obs = data$obs, fcst = data$fcst, hcst.full_val = data$hcst, obs.full_val = data$obs, + cat_lims = list(obs = lims), probs = list(hcst_ev = hcst_probs_ev, obs_ev = obs_probs_ev, fcst = fcst_probs), @@ -500,6 +501,6 @@ source("modules/Crossval/R/tmp/MergeDims.R") # .full_val are the original data to calculate bias (although it may be 0 in calibrated data) # hcst_ev and obs_ev probabilities are needed for the RPS # fcst_probs for an operaitonal forecast -#ref_obs_tr is needed for the CRPSS +# ref_obs_tr is needed for the CRPSS # observerd category limits for the operational forecast diff --git a/modules/Crossval/Crossval_metrics.R b/modules/Crossval/Crossval_metrics.R index 5a3a2239..af1a4425 100644 --- a/modules/Crossval/Crossval_metrics.R +++ b/modules/Crossval/Crossval_metrics.R @@ -128,7 +128,6 @@ Crossval_metrics <- function(recipe, data_crossval, data_crossval$obs <- CST_MergeDims(data_crossval$obs, merge_dims = c('sweek', 'syear'), rename_dim = 'syear', na.rm = FALSE) -browser() data_crossval$ref_obs_tr <- MergeDims(data_crossval$ref_obs_tr, merge_dims = c('sweek', 'syear'), rename_dim = 'syear', na.rm = FALSE) -- GitLab From 33dca08f2635ae4d4e43f511176ea4bcc81504b1 Mon Sep 17 00:00:00 2001 From: vagudets Date: Mon, 16 Jun 2025 10:02:44 +0200 Subject: [PATCH 4/6] Lint fix: rewrite path to outputs/Calibration/ --- modules/Crossval/Crossval_calibration.R | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/modules/Crossval/Crossval_calibration.R b/modules/Crossval/Crossval_calibration.R index 5b0694c3..9e6116a4 100644 --- a/modules/Crossval/Crossval_calibration.R +++ b/modules/Crossval/Crossval_calibration.R @@ -2,6 +2,7 @@ source("modules/Crossval/R/tmp/GetProbs.R") source("modules/Crossval/R/tmp/CST_MergeDims.R") + Crossval_calibration <- function(recipe, data, correct_negative = FALSE) { cal_method <- recipe$Analysis$Workflow$Calibration$method @@ -208,6 +209,7 @@ Crossval_calibration <- function(recipe, data, correct_negative = FALSE) { lims_cal_hcst_tr = lims_cal_hcst_tr, lims_obs_tr = lims_obs_tr)) } + indices <- array(1:length(cross), dim = c(syear = length(cross))) res <- Apply(data = list(indices = indices, exp = data$hcst$data, @@ -340,8 +342,8 @@ Crossval_calibration <- function(recipe, data, correct_negative = FALSE) { # SAVING 'lims' which are the observed category limits: # TODO saving: - recipe$Run$output_dir <- paste0(recipe$Run$output_dir, - "/outputs/Calibration/") + recipe$Run$output_dir <- paste0(recipe$Run$output_dir, "/" + "outputs/Calibration/") # Make categories rounded number to use as names: categories <- recipe$Analysis$Workflow$Probabilities$percentiles -- GitLab From 92380d67051dbf6879708d6a96dc583d2f83e5f6 Mon Sep 17 00:00:00 2001 From: vagudets Date: Mon, 16 Jun 2025 11:41:46 +0200 Subject: [PATCH 5/6] Fix pipeline --- modules/Crossval/Crossval_calibration.R | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/modules/Crossval/Crossval_calibration.R b/modules/Crossval/Crossval_calibration.R index 9e6116a4..394a3de5 100644 --- a/modules/Crossval/Crossval_calibration.R +++ b/modules/Crossval/Crossval_calibration.R @@ -245,10 +245,6 @@ Crossval_calibration <- function(recipe, data, correct_negative = FALSE) { rm(lims_tmp_hcst, lims_tmp_obs) gc() - - recipe$Run$output_dir <- paste0(recipe$Run$output_dir, - "/", "outputs/Calibration/") - # Forecast calibration: lims <- list() if (!is.null(data$fcst)) { @@ -342,7 +338,7 @@ Crossval_calibration <- function(recipe, data, correct_negative = FALSE) { # SAVING 'lims' which are the observed category limits: # TODO saving: - recipe$Run$output_dir <- paste0(recipe$Run$output_dir, "/" + recipe$Run$output_dir <- paste0(recipe$Run$output_dir, "/", "outputs/Calibration/") # Make categories rounded number to use as names: -- GitLab From 80ea9755119a16427afa61843faf2ab30486141c Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 16 Jun 2025 15:37:10 +0200 Subject: [PATCH 6/6] use subset and add z500 --- modules/Visualization/R/get_plot_time_labels.R | 3 ++- modules/Visualization/var_long_names.yml | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/modules/Visualization/R/get_plot_time_labels.R b/modules/Visualization/R/get_plot_time_labels.R index 2dfc4202..4cdc23f5 100644 --- a/modules/Visualization/R/get_plot_time_labels.R +++ b/modules/Visualization/R/get_plot_time_labels.R @@ -11,7 +11,8 @@ .get_plot_time_labels <- function(dates, horizon, init_date, start_date, i_syear, time_bounds = NULL) { - years <- lubridate::year(dates[1, 1, which(start_date == i_syear), ]) + years <- lubridate::year(Subset(dates, along = 'syear', + indices = which(start_date == i_syear))) if (is.null(attributes(time_bounds))) { if (horizon == "subseasonal") { time_labels <- dates[1, 1, which(start_date == i_syear), ] diff --git a/modules/Visualization/var_long_names.yml b/modules/Visualization/var_long_names.yml index a6293406..9f670abc 100644 --- a/modules/Visualization/var_long_names.yml +++ b/modules/Visualization/var_long_names.yml @@ -5,4 +5,4 @@ Mycase: prlr: "Total precipitation" sfcWind: "Wind speed" rsds: "Solar radiation" - + z500: "Geopotential height at 500 hPa" -- GitLab