Newer
Older
#'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.
#'
#'@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', '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'.
#'@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 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}}
#'Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31.\url{https://link.springer.com/10.1007/s00382-018-4404-z}
#'@importFrom s2dv MeanDims Reorder Corr RMS RMSSS InsertDim
#'@import abind
#'@importFrom easyVerification climFairRpss veriApply
#'mod <- 1 : (2 * 3 * 4 * 5 * 6 * 7)
#'dim(mod) <- c(dataset = 2, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
#'obs <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
#'dim(obs) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
#'lon <- seq(0, 30, 5)
#'lat <- seq(0, 25, 5)
#'exp <- list(data = mod, lat = lat, lon = lon)
#'obs <- list(data = obs, lat = lat, lon = lon)
#'attr(exp, 'class') <- 's2dv_cube'
#'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)
#'exp <- lonlat_data$exp
#'obs <- lonlat_data$obs
#'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')
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)
}
#'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{https://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)))) {
if (all(names(dim(exp)) %in% names(dim(obs)))) {
dimnames <- names(dim(exp))
stop("Dimension names of element 'data' from parameters 'exp'",
" and 'obs' should have the same name dimmension.")
stop("Element 'data' from parameters 'exp' and 'obs'",
" should have dimmension names.")
}
if (!is.logical(multimodel)) {
stop("Parameter 'multimodel' must be a logical value.")
}
if (length(multimodel) > 1) {
multimodel <- multimodel[1]
warning("Parameter 'multimodel' has length > 1 and only the first ",
"element will be used.")
}
if (length(metric) > 1) {
metric <- metric[1]
warning("Parameter 'multimodel' has length > 1 and only the first ",
"element will be used.")
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
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.")
}
stop("Element 'data' from parameters 'exp' and 'obs'",
" should have dimmension names.")
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)
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)
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'),
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)
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 <- 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))
}
corr <- s2dv::Corr(AvgExp, AvgObs, dat_dim = dataset_dim, time_dim = sdate_dim)
corr <- s2dv::RMS(AvgExp, AvgObs, dat_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 ",
"one of the options: 'correlation', 'rms', 'rmsss' or 'rpss'.")