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' 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}}
#'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
#'@import easyVerification
#'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')
#'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)
}
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.")
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
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.")
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
161
162
163
164
165
166
167
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))
}
corr <- s2dv::Corr(AvgExp, AvgObs, memb_dim = dataset_dim, time_dim = sdate_dim)
corr <- s2dv::RMS(AvgExp, AvgObs, memb_dim = dataset_dim, time_dim = sdate_dim)
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'.")