From 3504f22e9ec5bed1c56f32058b4e53077643f518 Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 3 Apr 2019 19:31:33 +0200 Subject: [PATCH 1/2] New anomaly function --- NAMESPACE | 1 + R/CST_Anomaly.R | 188 ++++++++++++++++++++++++++++++++++++++++ R/CST_MultiMetric.R | 13 +-- R/CST_MultivarRMSE.R | 22 ++--- man/CST_Anomaly.Rd | 64 ++++++++++++++ man/CST_MultiMetric.Rd | 8 +- man/CST_MultivarRMSE.Rd | 16 ++-- 7 files changed, 287 insertions(+), 25 deletions(-) create mode 100644 R/CST_Anomaly.R create mode 100644 man/CST_Anomaly.Rd diff --git a/NAMESPACE b/NAMESPACE index a2a12c10..872506f7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,6 @@ # Generated by roxygen2: do not edit by hand +export(CST_Anomaly) export(CST_BiasCorrection) export(CST_Calibration) export(CST_Load) diff --git a/R/CST_Anomaly.R b/R/CST_Anomaly.R new file mode 100644 index 00000000..e76c46ab --- /dev/null +++ b/R/CST_Anomaly.R @@ -0,0 +1,188 @@ +#'Anomalies relative to a climatology along selected dimension with or without cross-validation +#' +#'@author Perez-Zanon Nuria, \email{nuria.perez@bsc.es} +#'@author Pena Jesus, \email{jesus.pena@bsc.es} +#'@description This function computes the anomalies relative to a climatology computed along the +#'selected dimension (usually starting dates or forecast time) allowing the application or not of +#'crossvalidated climatologies. The computation is carried out independently for experimental and +#'observational data products. +#' +#'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}. +#'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed data in the element named \code{$data}.' +#'@param cross A logical value indicating whether cross-validation should be applied or not. Default = FALSE. +#'@param memb A logical value indicating whether Clim() computes one climatology for each experimental data +#'product member(TRUE) or it computes one sole climatology for all members (FALSE). Default = TRUE. +#' +#'@param dim_anom An integer indicating the dimension along which the climatology will be computed. It +#'usually corresponds to 3 (sdates) or 4 (ftime). Default = 3. +#' +#' @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 +#' +#'@examples +#'# Example 1: +#'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' +#' +#'anom1 <- CST_Anomaly(exp = exp, obs = obs, cross = FALSE, memb = TRUE) +#'str(anom1) +#'anom2 <- CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE) +#'str(anom2) +#' +#'anom3 <- CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = FALSE) +#'str(anom3) +#' +#'anom4 <- CST_Anomaly(exp = exp, obs = obs, cross = FALSE, memb = FALSE) +#'str(anom4) +#' +#'@seealso \code{\link[s2dverification]{Ano_CrossValid}}, \code{\link[s2dverification]{Clim}} and \code{\link{CST_Load}} +#' +#' +#'@export +CST_Anomaly <- function(exp = NULL, obs = NULL, cross = FALSE, memb = TRUE, dim_anom = 3) { + + 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.") + } + + if (!is.null(obs) & 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.") + } + + case_exp = case_obs = 0 + if (is.null(exp) & is.null(obs)) { + stop("One of the parameter 'exp' or 'obs' cannot be NULL.") + } + if (is.null(exp)) { + exp <- obs + case_obs = 1 + warning("Parameter 'exp' is not provided and will be recycled.") + } + if (is.null(obs)) { + obs <- exp + case_exp = 1 + warning("Parameter 'obs' is not provided and will be recycled.") + } + + 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)) + } 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.") + } + dim_exp <- dim(exp$data) + dim_obs <- dim(obs$data) + dimnames_data <- names(dim_exp) + + if (dim_exp[dim_anom] == 1 | dim_obs[dim_anom] == 1) { + stop("The length of dimension 'dim_anom' in label 'data' of the parameter ", + "'exp' and 'obs' must be greater than 1.") + } + if (!any(names(dim_exp)[dim_anom] == c('sdate', 'time', 'ftime'))) { + warning("Parameter 'dim_anom' correspond to a position name different ", + "than 'sdate', 'time' or 'ftime'.") + } + if (!is.logical(cross) | !is.logical(memb) ) { + stop("Parameters 'cross' and 'memb' must be logical.") + } + if (length(cross) > 1 | length(memb) > 1 ) { + cross <- cross[1] + warning("Parameter 'cross' has length greater than 1 and only the first element", + "will be used.") + } + if (length(memb) > 1) { + memb <- memb[1] + warning("Parameter 'memb' has length greater than 1 and only the first element", + "will be used.") + } + + # Selecting time dimension through dimensions permutation + if (dim_anom != 3) { + dimperm <- 1 : length(dim_exp) + dimperm[3] <- dim_anom + dimperm[dim_anom] <- 3 + + var_exp <- aperm(exp$data, perm = dimperm) + var_obs <- aperm(obs$data, perm = dimperm) + + #Updating permuted dimensions + dim_exp <- dim(exp$data) + dim_obs <- dim(obs$data) + } + + + # Computating anomalies + #---------------------- + + # With cross-validation + if (cross) { + ano <- Ano_CrossValid(var_exp = exp$data, var_obs = obs$data, memb = memb) + + # Without cross-validation + } else { + tmp <- Clim(var_exp = exp$data, var_obs = obs$data, memb = memb) + + if (memb) { + 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 <- InsertDim(clim_exp, 3, dim_exp[3]) + clim_obs <- InsertDim(clim_obs, 3, dim_obs[3]) + ano <- NULL + ano$ano_exp <- exp$data - clim_exp + ano$ano_obs <- obs$data - clim_obs + } + + # Permuting back dimensions to original order + if (dim_anom != 3) { + + if (case_obs == 0) { + ano$ano_exp <- aperm(ano$ano_exp, perm = dimperm) + } + if (case_exp == 0) { + ano$ano_obs <- aperm(ano$ano_obs, perm = dimperm) + } + + #Updating back permuted dimensions + dim_exp <- dim(exp$data) + dim_obs <- dim(obs$data) + } + + # Adding dimensions names + attr(ano$ano_exp, 'dimensions') <- dimnames_data + attr(ano$ano_obs, 'dimensions') <- dimnames_data + exp$data <- ano$ano_exp + obs$data <- ano$ano_obs + + # Outputs + # ~~~~~~~~~ + if (case_obs == 1) { + return(obs) + } + else if (case_exp == 1) { + return(exp) + } + else { + return(list(exp = exp, obs = obs)) + } +} diff --git a/R/CST_MultiMetric.R b/R/CST_MultiMetric.R index 596b6152..e7bfc729 100644 --- a/R/CST_MultiMetric.R +++ b/R/CST_MultiMetric.R @@ -4,8 +4,8 @@ #'@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_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}. -#'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed 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 experiment 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. #' @@ -17,6 +17,7 @@ #'@import s2dverification #'@import stats #'@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) @@ -27,7 +28,8 @@ #'obs <- list(data = obs, lat = lat, lon = lon) #'attr(exp, 'class') <- 's2dv_cube' #'attr(obs, 'class') <- 's2dv_cube' -#'a <- CST_MultiMetric(exp = exp, obs = obs) +#'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) #'@export CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE) { @@ -66,11 +68,10 @@ CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE) warning("Parameter 'multimodel' has length > 1 and only the first ", "element will be used.") } - ano <- Ano_CrossValid(var_exp = exp$data, var_obs = obs$data) # seasonal average of anomalies per model - AvgExp <- MeanListDim(ano$ano_exp, narm = T, c(2, 4)) - AvgObs <- MeanListDim(ano$ano_obs, narm = T, c(2, 4)) + AvgExp <- MeanListDim(exp$data, narm = T, c(2, 4)) + AvgObs <- MeanListDim(obs$data, narm = T, c(2, 4)) # indv model correlation if (metric == 'correlation') { diff --git a/R/CST_MultivarRMSE.R b/R/CST_MultivarRMSE.R index d1fc4f29..7c623a0f 100644 --- a/R/CST_MultivarRMSE.R +++ b/R/CST_MultivarRMSE.R @@ -3,8 +3,8 @@ #'@author Deborah Verfaillie, \email{deborah.verfaillie@bsc.es} #'@description This function calculates the RMSE from multiple variables, as the mean of each variable's RMSE scaled by its observed standard deviation. Variables can be weighted based on their relative importance (defined by the user). #' -#'@param exp a list of objects, one for each variable, of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}. -#'@param obs a list of objects, one for each variable (in the same order than the input in 'exp') of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed data in the element named \code{$data}. +#'@param exp a list of objects, one for each variable, 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 obs a list of objects, one for each variable (in the same order than the input in 'exp') of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the observed anomaly data in the element named \code{$data}. #' #'@param weight (optional) a vector of weight values to assign to each variable. If no weights are defined, a value of 1 is assigned to every variable. #' @@ -15,7 +15,8 @@ #'@examples #'# Creation of sample s2dverification objects. These are not complete #'# s2dverification objects though. The Load function returns complete objects. -#' +#'# using package zeallot is optional: +#' library(zeallot) #'# Example with 2 variables #'mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7) #'mod2 <- 1 : (1 * 3 * 4 * 5 * 6 * 7) @@ -33,17 +34,19 @@ #'exp2 <- list(data = mod2, lat = lat, lon = lon, Datasets = "EXP2", #' source_files = "file2", Variable = list('tas')) #'attr(exp2, 'class') <- 's2dv_cube' -#'exp <- list(exp1, exp2) #'obs1 <- list(data = obs1, lat = lat, lon = lon, Datasets = "OBS1", #' source_files = "file1", Variable = list('pre')) #'attr(obs1, 'class') <- 's2dv_cube' #'obs2 <- list(data = obs2, lat = lat, lon = lon, Datasets = "OBS2", #' source_files = "file2", Variable = list('tas')) #'attr(obs2, 'class') <- 's2dv_cube' - -#'obs <- list(obs1, obs2) +#' +#'c(ano_exp1, ano_obs1) %<-% CST_Anomaly(exp1, obs1, cross = TRUE, memb = TRUE) +#'c(ano_exp2, ano_obs2) %<-% CST_Anomaly(exp2, obs2, cross = TRUE, memb = TRUE) +#'ano_exp <- list(exp1, exp2) +#'ano_obs <- list(ano_obs1, ano_obs2) #'weight <- c(1, 2) -#'a <- CST_MultivarRMSE(exp = exp, obs = obs, weight = weight) +#'a <- CST_MultivarRMSE(exp = ano_exp, obs = ano_obs, weight = weight) #'str(a) #'@export CST_MultivarRMSE <- function(exp, obs, weight = NULL) { @@ -104,10 +107,9 @@ CST_MultivarRMSE <- function(exp, obs, weight = NULL) { mvrmse <- 0 sumweights <- 0 for (j in 1 : nvar) { - ano <- Ano_CrossValid(exp[[j]]$data, obs[[j]]$data) # seasonal average of anomalies - AvgExp <- MeanListDim(ano$ano_exp, narm = TRUE, c(2, 4)) - AvgObs <- MeanListDim(ano$ano_obs, narm = TRUE, c(2, 4)) + AvgExp <- MeanListDim(exp[[j]]$data, narm = TRUE, c(2, 4)) + AvgObs <- MeanListDim(obs[[j]]$data, narm = TRUE, c(2, 4)) # multivariate RMSE (weighted) rmse <- s2dverification::RMS(AvgExp, AvgObs, posloop = 1, posRMS = 2, conf = FALSE) diff --git a/man/CST_Anomaly.Rd b/man/CST_Anomaly.Rd new file mode 100644 index 00000000..b214a31a --- /dev/null +++ b/man/CST_Anomaly.Rd @@ -0,0 +1,64 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_Anomaly.R +\name{CST_Anomaly} +\alias{CST_Anomaly} +\title{Anomalies relative to a climatology along selected dimension with or without cross-validation} +\usage{ +CST_Anomaly(exp = NULL, obs = NULL, cross = FALSE, memb = TRUE, + dim_anom = 3) +} +\arguments{ +\item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}.} + +\item{obs}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed data in the element named \code{$data}.'} + +\item{cross}{A logical value indicating whether cross-validation should be applied or not. Default = FALSE.} + +\item{memb}{A logical value indicating whether Clim() computes one climatology for each experimental data +product member(TRUE) or it computes one sole climatology for all members (FALSE). Default = TRUE.} + +\item{dim_anom}{An integer indicating the dimension along which the climatology will be computed. It +usually corresponds to 3 (sdates) or 4 (ftime). Default = 3.} +} +\value{ +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. +} +\description{ +This function computes the anomalies relative to a climatology computed along the +selected dimension (usually starting dates or forecast time) allowing the application or not of +crossvalidated climatologies. The computation is carried out independently for experimental and +observational data products. +} +\examples{ +# Example 1: +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' + +anom1 <- CST_Anomaly(exp = exp, obs = obs, cross = FALSE, memb = TRUE) +str(anom1) +anom2 <- CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE) +str(anom2) + +anom3 <- CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = FALSE) +str(anom3) + +anom4 <- CST_Anomaly(exp = exp, obs = obs, cross = FALSE, memb = FALSE) +str(anom4) + +} +\seealso{ +\code{\link[s2dverification]{Ano_CrossValid}}, \code{\link[s2dverification]{Clim}} and \code{\link{CST_Load}} +} +\author{ +Perez-Zanon Nuria, \email{nuria.perez@bsc.es} + +Pena Jesus, \email{jesus.pena@bsc.es} +} diff --git a/man/CST_MultiMetric.Rd b/man/CST_MultiMetric.Rd index 0f5f2acf..99fdb400 100644 --- a/man/CST_MultiMetric.Rd +++ b/man/CST_MultiMetric.Rd @@ -7,9 +7,9 @@ CST_MultiMetric(exp, obs, metric = "correlation", multimodel = TRUE) } \arguments{ -\item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing 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 experiment data in the element named \code{$data}.} -\item{obs}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed 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.} @@ -22,6 +22,7 @@ an object of class \code{s2dv_cube} containing the statistics of the selected me 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. } \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) @@ -32,7 +33,8 @@ 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' -a <- CST_MultiMetric(exp = exp, obs = obs) +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) } \references{ diff --git a/man/CST_MultivarRMSE.Rd b/man/CST_MultivarRMSE.Rd index 5db63421..b81d0a49 100644 --- a/man/CST_MultivarRMSE.Rd +++ b/man/CST_MultivarRMSE.Rd @@ -7,9 +7,9 @@ CST_MultivarRMSE(exp, obs, weight = NULL) } \arguments{ -\item{exp}{a list of objects, one for each variable, of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}.} +\item{exp}{a list of objects, one for each variable, 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{obs}{a list of objects, one for each variable (in the same order than the input in 'exp') of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed data in the element named \code{$data}.} +\item{obs}{a list of objects, one for each variable (in the same order than the input in 'exp') of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the observed anomaly data in the element named \code{$data}.} \item{weight}{(optional) a vector of weight values to assign to each variable. If no weights are defined, a value of 1 is assigned to every variable.} } @@ -22,7 +22,8 @@ This function calculates the RMSE from multiple variables, as the mean of each v \examples{ # Creation of sample s2dverification objects. These are not complete # s2dverification objects though. The Load function returns complete objects. - +# using package zeallot is optional: +library(zeallot) # Example with 2 variables mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7) mod2 <- 1 : (1 * 3 * 4 * 5 * 6 * 7) @@ -40,16 +41,19 @@ attr(exp1, 'class') <- 's2dv_cube' exp2 <- list(data = mod2, lat = lat, lon = lon, Datasets = "EXP2", source_files = "file2", Variable = list('tas')) attr(exp2, 'class') <- 's2dv_cube' -exp <- list(exp1, exp2) obs1 <- list(data = obs1, lat = lat, lon = lon, Datasets = "OBS1", source_files = "file1", Variable = list('pre')) attr(obs1, 'class') <- 's2dv_cube' obs2 <- list(data = obs2, lat = lat, lon = lon, Datasets = "OBS2", source_files = "file2", Variable = list('tas')) attr(obs2, 'class') <- 's2dv_cube' -obs <- list(obs1, obs2) + +c(ano_exp1, ano_obs1) \%<-\% CST_Anomaly(exp1, obs1, cross = TRUE, memb = TRUE) +c(ano_exp2, ano_obs2) \%<-\% CST_Anomaly(exp2, obs2, cross = TRUE, memb = TRUE) +ano_exp <- list(exp1, exp2) +ano_obs <- list(ano_obs1, ano_obs2) weight <- c(1, 2) -a <- CST_MultivarRMSE(exp = exp, obs = obs, weight = weight) +a <- CST_MultivarRMSE(exp = ano_exp, obs = ano_obs, weight = weight) str(a) } \seealso{ -- GitLab From c9e74765bb4d620325c0ab25903f6b9e9e32aca6 Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 4 Apr 2019 12:10:22 +0200 Subject: [PATCH 2/2] Vignettes updated with CST_Anomaly function --- R/CST_MultivarRMSE.R | 4 +-- man/CST_MultivarRMSE.Rd | 4 +-- vignettes/MultiModelSkill_vignette.md | 37 ++++++++++++++++++++++----- vignettes/MultivarRMSE_vignette.md | 34 ++++++++++++++++++------ 4 files changed, 60 insertions(+), 19 deletions(-) diff --git a/R/CST_MultivarRMSE.R b/R/CST_MultivarRMSE.R index 7c623a0f..88967a13 100644 --- a/R/CST_MultivarRMSE.R +++ b/R/CST_MultivarRMSE.R @@ -26,8 +26,8 @@ #'obs2 <- 1 : (1 * 1 * 4 * 5 * 6 * 7) #'dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7) #'dim(obs2) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7) -#'lon <- seq(0, 25, 5) -#'lat <- seq(0, 30, 5) +#'lon <- seq(0, 30, 5) +#'lat <- seq(0, 25, 5) #'exp1 <- list(data = mod1, lat = lat, lon = lon, Datasets = "EXP1", #' source_files = "file1", Variable = list('pre')) #'attr(exp1, 'class') <- 's2dv_cube' diff --git a/man/CST_MultivarRMSE.Rd b/man/CST_MultivarRMSE.Rd index b81d0a49..24af608c 100644 --- a/man/CST_MultivarRMSE.Rd +++ b/man/CST_MultivarRMSE.Rd @@ -33,8 +33,8 @@ obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7) obs2 <- 1 : (1 * 1 * 4 * 5 * 6 * 7) dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7) dim(obs2) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7) -lon <- seq(0, 25, 5) -lat <- seq(0, 30, 5) +lon <- seq(0, 30, 5) +lat <- seq(0, 25, 5) exp1 <- list(data = mod1, lat = lat, lon = lon, Datasets = "EXP1", source_files = "file1", Variable = list('pre')) attr(exp1, 'class') <- 's2dv_cube' diff --git a/vignettes/MultiModelSkill_vignette.md b/vignettes/MultiModelSkill_vignette.md index 6bebe138..6fc67610 100644 --- a/vignettes/MultiModelSkill_vignette.md +++ b/vignettes/MultiModelSkill_vignette.md @@ -71,16 +71,19 @@ Using the `CST_Load` function, the data available in our data store can be loade Ask nuria.perez at bsc.es to achieve the data to run the recipe. ```r -#glosea5 <- '/esnas/exp/glosea5/specs-seasonal_i1p1/$STORE_FREQ$_mean/$VAR_NAME$-allmemb/$VAR_NAME$_$START_DATE$.nc' +glosea5 <- '/esnas/exp/glosea5/specs-seasonal_i1p1/$STORE_FREQ$_mean/$VAR_NAME$-allmemb/$VAR_NAME$_$START_DATE$.nc' c(exp, obs) %<-% - CST_Load(var = clim_var, exp = list(list(name = 'glosea5', path = glosea5), list(name = 'ecmwf/system4_m1'), + CST_Load(var = clim_var, exp = list(list(name = 'glosea5', path = glosea5), + list(name = 'ecmwf/system4_m1'), list(name = 'meteofrance/system5_m1')), obs = obs, sdates = dateseq, leadtimemin = 2, leadtimemax = 4, lonmin = -20, lonmax = 70, latmin = 25, latmax = 75, storefreq = "monthly", sampleperiod = 1, nmember = 9, output = "lonlat", method = "bilinear", grid = paste("r", grid, sep = "")) -#save(exp, obs, file = "../tas_toydata.RData") +#save(exp, obs, file = "../tas_toydata.RData") + +# Or use the following line to load the file provided in .RData format: load(file = "./tas_toydata.RData") ``` @@ -108,11 +111,31 @@ Lon <- exp$lon ### 2.- Computing and plotting Anomaly Correlation Coefficient -The Anomaly Correlation Coefficient (ACC) is the most widely used skill metric for Seasonal Climate Forecast quality (Mishra et al., 2018). The ACC is obtained by running the `CST_MultiMetric` function defining the parameter 'metric' as correlation. Furthermore, the function includes the option of computing the Multi-Model Mean ensemble (MMM). +The Anomaly Correlation Coefficient (ACC) is the most widely used skill metric for Seasonal Climate Forecast quality (Mishra et al., 2018). + + +First step is to compute the anomalies over the loaded data applying cross validation technique on individual members by running: + +``` +c(ano_exp, ano_obs) %<-% CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE) +``` + +The dimensions are preserved: + +``` +> str(ano_exp$data) + num [1:3, 1:9, 1:21, 1:3, 1:35, 1:64] -1.647 -0.478 -0.096 1.575 1.086 ... + - attr(*, "dimensions")= chr [1:6] "dataset" "member" "sdate" "ftime" ... +> str(ano_obs$data) + num [1, 1, 1:21, 1:3, 1:35, 1:64] 0.0235 1.546 1.3885 -0.344 -5.972 ... + - attr(*, "dimensions")= chr [1:6] "dataset" "member" "sdate" "ftime" ... +``` + +The ACC is obtained by running the `CST_MultiMetric` function defining the parameter 'metric' as correlation. The function also includes the option of computing the Multi-Model Mean ensemble (MMM). ```r -AnomDJF <- CST_MultiMetric(exp = exp, obs = obs, metric = 'correlation', +AnomDJF <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'correlation', multimodel = TRUE) ``` @@ -168,7 +191,7 @@ The next figure is the map of the maximum positive Anomaly Correlation Coefficie The same function can be used to compute the RMS error by defining the parameter `metric` as 'rms'. ```r -AnomDJF <- CST_MultiMetric(exp = exp, obs = obs, metric = 'rms', +AnomDJF <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'rms', multimodel = TRUE) RMS <- AnomDJF$data[ , , 2, , ] ``` @@ -201,7 +224,7 @@ Notice that the perfect RMSSS is 1 and the parameter `map_select_fun` from `Plo ```r -AnomDJF <- CST_MultiMetric(exp = exp, obs = obs, metric = 'rmsss', +AnomDJF <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'rmsss', multimodel = TRUE) RMSSS <- AnomDJF$data[ , , 1, , ] names(dim(RMSSS)) <- c("maps", "lat", "lon") diff --git a/vignettes/MultivarRMSE_vignette.md b/vignettes/MultivarRMSE_vignette.md index fc18b99f..78af6080 100644 --- a/vignettes/MultivarRMSE_vignette.md +++ b/vignettes/MultivarRMSE_vignette.md @@ -17,7 +17,6 @@ To run this vignette, the next R packages should be installed and loaded: ```r library(s2dverification) -library(startR) library(RColorBrewer) library(zeallot) ``` @@ -70,14 +69,13 @@ Using the `CST_Load` function from **CSTool package**, the data available in our Ask nuria.perez at bsc.es for the data to run the recipe. ```r -#glosea5 <- list(path = '/esnas/exp/glosea5/specs-seasonal_i1p1/$STORE_FREQ$_mean/$VAR_NAME$-allmemb/$VAR_NAME$_$START_DATE$.nc') +glosea5 <- list(path = '/esnas/exp/glosea5/specs-seasonal_i1p1/$STORE_FREQ$_mean/$VAR_NAME$-allmemb/$VAR_NAME$_$START_DATE$.nc') c(exp_T, obs_T) %<-% CST_Load(var = temp, exp = list(glosea5), obs = obs, sdates = dateseq, leadtimemin = 2, leadtimemax = 4, latmin = 25, latmax = 75, lonmin = -20, lonmax = 70, output = 'lonlat', nprocs = 1, storefreq = "monthly", sampleperiod = 1, nmember = 9, method = "bilinear", grid = paste("r", grid, sep = "")) -#ecmwf <- list(path = '/esnas/exp/ecmwf/system4_m1/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc') c(exp_P, obs_P) %<-% CST_Load(var = precip, exp = list(glosea5), obs = obs, sdates = dateseq, leadtimemin = 2, leadtimemax = 4, @@ -85,7 +83,9 @@ c(exp_P, obs_P) %<-% nprocs = 1, storefreq = "monthly", sampleperiod = 1, nmember = 9, method = "bilinear", grid = paste("r", grid, sep = "")) #save(exp_T, obs_T, exp_P, obs_P, file = "./tas_prlr_toydata.RData") -load(file = "../tas_prlr_toydata.RData") + +# Or use the following line to load the file provided in .RData format: +load(file = "./tas_prlr_toydata.RData") ``` There should be four new elements loaded in the R working environment: `exp_T`, `obs_T`, `exp_P` and `obs_P`. The first two elements correspond to the experimental and observed data for temperature and the other are the equivalent for the precipitation data. It's possible to check that they are of class `sd2v_cube` by running: @@ -118,14 +118,32 @@ Lat <- exp_T$lat Lon <- exp_T$lon ``` +The next step is to compute the anomalies of the experimental and observational data using `CST_Anomaly` function, which could be applied over data from each variable, and in this case it's compute applying cross validation technique over individual members: + +``` +c(ano_exp_T, ano_obs_T) %<-% CST_Anomaly(exp = exp_T, obs = obs_T, cross = TRUE, memb = TRUE) +c(ano_exp_P, ano_obs_P) %<-% CST_Anomaly(exp = exp_P, obs = obs_P, cross = TRUE, memb = TRUE) +``` + +The original dimensions are preserved and the anomalies are stored in the `data` element of the correspondent object: + +``` +> str(ano_exp_T$data) + num [1, 1:9, 1:21, 1:3, 1:35, 1:64] -1.647 1.575 2.77 0.048 -1.886 ... + - attr(*, "dimensions")= chr [1:6] "dataset" "member" "sdate" "ftime" ... +> str(ano_obs_T$data) + num [1, 1, 1:21, 1:3, 1:35, 1:64] 0.0235 1.546 1.3885 -0.344 -5.972 ... + - attr(*, "dimensions")= chr [1:6] "dataset" "member" "sdate" "ftime" ... +``` + +Two lists containing the experiment ,`ano_exp`, and the observation, `ano_obs`, lists should be put together to serve as input of the function to compute multivariate RMSEs. -Two lists containing the `exp` and `obs` lists should be put together to serve as input of the function to compute multivariate RMSEs. Furthermore, some weights can be applied to the difference variables based on their relative importance (if no weights are given, a value of 1 is automatically assigned to each variable). For this example, we'll give a weight of 2 to the temperature dataset and a weight of 1 to the precipitation dataset: ```r -exp <- list(exp_T, exp_P) -obs <- list(obs_T, obs_P) +ano_exp <- list(ano_exp_T, ano_exp_P) +ano_obs <- list(ano_obs_T, ano_obs_P) weight <- c(2, 1) ``` @@ -136,7 +154,7 @@ It is obtained by running the `CST_MultivarRMSE` function: ```r -mvrmse <- CST_MultivarRMSE(exp = exp, obs = obs, weight) +mvrmse <- CST_MultivarRMSE(exp = ano_exp, obs = ano_obs, weight) ``` -- GitLab