CST_MultiMetric.R

#'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}}
#'@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
#'@import multiApply
#'@examples
#'library(zeallot)
#'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)
#'str(a)
#'\donttest{
#'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')
#'}
#'@export
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{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
#'@import multiApply
#'@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))
    } else {
      stop("Dimension names of element 'data' from parameters 'exp'",
           " and 'obs' should have the same name dimmension.")
    }
  } else {
    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.")
  }
  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("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)
    corr <- Apply(list(AvgExp, AvgObs),
       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))
    }
    if (metric == 'correlation') {
      corr <- s2dv::Corr(AvgExp, AvgObs, dat_dim = dataset_dim, time_dim = sdate_dim)
    } else if (metric == 'rms') {
      corr <- s2dv::RMS(AvgExp, AvgObs, dat_dim = dataset_dim, time_dim = sdate_dim)
    } else if (metric == 'rmsss') {
      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'.")
  }
  
  return(corr)
}