From c49e5730e9c096121c6904eb15045a4d8711e9c8 Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 6 Nov 2020 19:15:49 +0100 Subject: [PATCH 01/10] RPSS option into CSST_MultiMetric --- DESCRIPTION | 1 + NAMESPACE | 6 + R/CST_Calibration.R | 10 +- R/CST_CategoricalEnsCombination.R | 10 +- R/CST_MultiMetric.R | 165 ++++++++++++++++++-------- man/CST_MultiMetric.Rd | 24 +++- tests/testthat/test-CST_MultiMetric.R | 50 +++++--- 7 files changed, 190 insertions(+), 76 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index d7627660..9020da9e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -62,6 +62,7 @@ Imports: reshape2, ggplot2, RColorBrewer, + easyVerification, graphics, grDevices, stats, diff --git a/NAMESPACE b/NAMESPACE index 983a8926..6996a8b2 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -44,6 +44,7 @@ export(WeatherRegime) export(as.s2dv_cube) export(s2dv_cube) import(abind) +import(easyVerification) import(ggplot2) import(multiApply) import(ncdf4) @@ -86,6 +87,11 @@ importFrom(maps,map) importFrom(plyr,.) importFrom(plyr,dlply) importFrom(reshape2,melt) +importFrom(s2dv,Corr) +importFrom(s2dv,InsertDim) +importFrom(s2dv,MeanDims) +importFrom(s2dv,RMS) +importFrom(s2dv,RMSSS) importFrom(s2dv,Reorder) importFrom(s2dverification,Subset) importFrom(utils,glob2rx) diff --git a/R/CST_Calibration.R b/R/CST_Calibration.R index ca290397..83d61aba 100644 --- a/R/CST_Calibration.R +++ b/R/CST_Calibration.R @@ -251,7 +251,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o .calc.obs.fc.quant <- function(obs, fc){ #function to calculate different quantities of a series of ensemble forecasts and corresponding observations amt.mbr <- dim(fc)[1] - obs.per.ens <- InsertDim(var = obs, posdim = 1, lendim = amt.mbr) + obs.per.ens <- InsertDim(obs, posdim = 1, lendim = amt.mbr) fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE) cor.obs.fc <- cor(fc.ens.av, obs, use = "complete.obs") obs.av <- mean(obs, na.rm = TRUE) @@ -271,7 +271,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o .calc.obs.fc.quant.ext <- function(obs, fc){ #extended function to calculate different quantities of a series of ensemble forecasts and corresponding observations amt.mbr <- dim(fc)[1] - obs.per.ens <- InsertDim(var = obs, posdim = 1, lendim = amt.mbr) + obs.per.ens <- InsertDim(obs, posdim = 1, lendim = amt.mbr) fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE) cor.obs.fc <- cor(fc.ens.av, obs, use = "complete.obs") obs.av <- mean(obs, na.rm = TRUE) @@ -296,7 +296,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE) fc.ens.av.av <- mean(fc.ens.av, na.rm = TRUE) fc.ens.av.sd <- sd(fc.ens.av, na.rm = TRUE) - fc.ens.av.per.ens <- InsertDim(var = fc.ens.av, posdim = 1, lendim = amt.mbr) + fc.ens.av.per.ens <- InsertDim(fc.ens.av, posdim = 1, lendim = amt.mbr) fc.ens.sd <- apply(fc, c(2), sd, na.rm = TRUE) fc.ens.var.av.sqrt <- sqrt(mean(fc.ens.sd^2, na.rm = TRUE)) fc.dev <- fc - fc.ens.av.per.ens @@ -322,10 +322,10 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o .calc.fc.quant.ext <- function(fc){ #extended function to calculate different quantities of a series of ensemble forecasts amt.mbr <- dim(fc)[1] - repmat1.tmp <- InsertDim(var = fc, posdim = 1, lendim = amt.mbr) + repmat1.tmp <- InsertDim(fc, posdim = 1, lendim = amt.mbr) repmat2.tmp <- aperm(repmat1.tmp, c(2, 1, 3)) spr.abs <- apply(abs(repmat1.tmp - repmat2.tmp), c(3), mean, na.rm = TRUE) - spr.abs.per.ens <- InsertDim(var = spr.abs, posdim = 1, lendim = amt.mbr) + spr.abs.per.ens <- InsertDim(spr.abs, posdim = 1, lendim = amt.mbr) return( append(.calc.fc.quant(fc), diff --git a/R/CST_CategoricalEnsCombination.R b/R/CST_CategoricalEnsCombination.R index 5b70d068..f72415e3 100644 --- a/R/CST_CategoricalEnsCombination.R +++ b/R/CST_CategoricalEnsCombination.R @@ -89,7 +89,7 @@ CST_CategoricalEnsCombination <- function(exp, obs, cat.method = "pool", eval.me exp$data <- .CategoricalEnsCombination.wrap(fc = exp$data, obs = obs$data, cat.method = cat.method, eval.method = eval.method, amt.cat = amt.cat, ...) names.dim.tmp[which(names.dim.tmp == "member")] <- "category" names(dim(exp$data)) <- names.dim.tmp - exp$data <- InsertDim(var = exp$data, lendim = 1, posdim = 2) + exp$data <- InsertDim(exp$data, lendim = 1, posdim = 2) names(dim(exp$data))[2] <- "member" exp$Datasets <- c(exp$Datasets, obs$Datasets) exp$source_files <- c(exp$source_files, obs$source_files) @@ -277,8 +277,8 @@ comb.dims <- function(arr.in, dims.to.combine){ ui = constr.mtrx, ci = constr.vec, freq.per.mdl.at.obs = freq.per.mdl.at.obs) init.par <- optim.tmp$par * (1 - abs(rnorm(amt.coeff, 0, 0.01))) - var.cat.fc[ , eval.dexes] <- apply(InsertDim(var = - InsertDim(var = optim.tmp$par, lendim = amt.cat, posdim = 2), + var.cat.fc[ , eval.dexes] <- apply(InsertDim( + InsertDim(optim.tmp$par, lendim = amt.cat, posdim = 2), lendim = amt.sdate.ev, posdim = 3) * freq.per.mdl.ev[ , , , drop = FALSE], c(2,3), sum, na.rm = TRUE) } else if (cat.method == "comb") { @@ -363,7 +363,7 @@ comb.dims <- function(arr.in, dims.to.combine){ .funct.optim.grad <- function(par, freq.per.mdl.at.obs){ amt.model <- dim(freq.per.mdl.at.obs)[1] - return(-apply(freq.per.mdl.at.obs/InsertDim(var = drop(par %*% freq.per.mdl.at.obs), + return(-apply(freq.per.mdl.at.obs/InsertDim(drop(par %*% freq.per.mdl.at.obs), lendim = amt.model, posdim = 1), c(1), mean, na.rm = TRUE)) } @@ -373,7 +373,7 @@ comb.dims <- function(arr.in, dims.to.combine){ amt.mdl <- mdl.feat$amt.mdl mdl.msk.tmp <- mdl.feat$mdl.msk amt.coeff <- amt.mdl + 1 - msk.fc.obs <- (cat.fc == InsertDim(var = cat.obs, posdim = 1, lendim = amt.mbr)) + msk.fc.obs <- (cat.fc == InsertDim(cat.obs, posdim = 1, lendim = amt.mbr)) freq.per.mdl.at.obs <- array(NA, c(amt.coeff, amt.sdate)) for (i.mdl in seq(1, amt.mdl)){ freq.per.mdl.at.obs[i.mdl, ] <- apply(msk.fc.obs[mdl.msk.tmp[i.mdl, ], , drop = FALSE], diff --git a/R/CST_MultiMetric.R b/R/CST_MultiMetric.R index e85c8b68..96e64f40 100644 --- a/R/CST_MultiMetric.R +++ b/R/CST_MultiMetric.R @@ -4,17 +4,22 @@ #'@author Perez-Zanon Nuria, \email{nuria.perez@bsc.es} #'@description This function calculates correlation (Anomaly Correlation Coefficient; ACC), root mean square error (RMS) and the root mean square error skill score (RMSSS) of individual anomaly models and multi-models mean (if desired) with the observations. #' -#'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the anomaly of the seasonal forecast experiment data in the element named \code{$data}. +#'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the anomaly of the seasonal forecast experiments data in the element named \code{$data}. #'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the anomaly of observed data in the element named \code{$data}. #'@param metric a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms' or 'rmsss. #'@param multimodel a logical value indicating whether a Multi-Model Mean should be computed. #' +#'@param time_dim_name name of the temporal dimension where a mean will be applied. It can be NULL, the default value is 'ftime'. +#'@param memb_dim_name name of the member dimension. It can be NULL, the default value is 'member'. +#'@param sdate_dim_name name of the start date dimension or a dimension name identifiying the different forecast. It can be NULL, the default value is 'sdate'. #'@return an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is an array with two datset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest first dimension correspons to the Multi-Model Mean. The third dimension contains the statistics selected. For metric \code{correlation}, the third dimension is of length four and they corresponds to the lower limit of the 95\% confidence interval, the statistics itselfs, the upper limit of the 95\% confidence interval and the 95\% significance level. For metric \code{rms}, the third dimension is length three and they corresponds to the lower limit of the 95\% confidence interval, the RMSE and the upper limit of the 95\% confidence interval. For metric \code{rmsss}, the third dimension is length two and they corresponds to the statistics itselfs and the p-value of the one-sided Fisher test with Ho: RMSSS = 0. #'@seealso \code{\link[s2dverification]{Corr}}, \code{\link[s2dverification]{RMS}}, \code{\link[s2dverification]{RMSSS}} and \code{\link{CST_Load}} #'@references #'Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31.\url{http://link.springer.com/10.1007/s00382-018-4404-z} #' -#'@import s2dverification +#'@importFrom s2dv MeanDims Reorder Corr RMS RMSSS InsertDim +#'@import abind +#'@import easyVerification #'@import stats #'@examples #'library(zeallot) @@ -31,21 +36,31 @@ #'c(ano_exp, ano_obs) %<-% CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE) #'a <- CST_MultiMetric(exp = ano_exp, obs = ano_obs) #'str(a) +#'exp <- lonlat_data$exp +#'obs <- lonlat_data$obs +#'a <- CST_MultiMetric(exp, obs, metric = 'rpss') +#'a <- CST_MultiMetric(exp, obs, metric = 'correlation') +#'a <- CST_MultiMetric(exp, obs, metric = 'rms') +#'a <- CST_MultiMetric(exp, obs, metric = 'rmsss') #'@export -CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE) { +CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE, + time_dim = 'ftime', memb_dim = 'member', + sdate_dim = 'sdate') { if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) { stop("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ", "as output by CSTools::CST_Load.") } + result <- MultiMetric(exp$data, obs$data, metric = metric, multimodel = multimodel, + time_dim = time_dim, memb_dim = memb_dim, sdate_dim = sdate_dim) + exp$data <- result + return(exp) +} - if (dim(obs$data)['member'] != 1) { - stop("The length of the dimension 'member' in the component 'data' ", - "of the parameter 'obs' must be equal to 1.") - } - - if (!is.null(names(dim(exp$data))) & !is.null(names(dim(obs$data)))) { - if (all(names(dim(exp$data)) %in% names(dim(obs$data)))) { - dimnames <- names(dim(exp$data)) +MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE, + time_dim = 'ftime', memb_dim = 'member', sdate_dim = 'sdate') { + if (!is.null(names(dim(exp))) & !is.null(names(dim(obs)))) { + if (all(names(dim(exp)) %in% names(dim(obs)))) { + dimnames <- names(dim(exp)) } else { stop("Dimension names of element 'data' from parameters 'exp'", " and 'obs' should have the same name dimmension.") @@ -54,7 +69,6 @@ CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE) stop("Element 'data' from parameters 'exp' and 'obs'", " should have dimmension names.") } - if (!is.logical(multimodel)) { stop("Parameter 'multimodel' must be a logical value.") } @@ -68,49 +82,100 @@ CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE) warning("Parameter 'multimodel' has length > 1 and only the first ", "element will be used.") } - - # seasonal average of anomalies per model - AvgExp <- MeanListDim(exp$data, narm = T, c(2, 4)) - AvgObs <- MeanListDim(obs$data, narm = T, c(2, 4)) - - # indv model correlation - if (metric == 'correlation') { - corr <- Corr(AvgExp, AvgObs, posloop = 1, poscor = 2) - } else if (metric == 'rms') { - corr <- RMS(AvgExp, AvgObs, posloop = 1, posRMS = 2) - } else if (metric == 'rmsss') { - corr <- RMSSS(AvgExp, AvgObs, posloop = 1, posRMS = 2) + if (is.null(time_dim) | !is.character(time_dim)) { + time_dim <- 'time' + } + if (is.null(memb_dim) | !is.character(memb_dim)) { + memb_dim <- 'memb' + } + if( is.null(sdate_dim) | !is.character(sdate_dim)) { + sdate_dim <- 'sdate' + } + exp_dims <- dim(exp) + obs_dims <- dim(obs) + if (!is.null(names(exp_dims)) & !is.null(names(obs_dims))) { + if (all(names(exp_dims) == names(obs_dims))) { + if (!(time_dim %in% names(exp_dims))) { + warning("Parameter 'time_dim' does not match with a dimension name in 'exp'", + " and 'obs'. A 'time_dim' of length 1 is added.") + dim(exp) <- c(exp_dims, time_dim = 1) + names(dim(exp))[length(dim(exp))] <- time_dim + dim(obs) <- c(obs_dims, time_dim = 1) + names(dim(obs))[length(dim(obs))] <- time_dim + exp_dims <- dim(exp) + obs_dims <- dim(obs) + } + if (!(memb_dim %in% names(exp_dims))) { + warning("Parameter 'memb_dim' does not match with a dimension name in ", + "'exp' and 'obs'. A 'memb_dim' of length 1 is added.") + dim(exp) <- c(exp_dims, memb_dim = 1) + names(dim(exp))[length(dim(exp))] <- memb_dim + dim(obs) <- c(obs_dims, memb_dim = 1) + names(dim(obs))[length(dim(obs))] <- memb_dim + exp_dims <- dim(exp) + obs_dims <- dim(obs) + } + if (!(sdate_dim %in% names(exp_dims))) { + warning("Parameter 'sdate_dim' does not match with a dimension name in ", + "'exp' and 'obs'. A 'sdate_dim' of length 1 is added.") + dim(exp) <- c(exp_dims, sdate_dim = 1) + names(dim(exp))[length(dim(exp))] <- sdate_dim + dim(obs) <- c(obs_dims, sdate_dim = 1) + names(dim(obs))[length(dim(obs))] <- sdate_dim + exp_dims <- dim(exp) + obs_dims <- dim(obs) + } + } else { + stop("Dimension names of element 'data' from parameters 'exp'", + " and 'obs' should be the same and in the same order.") + } } else { - stop("Parameter 'metric' must be a character string indicating ", - "one of the options: 'correlation', 'rms' or 'rmse'.") + stop("Element 'data' from parameters 'exp' and 'obs'", + " should have dimmension names.") } - if (multimodel == TRUE) { - # seasonal avg of anomalies for multi-model - AvgExp_MMM <- MeanListDim(AvgExp, narm = TRUE, 1) - AvgObs_MMM <- MeanListDim(AvgObs, narm = TRUE, 1) - # multi model correlation + if (metric == 'rpss') { + if (multimodel == TRUE) { + warning("A probabilistic metric cannot be use to evaluate a multimodel mean.") + } + AvgExp <- MeanDims(exp, time_dim, na.rm = TRUE) + AvgObs <- MeanDims(obs, time_dim, na.rm = TRUE) + pos_memb <- which(names(dim(AvgExp)) == memb_dim) + dim(AvgObs) <- dim(AvgObs)[-pos_memb] + AvgExp <- Reorder(AvgExp, c(names(dim(AvgExp))[-pos_memb], memb_dim)) + pos_memb <- which(names(dim(AvgExp)) == memb_dim) + pos_sdate <- which(names(dim(AvgExp)) == sdate_dim) + corr <- veriApply('FairRpss', fcst = AvgExp, obs = AvgObs, ensdim = pos_memb, + tdim = pos_sdate, prob = c(1/3, 2/3)) + + } else if (metric %in% c('correlation', 'rms', 'rmsss')) { + AvgExp <- MeanDims(exp, c(memb_dim, time_dim), na.rm = TRUE) + AvgObs <- MeanDims(obs, c(memb_dim, time_dim), na.rm = TRUE) + dataset_dim <- c('data', 'dataset', 'datsets', 'models') + if (any(dataset_dim %in% names(exp_dims))) { + dataset_dim <- dataset_dim[dataset_dim %in% names(dim(AvgExp))] + } else { + warning("Parameter 'exp' and 'obs' does not have a dimension 'dataset'.") + } + if (multimodel == TRUE) { + # seasonal avg of anomalies for multi-model + AvgExp_MMM <- MeanDims(AvgExp, c(dataset_dim), na.rm = TRUE) + pos_dataset <- which(names(dim(AvgExp)) == dataset_dim) + AvgExp_MMM <- InsertDim(AvgExp_MMM, posdim = pos_dataset, lendim = 1, + name = dataset_dim) + AvgExp <- abind(AvgExp_MMM, AvgExp, along = pos_dataset) + names(dim(AvgExp)) <- names(dim(AvgExp_MMM)) + } if (metric == 'correlation') { - corr_MMM <- Corr(var_exp = InsertDim(AvgExp_MMM, 1, 1), - var_obs = InsertDim(AvgObs_MMM, 1, 1), - posloop = 1, poscor = 2) + corr <- s2dv::Corr(AvgExp, AvgObs, memb_dim = dataset_dim, time_dim = sdate_dim) } else if (metric == 'rms') { - corr_MMM <- RMS(var_exp = InsertDim(AvgExp_MMM, 1, 1), - var_obs = InsertDim(AvgObs_MMM, 1, 1), - posloop = 1, posRMS = 2) + corr <- s2dv::RMS(AvgExp, AvgObs, memb_dim = dataset_dim, time_dim = sdate_dim) } else if (metric == 'rmsss') { - corr_MMM <- RMSSS(var_exp = InsertDim(AvgExp_MMM, 1, 1), - var_obs = InsertDim(AvgObs_MMM, 1, 1), - posloop = 1, posRMS = 2) - } - corr <- abind::abind(corr, corr_MMM, along = 1) + corr <- s2dv::RMSSS(AvgExp, AvgObs, memb_dim = dataset_dim, time_dim = sdate_dim) + } + } else { + stop("Parameter 'metric' must be a character string indicating ", + "one of the options: 'correlation', 'rms', 'rmsss' or 'rpss'.") } - names(dim(corr)) <- c(dimnames[1], dimnames[1], 'statistics', dimnames[5 : 6]) - - #exp$data <- ano$ano_exp - #obs$data <- ano$ano_obs - exp$data <- corr - exp$Datasets <- c(exp$Datasets, obs$Datasets) - exp$source_files <- c(exp$source_files, obs$source_files) - return(exp) + return(corr) } diff --git a/man/CST_MultiMetric.Rd b/man/CST_MultiMetric.Rd index 8e3ce593..4526fa0e 100644 --- a/man/CST_MultiMetric.Rd +++ b/man/CST_MultiMetric.Rd @@ -4,16 +4,30 @@ \alias{CST_MultiMetric} \title{Multiple Metrics applied in Multiple Model Anomalies} \usage{ -CST_MultiMetric(exp, obs, metric = "correlation", multimodel = TRUE) +CST_MultiMetric( + exp, + obs, + metric = "correlation", + multimodel = TRUE, + time_dim = "ftime", + memb_dim = "member", + sdate_dim = "sdate" +) } \arguments{ -\item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the anomaly of the seasonal forecast experiment data in the element named \code{$data}.} +\item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the anomaly of the seasonal forecast experiments data in the element named \code{$data}.} \item{obs}{an object of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the anomaly of observed data in the element named \code{$data}.} \item{metric}{a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms' or 'rmsss.} \item{multimodel}{a logical value indicating whether a Multi-Model Mean should be computed.} + +\item{time_dim_name}{name of the temporal dimension where a mean will be applied. It can be NULL, the default value is 'ftime'.} + +\item{memb_dim_name}{name of the member dimension. It can be NULL, the default value is 'member'.} + +\item{sdate_dim_name}{name of the start date dimension or a dimension name identifiying the different forecast. It can be NULL, the default value is 'sdate'.} } \value{ an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is an array with two datset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest first dimension correspons to the Multi-Model Mean. The third dimension contains the statistics selected. For metric \code{correlation}, the third dimension is of length four and they corresponds to the lower limit of the 95\% confidence interval, the statistics itselfs, the upper limit of the 95\% confidence interval and the 95\% significance level. For metric \code{rms}, the third dimension is length three and they corresponds to the lower limit of the 95\% confidence interval, the RMSE and the upper limit of the 95\% confidence interval. For metric \code{rmsss}, the third dimension is length two and they corresponds to the statistics itselfs and the p-value of the one-sided Fisher test with Ho: RMSSS = 0. @@ -36,6 +50,12 @@ attr(obs, 'class') <- 's2dv_cube' c(ano_exp, ano_obs) \%<-\% CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE) a <- CST_MultiMetric(exp = ano_exp, obs = ano_obs) str(a) +exp <- lonlat_data$exp +obs <- lonlat_data$obs +a <- CST_MultiMetric(exp, obs, metric = 'rpss') +a <- CST_MultiMetric(exp, obs, metric = 'correlation') +a <- CST_MultiMetric(exp, obs, metric = 'rms') +a <- CST_MultiMetric(exp, obs, metric = 'rmsss') } \references{ Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31.\url{http://link.springer.com/10.1007/s00382-018-4404-z} diff --git a/tests/testthat/test-CST_MultiMetric.R b/tests/testthat/test-CST_MultiMetric.R index 1a41a83a..acbddfdc 100644 --- a/tests/testthat/test-CST_MultiMetric.R +++ b/tests/testthat/test-CST_MultiMetric.R @@ -11,30 +11,41 @@ test_that("basic use case", { attr(exp, 'class') <- 's2dv_cube' attr(obs, 'class') <- 's2dv_cube' - result <- list(data = array(rep(c(rep(1, 9), rep(0, 3)), 48), - dim = c(dataset = 3, dataset = 1, statistics = 4, - lat = 6, lon = 8)), - lat = lat, lon = lon) + result <- list(data = list(corr = array(rep(1, 3* 48), + #dim = c(nexp = 3, nobs = 1, + dim = c(n_exp = 3, n_obs = 1, + lat = 6, lon = 8)), + p.val = array(rep(0, 3 * 48), dim = c(n_exp = 3, n_obs = 1, + lat = 6, lon = 8)), + conf.lower = array(rep(1, 3* 48), + #dim = c(nexp = 3, nobs = 1, + dim = c(n_exp = 3, n_obs = 1, + lat = 6, lon = 8)), + conf.upper = array(rep(1, 3* 48), + #dim = c(nexp = 3, nobs = 1, + dim = c(n_exp = 3, n_obs = 1, + lat = 6, lon = 8))), + lat = lat, lon = lon) attr(result, 'class') <- 's2dv_cube' expect_equal(CST_MultiMetric(exp = exp, obs = obs), result) exp2 <- exp exp2$data[1, 1, 1, 2, 1, 1] = NA - CST_MultiMetric(exp = exp2, obs = obs) res <- CST_MultiMetric(exp = exp, obs = obs, metric = 'rms') expect_equal(length(res), 3) - expect_equal(dim(res$data), - c(dataset = 3, dataset = 1, statistics = 3, lat = 6, lon = 8)) + expect_equal(dim(res$data$rms), + c(n_exp = 3, n_obs = 1, lat = 6, lon = 8)) res <- CST_MultiMetric(exp = exp, obs = obs, metric = 'rms', multimodel = FALSE) - expect_equal(dim(res$data), - c(dataset = 2, dataset = 1, statistics = 3, lat = 6, lon = 8)) + expect_equal(dim(res$data$rms), + c(n_exp = 2, n_obs = 1, lat = 6, lon = 8)) + expect_equal(length(res$data), 3) res <- CST_MultiMetric(exp = exp, obs = obs, metric = 'rmsss') - expect_equal(dim(res$data), - c(dataset = 3, dataset = 1, statistics = 2, lat = 6, lon = 8)) + expect_equal(dim(res$data$rmsss), + c(nexp = 3, nobs = 1, lat = 6, lon = 8)) res <- CST_MultiMetric(exp = exp, obs = obs, metric = 'rmsss', multimodel = FALSE) - expect_equal(dim(res$data), - c(dataset = 2, dataset = 1, statistics = 2, lat = 6, lon =8)) + expect_equal(dim(res$data$rmsss), + c(nexp = 2, nobs = 1, lat = 6, lon = 8)) }) @@ -58,7 +69,7 @@ test_that("Sanity checks", { expect_error( CST_MultiMetric(exp = exp, obs = obs, metric = 1), paste0("Parameter 'metric' must be a character string indicating one ", - "of the options: 'correlation', 'rms' or 'rmse'")) + "of the options: 'correlation', 'rms', 'rmsss' or 'rpss'")) expect_error( CST_MultiMetric(exp = exp, obs = obs, metric = NA), "missing value where TRUE/FALSE needed") @@ -69,4 +80,15 @@ test_that("Sanity checks", { CST_MultiMetric(exp = exp, obs = obs, metric = "correlation", multimodel = NULL), "Parameter 'multimodel' must be a logical value.") + expect_error( + MultiMetric(exp = lonlat_data$exp, obs = lonlat_data$obs, metric = "rpss", + multimodel = TRUE), + "Element 'data' from parameters 'exp' and 'obs' should have dimmension names.") +exp <- lonlat_data$exp$data[1,,,,,] +obs <- lonlat_data$obs$data[1,,,,,] + expect_error( + MultiMetric(exp = exp, obs = obs, metric = "rpss", + multimodel = TRUE), + paste0("Dimension names of element 'data' from parameters 'exp' and ", + "'obs' should have the same name dimmension.")) }) -- GitLab From f743ed2775bf7fe62547ba303dcdb58c3d4bc173 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 9 Nov 2020 12:43:08 +0100 Subject: [PATCH 02/10] Fixing dependencies between s2dv and s2dverification --- DESCRIPTION | 4 +-- NAMESPACE | 13 +++++++- R/CST_Anomaly.R | 3 +- R/CST_BiasCorrection.R | 1 - R/CST_Calibration.R | 4 +-- R/CST_CategoricalEnsCombination.R | 2 +- R/CST_Load.R | 2 +- R/CST_MergeDims.R | 4 +-- R/CST_MultiMetric.R | 41 ++++++++++++++++++++----- R/CST_MultivarRMSE.R | 13 ++++---- R/CST_RegimesAssign.R | 6 ++-- R/CST_SaveExp.R | 6 ++-- R/CST_SplitDim.R | 4 +-- R/CST_WeatherRegimes.R | 4 +-- R/PlotCombinedMap.R | 2 +- R/PlotForecastPDF.R | 1 - R/PlotMostLikelyQuantileMap.R | 1 - R/PlotTriangles4Categories.R | 1 + man/CST_MultiMetric.Rd | 14 +++++---- man/CST_MultivarRMSE.Rd | 2 +- man/MultiMetric.Rd | 51 +++++++++++++++++++++++++++++++ 21 files changed, 136 insertions(+), 43 deletions(-) create mode 100644 man/MultiMetric.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 9020da9e..7fd6c5d4 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -47,7 +47,8 @@ Description: Exploits dynamical seasonal forecasts in order to provide Yiou et al. (2013) . Depends: R (>= 3.4.0), - maps + maps, + easyVerification Imports: s2dverification, s2dv, @@ -62,7 +63,6 @@ Imports: reshape2, ggplot2, RColorBrewer, - easyVerification, graphics, grDevices, stats, diff --git a/NAMESPACE b/NAMESPACE index 6996a8b2..98e6d022 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -28,6 +28,7 @@ export(Calibration) export(EnsClustering) export(MergeDims) export(MultiEOF) +export(MultiMetric) export(PlotCombinedMap) export(PlotForecastPDF) export(PlotMostLikelyQuantileMap) @@ -50,7 +51,6 @@ import(multiApply) import(ncdf4) import(qmap) import(rainfarmr) -import(s2dverification) import(stats) importFrom(ClimProjDiags,SelBox) importFrom(RColorBrewer,brewer.pal) @@ -58,6 +58,8 @@ importFrom(abind,abind) importFrom(data.table,CJ) importFrom(data.table,data.table) importFrom(data.table,setkey) +importFrom(easyVerification,climFairRpss) +importFrom(easyVerification,veriApply) importFrom(grDevices,adjustcolor) importFrom(grDevices,bmp) importFrom(grDevices,colorRampPalette) @@ -87,12 +89,21 @@ importFrom(maps,map) importFrom(plyr,.) importFrom(plyr,dlply) importFrom(reshape2,melt) +importFrom(s2dv,ColorBar) importFrom(s2dv,Corr) importFrom(s2dv,InsertDim) importFrom(s2dv,MeanDims) +importFrom(s2dv,PlotEquiMap) importFrom(s2dv,RMS) importFrom(s2dv,RMSSS) importFrom(s2dv,Reorder) +importFrom(s2dverification,ACC) +importFrom(s2dverification,Ano_CrossValid) +importFrom(s2dverification,Clim) +importFrom(s2dverification,EOF) +importFrom(s2dverification,Eno) +importFrom(s2dverification,Load) +importFrom(s2dverification,Mean1Dim) importFrom(s2dverification,Subset) importFrom(utils,glob2rx) importFrom(utils,head) diff --git a/R/CST_Anomaly.R b/R/CST_Anomaly.R index 6e33c341..52a786fb 100644 --- a/R/CST_Anomaly.R +++ b/R/CST_Anomaly.R @@ -18,7 +18,8 @@ #' #' @return A list with two S3 objects, 'exp' and 'obs', of the class 's2dv_cube', containing experimental and date-corresponding observational anomalies, respectively. These 's2dv_cube's can be ingested by other functions in CSTools. #' -#'@import s2dverification +#'@importFrom s2dverification Clim Ano_CrossValid +#'@importFrom s2dv InsertDim #' #'@examples #'# Example 1: diff --git a/R/CST_BiasCorrection.R b/R/CST_BiasCorrection.R index 1da7fb5b..9baf897b 100644 --- a/R/CST_BiasCorrection.R +++ b/R/CST_BiasCorrection.R @@ -11,7 +11,6 @@ #' #'@references Torralba, V., F.J. Doblas-Reyes, D. MacLeod, I. Christel and M. Davis (2017). Seasonal climate prediction: a new source of information for the management of wind energy resources. Journal of Applied Meteorology and Climatology, 56, 1231-1247, doi:10.1175/JAMC-D-16-0204.1. (CLIM4ENERGY, EUPORIAS, NEWA, RESILIENCE, SPECS) #' -#'@import s2dverification #'@import multiApply #'@examples #' diff --git a/R/CST_Calibration.R b/R/CST_Calibration.R index 83d61aba..874c351e 100644 --- a/R/CST_Calibration.R +++ b/R/CST_Calibration.R @@ -13,7 +13,7 @@ #'@param ncores is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one. #'@return an object of class \code{s2dv_cube} containing the calibrated forecasts in the element \code{$data} with the same dimensions as the one in the exp object. #' -#'@import s2dverification +#'@importFrom s2dv InsertDim #'@import abind #' #'@seealso \code{\link{CST_Load}} @@ -84,7 +84,7 @@ CST_Calibration <- function(exp, obs, cal.method = "mse_min", #'@param ncores is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one. #'@return an array containing the calibrated forecasts with the same dimensions as the \code{exp} array. #' -#'@import s2dverification +#'@importFrom s2dv InsertDim #'@import abind #' #'@seealso \code{\link{CST_Load}} diff --git a/R/CST_CategoricalEnsCombination.R b/R/CST_CategoricalEnsCombination.R index f72415e3..148858ee 100644 --- a/R/CST_CategoricalEnsCombination.R +++ b/R/CST_CategoricalEnsCombination.R @@ -55,7 +55,7 @@ #'@references Robertson, A. W., Lall, U., Zebiak, S. E., & Goddard, L. (2004). Improved combination of multiple atmospheric GCM ensembles for seasonal prediction. Monthly Weather Review, 132(12), 2732-2744. #'@references Van Schaeybroeck, B., & Vannitsem, S. (2019). Postprocessing of Long-Range Forecasts. In Statistical Postprocessing of Ensemble Forecasts (pp. 267-290). #' -#'@import s2dverification +#'@importFrom s2dv InsertDim #'@import abind #'@examples #' diff --git a/R/CST_Load.R b/R/CST_Load.R index 76d5deb3..65b695cd 100644 --- a/R/CST_Load.R +++ b/R/CST_Load.R @@ -9,7 +9,7 @@ #' @param ... Parameters that are automatically forwarded to the `s2dverification::Load` function. See details in `?s2dverification::Load`. #' @return A list with one or two S3 objects, named 'exp' and 'obs', of the class 's2dv_cube', containing experimental and date-corresponding observational data, respectively. These 's2dv_cube's can be ingested by other functions in CSTools. If the parameter `exp` in the call to `CST_Load` is set to `NULL`, then only the 'obs' component is returned, and viceversa. #' @author Nicolau Manubens, \email{nicolau.manubens@bsc.es} -#' @import s2dverification +#' @importFrom s2dverification Load #' @importFrom utils glob2rx #' @export #' @examples diff --git a/R/CST_MergeDims.R b/R/CST_MergeDims.R index a9923c58..720448c7 100644 --- a/R/CST_MergeDims.R +++ b/R/CST_MergeDims.R @@ -10,7 +10,7 @@ #'@param na.rm a logical indicating if the NA values should be removed or not. #' #'@import abind -#'@import s2dverification +#'@importFrom s2dverification Subset #'@examples #' #'data <- 1 : c(2 * 3 * 4 * 5 * 6 * 7) @@ -49,7 +49,7 @@ CST_MergeDims <- function(data, merge_dims = c('ftime', 'monthly'), rename_dim = #'@param na.rm a logical indicating if the NA values should be removed or not. #' #'@import abind -#'@import s2dverification +#'@importFrom s2dverification Subset #'@examples #' #'data <- 1 : 20 diff --git a/R/CST_MultiMetric.R b/R/CST_MultiMetric.R index 96e64f40..5723e142 100644 --- a/R/CST_MultiMetric.R +++ b/R/CST_MultiMetric.R @@ -9,11 +9,11 @@ #'@param metric a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms' or 'rmsss. #'@param multimodel a logical value indicating whether a Multi-Model Mean should be computed. #' -#'@param time_dim_name name of the temporal dimension where a mean will be applied. It can be NULL, the default value is 'ftime'. -#'@param memb_dim_name name of the member dimension. It can be NULL, the default value is 'member'. -#'@param sdate_dim_name name of the start date dimension or a dimension name identifiying the different forecast. It can be NULL, the default value is 'sdate'. -#'@return an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is an array with two datset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest first dimension correspons to the Multi-Model Mean. The third dimension contains the statistics selected. For metric \code{correlation}, the third dimension is of length four and they corresponds to the lower limit of the 95\% confidence interval, the statistics itselfs, the upper limit of the 95\% confidence interval and the 95\% significance level. For metric \code{rms}, the third dimension is length three and they corresponds to the lower limit of the 95\% confidence interval, the RMSE and the upper limit of the 95\% confidence interval. For metric \code{rmsss}, the third dimension is length two and they corresponds to the statistics itselfs and the p-value of the one-sided Fisher test with Ho: RMSSS = 0. -#'@seealso \code{\link[s2dverification]{Corr}}, \code{\link[s2dverification]{RMS}}, \code{\link[s2dverification]{RMSSS}} and \code{\link{CST_Load}} +#'@param time_dim name of the temporal dimension where a mean will be applied. It can be NULL, the default value is 'ftime'. +#'@param memb_dim name of the member dimension. It can be NULL, the default value is 'member'. +#'@param sdate_dim name of the start date dimension or a dimension name identifiying the different forecast. It can be NULL, the default value is 'sdate'. +#'@return an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest position in the first dimension correspons to the Multi-Model Mean. +#'@seealso \code{\link[s2dv]{Corr}}, \code{\link[s2dv]{RMS}}, \code{\link[s2dv]{RMSSS}} and \code{\link{CST_Load}} #'@references #'Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31.\url{http://link.springer.com/10.1007/s00382-018-4404-z} #' @@ -36,12 +36,14 @@ #'c(ano_exp, ano_obs) %<-% CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE) #'a <- CST_MultiMetric(exp = ano_exp, obs = ano_obs) #'str(a) +#'\donttest{ #'exp <- lonlat_data$exp #'obs <- lonlat_data$obs -#'a <- CST_MultiMetric(exp, obs, metric = 'rpss') +#'a <- CST_MultiMetric(exp, obs, metric = 'rpss', multimodel = FALSE) #'a <- CST_MultiMetric(exp, obs, metric = 'correlation') #'a <- CST_MultiMetric(exp, obs, metric = 'rms') #'a <- CST_MultiMetric(exp, obs, metric = 'rmsss') +#'} #'@export CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE, time_dim = 'ftime', memb_dim = 'member', @@ -55,7 +57,32 @@ CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE, exp$data <- result return(exp) } - +#'Multiple Metrics applied in Multiple Model Anomalies +#' +#'@author Mishra Niti, \email{niti.mishra@bsc.es} +#'@author Perez-Zanon Nuria, \email{nuria.perez@bsc.es} +#'@description This function calculates correlation (Anomaly Correlation Coefficient; ACC), root mean square error (RMS) and the root mean square error skill score (RMSSS) of individual anomaly models and multi-models mean (if desired) with the observations on arrays with named dimensions. +#' +#'@param exp a multidimensional array with named dimensions. +#'@param obs a multidimensional array with named dimensions. +#'@param metric a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms' or 'rmsss. +#'@param multimodel a logical value indicating whether a Multi-Model Mean should be computed. +#' +#'@param time_dim name of the temporal dimension where a mean will be applied. It can be NULL, the default value is 'ftime'. +#'@param memb_dim name of the member dimension. It can be NULL, the default value is 'member'. +#'@param sdate_dim name of the start date dimension or a dimension name identifiying the different forecast. It can be NULL, the default value is 'sdate'. +#'@return a list of arrays containing the statistics of the selected metric in the element \code{$data} which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest position in the first dimension correspons to the Multi-Model Mean. +#'@seealso \code{\link[s2dv]{Corr}}, \code{\link[s2dv]{RMS}}, \code{\link[s2dv]{RMSSS}} and \code{\link{CST_Load}} +#'@references +#'Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31.\url{http://link.springer.com/10.1007/s00382-018-4404-z} +#' +#'@importFrom s2dv MeanDims Reorder Corr RMS RMSSS InsertDim +#'@import abind +#'@importFrom easyVerification climFairRpss veriApply +#'@import stats +#'@examples +#'res <- MultiMetric(lonlat_data$exp$data, lonlat_data$obs$data) +#'@export MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE, time_dim = 'ftime', memb_dim = 'member', sdate_dim = 'sdate') { if (!is.null(names(dim(exp))) & !is.null(names(dim(obs)))) { diff --git a/R/CST_MultivarRMSE.R b/R/CST_MultivarRMSE.R index c5ab53b5..a82dde24 100644 --- a/R/CST_MultivarRMSE.R +++ b/R/CST_MultivarRMSE.R @@ -10,8 +10,8 @@ #' #'@return an object of class \code{s2dv_cube} containing the RMSE in the element \code{$data} which is an array with two datset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. An array with dimensions: c(number of exp, number of obs, 1 (the multivariate RMSE value), number of lat, number of lon) #' -#'@seealso \code{\link[s2dverification]{RMS}} and \code{\link{CST_Load}} -#'@import s2dverification +#'@seealso \code{\link[s2dv]{RMS}} and \code{\link{CST_Load}} +#'@importFrom s2dv RMS MeanDims #'@examples #'# Creation of sample s2dverification objects. These are not complete #'# s2dverification objects though. The Load function returns complete objects. @@ -108,17 +108,18 @@ CST_MultivarRMSE <- function(exp, obs, weight = NULL) { sumweights <- 0 for (j in 1 : nvar) { # seasonal average of anomalies - AvgExp <- MeanListDim(exp[[j]]$data, narm = TRUE, c(2, 4)) - AvgObs <- MeanListDim(obs[[j]]$data, narm = TRUE, c(2, 4)) + AvgExp <- MeanDims(exp[[j]]$data, c('member', 'ftime'), na.rm = TRUE) + AvgObs <- MeanDims(obs[[j]]$data, c('member', 'ftime'), na.rm = TRUE) # multivariate RMSE (weighted) - rmse <- RMS(AvgExp, AvgObs, posloop = 1, posRMS = 2, conf = FALSE) + rmse <- RMS(AvgExp, AvgObs, memb_dim = 'dataset', time_dim = 'sdate', + conf = FALSE)$rms stdev <- sd(AvgObs) mvrmse <- mvrmse + (rmse / stdev * as.numeric(weight[j])) sumweights <- sumweights + as.numeric(weight[j]) } mvrmse <- mvrmse / sumweights - names(dim(mvrmse)) <- c(dimnames[1], dimnames[1], 'statistics', dimnames[5 : 6]) + # names(dim(mvrmse)) <- c(dimnames[1], dimnames[1], 'statistics', dimnames[5 : 6]) exp_Datasets <- unlist(lapply(exp, function(x) { x[[which(names(x) == 'Datasets')]]})) exp_source_files <- unlist(lapply(exp, function(x) { diff --git a/R/CST_RegimesAssign.R b/R/CST_RegimesAssign.R index 60f0d704..fa6e27e8 100644 --- a/R/CST_RegimesAssign.R +++ b/R/CST_RegimesAssign.R @@ -27,7 +27,8 @@ #' that accounts for the serial dependence of the data with the same structure as Composite.)(only when composite = 'TRUE'), #' \code{$cluster} (array with the same dimensions as data (except latitude and longitude which are removed) indicating the ref_maps to which each point is allocated.) , #' \code{$frequency} (A vector of integers (from k=1,...k n reference maps) indicating the percentage of assignations corresponding to each map.), -#'@import s2dverification +#'@importFrom s2dverification ACC Mean1Dim +#'@importFrom s2dv InsertDim #'@import multiApply #'@examples #'\dontrun{ @@ -102,7 +103,8 @@ CST_RegimesAssign <- function(data, ref_maps, #' \code{$cluster} (array with the same dimensions as data (except latitude and longitude which are removed) indicating the ref_maps to which each point is allocated.) , #' \code{$frequency} (A vector of integers (from k = 1, ... k n reference maps) indicating the percentage of assignations corresponding to each map.), #' -#'@import s2dverification +#'@importFrom s2dverification ACC Mean1Dim Eno +#'@importFrom s2dv InsertDim #'@import multiApply #'@examples #'\dontrun{ diff --git a/R/CST_SaveExp.R b/R/CST_SaveExp.R index 1bdd4779..ac377101 100644 --- a/R/CST_SaveExp.R +++ b/R/CST_SaveExp.R @@ -17,7 +17,7 @@ #'@seealso \code{\link{CST_Load}}, \code{\link{as.s2dv_cube}} and \code{\link{s2dv_cube}} #' #'@import ncdf4 -#'@importFrom s2dv Reorder +#'@importFrom s2dv Reorder InsertDim #'@import multiApply #' #'@examples @@ -78,7 +78,7 @@ CST_SaveExp <- function(data, destination = "./CST_Data") { #' The path will be created with the name of the variable and each Datasets. #' #'@import ncdf4 -#'@importFrom s2dv Reorder +#'@importFrom s2dv Reorder InsertDim #'@import multiApply #' #'@examples @@ -142,7 +142,7 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates, if (length(dataset_pos) == 0) { warning("Element 'data' in parameter 'data' hasn't 'dataset' dimension. ", "All data is stored in the same 'dataset' folder.") - data$data <- InsertDim(var = data, posdim = 1, lendim = 1) + data$data <- InsertDim(data, posdim = 1, lendim = 1) names(dim(data))[1] <- "dataset" dimname <- c("dataset", dimname) dataset_pos = 1 diff --git a/R/CST_SplitDim.R b/R/CST_SplitDim.R index 3326d435..6332684f 100644 --- a/R/CST_SplitDim.R +++ b/R/CST_SplitDim.R @@ -11,7 +11,7 @@ #'@param new_dim_name a character string indicating the name of the new dimension. #' #'@import abind -#'@import s2dverification +#'@importFrom s2dverification Subset #'@examples #' #'data <- 1 : 20 @@ -73,7 +73,7 @@ CST_SplitDim <- function(data, split_dim = 'time', indices = NULL, #'@param freq a character string indicating the frequency: by 'day', 'month' and 'year' or 'monthly' (by default). 'month' identifies months between 1 and 12 independetly of the year they belong to, while 'monthly' differenciates months from different years. Parameter 'freq' can also be numeric indicating the length in which to subset the dimension. #'@param new_dim_name a character string indicating the name of the new dimension. #'@import abind -#'@import s2dverification +#'@importFrom s2dverification Subset #'@examples #' #'data <- 1 : 20 diff --git a/R/CST_WeatherRegimes.R b/R/CST_WeatherRegimes.R index 72ab3987..e7cc925c 100644 --- a/R/CST_WeatherRegimes.R +++ b/R/CST_WeatherRegimes.R @@ -34,7 +34,7 @@ #' \code{cluster} (A matrix or vector with integers (from 1:k) indicating the cluster to which each time step is allocated.), #' \code{persistence} (Percentage of days in a month/season before a cluster is replaced for a new one (only if method=’kmeans’ has been selected.)), #' \code{frequency} (Percentage of days in a month/season belonging to each cluster (only if method=’kmeans’ has been selected).), -#'@import s2dverification +#'@importFrom s2dverification EOF #'@import multiApply #'@examples #'\dontrun{ @@ -109,7 +109,7 @@ CST_WeatherRegimes <- function(data, ncenters = NULL, #' \code{cluster} (A matrix or vector with integers (from 1:k) indicating the cluster to which each time step is allocated.), #' \code{persistence} (Percentage of days in a month/season before a cluster is replaced for a new one (only if method=’kmeans’ has been selected.)), #' \code{frequency} (Percentage of days in a month/season belonging to each cluster (only if method=’kmeans’ has been selected).), -#'@import s2dverification +#'@importFrom s2dverification EOF #'@import multiApply #'@examples #'\dontrun{ diff --git a/R/PlotCombinedMap.R b/R/PlotCombinedMap.R index 7f8f5d64..2aee71b8 100644 --- a/R/PlotCombinedMap.R +++ b/R/PlotCombinedMap.R @@ -25,7 +25,7 @@ #'@seealso \code{PlotCombinedMap} and \code{PlotEquiMap} #' -#'@import s2dverification +#'@importFrom s2dv PlotEquiMap ColorBar #'@importFrom maps map #'@importFrom graphics box image layout mtext par plot.new #'@importFrom grDevices adjustcolor bmp colorRampPalette dev.cur dev.new dev.off hcl jpeg pdf png postscript svg tiff diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 9cd4dfeb..5afdcbe4 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -23,7 +23,6 @@ #'@importFrom reshape2 melt #'@importFrom plyr . #'@importFrom plyr dlply -#'@import s2dverification #' #'@examples #'fcsts <- data.frame(fcst1 = rnorm(10), fcst2 = rnorm(10, 0.5, 1.2), diff --git a/R/PlotMostLikelyQuantileMap.R b/R/PlotMostLikelyQuantileMap.R index 57be8643..b4e974a8 100644 --- a/R/PlotMostLikelyQuantileMap.R +++ b/R/PlotMostLikelyQuantileMap.R @@ -12,7 +12,6 @@ #'@param ... additional parameters to be sent to \code{PlotCombinedMap} and \code{PlotEquiMap}. #'@seealso \code{PlotCombinedMap} and \code{PlotEquiMap} #' -#'@import s2dverification #'@importFrom maps map #'@importFrom graphics box image layout mtext par plot.new #'@importFrom grDevices adjustcolor bmp colorRampPalette dev.cur dev.new dev.off hcl jpeg pdf png postscript svg tiff diff --git a/R/PlotTriangles4Categories.R b/R/PlotTriangles4Categories.R index fcfb36bb..14d578a5 100644 --- a/R/PlotTriangles4Categories.R +++ b/R/PlotTriangles4Categories.R @@ -72,6 +72,7 @@ #'@importFrom grDevices dev.new dev.off dev.cur #'@importFrom graphics plot points polygon text title axis #'@importFrom RColorBrewer brewer.pal +#'@importFrom s2dv ColorBar #'@export PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL, toptitle = NULL, sig_data = NULL, diff --git a/man/CST_MultiMetric.Rd b/man/CST_MultiMetric.Rd index 4526fa0e..c9aacd26 100644 --- a/man/CST_MultiMetric.Rd +++ b/man/CST_MultiMetric.Rd @@ -23,14 +23,14 @@ CST_MultiMetric( \item{multimodel}{a logical value indicating whether a Multi-Model Mean should be computed.} -\item{time_dim_name}{name of the temporal dimension where a mean will be applied. It can be NULL, the default value is 'ftime'.} +\item{time_dim}{name of the temporal dimension where a mean will be applied. It can be NULL, the default value is 'ftime'.} -\item{memb_dim_name}{name of the member dimension. It can be NULL, the default value is 'member'.} +\item{memb_dim}{name of the member dimension. It can be NULL, the default value is 'member'.} -\item{sdate_dim_name}{name of the start date dimension or a dimension name identifiying the different forecast. It can be NULL, the default value is 'sdate'.} +\item{sdate_dim}{name of the start date dimension or a dimension name identifiying the different forecast. It can be NULL, the default value is 'sdate'.} } \value{ -an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is an array with two datset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest first dimension correspons to the Multi-Model Mean. The third dimension contains the statistics selected. For metric \code{correlation}, the third dimension is of length four and they corresponds to the lower limit of the 95\% confidence interval, the statistics itselfs, the upper limit of the 95\% confidence interval and the 95\% significance level. For metric \code{rms}, the third dimension is length three and they corresponds to the lower limit of the 95\% confidence interval, the RMSE and the upper limit of the 95\% confidence interval. For metric \code{rmsss}, the third dimension is length two and they corresponds to the statistics itselfs and the p-value of the one-sided Fisher test with Ho: RMSSS = 0. +an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest position in the first dimension correspons to the Multi-Model Mean. } \description{ This function calculates correlation (Anomaly Correlation Coefficient; ACC), root mean square error (RMS) and the root mean square error skill score (RMSSS) of individual anomaly models and multi-models mean (if desired) with the observations. @@ -50,18 +50,20 @@ attr(obs, 'class') <- 's2dv_cube' c(ano_exp, ano_obs) \%<-\% CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE) a <- CST_MultiMetric(exp = ano_exp, obs = ano_obs) str(a) +\donttest{ exp <- lonlat_data$exp obs <- lonlat_data$obs -a <- CST_MultiMetric(exp, obs, metric = 'rpss') +a <- CST_MultiMetric(exp, obs, metric = 'rpss', multimodel = FALSE) a <- CST_MultiMetric(exp, obs, metric = 'correlation') a <- CST_MultiMetric(exp, obs, metric = 'rms') a <- CST_MultiMetric(exp, obs, metric = 'rmsss') } +} \references{ Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31.\url{http://link.springer.com/10.1007/s00382-018-4404-z} } \seealso{ -\code{\link[s2dverification]{Corr}}, \code{\link[s2dverification]{RMS}}, \code{\link[s2dverification]{RMSSS}} and \code{\link{CST_Load}} +\code{\link[s2dv]{Corr}}, \code{\link[s2dv]{RMS}}, \code{\link[s2dv]{RMSSS}} and \code{\link{CST_Load}} } \author{ Mishra Niti, \email{niti.mishra@bsc.es} diff --git a/man/CST_MultivarRMSE.Rd b/man/CST_MultivarRMSE.Rd index 24af608c..c318c105 100644 --- a/man/CST_MultivarRMSE.Rd +++ b/man/CST_MultivarRMSE.Rd @@ -57,7 +57,7 @@ a <- CST_MultivarRMSE(exp = ano_exp, obs = ano_obs, weight = weight) str(a) } \seealso{ -\code{\link[s2dverification]{RMS}} and \code{\link{CST_Load}} +\code{\link[s2dv]{RMS}} and \code{\link{CST_Load}} } \author{ Deborah Verfaillie, \email{deborah.verfaillie@bsc.es} diff --git a/man/MultiMetric.Rd b/man/MultiMetric.Rd new file mode 100644 index 00000000..f688f0a3 --- /dev/null +++ b/man/MultiMetric.Rd @@ -0,0 +1,51 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_MultiMetric.R +\name{MultiMetric} +\alias{MultiMetric} +\title{Multiple Metrics applied in Multiple Model Anomalies} +\usage{ +MultiMetric( + exp, + obs, + metric = "correlation", + multimodel = TRUE, + time_dim = "ftime", + memb_dim = "member", + sdate_dim = "sdate" +) +} +\arguments{ +\item{exp}{a multidimensional array with named dimensions.} + +\item{obs}{a multidimensional array with named dimensions.} + +\item{metric}{a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms' or 'rmsss.} + +\item{multimodel}{a logical value indicating whether a Multi-Model Mean should be computed.} + +\item{time_dim}{name of the temporal dimension where a mean will be applied. It can be NULL, the default value is 'ftime'.} + +\item{memb_dim}{name of the member dimension. It can be NULL, the default value is 'member'.} + +\item{sdate_dim}{name of the start date dimension or a dimension name identifiying the different forecast. It can be NULL, the default value is 'sdate'.} +} +\value{ +a list of arrays containing the statistics of the selected metric in the element \code{$data} which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest position in the first dimension correspons to the Multi-Model Mean. +} +\description{ +This function calculates correlation (Anomaly Correlation Coefficient; ACC), root mean square error (RMS) and the root mean square error skill score (RMSSS) of individual anomaly models and multi-models mean (if desired) with the observations on arrays with named dimensions. +} +\examples{ +res <- MultiMetric(lonlat_data$exp$data, lonlat_data$obs$data) +} +\references{ +Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31.\url{http://link.springer.com/10.1007/s00382-018-4404-z} +} +\seealso{ +\code{\link[s2dv]{Corr}}, \code{\link[s2dv]{RMS}}, \code{\link[s2dv]{RMSSS}} and \code{\link{CST_Load}} +} +\author{ +Mishra Niti, \email{niti.mishra@bsc.es} + +Perez-Zanon Nuria, \email{nuria.perez@bsc.es} +} -- GitLab From 40aaaa0ab70d3a0661ccf1b17eee4febfb08e54a Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 10 Nov 2020 10:26:28 +0100 Subject: [PATCH 03/10] Improve documentation --- R/CST_MultiMetric.R | 2 +- man/CST_MultiMetric.Rd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/CST_MultiMetric.R b/R/CST_MultiMetric.R index 5723e142..689c7ccc 100644 --- a/R/CST_MultiMetric.R +++ b/R/CST_MultiMetric.R @@ -6,7 +6,7 @@ #' #'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the anomaly of the seasonal forecast experiments data in the element named \code{$data}. #'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the anomaly of observed data in the element named \code{$data}. -#'@param metric a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms' or 'rmsss. +#'@param metric a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms', 'rmsss' and 'rpss'. If 'rpss' is chossen the terciles probabilities are evaluated. #'@param multimodel a logical value indicating whether a Multi-Model Mean should be computed. #' #'@param time_dim name of the temporal dimension where a mean will be applied. It can be NULL, the default value is 'ftime'. diff --git a/man/CST_MultiMetric.Rd b/man/CST_MultiMetric.Rd index c9aacd26..e8e047dc 100644 --- a/man/CST_MultiMetric.Rd +++ b/man/CST_MultiMetric.Rd @@ -19,7 +19,7 @@ CST_MultiMetric( \item{obs}{an object of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the anomaly of observed data in the element named \code{$data}.} -\item{metric}{a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms' or 'rmsss.} +\item{metric}{a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms', 'rmsss' and 'rpss'. If 'rpss' is chossen the terciles probabilities are evaluated.} \item{multimodel}{a logical value indicating whether a Multi-Model Mean should be computed.} -- GitLab From 3abd1199e9cf41034fdb23cf61e17a771fd0ad10 Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 10 Nov 2020 11:18:44 +0100 Subject: [PATCH 04/10] Fixing InsertDim use from 2sdv and s2dverification --- NAMESPACE | 1 + R/CST_CategoricalEnsCombination.R | 10 +++++----- R/CST_MultiMetric.R | 2 +- 3 files changed, 7 insertions(+), 6 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 98e6d022..c1e587ea 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -102,6 +102,7 @@ importFrom(s2dverification,Ano_CrossValid) importFrom(s2dverification,Clim) importFrom(s2dverification,EOF) importFrom(s2dverification,Eno) +importFrom(s2dverification,InsertDim) importFrom(s2dverification,Load) importFrom(s2dverification,Mean1Dim) importFrom(s2dverification,Subset) diff --git a/R/CST_CategoricalEnsCombination.R b/R/CST_CategoricalEnsCombination.R index 148858ee..c682c534 100644 --- a/R/CST_CategoricalEnsCombination.R +++ b/R/CST_CategoricalEnsCombination.R @@ -55,7 +55,7 @@ #'@references Robertson, A. W., Lall, U., Zebiak, S. E., & Goddard, L. (2004). Improved combination of multiple atmospheric GCM ensembles for seasonal prediction. Monthly Weather Review, 132(12), 2732-2744. #'@references Van Schaeybroeck, B., & Vannitsem, S. (2019). Postprocessing of Long-Range Forecasts. In Statistical Postprocessing of Ensemble Forecasts (pp. 267-290). #' -#'@importFrom s2dv InsertDim +#'@importFrom s2dverification InsertDim #'@import abind #'@examples #' @@ -277,8 +277,8 @@ comb.dims <- function(arr.in, dims.to.combine){ ui = constr.mtrx, ci = constr.vec, freq.per.mdl.at.obs = freq.per.mdl.at.obs) init.par <- optim.tmp$par * (1 - abs(rnorm(amt.coeff, 0, 0.01))) - var.cat.fc[ , eval.dexes] <- apply(InsertDim( - InsertDim(optim.tmp$par, lendim = amt.cat, posdim = 2), + var.cat.fc[ , eval.dexes] <- apply(s2dverification::InsertDim( + s2dverification::InsertDim(optim.tmp$par, lendim = amt.cat, posdim = 2), lendim = amt.sdate.ev, posdim = 3) * freq.per.mdl.ev[ , , , drop = FALSE], c(2,3), sum, na.rm = TRUE) } else if (cat.method == "comb") { @@ -363,7 +363,7 @@ comb.dims <- function(arr.in, dims.to.combine){ .funct.optim.grad <- function(par, freq.per.mdl.at.obs){ amt.model <- dim(freq.per.mdl.at.obs)[1] - return(-apply(freq.per.mdl.at.obs/InsertDim(drop(par %*% freq.per.mdl.at.obs), + return(-apply(freq.per.mdl.at.obs/s2dverification::InsertDim(drop(par %*% freq.per.mdl.at.obs), lendim = amt.model, posdim = 1), c(1), mean, na.rm = TRUE)) } @@ -373,7 +373,7 @@ comb.dims <- function(arr.in, dims.to.combine){ amt.mdl <- mdl.feat$amt.mdl mdl.msk.tmp <- mdl.feat$mdl.msk amt.coeff <- amt.mdl + 1 - msk.fc.obs <- (cat.fc == InsertDim(cat.obs, posdim = 1, lendim = amt.mbr)) + msk.fc.obs <- (cat.fc == s2dverification::InsertDim(cat.obs, posdim = 1, lendim = amt.mbr)) freq.per.mdl.at.obs <- array(NA, c(amt.coeff, amt.sdate)) for (i.mdl in seq(1, amt.mdl)){ freq.per.mdl.at.obs[i.mdl, ] <- apply(msk.fc.obs[mdl.msk.tmp[i.mdl, ], , drop = FALSE], diff --git a/R/CST_MultiMetric.R b/R/CST_MultiMetric.R index 689c7ccc..41a3d535 100644 --- a/R/CST_MultiMetric.R +++ b/R/CST_MultiMetric.R @@ -187,7 +187,7 @@ MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE, # seasonal avg of anomalies for multi-model AvgExp_MMM <- MeanDims(AvgExp, c(dataset_dim), na.rm = TRUE) pos_dataset <- which(names(dim(AvgExp)) == dataset_dim) - AvgExp_MMM <- InsertDim(AvgExp_MMM, posdim = pos_dataset, lendim = 1, + AvgExp_MMM <- s2dv::InsertDim(AvgExp_MMM, posdim = pos_dataset, lendim = 1, name = dataset_dim) AvgExp <- abind(AvgExp_MMM, AvgExp, along = pos_dataset) names(dim(AvgExp)) <- names(dim(AvgExp_MMM)) -- GitLab From d7bc662368c2dbf8d069c0ee9009dc6054e7d1a3 Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 11 Nov 2020 13:10:13 +0100 Subject: [PATCH 05/10] InsertDim s2dv and s2dverification --- R/CST_Anomaly.R | 11 +++++------ R/CST_Calibration.R | 14 +++++++------- R/CST_CategoricalEnsCombination.R | 5 +++-- R/CST_RegimesAssign.R | 20 +++++++++----------- R/CST_SaveExp.R | 8 +++++--- R/PlotForecastPDF.R | 6 +++--- 6 files changed, 32 insertions(+), 32 deletions(-) diff --git a/R/CST_Anomaly.R b/R/CST_Anomaly.R index 52a786fb..ef616225 100644 --- a/R/CST_Anomaly.R +++ b/R/CST_Anomaly.R @@ -18,8 +18,7 @@ #' #' @return A list with two S3 objects, 'exp' and 'obs', of the class 's2dv_cube', containing experimental and date-corresponding observational anomalies, respectively. These 's2dv_cube's can be ingested by other functions in CSTools. #' -#'@importFrom s2dverification Clim Ano_CrossValid -#'@importFrom s2dv InsertDim +#'@importFrom s2dverification Clim Ano_CrossValid InsertDim #' #'@examples #'# Example 1: @@ -166,12 +165,12 @@ CST_Anomaly <- function(exp = NULL, obs = NULL, cross = FALSE, memb = TRUE, clim_exp <- tmp$clim_exp clim_obs <- tmp$clim_obs } else { - clim_exp <- InsertDim(tmp$clim_exp, 2, dim_exp[2]) - clim_obs <- InsertDim(tmp$clim_obs, 2, dim_obs[2]) + clim_exp <- s2dverification::InsertDim(tmp$clim_exp, 2, dim_exp[2]) + clim_obs <- s2dverification::InsertDim(tmp$clim_obs, 2, dim_obs[2]) } - clim_exp <- InsertDim(clim_exp, 3, dim_exp[3]) - clim_obs <- InsertDim(clim_obs, 3, dim_obs[3]) + clim_exp <- s2dverification::InsertDim(clim_exp, 3, dim_exp[3]) + clim_obs <- s2dverification::InsertDim(clim_obs, 3, dim_obs[3]) ano <- NULL ano$ano_exp <- exp$data - clim_exp ano$ano_obs <- obs$data - clim_obs diff --git a/R/CST_Calibration.R b/R/CST_Calibration.R index 874c351e..fd381790 100644 --- a/R/CST_Calibration.R +++ b/R/CST_Calibration.R @@ -13,7 +13,7 @@ #'@param ncores is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one. #'@return an object of class \code{s2dv_cube} containing the calibrated forecasts in the element \code{$data} with the same dimensions as the one in the exp object. #' -#'@importFrom s2dv InsertDim +#'@importFrom s2dverification InsertDim #'@import abind #' #'@seealso \code{\link{CST_Load}} @@ -84,7 +84,7 @@ CST_Calibration <- function(exp, obs, cal.method = "mse_min", #'@param ncores is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one. #'@return an array containing the calibrated forecasts with the same dimensions as the \code{exp} array. #' -#'@importFrom s2dv InsertDim +#'@importFrom s2dverification InsertDim #'@import abind #' #'@seealso \code{\link{CST_Load}} @@ -251,7 +251,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o .calc.obs.fc.quant <- function(obs, fc){ #function to calculate different quantities of a series of ensemble forecasts and corresponding observations amt.mbr <- dim(fc)[1] - obs.per.ens <- InsertDim(obs, posdim = 1, lendim = amt.mbr) + obs.per.ens <- s2dverification::InsertDim(obs, posdim = 1, lendim = amt.mbr) fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE) cor.obs.fc <- cor(fc.ens.av, obs, use = "complete.obs") obs.av <- mean(obs, na.rm = TRUE) @@ -271,7 +271,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o .calc.obs.fc.quant.ext <- function(obs, fc){ #extended function to calculate different quantities of a series of ensemble forecasts and corresponding observations amt.mbr <- dim(fc)[1] - obs.per.ens <- InsertDim(obs, posdim = 1, lendim = amt.mbr) + obs.per.ens <- s2dverification::InsertDim(obs, posdim = 1, lendim = amt.mbr) fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE) cor.obs.fc <- cor(fc.ens.av, obs, use = "complete.obs") obs.av <- mean(obs, na.rm = TRUE) @@ -296,7 +296,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE) fc.ens.av.av <- mean(fc.ens.av, na.rm = TRUE) fc.ens.av.sd <- sd(fc.ens.av, na.rm = TRUE) - fc.ens.av.per.ens <- InsertDim(fc.ens.av, posdim = 1, lendim = amt.mbr) + fc.ens.av.per.ens <- s2dverification::InsertDim(fc.ens.av, posdim = 1, lendim = amt.mbr) fc.ens.sd <- apply(fc, c(2), sd, na.rm = TRUE) fc.ens.var.av.sqrt <- sqrt(mean(fc.ens.sd^2, na.rm = TRUE)) fc.dev <- fc - fc.ens.av.per.ens @@ -322,10 +322,10 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o .calc.fc.quant.ext <- function(fc){ #extended function to calculate different quantities of a series of ensemble forecasts amt.mbr <- dim(fc)[1] - repmat1.tmp <- InsertDim(fc, posdim = 1, lendim = amt.mbr) + repmat1.tmp <- s2dverification::InsertDim(fc, posdim = 1, lendim = amt.mbr) repmat2.tmp <- aperm(repmat1.tmp, c(2, 1, 3)) spr.abs <- apply(abs(repmat1.tmp - repmat2.tmp), c(3), mean, na.rm = TRUE) - spr.abs.per.ens <- InsertDim(spr.abs, posdim = 1, lendim = amt.mbr) + spr.abs.per.ens <- s2dverification::InsertDim(spr.abs, posdim = 1, lendim = amt.mbr) return( append(.calc.fc.quant(fc), diff --git a/R/CST_CategoricalEnsCombination.R b/R/CST_CategoricalEnsCombination.R index c682c534..26ac5449 100644 --- a/R/CST_CategoricalEnsCombination.R +++ b/R/CST_CategoricalEnsCombination.R @@ -89,7 +89,7 @@ CST_CategoricalEnsCombination <- function(exp, obs, cat.method = "pool", eval.me exp$data <- .CategoricalEnsCombination.wrap(fc = exp$data, obs = obs$data, cat.method = cat.method, eval.method = eval.method, amt.cat = amt.cat, ...) names.dim.tmp[which(names.dim.tmp == "member")] <- "category" names(dim(exp$data)) <- names.dim.tmp - exp$data <- InsertDim(exp$data, lendim = 1, posdim = 2) + exp$data <- s2dverification::InsertDim(exp$data, lendim = 1, posdim = 2) names(dim(exp$data))[2] <- "member" exp$Datasets <- c(exp$Datasets, obs$Datasets) exp$source_files <- c(exp$source_files, obs$source_files) @@ -363,7 +363,8 @@ comb.dims <- function(arr.in, dims.to.combine){ .funct.optim.grad <- function(par, freq.per.mdl.at.obs){ amt.model <- dim(freq.per.mdl.at.obs)[1] - return(-apply(freq.per.mdl.at.obs/s2dverification::InsertDim(drop(par %*% freq.per.mdl.at.obs), + return(-apply(freq.per.mdl.at.obs/s2dverification::InsertDim(drop(par %*% +freq.per.mdl.at.obs), lendim = amt.model, posdim = 1), c(1), mean, na.rm = TRUE)) } diff --git a/R/CST_RegimesAssign.R b/R/CST_RegimesAssign.R index fa6e27e8..72ae69f1 100644 --- a/R/CST_RegimesAssign.R +++ b/R/CST_RegimesAssign.R @@ -27,8 +27,7 @@ #' that accounts for the serial dependence of the data with the same structure as Composite.)(only when composite = 'TRUE'), #' \code{$cluster} (array with the same dimensions as data (except latitude and longitude which are removed) indicating the ref_maps to which each point is allocated.) , #' \code{$frequency} (A vector of integers (from k=1,...k n reference maps) indicating the percentage of assignations corresponding to each map.), -#'@importFrom s2dverification ACC Mean1Dim -#'@importFrom s2dv InsertDim +#'@importFrom s2dverification ACC Mean1Dim InsertDim #'@import multiApply #'@examples #'\dontrun{ @@ -103,8 +102,7 @@ CST_RegimesAssign <- function(data, ref_maps, #' \code{$cluster} (array with the same dimensions as data (except latitude and longitude which are removed) indicating the ref_maps to which each point is allocated.) , #' \code{$frequency} (A vector of integers (from k = 1, ... k n reference maps) indicating the percentage of assignations corresponding to each map.), #' -#'@importFrom s2dverification ACC Mean1Dim Eno -#'@importFrom s2dv InsertDim +#'@importFrom s2dverification ACC Mean1Dim Eno InsertDim #'@import multiApply #'@examples #'\dontrun{ @@ -192,9 +190,9 @@ RegimesAssign <- function(data, ref_maps, lat, method = "distance", composite = if (any(is.na(index))) { recon <-list( - composite = InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), + composite = s2dverification::InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), postime, dim(ref_maps)['composite.cluster']), - pvalue = InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), + pvalue = s2dverification::InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), postime, dim(ref_maps)['composite.cluster'])) } else { if (memb) { @@ -257,7 +255,7 @@ RegimesAssign <- function(data, ref_maps, lat, method = "distance", composite = which(names(dim(target)) == 'lon'))) # weights are defined - latWeights <- InsertDim(sqrt(cos(lat * pi / 180)), 2, dim(ref)[3]) + latWeights <- s2dverification::InsertDim(sqrt(cos(lat * pi / 180)), 2, dim(ref)[3]) rmsdiff <- function(x, y) { @@ -280,11 +278,11 @@ RegimesAssign <- function(data, ref_maps, lat, method = "distance", composite = corr <- rep(NA, nclust) for (i in 1:nclust) { corr[i] <- - ACC(InsertDim(InsertDim( - InsertDim(ref[i, , ] * latWeights, 1, 1), 2, 1 + ACC(s2dverification::InsertDim(InsertDim( + s2dverification::InsertDim(ref[i, , ] * latWeights, 1, 1), 2, 1 ), 3, 1), - InsertDim(InsertDim( - InsertDim(target * latWeights, 1, 1), 2, 1 + s2dverification::InsertDim(InsertDim( + s2dverification::InsertDim(target * latWeights, 1, 1), 2, 1 ), 3, 1))$ACC[2] } assign <- which(corr == max(corr)) diff --git a/R/CST_SaveExp.R b/R/CST_SaveExp.R index ac377101..4ae560bf 100644 --- a/R/CST_SaveExp.R +++ b/R/CST_SaveExp.R @@ -17,7 +17,8 @@ #'@seealso \code{\link{CST_Load}}, \code{\link{as.s2dv_cube}} and \code{\link{s2dv_cube}} #' #'@import ncdf4 -#'@importFrom s2dv Reorder InsertDim +#'@importFrom s2dv Reorder +#'@importFrom s2dverification InsertDim #'@import multiApply #' #'@examples @@ -78,7 +79,8 @@ CST_SaveExp <- function(data, destination = "./CST_Data") { #' The path will be created with the name of the variable and each Datasets. #' #'@import ncdf4 -#'@importFrom s2dv Reorder InsertDim +#'@importFrom s2dv Reorder +#'@importFrom s2dverification InsertDim #'@import multiApply #' #'@examples @@ -142,7 +144,7 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates, if (length(dataset_pos) == 0) { warning("Element 'data' in parameter 'data' hasn't 'dataset' dimension. ", "All data is stored in the same 'dataset' folder.") - data$data <- InsertDim(data, posdim = 1, lendim = 1) + data$data <- s2dverification::InsertDim(data, posdim = 1, lendim = 1) names(dim(data))[1] <- "dataset" dimname <- c("dataset", dimname) dataset_pos = 1 diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 5afdcbe4..8c4b1068 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -23,7 +23,7 @@ #'@importFrom reshape2 melt #'@importFrom plyr . #'@importFrom plyr dlply -#' +#'@importFrom s2dverification InsertDim #'@examples #'fcsts <- data.frame(fcst1 = rnorm(10), fcst2 = rnorm(10, 0.5, 1.2), #' fcst3 = rnorm(10, -0.5, 0.9)) @@ -108,7 +108,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N if (length(tercile.limits) != 2) { stop("Provide two tercile limits") } - tercile.limits <- InsertDim(tercile.limits, 1, npanels) + tercile.limits <- s2dverification::InsertDim(tercile.limits, 1, npanels) } else if (is.array(tercile.limits)) { if (length(dim(tercile.limits)) == 2) { if (dim(tercile.limits)[2] != 2) { @@ -135,7 +135,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N if (length(extreme.limits) != 2) { stop("Provide two extreme limits") } - extreme.limits <- InsertDim(extreme.limits, 1, npanels) + extreme.limits <- s2dverification::InsertDim(extreme.limits, 1, npanels) } else if (is.array(extreme.limits)) { if (length(dim(extreme.limits)) == 2) { if (dim(extreme.limits)[2] != 2) { -- GitLab From b5bfffe2f04c36e7b5bce35dd086867c693f58b7 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 30 Nov 2020 14:39:06 +0100 Subject: [PATCH 06/10] Functions addapted to new s2dv version --- NAMESPACE | 1 - R/CST_MultiMetric.R | 24 ++++++++++++++++++------ R/CST_MultivarRMSE.R | 2 +- tests/testthat/test-CST_MultiMetric.R | 15 ++++++--------- 4 files changed, 25 insertions(+), 17 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index c1e587ea..29dc0811 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -45,7 +45,6 @@ export(WeatherRegime) export(as.s2dv_cube) export(s2dv_cube) import(abind) -import(easyVerification) import(ggplot2) import(multiApply) import(ncdf4) diff --git a/R/CST_MultiMetric.R b/R/CST_MultiMetric.R index 41a3d535..a932b844 100644 --- a/R/CST_MultiMetric.R +++ b/R/CST_MultiMetric.R @@ -19,8 +19,9 @@ #' #'@importFrom s2dv MeanDims Reorder Corr RMS RMSSS InsertDim #'@import abind -#'@import easyVerification +#'@importFrom easyVerification climFairRpss veriApply #'@import stats +#'@import multiApply #'@examples #'library(zeallot) #'mod <- 1 : (2 * 3 * 4 * 5 * 6 * 7) @@ -80,6 +81,7 @@ CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE, #'@import abind #'@importFrom easyVerification climFairRpss veriApply #'@import stats +#'@import multiApply #'@examples #'res <- MultiMetric(lonlat_data$exp$data, lonlat_data$obs$data) #'@export @@ -166,13 +168,23 @@ MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE, } AvgExp <- MeanDims(exp, time_dim, na.rm = TRUE) AvgObs <- MeanDims(obs, time_dim, na.rm = TRUE) + dif_dims <- which(dim(AvgExp) != dim(AvgObs)) + dif_dims <- names(dif_dims[-which(names(dif_dims) == memb_dim)]) + lapply(dif_dims, function(x) { + names(dim(AvgExp))[which(names(dim(AvgExp)) == x)] <<- paste0(dif_dims, '_exp')}) + pos_memb <- which(names(dim(AvgExp)) == memb_dim) dim(AvgObs) <- dim(AvgObs)[-pos_memb] AvgExp <- Reorder(AvgExp, c(names(dim(AvgExp))[-pos_memb], memb_dim)) pos_memb <- which(names(dim(AvgExp)) == memb_dim) pos_sdate <- which(names(dim(AvgExp)) == sdate_dim) - corr <- veriApply('FairRpss', fcst = AvgExp, obs = AvgObs, ensdim = pos_memb, - tdim = pos_sdate, prob = c(1/3, 2/3)) + corr <- Apply(list(AvgExp, AvgObs), + target_dims = list(c('sdate', 'lat', 'lon', 'member'), c('sdate', 'lat', 'lon')), + fun = function(x, y) { + veriApply('FairRpss', fcst = x, obs = y, + ensdim = which(names(dim(x)) == 'member'), + tdim = which(names(dim(x)) == 'sdate'), + prob = c(1/3, 2/3))}) } else if (metric %in% c('correlation', 'rms', 'rmsss')) { AvgExp <- MeanDims(exp, c(memb_dim, time_dim), na.rm = TRUE) @@ -193,11 +205,11 @@ MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE, names(dim(AvgExp)) <- names(dim(AvgExp_MMM)) } if (metric == 'correlation') { - corr <- s2dv::Corr(AvgExp, AvgObs, memb_dim = dataset_dim, time_dim = sdate_dim) + corr <- s2dv::Corr(AvgExp, AvgObs, dat_dim = dataset_dim, time_dim = sdate_dim) } else if (metric == 'rms') { - corr <- s2dv::RMS(AvgExp, AvgObs, memb_dim = dataset_dim, time_dim = sdate_dim) + corr <- s2dv::RMS(AvgExp, AvgObs, dat_dim = dataset_dim, time_dim = sdate_dim) } else if (metric == 'rmsss') { - corr <- s2dv::RMSSS(AvgExp, AvgObs, memb_dim = dataset_dim, time_dim = sdate_dim) + corr <- s2dv::RMSSS(AvgExp, AvgObs, dat_dim = dataset_dim, time_dim = sdate_dim) } } else { stop("Parameter 'metric' must be a character string indicating ", diff --git a/R/CST_MultivarRMSE.R b/R/CST_MultivarRMSE.R index a82dde24..59d0a759 100644 --- a/R/CST_MultivarRMSE.R +++ b/R/CST_MultivarRMSE.R @@ -111,7 +111,7 @@ CST_MultivarRMSE <- function(exp, obs, weight = NULL) { AvgExp <- MeanDims(exp[[j]]$data, c('member', 'ftime'), na.rm = TRUE) AvgObs <- MeanDims(obs[[j]]$data, c('member', 'ftime'), na.rm = TRUE) # multivariate RMSE (weighted) - rmse <- RMS(AvgExp, AvgObs, memb_dim = 'dataset', time_dim = 'sdate', + rmse <- RMS(AvgExp, AvgObs, dat_dim = 'dataset', time_dim = 'sdate', conf = FALSE)$rms stdev <- sd(AvgObs) mvrmse <- mvrmse + (rmse / stdev * as.numeric(weight[j])) diff --git a/tests/testthat/test-CST_MultiMetric.R b/tests/testthat/test-CST_MultiMetric.R index acbddfdc..9d048ab3 100644 --- a/tests/testthat/test-CST_MultiMetric.R +++ b/tests/testthat/test-CST_MultiMetric.R @@ -12,18 +12,15 @@ test_that("basic use case", { attr(obs, 'class') <- 's2dv_cube' result <- list(data = list(corr = array(rep(1, 3* 48), - #dim = c(nexp = 3, nobs = 1, - dim = c(n_exp = 3, n_obs = 1, + dim = c(nexp = 3, nobs = 1, lat = 6, lon = 8)), - p.val = array(rep(0, 3 * 48), dim = c(n_exp = 3, n_obs = 1, + p.val = array(rep(0, 3 * 48), dim = c(nexp = 3, nobs = 1, lat = 6, lon = 8)), conf.lower = array(rep(1, 3* 48), - #dim = c(nexp = 3, nobs = 1, - dim = c(n_exp = 3, n_obs = 1, + dim = c(nexp = 3, nobs = 1, lat = 6, lon = 8)), conf.upper = array(rep(1, 3* 48), - #dim = c(nexp = 3, nobs = 1, - dim = c(n_exp = 3, n_obs = 1, + dim = c(nexp = 3, nobs = 1, lat = 6, lon = 8))), lat = lat, lon = lon) attr(result, 'class') <- 's2dv_cube' @@ -34,11 +31,11 @@ test_that("basic use case", { res <- CST_MultiMetric(exp = exp, obs = obs, metric = 'rms') expect_equal(length(res), 3) expect_equal(dim(res$data$rms), - c(n_exp = 3, n_obs = 1, lat = 6, lon = 8)) + c(nexp = 3, nobs = 1, lat = 6, lon = 8)) res <- CST_MultiMetric(exp = exp, obs = obs, metric = 'rms', multimodel = FALSE) expect_equal(dim(res$data$rms), - c(n_exp = 2, n_obs = 1, lat = 6, lon = 8)) + c(nexp = 2, nobs = 1, lat = 6, lon = 8)) expect_equal(length(res$data), 3) res <- CST_MultiMetric(exp = exp, obs = obs, metric = 'rmsss') expect_equal(dim(res$data$rmsss), -- GitLab From 1ef7b93c86b66ecdbf21518763805a7b7b55def5 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 30 Nov 2020 16:27:50 +0100 Subject: [PATCH 07/10] fix dimnames --- R/CST_MultiMetric.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/CST_MultiMetric.R b/R/CST_MultiMetric.R index a932b844..3781db6c 100644 --- a/R/CST_MultiMetric.R +++ b/R/CST_MultiMetric.R @@ -179,7 +179,7 @@ MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE, pos_memb <- which(names(dim(AvgExp)) == memb_dim) pos_sdate <- which(names(dim(AvgExp)) == sdate_dim) corr <- Apply(list(AvgExp, AvgObs), - target_dims = list(c('sdate', 'lat', 'lon', 'member'), c('sdate', 'lat', 'lon')), + target_dims = list(c(sdate_dim, 'lat', 'lon', memb_dim), c(sdate_dim, 'lat', 'lon')), fun = function(x, y) { veriApply('FairRpss', fcst = x, obs = y, ensdim = which(names(dim(x)) == 'member'), -- GitLab From cd8d03fede359587d64522971ecaddd467018b29 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 30 Nov 2020 17:34:28 +0100 Subject: [PATCH 08/10] Fix vignette --- R/CST_MultiMetric.R | 2 +- vignettes/MultiModelSkill_vignette.Rmd | 56 ++++++++++++-------------- 2 files changed, 26 insertions(+), 32 deletions(-) diff --git a/R/CST_MultiMetric.R b/R/CST_MultiMetric.R index 3781db6c..d25ed5e4 100644 --- a/R/CST_MultiMetric.R +++ b/R/CST_MultiMetric.R @@ -12,7 +12,7 @@ #'@param time_dim name of the temporal dimension where a mean will be applied. It can be NULL, the default value is 'ftime'. #'@param memb_dim name of the member dimension. It can be NULL, the default value is 'member'. #'@param sdate_dim name of the start date dimension or a dimension name identifiying the different forecast. It can be NULL, the default value is 'sdate'. -#'@return an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest position in the first dimension correspons to the Multi-Model Mean. +#'@return an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the first position in the first 'nexp' dimension correspons to the Multi-Model Mean. #'@seealso \code{\link[s2dv]{Corr}}, \code{\link[s2dv]{RMS}}, \code{\link[s2dv]{RMSSS}} and \code{\link{CST_Load}} #'@references #'Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31.\url{http://link.springer.com/10.1007/s00382-018-4404-z} diff --git a/vignettes/MultiModelSkill_vignette.Rmd b/vignettes/MultiModelSkill_vignette.Rmd index 734652fe..e44be852 100644 --- a/vignettes/MultiModelSkill_vignette.Rmd +++ b/vignettes/MultiModelSkill_vignette.Rmd @@ -145,9 +145,12 @@ While other relevant data is being stored in the corresponding element of the ob ```r -> dim(AnomDJF$data) - dataset dataset statistics lat lon - 4 1 4 35 43 +> str(AnomDJF$data) +List of 4 + $ corr : num [1:4, 1, 1:35, 1:64] 0.586 0.614 0.143 0.501 0.419 ... + $ p.val : num [1:4, 1, 1:35, 1:64] 0.0026 0.00153 0.26805 0.01036 0.02931 ... + $ conf.lower: num [1:4, 1, 1:35, 1:64] 0.2073 0.2485 -0.3076 0.0883 -0.0154 ... + $ conf.upper: num [1:4, 1, 1:35, 1:64] 0.812 0.827 0.541 0.767 0.72 ... > names(AnomDJF) [1] "data" "lon" "lat" "Variable" "Datasets" "Dates" [7] "when" "source_files" "load_parameters" @@ -155,26 +158,19 @@ While other relevant data is being stored in the corresponding element of the ob [1] "glosea5" "ecmwf/system4_m1" "meteofrance/system5_m1" "erainterim" ``` -In the element $data of the `AnomDJF` object, the third dimension contains the lower limit of the 95% confidence interval, the correlation, the upper limit of the 95% confidence interval and the 95% significance level given by a one-sided T-test. - - -```r -corre <- AnomDJF$data[ , , 2, , ] -names(dim(corre)) <- c("maps", "lat", "lon") -``` - +In the element $data of the `AnomDJF` object is a list of object for the metric and its statistics: correlation, p-value, the lower limit of the 95% confidence interval and the upper limit of the 95% confidence interval and the 95% significance level given by a one-sided T-test. To obtain a spatial plot with a scale from -1 to 1 value of correlation for the model with the highest correlation for each grid point, the following lines should be run: ```r -PlotCombinedMap(corre, lon = Lon, lat = Lat, map_select_fun = max, - display_range = c(0, 1), map_dim = 'maps', +PlotCombinedMap(AnomDJF$data$corr[,1,,], lon = Lon, lat = Lat, map_select_fun = max, + display_range = c(0, 1), map_dim = 'nexp', legend_scale = 0.5, brks = 11, - cols = list(c('white', 'darkblue'), + cols = list(c('white', 'black'), + c('white', 'darkblue'), c('white', 'darkred'), - c('white', 'darkorange'), - c('white', 'black')), - bar_titles = c(names(AnomDJF$Datasets)[-4], "MMM"), + c('white', 'darkorange')), + bar_titles = c("MMM", names(AnomDJF$Datasets)), fileout = "./vignettes/Figures/MultiModelSkill_cor_tas_1992-2012.png", width = 14, height = 8) ``` @@ -200,14 +196,14 @@ The following lines are necessary to obtain the plot which visualizes the best m ```r names(dim(RMS)) <- c("maps", "lat", "lon") -PlotCombinedMap(RMS, lon = Lon, lat = Lat, map_select_fun = min, - display_range = c(0, ceiling(max(abs(RMS)))), map_dim = 'maps', +PlotCombinedMap(AnomDJF$data$rms[,1,,], lon = Lon, lat = Lat, map_select_fun = min, + display_range = c(0, ceiling(max(abs(AnomDJF$data$rms)))), map_dim = 'nexp', legend_scale = 0.5, brks = 11, - cols = list(c('darkblue', 'white'), + cols = list(c('black', 'white'), + c('darkblue', 'white'), c('darkred', 'white'), - c('darkorange', 'white'), - c('black', 'white')), - bar_titles = c(names(AnomDJF$Datasets)[-4], "MMM"), + c('darkorange', 'white')), + bar_titles = c("MMM", names(AnomDJF$Datasets)), fileout = "./vignettes/Figures/MultiModelSkill_rms_tas_1992-2012.png", width = 14, height = 8) ``` @@ -226,19 +222,17 @@ Notice that the perfect RMSSS is 1 and the parameter `map_select_fun` from `Plo ```r AnomDJF <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'rmsss', multimodel = TRUE) -RMSSS <- AnomDJF$data[ , , 1, , ] -names(dim(RMSSS)) <- c("maps", "lat", "lon") -PlotCombinedMap(RMSSS, lon = Lon, lat = Lat, +PlotCombinedMap(AnomDJF$data$rmsss[,1,,], lon = Lon, lat = Lat, map_select_fun = function(x) {x[which.min(abs(x - 1))]}, display_range = c(0, - ceiling(max(abs(RMSSS)))), map_dim = 'maps', + ceiling(max(abs(AnomDJF$data$rmsss)))), map_dim = 'nexp', legend_scale = 0.5, brks = 11, - cols = list(c('white', 'darkblue'), + cols = list(c('white', 'black'), + c('white', 'darkblue'), c('white', 'darkred'), - c('white', 'darkorange'), - c('white', 'black')), - bar_titles = c(names(AnomDJF$Datasets)[-4], "MMM"), + c('white', 'darkorange')), + bar_titles = c("MMM", names(AnomDJF$Datasets)), fileout = "./vignettes/Figures/MultiModelSkill_rmsss_tas_1992-2012.png", width = 14, height = 8) ``` -- GitLab From 58beb349513ee20d913d24a7a86dd747ab003fad Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 30 Nov 2020 17:39:24 +0100 Subject: [PATCH 09/10] package specification --- R/CST_MultivarRMSE.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/CST_MultivarRMSE.R b/R/CST_MultivarRMSE.R index 59d0a759..70c88c09 100644 --- a/R/CST_MultivarRMSE.R +++ b/R/CST_MultivarRMSE.R @@ -111,7 +111,7 @@ CST_MultivarRMSE <- function(exp, obs, weight = NULL) { AvgExp <- MeanDims(exp[[j]]$data, c('member', 'ftime'), na.rm = TRUE) AvgObs <- MeanDims(obs[[j]]$data, c('member', 'ftime'), na.rm = TRUE) # multivariate RMSE (weighted) - rmse <- RMS(AvgExp, AvgObs, dat_dim = 'dataset', time_dim = 'sdate', + rmse <- s2dv::RMS(AvgExp, AvgObs, dat_dim = 'dataset', time_dim = 'sdate', conf = FALSE)$rms stdev <- sd(AvgObs) mvrmse <- mvrmse + (rmse / stdev * as.numeric(weight[j])) -- GitLab From ada1b9fb02d51ef5f2bf8174970b6bd8289176a1 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 30 Nov 2020 17:48:32 +0100 Subject: [PATCH 10/10] NEWs updated --- NEWS.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/NEWS.md b/NEWS.md index e57b8358..8889b807 100644 --- a/NEWS.md +++ b/NEWS.md @@ -6,11 +6,13 @@ + CST_RFWeights accepts s2dv_cube objects as input and new 'ncores' paramenter allows parallel computation + RFWeights is exposed to users + CST_RainFARM accepts multi-dimensional slopes and weights and handless missing values in sample dimensions. + + CST_MultiMetric includes 'rpss' metric and it is addapted to s2dv. - Fixes: + PlotForecastPDF correctly displays terciles labels + CST_SaveExp correctly save time units + CST_SplitDims returns ordered output following ascending order provided in indices when it is numeric + ### CSTools 3.1.0 **Submission date to CRAN: 02-07-2020** -- GitLab