diff --git a/modules/Skill/R/CRPS_clim.R b/modules/Skill/R/CRPS_clim.R new file mode 100644 index 0000000000000000000000000000000000000000..36db4e945fcfc7256ed66c05b98bf9db8c6892e3 --- /dev/null +++ b/modules/Skill/R/CRPS_clim.R @@ -0,0 +1,13 @@ +CRPS_clim <- function(obs, memb_dim ='ensemble'){ + time_dim <- names(dim(obs)) + obs_time_len <- dim(obs)[time_dim] + + ref <- array(data = rep(obs, each = obs_time_len), dim = c(obs_time_len, obs_time_len)) + names(dim(ref)) <- c(time_dim, memb_dim) + # ref: [sdate, memb] + # obs: [sdate] + crps_ref <- s2dv:::.CRPS(exp = ref, obs = obs, time_dim = time_dim, memb_dim = memb_dim, + dat_dim = NULL, Fair = FALSE) + # crps_ref should be [sdate] + return(mean(crps_ref)) +} diff --git a/modules/Skill/R/RPS_clim.R b/modules/Skill/R/RPS_clim.R new file mode 100644 index 0000000000000000000000000000000000000000..e8b6452d14f4bfa2972af651b16f56011cf20b13 --- /dev/null +++ b/modules/Skill/R/RPS_clim.R @@ -0,0 +1,20 @@ +# RPS version for climatology +RPS_clim <- function(obs, indices_for_clim = NULL, prob_thresholds = c(1/3, 2/3)) { + + if (is.null(indices_for_clim)){ + indices_for_clim <- 1:length(obs) + } + + obs_probs <- .GetProbs(data = obs, indices_for_quantiles = indices_for_clim, ## temporarily removed s2dv::: + prob_thresholds = prob_thresholds, weights = NULL, cross.val = T) ## here! + # clim_probs: [bin, sdate] + clim_probs <- c(prob_thresholds[1], diff(prob_thresholds), 1 - prob_thresholds[length(prob_thresholds)]) + clim_probs <- array(clim_probs, dim = dim(obs_probs)) + + # Calculate RPS for each time step + probs_clim_cumsum <- apply(clim_probs, 2, cumsum) + probs_obs_cumsum <- apply(obs_probs, 2, cumsum) + rps_ref <- apply((probs_clim_cumsum - probs_obs_cumsum)^2, 2, sum) + + return(mean(rps_ref)) +} diff --git a/modules/Skill/Skill.R b/modules/Skill/Skill.R index d3f74b1374e4fc8607abe7a73a5cf4db48d17f99..7eb443215af8bdf26947dc5d5be9825435cd3aaf 100644 --- a/modules/Skill/Skill.R +++ b/modules/Skill/Skill.R @@ -15,6 +15,8 @@ source("modules/Skill/R/tmp/RPS.R") source("modules/Skill/R/tmp/RPSS.R") source("modules/Skill/R/tmp/GetProbs.R") source("modules/Skill/R/tmp/RandomWalkTest.R") +source("modules/Skill/R/RPS_clim.R") +source("modules/Skill/R/CRPS_clim.R") ## TODO: Implement this in the future ## Which parameter are required? # if (!("obs" %in% ls()) || is.null(obs)) { @@ -90,6 +92,13 @@ compute_skill_metrics <- function(recipe, data) { } else { na.rm = recipe$Analysis$remove_NAs } + if (is.null(recipe$Analysis$Workflow$Skill$cross_validation)) { + warn(recipe$Run$logger, + "cross_validation parameter not defined, setting it to FALSE.") + cross.val <- FALSE + } else { + cross.val <- recipe$Analysis$Workflow$Skill$cross_validation + } skill_metrics <- list() for (metric in strsplit(metrics, ", | |,")[[1]]) { # Whether the fair version of the metric is to be computed @@ -111,15 +120,22 @@ compute_skill_metrics <- function(recipe, data) { time_dim = time_dim, memb_dim = memb_dim, Fair = Fair, + cross.val = cross.val, ncores = ncores) skill <- .drop_dims(skill) skill_metrics[[ metric ]] <- skill + rps_clim <- Apply(list(data$obs$data), + target_dims = c(time_dim, memb_dim), + RPS_clim)$output1 + rps_clim <- .drop_dims(rps_clim) + skill_metrics[[paste0(metric, "_clim")]] <- rps_clim # Ranked Probability Skill Score and Fair version } else if (metric %in% c('rpss', 'frpss')) { skill <- RPSS(data$hcst$data, data$obs$data, time_dim = time_dim, memb_dim = memb_dim, Fair = Fair, + cross.val = cross.val, ncores = ncores) skill <- lapply(skill, function(x) { .drop_dims(x)}) @@ -131,6 +147,7 @@ compute_skill_metrics <- function(recipe, data) { time_dim = time_dim, memb_dim = memb_dim, prob_thresholds = 0.1, + cross.val = cross.val, Fair = Fair, ncores = ncores) skill <- lapply(skill, function(x) { @@ -143,6 +160,7 @@ compute_skill_metrics <- function(recipe, data) { time_dim = time_dim, memb_dim = memb_dim, prob_thresholds = 0.9, + cross.val = cross.val, Fair = Fair, ncores = ncores) skill <- lapply(skill, function(x) { @@ -158,6 +176,10 @@ compute_skill_metrics <- function(recipe, data) { ncores = ncores) skill <- .drop_dims(skill) skill_metrics[[ metric ]] <- skill + crps_clim <- Apply(list(data$obs$data), target_dims = time_dim, + fun = CRPS_clim, memb_dim = memb_dim)$output1 + crps_clim <- .drop_dims(crps_clim) + skill_metrics[['crps_clim']] <- crps_clim # CRPSS and FCRPSS } else if (metric %in% c('crpss', 'fcrpss')) { skill <- CRPSS(data$hcst$data, data$obs$data, @@ -169,6 +191,18 @@ compute_skill_metrics <- function(recipe, data) { .drop_dims(x)}) skill_metrics[[ metric ]] <- skill$crpss skill_metrics[[ paste0(metric, "_significance") ]] <- skill$sign + } else if (metric == 'rms') { + source("https://earth.bsc.es/gitlab/es/s2dv/-/raw/master/R/RMS.R") + hcst_mean <- Apply(list(data$hcst$data), target_dims = memb_dim, + fun = mean, na.rm = na.rm, ncores = ncores)$output1 + hcst_mean <- InsertDim(hcst_mean, pos = 1, len = 1, name = memb_dim) + skill <- RMS(exp = hcst_mean, + obs = data$obs$data, + time_dim = time_dim, dat_dim = NULL, comp_dim = NULL, + limits = NULL, conf = FALSE, alpha = 0.05, ncores = ncores) + skill <- lapply(skill, function(x) { + .drop_dims(x)}) + skill_metrics[[ metric ]] <- skill$rms # Mean bias (climatology) } else if (metric == 'mean_bias') { ## TODO: Eliminate option to compute from anomalies diff --git a/tests/testthat/test-decadal_monthly_2.R b/tests/testthat/test-decadal_monthly_2.R index a6ca9254464918bd99295501dcb9599117a7b2c6..01fe0440c3f218b8c78389c48ac3b489e3e51793 100644 --- a/tests/testthat/test-decadal_monthly_2.R +++ b/tests/testthat/test-decadal_monthly_2.R @@ -164,7 +164,7 @@ TRUE ) expect_equal( names(skill_metrics), -c("rpss_specs", "enscorr_specs", "frps_specs", "frpss_specs", "bss10_specs", "frps") +c("rpss_specs", "enscorr_specs", "frps_specs", "frpss_specs", "bss10_specs", "frps", "frps_clim") ) expect_equal( class(skill_metrics$rpss_specs),