diff --git a/DESCRIPTION b/DESCRIPTION index 3e9a394a2917c44f1b4f2dcc11c51667c392016a..59fad67a128540b694134a4d5ea64e9439aacabd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -14,7 +14,8 @@ Authors@R: c( person("Deborah", "Verfaillie", , "deborah.verfaillie@bsc.es", role = "aut"), person("Lauriane", "Batte", , "lauriane.batte@meteo.fr", role = "ctb"), person("Jesus", "Peña", , "jesus.pena@bsc.es", role = "ctb"), - person("Bert", "van Schaeybroeck", , "bertvs@meteo.be", role = "aut")) + person("Bert", "van Schaeybroeck", , "bertvs@meteo.be", role = "aut"), + person("Eroteida", "Sanchez-Garcia", , "esanchezg@aemet.es", role = "ctb")) Description: Exploits dynamical seasonal forecasts in order to provide information relevant to stakeholders at the seasonal timescale. The package contains process-based methods for forecast calibration, bias correction, @@ -44,7 +45,8 @@ Imports: graphics, grDevices, stats, - utils + utils, + verification Suggests: zeallot, testthat, diff --git a/NAMESPACE b/NAMESPACE index f57c5ac4083998f071d9b2d6f96da1c9df30d00e..06463e5291fae05467cfb94047c2d3cbf9aea5f3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,9 @@ # Generated by roxygen2: do not edit by hand +export(BEI_PDFBest) +export(BEI_Weights) export(CST_Anomaly) +export(CST_BEI_Weighting) export(CST_BiasCorrection) export(CST_Calibration) export(CST_Load) @@ -50,3 +53,4 @@ importFrom(plyr,.) importFrom(plyr,dlply) importFrom(reshape2,melt) importFrom(utils,glob2rx) +importFrom(verification,verify) diff --git a/R/BEI_PDFBest.R b/R/BEI_PDFBest.R new file mode 100644 index 0000000000000000000000000000000000000000..46833b2c57bc716f8459c90bfebb1fac9a189350 --- /dev/null +++ b/R/BEI_PDFBest.R @@ -0,0 +1,948 @@ +#' Computing the Best Index PDFs combining Index PDFs from two SFSs +#' +#'@author Eroteida Sanchez-Garcia - AEMET, \email{esanchezg@aemet.es} +#' +#'@description This function implements the computation to obtain the index +#'Probability Density Functions (PDFs) (e.g. NAO index) obtained to combining +#'the Index PDFs for two Seasonal Forecast Systems (SFSs), the Best Index +#'estimation (see Sanchez-Garcia, E. et al (2019), +#'https://doi.org/10.5194/asr-16-165-2019 for more details about the +#'methodology applied to estimate the Best Index). +#' +#'@references Regionally improved seasonal forecast of precipitation through +#'Best estimation of winter NAO, Sanchez-Garcia, E. et al., +#' Adv. Sci. Res., 16, 165174, 2019, https://doi.org/10.5194/asr-16-165-2019 +#' +#'@param index_obs Index (e.g. NAO index) array from an observational database +#' or reanalysis with at least a temporal dimension (by default 'time'), +#' which must be greater than 2. +#'@param index_hind1 Index (e.g. NAO index) array from a SFS (named SFS1) +#' with at least two dimensions (time , member) or (time, statistic). +#' The temporal dimension, by default 'time', must be greater than 2. +#' The dimension 'member' must be greater than 1. +#' The dimension 'statistic' must be equal to 2, for containing the two paramenters of +#' a normal distribution (mean and sd) representing the ensemble of a SFS. +#' It is not possible to have the dimension 'member' and 'statistic' at the same time. +#'@param index_hind2 Index (e.g. NAO index) array from a SFS (named SFS2) +#' with at least two dimensions (time , member) or (time, statistic). +#' The temporal dimension, by default 'time', must be greater than 2. +#' The dimension 'member' must be greater than 1. +#' The dimension 'statistic' must be equal to 2, for containing the two paramenters of +#' a normal distribution (mean and sd) representing the ensemble of a SFS. +#' It is not possible to have the dimension 'member' and 'statistic' together. +#'@param index_fcst1 (optional, default = NULL) Index (e.g. NAO index) array from forescating of SFS1 +#' with at least two dimensions (time , member) or (time, statistic). +#' The temporal dimension, by default 'time', must be equal to 1, the forecast year target. +#' The dimension 'member' must be greater than 1. +#' The dimension 'statistic' must be equal to 2, for containing the two paramenters of +#' a normal distribution (mean and sd) representing the ensemble of a SFS. +#' It is not possible to have the dimension 'member' and 'statistic' together. +#'@param index_fcst2 (optional, default = NULL) Index (e.g. NAO index) array from forescating of SFS2 +#' with at least two dimensions (time , member) or (time, statistic). +#' The temporal dimension, by default 'time', must be equal to 1, the forecast year target. +#' The dimension 'member' must be greater than 1. +#' The dimension 'statistic' must be equal to 2, for containing the two paramenters of +#' a normal distribution (mean and sd) representing the ensemble of a SFS. +#' It is not possible to have the dimension 'member' and 'statistic' together. +#'@param method_BC A character vector of maximun length 2 indicating the bias +#'correction methodology to be applied on each SFS. If it is 'none' or any of +#'its elements is 'none', the bias correction won't be applied. +#'Available methods developped are "ME" (a bias correction scheme based on the +#'mean error or bias between observation and predictions to correct the +#'predicted values), and "LMEV" (a bias correction scheme based on a linear +#'model using ensemble variance of index as predictor). (see Sanchez-Garcia, +#'E. et al (2019), https://doi.org/10.5194/asr-16-165-2019 for more details). +#'@param time_dim_name A character string indicating the name of the temporal +#'dimension, by default 'time'. +#'@param na.rm Logical (default = FALSE). Should missing values be removed? +#' +#' @return BEI_PDFBest() returns an array with the parameters that caracterize +#' the PDFs, with at least a temporal dimension, by default 'time' and dimension +#' 'statistic' equal to 2. +#' The firt statistic is the parameter 'mean' of the PDF for the best estimation +#' combining the two SFSs PDFs. +#' The second statistic is the parameter 'standard deviation' of the PDF for +#' the best estimation combining the two SFSs PDFs. +#' If index_fcst1 and/or index_fcst2 are null, returns the values for hindcast period. +#' Otherwise, it returns the values for a forecast year. +#'@import multiApply +#' +#'@examples +#' # Example 1 for the BEI_PDFBest function +#' index_obs<- rnorm(10, sd = 3) +#' dim(index_obs) <- c(time = 5, season = 2) +#' index_hind1 <- rnorm(40, mean = 0.2, sd = 3) +#' dim(index_hind1) <- c(time = 5, member = 4, season = 2) +#' index_hind2 <- rnorm(60, mean = -0.5, sd = 4) +#' dim(index_hind2) <- c(time = 5, member = 6, season = 2) +#' index_fcst1 <- rnorm(16, mean = 0.2, sd = 3) +#' dim(index_fcst1) <- c(time = 1, member = 8, season = 2) +#' index_fcst2 <- rnorm(18, mean = -0.5, sd = 4) +#' dim(index_fcst2) <- c(time = 1, member = 9, season = 2) +#' method_BC <- 'ME' +#' res <- BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, +#' index_fcst2, method_BC) +#' dim(res) +#' # time statistic season +#' # 1 2 2 +#' # Example 2 for the BEI_PDFBest function +#' index_obs<- rnorm(10, sd = 3) +#' dim(index_obs) <- c(time = 5, season = 2) +#' index_hind1 <- rnorm(40, mean = 0.2, sd = 3) +#' dim(index_hind1) <- c(time = 5, member = 4, season = 2) +#' index_hind2 <- rnorm(60, mean = -0.5, sd = 4) +#' dim(index_hind2) <- c(time = 5, member = 6, season = 2) +#' index_fcst1 <- rnorm(16, mean = 0.2, sd = 3) +#' dim(index_fcst1) <- c(time = 1, member = 8, season = 2) +#' index_fcst2 <- rnorm(18, mean = -0.5, sd = 4) +#' dim(index_fcst2) <- c(time = 1, member = 9, season = 2) +#' method_BC <- c('LMEV', 'ME') +#' res <- BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, index_fcst2, method_BC) +#' dim(res) +#' # time statistic season +#' # 1 2 2 +#'@export + +BEI_PDFBest <- function(index_obs, index_hind1, index_hind2, + index_fcst1 = NULL, index_fcst2 = NULL, method_BC = 'none', + time_dim_name = 'time', na.rm = FALSE) { + + if (!is.logical(na.rm)) { + stop("Parameter 'na.rm' must be a logical value.") + } + if (!is.character(time_dim_name)) { + stop("Parameter 'time_dim_name' must be a character string ", + "indicating the name of the temporal dimension.") + } + if (length(time_dim_name) > 1) { + warning("Parameter 'time_dim_name' has length greater than 1 and ", + "only the first element will be used.") + time_dim_name <- time_dim_name[1] + } + if (!is.character(method_BC) || !is.vector(method_BC)){ + stop("Parameter 'method_BC' must be a character vector.") + } + if (!(length(method_BC) == 1 || length(method_BC) == 2)){ + stop("Length of parameter 'method_BC' must be 1 or 2.") + } + if(!all(method_BC %in% c('ME', 'LMEV', 'none'))){ + stop("Elements of parameter 'method_BC' must be equals to ", + "'none, 'ME' or 'LMEV'") + } + if (!is.array(index_obs)) { + stop("Parameter 'index_obs' must be an array.") + } + if (!is.array(index_hind1)) { + stop("Parameter 'index_hind1' must be an array.") + } + if (!is.array(index_hind2)) { + stop("Parameter 'index_hind2' must be an array.") + } + if (is.null(names(dim(index_hind1))) || + is.null(names(dim(index_obs))) || + is.null(names(dim(index_hind2)))) { + stop("Parameters 'index_obs', 'index_hind1' and 'index_hind2' ", + "should have dimmension names.") + } + if(!(time_dim_name %in% names(dim(index_obs)))) { + stop("Parameter 'index_obs' must have temporal dimension.") + } + if(!(time_dim_name %in% names(dim(index_hind1)))) { + stop("Parameter 'index_hind1' must have temporal dimension.") + } + if(!(time_dim_name %in% names(dim(index_hind2)))) { + stop("Parameter 'index_hind2' must have temporal dimension.") + } + if (dim(index_obs)[time_dim_name] <= 2) { + stop("Length of temporal dimension ", + "of parameter 'index_obs', 'index_hind1' and 'index_hind2' ", + "must be greater than 2.") + } + if (dim(index_obs)[time_dim_name] != dim(index_hind1)[time_dim_name]) { + stop("Length of temporal dimensions ", + "of parameter 'index_obs' and 'index_hind1' must be equals.") + } + if (dim(index_hind1)[time_dim_name] != dim(index_hind2)[time_dim_name]) { + stop("Length of temporal dimensions ", + "of parameter 'index_hind1' and 'index_hind2' must be equals.") + } + if('member' %in% names(dim(index_hind1)) & + 'statistic' %in% names(dim(index_hind1))) { + stop("Parameter 'index_hind1' must have at least ", + "dimension 'member' or 'statistic', not 'member' and 'statistic' ", + "together.") + } + if('member' %in% names(dim(index_hind2)) & + 'statistic' %in% names(dim(index_hind2))) { + stop("Parameter 'index_hind2' must have at least ", + "dimension 'member' or 'statistic', not 'member' and 'statistic' ", + "together.") + } + if(!('member' %in% names(dim(index_hind1))) & + !('statistic' %in% names(dim(index_hind1)))) { + stop("Parameter 'index_hind1' must have dimension ", + "'member' or 'statistic'") + } + if(!('member' %in% names(dim(index_hind2))) & + !('statistic' %in% names(dim(index_hind2)))) { + stop("Parameter 'index_hind2' must have dimension ", + "'member' or 'statistic'") + } + if ('member' %in% names(dim(index_hind1))){ + if (dim(index_hind1)['member'] == 1) { + stop("Length of dimension 'member' ", + "of parameter 'index_hind1' must be greater than 1.") + } + } + if ('member' %in% names(dim(index_hind2))){ + if (dim(index_hind2)['member'] == 1) { + stop("Length of dimension 'member' ", + "of parameter 'index_hind2' must be greater than 1.") + } + } + if ('statistic' %in% names(dim(index_hind1))){ + if (dim(index_hind1)['statistic'] != 2) { + stop("Length of dimension 'statistic' ", + "of parameter 'index_hind1' must be equal to 2.") + } + } + if ('statistic' %in% names(dim(index_hind2))){ + if (dim(index_hind2)['statistic'] != 2) { + stop("Length of dimension 'statistic' ", + "of parameter 'index_hind2' must be equal to 2.") + } + } + if (!is.null(index_fcst1)){ + if (!is.array(index_fcst1)) { + stop("Parameter 'index_fcst1' must be an array.") + } + if (is.null(names(dim(index_fcst1)))){ + stop("Parameters 'index_fcst1' should have dimmension names.") + } + if(!(time_dim_name %in% names(dim(index_fcst1)))) { + stop("Parameter 'index_fcst1' must have temporal dimension.") + } + if (dim(index_fcst1)[time_dim_name] != 1) { + stop("Length of temporal dimension ", + "of parameter 'index_fcst1' must be equal to 1.") + } + if((length(names(dim(index_hind1))) != length(names(dim(index_fcst1)))) + || (!all(names(dim(index_hind1)) == names(dim(index_fcst1))))){ + stop("Dimension names of parameter 'index_hind1' and 'index_fcst1' ", + "must be equals.") + } + if ('member' %in% names(dim(index_fcst1))){ + if (dim(index_fcst1)['member'] == 1) { + stop("Length of dimension 'member' ", + "of parameter 'index_fcst1' must be greater than 1.") + } + } + if ('statistic' %in% names(dim(index_fcst1))){ + if (dim(index_fcst1)['statistic'] != 2) { + stop("Length of dimension 'statistic' ", + "of parameter 'index_fcst1' must be equal to 2.") + } + } + } + if (!is.null(index_fcst2)){ + if (!is.array(index_fcst2)) { + stop("Parameter 'index_fcst2' must be an array.") + } + if (is.null(names(dim(index_fcst2)))){ + stop("Parameters 'index_fcst2' should have dimmension names.") + } + if(!(time_dim_name %in% names(dim(index_fcst2)))) { + stop("Parameter 'index_fcst2' must have temporal dimension.") + } + if (dim(index_fcst2)[time_dim_name] != 1) { + stop("Length of temporal dimension ", + "of parameter 'index_fcst2' must be equal to 1.") + } + if((length(names(dim(index_hind2))) != length(names(dim(index_fcst2)))) + || (!all(names(dim(index_hind2)) == names(dim(index_fcst2))))){ + stop("Dimension names of parameter 'index_hind2' and 'index_fcst2' ", + "must be equals.") + } + if ('member' %in% names(dim(index_fcst2))){ + if (dim(index_fcst2)['member'] == 1) { + stop("Length of dimension 'member' ", + "of parameter 'index_fcst2' must be greater than 1.") + } + } + if ('statistic' %in% names(dim(index_fcst1))){ + if (dim(index_fcst1)['statistic'] != 2) { + stop("Length of dimension 'statistic' ", + "of parameter 'index_fcst1' must be equal to 2.") + } + } + } + + if (all(method_BC == 'none')){ + method_BC1 <- 'ME' + method_BC2 <- 'ME' + bc_dataset1 <- FALSE + bc_dataset2 <- FALSE + } else if (length(method_BC) == 1){ + method_BC1 <- method_BC + method_BC2 <- method_BC + bc_dataset1 <- TRUE + bc_dataset2 <- TRUE + } else { + if(method_BC[1] == 'none'){ + method_BC1 <- 'ME' + bc_dataset1 <- FALSE + } else { + method_BC1 <- method_BC[1] + bc_dataset1 <- TRUE + } + if(method_BC[2] == 'none'){ + method_BC2 <- 'ME' + bc_dataset2 <- FALSE + } else { + method_BC2 <- method_BC[2] + bc_dataset2 <- TRUE + } + } + + pdf_hind1 <- PDFIndexHind(index_hind1, index_obs, method = method_BC1, + time_dim_name = time_dim_name, na.rm = na.rm) + pdf_hind2 <- PDFIndexHind(index_hind2, index_obs, method = method_BC2, + time_dim_name = time_dim_name, na.rm = na.rm) + if (is.null(index_fcst1) || is.null(index_fcst2)){ + pdf_best <- Apply(list(pdf_hind1, pdf_hind2), + target_dims = list('statistic', + 'statistic'), + fun = .BEI_PDFBest, + bc_dataset1, bc_dataset2) + } else { + pdf_fcst1 <- PDFIndexFcst(index_hind1, index_obs, index_fcst1, + method_BC1, time_dim_name = time_dim_name, + na.rm = na.rm) + pdf_fcst2 <- PDFIndexFcst(index_hind2, index_obs, index_fcst2, + method_BC2, time_dim_name = time_dim_name, + na.rm = na.rm) + pdf_best <- Apply(list(pdf_fcst1, pdf_fcst2), + target_dims = list('statistic', + 'statistic'), + fun = .BEI_PDFBest, + bc_dataset1, bc_dataset2) + } + Dim <- names(dim(index_hind1)) + if('member' %in% Dim){ + pos_aux <- match('member', Dim) + Dim[pos_aux] <- 'statistic' + } + pos <- match(Dim, names(dim(pdf_best$output1))) + res <- aperm(pdf_best$output1,pos) + names(dim(res)) <- Dim + return(res) +} + +#' Atomic BEI_PDFBest +#'@param pdf_1 Statistics array for the first SFS PDF with one dimension +#' 'statistic' equal to 4. +#' @param pdf_2 Statistics array for the second SFS PDF with one dimension +#' 'statistic' equal to 4. +#' @param bc_dataset1 Logical (default = TRUE). +#' If TRUE the Index PDFs for the first SFS has been computed +#' with bias corrected. +#' @param bc_dataset2 Logical (default = TRUE). +#' If TRUE the Index PDFs for the second SFS has been computed +#' with bias corrected. +#' +#' @return .BEI_PDFBest returns an array with dimensions (statistic = 2). +#' The firt statistic is the parameter 'mean' of the PDF for the best estimation +#' combining the two SFSs PDF. +#' The second statistic is the parameter 'standard deviation' of the PDF for +#' the best estimation combining the two SFSs PDF. +#' +#' @examples +#' # Example for the Atomic BEI_PDFBest function +#' pdf_1 <- c(1.1,0.6,1.6,0.9) +#' dim(pdf_1) <- c(statistic = 4) +#' pdf_2 <- c(1,0.5,1.5,0.8) +#' dim(pdf_2) <- c(statistic = 4) +#' res <- .BEI_PDFBest(pdf_1, pdf_2, bc_dataset1 = TRUE, bc_dataset2 = FALSE) +#' str(res) +#' dim(res) +#' # statistic +#' # 2 +#' @noRd +.BEI_PDFBest <- function(pdf_1, pdf_2, bc_dataset1 = TRUE, bc_dataset2 = TRUE) { + if(bc_dataset1){ + # apply bias correction to model 1 + mean_1 <- pdf_1[3] + sigma_1 <- pdf_1[4] + } else{ + # not bias correction to model 1 + mean_1 <- pdf_1[1] + sigma_1 <- pdf_1[2] + } + if(bc_dataset2){ + # apply bias correction to model 2 + mean_2 <- pdf_2[3] + sigma_2 <- pdf_2[4] + } else { + # not bias correction to model 2 + mean_2 <- pdf_2[1] + sigma_2 <- pdf_2[2] + } + a_1 <- (sigma_2^2)/((sigma_1^2)+(sigma_2^2)) + a_2 <- (sigma_1^2)/((sigma_1^2)+(sigma_2^2)) + pdf_mean <- a_1*mean_1 + a_2*mean_2 + pdf_sigma <- sqrt((sigma_1^2)*(sigma_2^2)/((sigma_1^2)+(sigma_2^2))) + data <- c(pdf_mean, pdf_sigma) + dim(data) <- c(statistic = 2) + return(data) +} + +#'Computing the Index PDFs for a dataset of SFSs for a hindcats period. +#' +#'@author Eroteida Sanchez-Garcia - AEMET, \email{esanchezg@aemet.es} +#' +#'@description This function implements the computation to obtain the index PDFs +#' (e.g. NAO index) to improve the index estimate from SFSs for a hindcast period. +#' +#'@references Regionally improved seasonal forecast of precipitation through Best +#' estimation of winter NAO, Sanchez-Garcia, E. et al., +#' Adv. Sci. Res., 16, 165174, 2019, https://doi.org/10.5194/asr-16-165-2019 +#' +#'@param index_hind Index (e.g. NAO index) array from SFSs +#' with at least two dimensions (time , member) or (time, statistic). +#' The temporal dimension, by default 'time', must be greater than 2. +#' The dimension 'member' must be greater than 1. +#' The dimension 'statistic' must be equal to 2, for containing the two paramenters of +#' a normal distribution (mean and sd) representing the ensemble of a SFS. +#' It is not possible to have the dimension 'member' and 'statistic' together. +#'@param index_obs Index (e.g. NAO index) array from an observational database +#' or reanalysis with at least a temporal dimension (by default 'time'), +#' which must be greater than 2. +#'@param method A character string indicating which methodology is applied +#' to compute the PDFs. One of "ME" (default) or "LMEV". +#'@param time_dim_name A character string indicating the name of the temporal dimension, by default 'time'. +#'@param na.rm Logical (default = FALSE). Should missing values be removed? +#' +#'@return an array with at least two dimensions (time, statistic = 4). +#' The firt statistic is the parameter 'mean' of the PDF with not bias corrected +#' The second statistic is the parameter 'standard deviation' of the PDF +#' with not bias corrected +#' The third statistic is the parameter 'mean' of the PDF with bias corrected +#' The fourth statistic is the parameter 'standard deviation' of the PDF +#' with bias corrected +#' @import multiApply +#' @examples +#' # Example for the PDFIndexHind function +#' # Example 1 +#' index_obs <- 1 : (5 * 3 ) +#' dim(index_obs) <- c(time = 5, season = 3) +#' index_hind <- 1 : (5 * 4 * 3) +#' dim(index_hind) <- c(time = 5, member = 4, season = 3) +#' res <- PDFIndexHind(index_hind, index_obs) +#' dim(res) +#' # time statistic season +#' # 5 4 3 +#' # Example 2 +#' index_obs <- 1 : (5 * 3) +#' dim(index_obs) <- c(time = 5, season = 3) +#' index_hind <- 1 : (5 * 2 * 3) +#' dim(index_hind) <- c(time = 5, statistic = 2, season = 3) +#' res <- PDFIndexHind(index_hind, index_obs) +#' dim(res) +#' # time statistic season +#' # 5 4 3 +#'@noRd +PDFIndexHind <- function(index_hind, index_obs, method ='ME', + time_dim_name = 'time', na.rm = FALSE) { + if (!is.character(time_dim_name)) { + stop("Parameter 'time_dim_name' must be a character string indicating", + " the name of the temporal dimension.") + } + if (length(time_dim_name) > 1) { + warning("Parameter 'time_dim_name' has length greater than 1 and ", + "only the first element will be used.") + time_dim_name <- time_dim_name[1] + } + if (!(method %in% c('ME', 'LMEV'))){ + stop("Parameter 'method' must be equal to 'ME' or 'LMEV'") + } + if (!is.logical(na.rm)) { + stop("Parameter 'na.rm' must be a logical value.") + } + if (!is.array(index_hind)) { + stop("Parameter 'index_hind' must be an array.") + } + if (!is.array(index_obs)) { + stop("Parameter 'index_obs' must be an array.") + } + if (is.null(names(dim(index_hind))) || is.null(names(dim(index_obs)))) { + stop("Parameters 'index_hind' and 'index_obs'", + " should have dimmension names.") + } + if(!(time_dim_name %in% names(dim(index_obs)))) { + stop("Parameter 'index_obs' must have temporal dimension.") + } + if(!(time_dim_name %in% names(dim(index_hind)))) { + stop("Parameter 'index_hind' must have temporal dimension.") + } + if (dim(index_obs)[time_dim_name] != dim(index_hind)[time_dim_name]) { + stop("Length of the temporal dimensions ", + "of the parameter 'index_obs' and 'index_hind' must be equals.") + } + if (dim(index_obs)[time_dim_name] <= 2) { + stop("Length of temporal dimension ", + "of parameter 'index_obs' and 'index_hind' must be greater than 2.") + } + if('member' %in% names(dim(index_hind)) & + 'statistic' %in% names(dim(index_hind))) { + stop("Parameter 'index_hind' must have at least dimension ", + "'member' or 'statistic', not 'member' and 'statistic' together.") + } + + if(!('member' %in% names(dim(index_hind))) & + !('statistic' %in% names(dim(index_hind)))) { + stop("Parameter 'index_hind' must have dimension ", + "'member' or 'statistic'") + } + + if ('member' %in% names(dim(index_hind))){ + if (dim(index_hind)['member'] == 1) { + stop("Length of dimension 'member' ", + "of parameter 'index_hind' must be greater than 1.") + } + } + if ('statistic' %in% names(dim(index_hind))){ + if (dim(index_hind)['statistic'] != 2) { + stop("Length of dimension 'statistic' ", + "of parameter 'index_hind' must be equal to 2.") + } + } + + if ('member' %in% names(dim(index_hind))) { + PDFstatistics <- Apply(list(index_hind, index_obs), + target_dims = list(c(time_dim_name, 'member'), time_dim_name), + fun = .PDFIndexHind, method, time_dim_name, na.rm) + } else if ('statistic' %in% names(dim(index_hind))){ + PDFstatistics <- Apply(list(index_hind, index_obs), + target_dims = list(c(time_dim_name, 'statistic'), time_dim_name), + fun = .PDFIndexHind, method, time_dim_name, na.rm) + } + return(PDFstatistics$output1) +} + + + +#' Atomic PDFIndexHind +#' @param index_hind Index (e.g. NAO index) array from a SFS with dimensions +#' (time, member) or (time, statistic) for a hindcast period. +#' The temporal dimension, by default 'time', must be greater +#' than 2. +#' @param index_obs Index (e.g. NAO index) array from an observational dataset +#' or reanalysis with dimension (time). The temporal dimension, +#' by default 'time', must be greater than 2. +#' @param method A character string indicating which methodology is applied +#' to compute the PDF. One of "ME" (default) or "LMEV". +#' @param time_dim_name A character string indicating the name of the temporal dimension, by default 'time'. +#' @param na.rm Logical. Should missing values be removed? +#' @return .PDFIndexHind returns an array with dimensions (time, statistic = 4). +#' The firt statistic is the parameter 'mean' of the PDF with not bias corrected +#' for the hindcast period. +#' The second statistic is the parameter 'standard deviation' of the PDF +#' with not bias corrected for the hindcast period. +#' The third statistic is the parameter 'mean' of the PDF with bias corrected +#' for the hindcast period. +#' The fourth statistic is the parameter 'standard deviation' of the PDF +#' with bias corrected for the hindcast period. +#' @import multiApply +#' @importFrom verification verify +#' @examples +#' # Example for the Atomic PDFIndexHind function +#' index_obs <- 1 : 10 +#' dim(index_obs) <- c(time = length(index_obs)) +#' index_hind <- 1 : (10 * 3) +#' dim(index_hind) <- c(time = 10, member = 3) +#' res <- .PDFIndexHind(index_hind, index_obs) +#' dim(res) +#' # time statistic +#' # 10 4 +#' @noRd +.PDFIndexHind <- function(index_hind, index_obs, method = 'ME', + time_dim_name = 'time', na.rm = FALSE) { + dimnameshind <- names(dim(index_hind)) + + if ('member' %in% dimnameshind) { + pos <- match(c(time_dim_name, 'member'), names(dim(index_hind))) + index_hind <- aperm(index_hind, pos) + names(dim(index_hind)) <- c(time_dim_name, 'member') + hind_ens_mean <- Apply(list(index_hind), target_dims = 'member', + fun = mean, na.rm = na.rm)$output1 + hind_ens_sd <- Apply(list(index_hind), target_dims = 'member', + fun = sd, na.rm = na.rm)$output1 + } else if ('statistic' %in% dimnameshind){ + pos <- match(c(time_dim_name, 'statistic'), names(dim(index_hind))) + index_hind <- aperm(index_hind,pos) + names(dim(index_hind)) <- c(time_dim_name, 'statistic') + hind_ens_mean <- index_hind[,1] + hind_ens_sd <- index_hind[,2] + } else { + stop("Wrong dimension of parameter 'index_hind'.") + } + error <- hind_ens_mean - index_obs + pdfnotbc_mean <- hind_ens_mean + ind <- 1 : length(index_obs) + pdfnotbc_sd <- unlist(lapply(ind, function(x) {Sigma_notBC(hind_ens_mean[-x], + index_obs[-x], + na.rm)})) + + if (method == 'ME') { + pdfbc <- unlist(lapply(ind, function(x) {MEBC(hind_ens_mean[-x], + index_obs[-x], na.rm)})) + pdfbc_mean <- hind_ens_mean - pdfbc + pdfbc_sd <- unlist(lapply(ind, function(x) {Sigma_BC(pdfbc_mean[-x], + index_obs[-x], + na.rm)})) + } else if (method == 'LMEV') { + # linear model error-variance + hind_ens_var <- hind_ens_sd^2 + dim_hind <- names(dim(index_hind)) + dim_obs <- names(dim(index_obs)) + coeff <- sapply(ind, function(x) {LMEV(index_hind[-x,], index_obs[-x], + dim_hind, dim_obs, na.rm)}) + pdfbc_mean <- hind_ens_mean - (unlist(coeff['coeff1',]) + + unlist(coeff['coeff2',])*hind_ens_var) + + pdfbc_sd <- unlist(lapply(ind, function(x) {Sigma_BC(pdfbc_mean[-x], + index_obs[-x], + na.rm)})) + } else { + stop("Error: Parameter 'method' is wrong") + } + data <- cbind(pdfnotbc_mean, pdfnotbc_sd, pdfbc_mean, pdfbc_sd) + names(dim(data)) <- c(names(dim(index_obs)), 'statistic') + return(data) +} + +#' Computing the Index PDFs for a dataset of SFSs for a forecast year. +#' +#'@author Eroteida Sanchez-Garcia - AEMET, \email{esanchezg@aemet.es} +#' +#'@description This function implements the computation to obtain the index PDFs +#' (e.g. NAO index) to improve the index estimate from SFSs for a forecast year. +#' +#'@references Regionally improved seasonal forecast of precipitation through Best +#' estimation of winter NAO, Sanchez-Garcia, E. et al., +#' Adv. Sci. Res., 16, 165174, 2019, https://doi.org/10.5194/asr-16-165-2019 +#' +#'@param index_hind Index (e.g. NAO index) array from SFSs +#' with at least two dimensions (time , member) or (time, statistic). +#' The temporal dimension, by default 'time', must be greater than 2. +#' The dimension 'member' must be greater than 1. +#' The dimension 'statistic' must be equal to 2, for containing the two paramenters of +#' a normal distribution (mean and sd) representing the ensemble of a SFS. +#' It is not possible to have the dimension 'member' and 'statistic' together. +#'@param index_obs Index (e.g. NAO index) array from an observational database +#' or reanalysis with at least a temporal dimension (by default 'time'), +#' which must be greater than 2. +#'@param index_fcst Index (e.g. NAO index) array from SFSs +#' with at least two dimensions (time , member) or (time, statistic). +#' The temporal dimension, by default 'time', must be equal to 1, the forecast year target. +#' The dimension 'member' must be greater than 1. +#' The dimension 'statistic' must be equal to 2, for containing the two paramenters of +#' a normal distribution (mean and sd) representing the ensemble of a SFS. +#' It is not possible to have the dimension 'member' and 'statistic' together. +#'@param method A character string indicating which methodology is applied +#' to compute the PDFs. One of "ME" (default) or "LMEV". +#'@param time_dim_name A character string indicating the name of the temporal dimension, by default 'time'. +#'@param na.rm Logical (default = FALSE). Should missing values be removed? +#' +#'@return an array with at least two dimensions (time = 1, statistic = 4). +#' The firt statistic is the parameter 'mean' of the PDF with not bias corrected +#' The second statistic is the parameter 'standard deviation' of the PDF +#' with not bias corrected +#' The third statistic is the parameter 'mean' of the PDF with bias corrected +#' The fourth statistic is the parameter 'standard deviation' of the PDF +#' with bias corrected +#' @import multiApply +#' @examples +#' # Example for the PDFIndexFcst function +#' index_fcst <- 1 : (8 * 4) +#' dim(index_fcst) <- c(time = 1, member = 8, season = 4) +#' index_obs<- 1 : (5 * 4) +#' dim(index_obs) <- c(time = 5, season = 4) +#' index_hind <- 1 : (5 * 4 * 4) +#' dim(index_hind) <- c(time = 5, member = 4, season = 4) +#' res <- PDFIndexFcst(index_hind, index_obs, index_fcst) +#' dim(res) +#' # time statistic season +#' # 1 4 4 +#' +#'@noRd +PDFIndexFcst <- function(index_hind, index_obs, index_fcst, + method ='ME', time_dim_name = 'time', + na.rm = FALSE) { + if (!is.character(time_dim_name)) { + stop("Parameter 'time_dim_name' must be a character string indicating", + " the name of the temporal dimension.") + } + if (length(time_dim_name) > 1) { + warning("Parameter 'time_dim_name' has length greater than 1 and ", + "only the first element will be used.") + time_dim_name <- time_dim_name[1] + } + if(!(method %in% c('ME', 'LMEV'))){ + stop("Parameter 'method' must be equal to 'ME' or 'LMEV'") + } + if (!is.logical(na.rm)) { + stop("Parameter 'na.rm' must be a logical value.") + } + if (!is.array(index_hind)) { + stop("Parameter 'index_hind' must be an array.") + } + if (!is.array(index_obs)) { + stop("Parameter 'index_obs' must be an array.") + } + if (!is.array(index_fcst)) { + stop("Parameter 'index_fcst' must be an array.") + } + if (is.null(names(dim(index_hind))) || + is.null(names(dim(index_obs))) || + is.null(names(dim(index_fcst)))) { + stop("Parameters 'index_hind', 'index_obs' and 'index_fcst' ", + "should have dimmension names.") + } + if(!(time_dim_name %in% names(dim(index_hind)))) { + stop("Parameter 'index_hind' must have temporal dimension.") + } + if(!(time_dim_name %in% names(dim(index_obs)))) { + stop("Parameter 'index_obs' must have temporal dimension.") + } + if(!(time_dim_name %in% names(dim(index_fcst)))) { + stop("Parameter 'index_fcst' must have temporal dimension.") + } + if (dim(index_obs)[time_dim_name] != dim(index_hind)[time_dim_name]) { + stop("Length of temporal dimensions ", + "of parameter 'index_obs' and 'index_hind' must be equals.") + } + if (dim(index_obs)[time_dim_name] <= 2) { + stop("Length of temporal dimension ", + "of parameter 'index_obs' and 'index_hind' must be greater than 2.") + } + if (dim(index_fcst)[time_dim_name] != 1) { + stop("Length of temporal dimension ", + "of parameter 'index_fcst' must be equal to 1.") + } + if((length(names(dim(index_hind))) != length(names(dim(index_fcst)))) + || (!all(names(dim(index_hind)) == names(dim(index_fcst))))){ + stop("Dimension names of parameter 'index_hind' and 'index_fcst' ", + "must be equals.") + } + if('member' %in% names(dim(index_hind)) & + 'statistic' %in% names(dim(index_hind))) { + stop("Parameter 'index_hind' and 'index_fcst' must have at least ", + "dimension 'member' or 'statistic', not 'member' and 'statistic' ", + "together.") + } + if(!('member' %in% names(dim(index_hind))) & + !('statistic' %in% names(dim(index_hind)))) { + stop("Parameter 'index_hind' and 'index_fcst' must have dimension ", + "'member' or 'statistic'") + } + if ('member' %in% names(dim(index_hind))){ + if (dim(index_hind)['member'] == 1) { + stop("Length of dimension 'member' ", + "of parameter 'index_hind' must be greater than 1.") + } + } + if ('member' %in% names(dim(index_fcst))){ + if (dim(index_fcst)['member'] == 1) { + stop("Length of dimension 'member' ", + "of parameter 'index_fcst' must be greater than 1.") + } + } + if ('statistic' %in% names(dim(index_hind))){ + if (dim(index_hind)['statistic'] != 2) { + stop("Length of dimension 'statistic' ", + "of parameter 'index_hind' must be equal to 2.") + } + } + if ('statistic' %in% names(dim(index_fcst))){ + if (dim(index_fcst)['statistic'] != 2) { + stop("Length of dimension 'statistic' ", + "of parameter 'index_fcst' must be equal to 2.") + } + } + + if ('member' %in% names(dim(index_hind))){ + PDFstatistics <- Apply(list(index_hind, index_obs, index_fcst), + target_dims = list(c(time_dim_name, 'member'), + time_dim_name, + c(time_dim_name,'member')), + fun = .PDFIndexFcst, method, time_dim_name, na.rm) + } else if ('statistic' %in% names(dim(index_hind))){ + PDFstatistics <- Apply(list(index_hind, index_obs, index_fcst), + target_dims = list(c(time_dim_name, 'statistic'), + time_dim_name, + c(time_dim_name,'statistic')), + fun = .PDFIndexFcst, method, time_dim_name, na.rm) + } + return(PDFstatistics$output1) +} + +#' Atomic PDFIndexFcst +#' @param index_hind Index (e.g. NAO index) array from a SFS with dimensions +#' (time, member) or (time, statistic) for a hindcast period. +#' The temporal dimension, by default 'time', must be greater +#' than 2. +#' @param index_obs Index (e.g. NAO index) array from an observational dataset +#' or reanalysis with dimension (time). The temporal dimension, +#' by default 'time', must be greater than 2. +#' @param index_fcst Index (e.g. NAO index) array from SFSs +#' with dimensions (time , member) or (time, statistic). +#' The temporal dimension, by default 'time', must be equal to 1, +#' the forecast year target. +#' The dimension 'member' must be greater than 1. +#' The dimension 'statistic' must be equal to 2, for containing the two paramenters of +#' a normal distribution (mean and sd) representing the ensemble of a SFS. +#' It is not possible to have the dimension 'member' and 'statistic' together. +#' @param method A character string indicating which methodology is applied +#' to compute the PDF. One of "ME" (default) or "LMEV". +#' @param time_dim_name A character string indicating the name of the temporal dimension, by default 'time'. +#' @param na.rm Logical. Should missing values be removed? +#' @return .PDFIndexFcst returns an array with dimensions (time = 1, statistic=4). +#' The firt statistic is the parameter 'mean' of the PDF with not bias corrected +#' for the forecast year. +#' The second statistic is the parameter 'standard deviation' of the PDF +#' with not bias corrected for the forecast year. +#' The third statistic is the parameter 'mean' of the PDF with bias corrected +#' for the forecast year. +#' The fourth statistic is the parameter 'standard deviation' of the PDF +#' with bias corrected for the forecast year. +#' @import multiApply +#' @importFrom verification verify +#' @examples +#' # Example 1 for the Atomic PDFIndexFcst function +#' index_fcst <- 1 : (1 * 6) +#' dim(index_fcst) <- c(time = 1, member = 6) +#' index_obs <- 1 : 10 +#' dim(index_obs) <- c(time = length(index_obs)) +#' index_hind <- 1 : (10 * 3) +#' dim(index_hind) <- c(time = 10, member = 3) +#' res <- .PDFIndexFcst(index_hind, index_obs, index_fcst) +#' dim(res) +#' # time statistic +#' # 1 4 +#' # Example 2 for the Atomic PDFIndexFcst function +#' index_fcst <- 1 : (1 * 2) +#' dim(index_fcst) <- c(time = 1, statistic = 2) +#' index_obs <- 1 : 10 +#' dim(index_obs) <- c(time = 10) +#' index_hind <- 1 : (10 * 2) +#' dim(index_hind) <- c(time = 10, statistic = 2) +#' res <- .PDFIndexFcst(index_hind, index_obs, index_fcst) +#' dim(res) +#' # time statistic +#' # 1 4 +#' @noRd +.PDFIndexFcst <- function(index_hind, index_obs, index_fcst, + method = 'ME', + time_dim_name = 'time', + na.rm = FALSE) { + dimnameshind <- names(dim(index_hind)) + if ('member' %in% dimnameshind) { + pos <- match(c(time_dim_name, 'member'), names(dim(index_hind))) + index_hind <- aperm(index_hind, pos) + names(dim(index_hind)) <- c(time_dim_name, 'member') + exp_fcst_ens_mean <- Apply(list(index_fcst), target_dims = 'member', + fun = mean)$output1 + exp_fcst_ens_sd <- Apply(list(index_fcst), target_dims = 'member', + fun = sd)$output1 + } else { + pos <- match(c(time_dim_name,'statistic'), names(dim(index_fcst))) + index_fcst <- aperm(index_fcst,pos) + names(dim(index_fcst)) <- c(time_dim_name, 'statistic') + exp_fcst_ens_mean <- index_fcst[,1] + exp_fcst_ens_sd <- index_fcst[,2] + } + data_hind <- .PDFIndexHind(index_hind, index_obs, method, time_dim_name, na.rm) + exp_hind_ens_mean <- data_hind[,1] + pdfnotbc_mean <- exp_fcst_ens_mean + pdfnotbc_sd <- Sigma_notBC(exp_hind_ens_mean, index_obs, na.rm) + if (method == 'ME') { + pdfbc <- MEBC(exp_hind_ens_mean, index_obs, na.rm) + pdfbc_mean <- exp_fcst_ens_mean - pdfbc + } else if (method == 'LMEV') { + # linear model error-variance + exp_fcst_ens_var <- exp_fcst_ens_sd^2 + dim_hind <- names(dim(index_hind)) + dim_obs <- names(dim(index_obs)) + coeff <- LMEV(index_hind, index_obs, dim_hind, dim_obs, time_dim_name, na.rm) + pdfbc_mean <- exp_fcst_ens_mean - (coeff$coeff1 + coeff$coeff2*exp_fcst_ens_var) + } else { + stop("Parameter 'method' must be 'ME' or 'LMEV'.") + } + pdfbc_sd <- Sigma_BC(data_hind[,3], index_obs, na.rm) + data <- cbind(pdfnotbc_mean, pdfnotbc_sd, pdfbc_mean, pdfbc_sd) + names(dim(data)) <- c(names(dim(index_obs)), 'statistic') + return(data) +} + +# Auxiliar function to compute the mean squared error between 'exp' and 'obs' +Sigma_notBC <- function(exp, obs, na.rm = FALSE) { + if (na.rm){ + indNA <- c(which(is.na(exp)),which(is.na(obs))) + if (length(indNA) > 0){ + exp <- exp[-indNA] + obs <- obs[-indNA] + } + } + return(sqrt(verify(obs, exp, frcst.type = "cont", + obs.type = "cont")$MSE)) +} + +# Auxiliar function to compute the bias between 'exp' and 'obs' +MEBC <- function(exp, obs, na.rm = FALSE) { + if (na.rm){ + indNA <- c(which(is.na(exp)),which(is.na(obs))) + if (length(indNA) > 0){ + exp <- exp[-indNA] + obs <- obs[-indNA] + } + } + return(verify(obs, exp, frcst.type = "cont", + obs.type = "cont")$ME) +} + +# Auxiliar function to compute the standard deviation of errors +# between 'obs' and 'exp' +Sigma_BC <- function(exp, obs, na.rm = FALSE) { + return(sd(obs - exp, na.rm = na.rm)) +} + +# Auxiliar function to compute the linear model used in the method +# called 'LMEV' to correct the bias in the Best NAO weighting +# methogology based on a linear model using ensemble variance of NAO index +# as predictor. +LMEV <- function(exp, obs, dim_exp, + dim_obs, time_dim_name = 'time', na.rm = FALSE) { + names(dim(exp)) <- dim_exp + names(dim(obs)) <- dim_obs + if (any(names(dim(exp)) == 'member')) { + exp_ens_mean <- Apply(list(exp), target_dims = 'member', + fun = mean, na.rm = na.rm)$output1 + exp_ens_sd <- Apply(list(exp), target_dims = 'member', + fun = sd, na.rm = na.rm)$output1 + } else { + if(na.rm & any(is.na(exp))){ + print("Some value are NA") + } + pos <- match(c(time_dim_name,'statistic'), names(dim(exp))) + exp <- aperm(exp,pos) + exp_ens_mean <- exp[,1] + exp_ens_sd <- exp[,2] + } + error <- exp_ens_mean - obs + exp_ens_var <- exp_ens_sd^2 + statModel <- lm(error ~ exp_ens_var) + return(list(coeff1 = statModel$coefficients[1], coeff2 = statModel$coefficients[2])) +} + + diff --git a/R/BEI_Weights.R b/R/BEI_Weights.R new file mode 100644 index 0000000000000000000000000000000000000000..63fb22c311618ecbe8e9259203ccc4796c7d9d9c --- /dev/null +++ b/R/BEI_Weights.R @@ -0,0 +1,133 @@ +#'Computing the weights for SFSs using the Best Index PDFs. +#' +#'@author Eroteida Sanchez-Garcia - AEMET, \email{esanchezg@aemet.es} +#' +#'@description This function implements the computation to obtain the +#' normalized weights for each member of each Seasonal Forecast Systems (SFS) +#' or dataset using the Probability Density Functions (PDFs) indicated by the +#' parameter 'pdf_weight' (for instance the Best Index estimation obtained +#' using the 'PDFBest' function). The weight of each member is proportional to +#' the probability of its index calculated with the PDF "pdf_weight". +#' +#'@references Regionally improved seasonal forecast of precipitation through Best +#' estimation of winter NAO, Sanchez-Garcia, E. et al., +#' Adv. Sci. Res., 16, 165174, 2019, https://doi.org/10.5194/asr-16-165-2019 +#' +#'@param index_weight Index (e.g. NAO index) array, from a dataset of SFSs +#' for a period of years, with at least dimensions 'member'. +#' Additional dimensions, for instance, a temporal dimension as 'time', +#' must have the same lenght in both parameters, +#' 'index_weight' and 'pdf_weight'. +#' @param pdf_weight Statistics array to define a Gaussian PDF with at least +#' dimensions 'statistic'. +#' The firt statistic is the parameter 'mean' of the PDF and +#' the second statistic is the parameter 'standard deviation' of the PDF. +#' @param time_dim_name A character string indicating the name of the temporal dimension, by default 'time'. +#' +#'@return BEI_Weights() returns a normalized weights array with the same +#' dimensions that index_weight. +#' +#'@import multiApply +#' +#'@examples +#' # Example for the BEI_Weights function +#' index_weight <- 1 : (10 * 3 * 5 * 1) +#' dim(index_weight) <- c(sdate = 10, dataset = 3, member = 5, season = 1) +#' pdf_weight <- 1 : (10 * 3 * 2 * 1) +#' dim(pdf_weight) <- c(sdate = 10, dataset = 3, statistic = 2, season = 1) +#' res <- BEI_Weights(index_weight, pdf_weight) +#' dim(res) +#' # sdate dataset member season +#' # 10 3 5 1 +#' +#'@export +BEI_Weights <- function(index_weight, pdf_weight, time_dim_name = 'time') { + + if (!is.character(time_dim_name)) { + stop("Parameter 'time_dim_name' must be a character string ", + "indicating the name of the temporal dimension.") + } + if (length(time_dim_name) > 1) { + warning("Parameter 'time_dim_name' has length greater than 1 and ", + "only the first element will be used.") + time_dim_name <- time_dim_name[1] + } + if (!is.array(index_weight)) { + stop("Parameter 'index_weight' must be an array.") + } + if (!is.array(pdf_weight)) { + stop("Parameter 'pdf_weight' must be an array.") + } + if (is.null(names(dim(index_weight))) || is.null(names(dim(pdf_weight)))) { + stop("Parameters 'index_weight' and 'pdf_weight'", + " should have dimmension names.") + } + if(!('member' %in% names(dim(index_weight)))) { + stop("Parameter 'index_weight' must have dimension 'member'.") + } + if(!('statistic' %in% names(dim(pdf_weight)))) { + stop("Parameter 'pdf_weight' must have dimension 'statistic'.") + } + if(time_dim_name %in% names(dim(index_weight)) & + !time_dim_name %in% names(dim(pdf_weight))) { + stop("Parameter 'pdf_weight' must have temporal dimension.") + } + if(!time_dim_name %in% names(dim(index_weight)) & + time_dim_name %in% names(dim(pdf_weight))) { + stop("Parameter 'index_weight' must have temporal dimension.") + } + if(time_dim_name %in% names(dim(index_weight)) & + time_dim_name %in% names(dim(pdf_weight)) & + dim(index_weight)[time_dim_name] != dim(pdf_weight)[time_dim_name]){ + stop("Length of temporal dimension of parameters must be equal") + } + if (dim(pdf_weight)['statistic'] != 2) { + stop("Length of dimension 'statistic' ", + "of the parameter 'pdf_weight' must be equal to 2.") + } + + + aweights <- Apply(list(index_weight, pdf_weight), + target_dims = list('member', 'statistic'), + fun = .BEI_Weights)$output1 + + dimnames <- names(dim(index_weight)) + pos <- match(dimnames, names(dim(aweights))) + aweights <- aperm(aweights, pos) + names(dim(aweights)) <- dimnames + return(aweights) +} + +#' Atomic BEI_Weights +#' @param index_weight Index (e.g. NAO index) array from a SFS with dimensions +#' (member) +#' @param pdf_weight Statistics array to define a Gaussian PDF with dimensions +#' (statistic = 2). +#' The firt statistic is the parameter 'mean' of the PDF and +#' the second statistic is the parameter 'standard deviation' of the PDF. +#' @return .BEI_Weights returns an array of with dimensions (member), +#' the normalized weights for each member of a SFS using a Best NAO PDF. +#' @examples +#' # Example for the Atomic BEI_Weights function +#' index_weight <- c(1.3,3,-1) +#' dim(index_weight) <- c(member = 3) +#' pdf_weight <- c(1.5,0.8) +#' dim(pdf_weight) <- c(statistic = 2) +#' res <- .BEI_Weights(index_weight, pdf_weight) +#' dim(res) +#' # member +#' # 3 +#' @noRd +.BEI_Weights <- function(index_weight, pdf_weight) { + aweights <- apply(index_weight, 1, dnorm, mean = pdf_weight[1], sd = pdf_weight[2]) + dim(aweights) <- dim(index_weight) + sumWeights <- sum(aweights) + aweightsNorm <- apply(aweights, 1, NormWeight, sumWeights) + dim(aweightsNorm) <- dim(index_weight) + return(aweightsNorm) +} + +# Auxiliar function to normalize a weight value +NormWeight <- function(weight, sumWeights) { + return(weight/sumWeights) +} diff --git a/R/CST_BEI_Weighting.R b/R/CST_BEI_Weighting.R new file mode 100644 index 0000000000000000000000000000000000000000..51c8ba32fdaede17917f610ecd685aa6991ebcda --- /dev/null +++ b/R/CST_BEI_Weighting.R @@ -0,0 +1,453 @@ +#' Weighting SFSs of a CSTools object. +#' +#' @author Eroteida Sanchez-Garcia - AEMET, \email{esanchezg@aemet.es} +#' +#' @description Function to apply weights to a 's2dv_cube' object. +#' It could return a weighted ensemble mean (deterministic output) or +#' the terciles probabilities (probabilistic output) for Seasonal Forecast +#' Systems (SFSs). +#' +#' @references Regionally improved seasonal forecast of precipitation through +#' Best estimation of winter NAO, Sanchez-Garcia, E. et al., +#' Adv. Sci. Res., 16, 165174, 2019, https://doi.org/10.5194/asr-16-165-2019 +#' +#' @param var_exp An object of the class 's2dv_cube' containing the variable +#' (e.g. precipitation, temperature, NAO index) array. +#' The var_exp object is expected to have an element named \code{$data} with +#' at least a temporal dimension and a dimension named 'member'. +#' @param aweights Normalized weights array with at least dimensions +#' (time, member), when 'time' is the temporal dimension as default. +#' When 'aweights' parameter has any other dimensions (as e.g. 'lat') and +#' 'var_exp' parameter has also the same dimension, they must be equals. +#' @param type A character string indicating the type of output. +#' If 'type' = 'probs', the function returns, in the element data from +#' 'var_exp' parameter, an array with at least two +#' or four dimensions depending if the variable is spatially aggregated variable +#' (as e.g. NAO index), dimension (time, tercil) or it is spatial variable +#' (as e.g. precipitation or temperature), dimension (time, tercile, lat, lon), +#' containing the terciles probabilities computing with weighted members. +#' The first tercil is the lower tercile, the second is the normal tercile and +#' the third is the upper tercile. +#' If 'type' = 'ensembleMean', the function returns, in the element data from +#' 'var_exp' parameter, an array with at least one or three dimensions +#' depending if the variable is a spatially aggregated variable +#' (as e.g. NAO index)(time) or it is spatial variable (as e.g. precipitation +#' or temperature) (time, lat, lon), containing the ensemble means computing +#' with weighted members. +#' @param time_dim_name A character string indicating the name of the +#' temporal dimension, by default 'time'. +#' +#' @return CST_BEI_Weighting() returns a CSTools object (i.e., of the +#' class 's2dv_cube'). +#' This object has at least an element named \code{$data} +#' with at least a temporal dimension (and dimension 'tercil' when the output +#' are tercile probabilities), containing the ensemble means computing with +#' weighted members or probabilities of terciles. +#' +#' @examples +#' var_exp <- 1 : (2 * 4 * 3 * 2) +#' dim(var_exp) <- c(time = 2, member = 4, lat = 3, lon = 2) +#' aweights <- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3, 0.1, 0.2, 0.4, 0.4, 0.1, 0.2, 0.4, 0.2) +#' dim(aweights) <- c(time = 2, member = 4, dataset = 2) +#' var_exp <- list(data = var_exp) +#' class(var_exp) <- 's2dv_cube' +#' res_CST <- CST_BEI_Weighting(var_exp, aweights) +#' dim(res_CST$data) +#' # time lat lon dataset +#' # 2 3 2 2 +#' @export +CST_BEI_Weighting <- function(var_exp, aweights, type = 'ensembleMean', + time_dim_name = 'time') { + if (!is.character(time_dim_name)) { + stop("Parameter 'time_dim_name' must be a character string indicating", + " the name of the temporal dimension.") + } + if (length(time_dim_name) > 1) { + warning("Parameter 'time_dim_name' has length greater than 1 and ", + "only the first element will be used.") + time_dim_name <- time_dim_name[1] + } + if (!is.character(type)) { + stop("Parameter 'type' must be a character string, 'probs' or ", + "'ensembleMean', indicating the type of output.") + } + if (length(type) > 1) { + warning("Parameter 'type' has length greater than 1 and ", + "only the first element will be used.") + type <- type[1] + } + if (!inherits(var_exp, 's2dv_cube')) { + stop("Parameter 'var_exp' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + if (!is.array(aweights)) { + stop("Parameter 'aweights' must be an array.") + } + if (is.null(names(dim(var_exp$data))) || is.null(names(dim(aweights)))) { + stop("Element 'data' from parameter 'var_exp' and parameter 'aweights'", + " should have dimmension names.") + } + if(!(time_dim_name %in% names(dim(var_exp$data)))) { + stop("Element 'data' from parameter 'var_exp' must have ", + "temporal dimension.") + } + if(!(time_dim_name %in% names(dim(aweights)))) { + stop("Parameter 'aweights' must have temporal dimension.") + } + if(!('member' %in% names(dim(var_exp$data)))) { + stop("Element 'data' from parameter 'var_exp' must have ", + "dimension 'member'.") + } + if(!('member' %in% names(dim(aweights)))) { + stop("Parameter 'aweights' must have dimension 'member'.") + } + if (dim(var_exp$data)[time_dim_name] != dim(aweights)[time_dim_name]) { + stop("Length of temporal dimensions ", + "of element 'data' from parameter 'var_exp' and parameter ", + "'aweights' must be equals.") + } + if (dim(var_exp$data)['member'] != dim(aweights)['member']) { + stop("Length of dimension 'member' of element 'data' from ", + "parameter 'var_exp' and parameter 'aweights' must be equals.") + } + + if (type == 'ensembleMean'){ + em <- BEI_EMWeighting(var_exp$data, aweights, time_dim_name) + var_exp$data <- em + } else if (type == 'probs'){ + probs <- BEI_ProbsWeighting(var_exp$data, aweights, time_dim_name) + var_exp$data <- probs + } else { + stop("Parameter 'type' must be a character string ('probs' or ", + "'ensembleMean'), indicating the type of output.") + } + return(var_exp) +} + + +#' @title Computing the weighted ensemble means for SFSs. +#' @author Eroteida Sanchez-Garcia - AEMET, \email{esanchezg@aemet.es} +#' @description This function implements the computation to obtain the weighted +#' ensemble means for SFSs using a normalized weights array, +#' +#' @references Regionally improved seasonal forecast of precipitation through Best +#' estimation of winter NAO, Sanchez-Garcia, E. et al., +#' Adv. Sci. Res., 16, 165174, 2019, https://doi.org/10.5194/asr-16-165-2019 +#' +#' @param var_exp Variable (e.g. precipitation, temperature, NAO index) +#' array from a SFS with at least dimensions (time, member) for a spatially +#' aggregated variable or dimensions (time, member, lat, lon) for a spatial +#' variable, as 'time' the spatial dimension by default. +#' @param aweights Normalized weights array with at least dimensions +#' (time, member), when 'time' is the temporal dimension as default. +#' @param time_dim_name A character string indicating the name of the +#' temporal dimension, by default 'time'. +#' +#' @return BEI_EMWeighting() returns an array with at least one or three +#' dimensions depending if the variable is spatially aggregated variable +#' (as e.g. NAO index)(time) or it is spatial variable (as e.g. precipitation +#' or temperature) (time, lat, lon), containing the ensemble means computing +#' with weighted members. +#' @import multiApply +#' +#' @examples +#' # Example 1 +#' var_exp <- 1 : (2 * 3 * 4) +#' dim(var_exp) <- c(time = 2, dataset = 3, member = 4) +#' aweights<- runif(24, min=0.001, max=0.999) +#' dim(aweights) <- c(time = 2, dataset = 3, member = 4) +#' res <- BEI_EMWeighting(var_exp, aweights) +#' dim(res) +#' # time dataset +#' # 2 3 +#' # Example 2 +#' var_exp <- 1 : (2 * 4 * 2 * 3) +#' dim(var_exp) <- c(time = 2, member = 4, lat = 2, lon = 3) +#' aweights<- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3) +#' dim(aweights) <- c(time = 2, member = 4) +#' res <- BEI_EMWeighting(var_exp, aweights) +#' dim(res) +#' # time lat lon +#' # 2 2 3 +#' +#' @noRd +BEI_EMWeighting <- function(var_exp, aweights, time_dim_name = 'time') { + if (!is.character(time_dim_name)) { + stop("Parameter 'time_dim_name' must be a character string indicating", + " the name of the temporal dimension.") + } + if (length(time_dim_name) > 1) { + warning("Parameter 'time_dim_name' has length greater than 1 and ", + "only the first element will be used.") + time_dim_name <- time_dim_name[1] + } + if (!is.array(var_exp)) { + stop("Parameter 'var_exp' must be an array.") + } + if (!is.array(aweights)) { + stop("Parameter 'aweights' must be an array.") + } + if (is.null(names(dim(var_exp))) || is.null(names(dim(aweights)))) { + stop("Parameters 'var_exp' and 'aweights'", + " should have dimmension names.") + } + if(!(time_dim_name %in% names(dim(var_exp)))) { + stop("Parameter 'var_exp' must have temporal dimension.") + } + if(!(time_dim_name %in% names(dim(aweights)))) { + stop("Parameter 'aweights' must have temporal dimension.") + } + if(!('member' %in% names(dim(var_exp)))) { + stop("Parameter 'var_exp' must have temporal dimension.") + } + if(!('member' %in% names(dim(aweights)))) { + stop("Parameter 'aweights' must have temporal dimension.") + } + if (dim(var_exp)[time_dim_name] != dim(aweights)[time_dim_name]) { + stop("Length of temporal dimensions ", + "of parameter 'var_exp' and 'aweights' must be equals.") + } + if (dim(var_exp)['member'] != dim(aweights)['member']) { + stop("Length of dimension 'member' ", + "of parameter 'var_exp' and 'aweights' must be equals.") + } + + res <- Apply(list(var_exp, aweights), + target_dims = list(c(time_dim_name,'member'), + c(time_dim_name,'member')), + fun = .BEI_EMWeighting, time_dim_name)$output1 + return(res) +} + +#' Atomic BEI_EMWeighting +#' @param var_exp Variable (e.g. precipitation, temperature, NAO index) +#' array from a SFS with a temporal dimension, +#' by default 'time', and dimension 'member'. +#' @param aweights Normalized weights array with a temporal dimension, +#' by default 'time', and dimension 'member' +#' @param time_dim_name A character string indicating the name of the +#' temporal dimension, by default 'time'. +#' @return .BEI_EMWeighting returns an array of with a temporal dimension, +#' by default 'time', containing the weighted ensemble means. +#' @examples +#' # Example for the Atomic BEI_EMWeighting function +#' var_exp <- 1 : 6 +#' dim(var_exp) <- c(time = 2, member = 3) +#' aweights <- c(0.28, 0.15, 0.69, 0.64, 0.42, 0.17) +#' dim(aweights) <- c(time = 2, member = 3) +#' res <- .BEI_EMWeighting(var_exp, aweights) +#' dim(res) +#' # time +#' # 2 +#' @noRd +.BEI_EMWeighting <- function(var_exp, aweights, time_dim_name = 'time') { + + posTime <- match(time_dim_name, names(dim(var_exp))) + var_exp_em <- as.array(apply(var_exp * aweights, posTime, sum)) + names(dim(var_exp_em)) <- time_dim_name + return(var_exp_em) +} + + +#' Computing the weighted tercile probabilities for SFSs. +#' @author Eroteida Sanchez-Garcia - AEMET, \email{esanchezg@aemet.es} +#' +#' @description This function implements the computation to obtain the tercile +#' probabilities for a weighted variable for SFSs using a normalized weights array, +#' +#' @references Regionally improved seasonal forecast of precipitation through Best +#' estimation of winter NAO, Sanchez-Garcia, E. et al., +#' Adv. Sci. Res., 16, 165174, 2019, https://doi.org/10.5194/asr-16-165-2019 +#' +#' @param var_exp Variable (e.g. precipitation, temperature, NAO index) +#' array from a SFS with at least dimensions (time, member) for a spatially +#' aggregated variable or dimensions (time, member, lat, lon) for a spatial +#' variable, as 'time' the spatial dimension by default. +#' @param aweights Normalized weights array with at least dimensions +#' (time, member), when 'time' is the temporal dimension as default. +#' @param time_dim_name A character string indicating the name of the +#' temporal dimension, by default 'time'. +#' +#' @return BEI_ProbsWeighting() returns an array with at least two or four +#' dimensions depending if the variable is a spatially aggregated variable +#' (as e.g. NAO index)(time, tercil) or it is spatial variable (as e.g. +#' precipitation or temperature)(time, tercile, lat, lon), containing the +#' terciles probabilities computing with weighted members. +#' The first tercil is the lower tercile, the second is the normal tercile and +#' the third is the upper tercile. +#' +#' @import multiApply +#' +#' @examples +#' # Example 1 +#' var_exp <- 1 : (2 * 4) +#' dim(var_exp) <- c(time = 2, member = 4) +#' aweights<- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3) +#' dim(aweights) <- c(time = 2, member = 4) +#' res <- BEI_ProbsWeighting(var_exp, aweights) +#' dim(res) +#' # time tercil +#' # 2 3 +#' # Example 2 +#' var_exp <- rnorm(48, 50, 9) +#' dim(var_exp) <- c(time = 2, member = 4, lat = 2, lon = 3) +#' aweights<- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3) +#' dim(aweights) <- c(time = 2, member = 4) +#' res <- BEI_ProbsWeighting(var_exp, aweights) +#' dim(res) +#' # time tercil lat lon +#' # 2 3 2 3 +#' @noRd +BEI_ProbsWeighting <- function(var_exp, aweights, time_dim_name = 'time') { + + if (!is.character(time_dim_name)) { + stop("Parameter 'time_dim_name' must be a character string indicating", + " the name of the temporal dimension.") + } + if (length(time_dim_name) > 1) { + warning("Parameter 'time_dim_name' has length greater than 1 and ", + "only the first element will be used.") + time_dim_name <- time_dim_name[1] + } + if (!is.array(var_exp)) { + stop("Parameter 'var_exp' must be an array.") + } + if (!is.array(aweights)) { + stop("Parameter 'aweights' must be an array.") + } + if (is.null(names(dim(var_exp))) || is.null(names(dim(aweights)))) { + stop("Parameters 'var_exp' and 'aweights'", + " should have dimmension names.") + } + if(!(time_dim_name %in% names(dim(var_exp)))) { + stop("Parameter 'var_exp' must have temporal dimension.") + } + if(!(time_dim_name %in% names(dim(aweights)))) { + stop("Parameter 'aweights' must have temporal dimension.") + } + if(!('member' %in% names(dim(var_exp)))) { + stop("Parameter 'var_exp' must have temporal dimension.") + } + if(!('member' %in% names(dim(aweights)))) { + stop("Parameter 'aweights' must have temporal dimension.") + } + if (dim(var_exp)[time_dim_name] != dim(aweights)[time_dim_name]) { + stop("Length of temporal dimensions ", + "of parameter 'var_exp' and 'aweights' must be equals.") + } + if (dim(var_exp)['member'] != dim(aweights)['member']) { + stop("Length of dimension 'member' ", + "of parameter 'var_exp' and 'aweights' must be equals.") + } + + res <- Apply(list(var_exp, aweights), + target_dims = list(c(time_dim_name,'member'), c(time_dim_name,'member')), + fun = .BEI_ProbsWeighting, time_dim_name)$output1 + return(res) +} + +#' Atomic BEI_ProbsWeighting +#' @param var_exp Variable (e.g. precipitation, temperature, NAO index) +#' array from a SFS with a temporal dimension, +#' by default 'time', and dimension 'member'. +#' @param aweights Normalized weights array with a temporal dimension, +#' by default 'time', and dimension 'member' +#' @param time_dim_name A character string indicating the name of the +#' temporal dimension, by default 'time'. +#' @return .BEI_ProbsWeighting returns an array of with a temporal dimension, +#' as default 'time', and 'tercil' dimension, containing the probabilities +#' for each tercile computing with weighted members. +#' The firt tercil is the lower tercile, the second is the normal tercile and +#' the third is the upper tercile. +#' @examples +#' # Example +#' var_exp <- 1 : 8 +#' dim(var_exp) <- c(stime = 2, member = 4) +#' aweights <- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.3, 0.3) +#' dim(aweights) <- c(stime = 2, member = 4) +#' res <- .BEI_ProbsWeighting(var_exp, aweights, time_dim_name = 'stime') +#' dim(res) +#' # stime tercil +#' # 2 3 +#' @noRd +.BEI_ProbsWeighting <- function(var_exp, aweights, time_dim_name = 'time') { + # computing terciles + terciles_exp <- WeightTerciles(var_exp, aweights, time_dim_name) + lowerTercile <- terciles_exp$lowerTercile + upperTercile <- terciles_exp$upperTercile + + # Probabilities + aTerciles <- Apply(list(var_exp), target_dims = list('member'), + fun = Data2Tercil, lowerTercile, upperTercile)$output1 + + pos <- match(names(dim(aTerciles)), c(time_dim_name,'member')) + aTerciles <- aperm(aTerciles,pos) + names(dim(aTerciles)) <- c(time_dim_name,'member') + + probTercile <- array(NA, dim = c(dim(var_exp)[time_dim_name], tercil = 3)) + for (idTercil in 1:3){ + probTercile[,idTercil] <- Apply(list(aTerciles, aweights), + target_dims = list('member','member'), + fun = WeightTercil2Prob, idTercil)$output1 + } + return(probTercile) +} + +# Auxiliar function to compute in which tercile is a data value +Data2Tercil <- function(x,lt,ut) { + y <- rep(2,length(x)) + y[x <= lt] <- 1 + y[x >= ut] <- 3 + if (lt == ut) { + warning("The upper and lower terciles are equals") + } + dim(y) <- c(member = length(x)) + return (y) +} + +# Auxiliar function to convers weighted terciles to probabilities +WeightTercil2Prob <- function(aTerciles, aWeights, idTercil) { + return(sum(aWeights[which(aTerciles == idTercil)])) +} + +# Auxiliar function to compute weighted terciles +WeightTerciles <- function(data, aweights, time_dim_name = 'time') { + namesdimdata <- names(dim(data)) + namesdimaweights <- names(dim(aweights)) + pos <- match(namesdimdata, namesdimaweights) + aweights <- aperm(aweights, pos) + names(dim(aweights)) <- namesdimdata + vectorData <- as.vector(data) + vectorWeights <- as.vector(aweights/dim(aweights)[time_dim_name]) # normalized + lSortData <- sort(vectorData,index.return=TRUE) + indSort <- lSortData$ix # index asociated to weight + # corresponding for this data + dataSort <- lSortData$x + # Adding normalized weights. When 1/3 is reached, the data value + # is lower tercile and when 2/3 is reached, it is the upper tercile. + sumWeights <- 0 + ilowerTercile <- 0 + while ((sumWeights < 1/3) & (ilowerTercile < length(aweights))){ + ilowerTercile<- ilowerTercile +1 + sumWeights <- sumWeights + vectorWeights[indSort[ilowerTercile]] + } + if (ilowerTercile == 1){ + lowerTercile <- dataSort[ilowerTercile] + } else { + lowerTercile <- (dataSort[ilowerTercile]+ + dataSort[ilowerTercile-1])/2 + } + sumWeights <- 0 + iupperTercile <- 0 + while ((sumWeights < 2/3) & (iupperTercile < length(aweights))){ + iupperTercile<- iupperTercile +1 + sumWeights <- sumWeights + vectorWeights[indSort[iupperTercile]] + } + if (iupperTercile == 1) { + upperTercile <- dataSort[iupperTercile] + } else { + upperTercile <- (dataSort[iupperTercile]+ + dataSort[iupperTercile-1])/2 + } + return(list(lowerTercile = lowerTercile, upperTercile = upperTercile)) +} diff --git a/R/zzz.R b/R/zzz.R index 662624d9cb94cff6436d71016703c558c974c579..1414e9490cd4f221d9d707bf3ff922f178369ba5 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -16,3 +16,4 @@ attributes(x) <- c(attributes(x), attr_bk) x } + diff --git a/man/BEI_PDFBest.Rd b/man/BEI_PDFBest.Rd new file mode 100644 index 0000000000000000000000000000000000000000..f836ab722211e359009ca03e7844a00a9b6d2162 --- /dev/null +++ b/man/BEI_PDFBest.Rd @@ -0,0 +1,124 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/BEI_PDFBest.R +\name{BEI_PDFBest} +\alias{BEI_PDFBest} +\title{Computing the Best Index PDFs combining Index PDFs from two SFSs} +\usage{ +BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1 = NULL, + index_fcst2 = NULL, method_BC = "none", time_dim_name = "time", + na.rm = FALSE) +} +\arguments{ +\item{index_obs}{Index (e.g. NAO index) array from an observational database +or reanalysis with at least a temporal dimension (by default 'time'), +which must be greater than 2.} + +\item{index_hind1}{Index (e.g. NAO index) array from a SFS (named SFS1) +with at least two dimensions (time , member) or (time, statistic). +The temporal dimension, by default 'time', must be greater than 2. +The dimension 'member' must be greater than 1. +The dimension 'statistic' must be equal to 2, for containing the two paramenters of +a normal distribution (mean and sd) representing the ensemble of a SFS. +It is not possible to have the dimension 'member' and 'statistic' at the same time.} + +\item{index_hind2}{Index (e.g. NAO index) array from a SFS (named SFS2) +with at least two dimensions (time , member) or (time, statistic). +The temporal dimension, by default 'time', must be greater than 2. +The dimension 'member' must be greater than 1. +The dimension 'statistic' must be equal to 2, for containing the two paramenters of +a normal distribution (mean and sd) representing the ensemble of a SFS. +It is not possible to have the dimension 'member' and 'statistic' together.} + +\item{index_fcst1}{(optional, default = NULL) Index (e.g. NAO index) array from forescating of SFS1 +with at least two dimensions (time , member) or (time, statistic). +The temporal dimension, by default 'time', must be equal to 1, the forecast year target. +The dimension 'member' must be greater than 1. +The dimension 'statistic' must be equal to 2, for containing the two paramenters of +a normal distribution (mean and sd) representing the ensemble of a SFS. +It is not possible to have the dimension 'member' and 'statistic' together.} + +\item{index_fcst2}{(optional, default = NULL) Index (e.g. NAO index) array from forescating of SFS2 +with at least two dimensions (time , member) or (time, statistic). +The temporal dimension, by default 'time', must be equal to 1, the forecast year target. +The dimension 'member' must be greater than 1. +The dimension 'statistic' must be equal to 2, for containing the two paramenters of +a normal distribution (mean and sd) representing the ensemble of a SFS. +It is not possible to have the dimension 'member' and 'statistic' together.} + +\item{method_BC}{A character vector of maximun length 2 indicating the bias +correction methodology to be applied on each SFS. If it is 'none' or any of +its elements is 'none', the bias correction won't be applied. +Available methods developped are "ME" (a bias correction scheme based on the +mean error or bias between observation and predictions to correct the +predicted values), and "LMEV" (a bias correction scheme based on a linear +model using ensemble variance of index as predictor). (see Sanchez-Garcia, +E. et al (2019), https://doi.org/10.5194/asr-16-165-2019 for more details).} + +\item{time_dim_name}{A character string indicating the name of the temporal +dimension, by default 'time'.} + +\item{na.rm}{Logical (default = FALSE). Should missing values be removed?} +} +\value{ +BEI_PDFBest() returns an array with the parameters that caracterize +the PDFs, with at least a temporal dimension, by default 'time' and dimension +'statistic' equal to 2. +The firt statistic is the parameter 'mean' of the PDF for the best estimation +combining the two SFSs PDFs. +The second statistic is the parameter 'standard deviation' of the PDF for +the best estimation combining the two SFSs PDFs. +If index_fcst1 and/or index_fcst2 are null, returns the values for hindcast period. +Otherwise, it returns the values for a forecast year. +} +\description{ +This function implements the computation to obtain the index +Probability Density Functions (PDFs) (e.g. NAO index) obtained to combining +the Index PDFs for two Seasonal Forecast Systems (SFSs), the Best Index +estimation (see Sanchez-Garcia, E. et al (2019), +https://doi.org/10.5194/asr-16-165-2019 for more details about the +methodology applied to estimate the Best Index). +} +\examples{ +# Example 1 for the BEI_PDFBest function +index_obs<- rnorm(10, sd = 3) +dim(index_obs) <- c(time = 5, season = 2) +index_hind1 <- rnorm(40, mean = 0.2, sd = 3) +dim(index_hind1) <- c(time = 5, member = 4, season = 2) +index_hind2 <- rnorm(60, mean = -0.5, sd = 4) +dim(index_hind2) <- c(time = 5, member = 6, season = 2) +index_fcst1 <- rnorm(16, mean = 0.2, sd = 3) +dim(index_fcst1) <- c(time = 1, member = 8, season = 2) +index_fcst2 <- rnorm(18, mean = -0.5, sd = 4) +dim(index_fcst2) <- c(time = 1, member = 9, season = 2) +method_BC <- 'ME' +res <- BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, +index_fcst2, method_BC) +dim(res) +# time statistic season +# 1 2 2 +# Example 2 for the BEI_PDFBest function +index_obs<- rnorm(10, sd = 3) +dim(index_obs) <- c(time = 5, season = 2) +index_hind1 <- rnorm(40, mean = 0.2, sd = 3) +dim(index_hind1) <- c(time = 5, member = 4, season = 2) +index_hind2 <- rnorm(60, mean = -0.5, sd = 4) +dim(index_hind2) <- c(time = 5, member = 6, season = 2) +index_fcst1 <- rnorm(16, mean = 0.2, sd = 3) +dim(index_fcst1) <- c(time = 1, member = 8, season = 2) +index_fcst2 <- rnorm(18, mean = -0.5, sd = 4) +dim(index_fcst2) <- c(time = 1, member = 9, season = 2) +method_BC <- c('LMEV', 'ME') +res <- BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, index_fcst2, method_BC) +dim(res) +# time statistic season +# 1 2 2 +} +\author{ +Eroteida Sanchez-Garcia - AEMET, \email{esanchezg@aemet.es} +} +\references{ +Regionally improved seasonal forecast of precipitation through +Best estimation of winter NAO, Sanchez-Garcia, E. et al., +Adv. Sci. Res., 16, 165174, 2019, https://doi.org/10.5194/asr-16-165-2019 +} + diff --git a/man/BEI_Weights.Rd b/man/BEI_Weights.Rd new file mode 100644 index 0000000000000000000000000000000000000000..61db33af3a5edcf6a140f35d7a5f3240320b663b --- /dev/null +++ b/man/BEI_Weights.Rd @@ -0,0 +1,55 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/BEI_Weights.R +\name{BEI_Weights} +\alias{BEI_Weights} +\title{Computing the weights for SFSs using the Best Index PDFs.} +\usage{ +BEI_Weights(index_weight, pdf_weight, time_dim_name = "time") +} +\arguments{ +\item{index_weight}{Index (e.g. NAO index) array, from a dataset of SFSs +for a period of years, with at least dimensions 'member'. +Additional dimensions, for instance, a temporal dimension as 'time', +must have the same lenght in both parameters, +'index_weight' and 'pdf_weight'.} + +\item{pdf_weight}{Statistics array to define a Gaussian PDF with at least +dimensions 'statistic'. +The firt statistic is the parameter 'mean' of the PDF and +the second statistic is the parameter 'standard deviation' of the PDF.} + +\item{time_dim_name}{A character string indicating the name of the temporal dimension, by default 'time'.} +} +\value{ +BEI_Weights() returns a normalized weights array with the same +dimensions that index_weight. +} +\description{ +This function implements the computation to obtain the +normalized weights for each member of each Seasonal Forecast Systems (SFS) +or dataset using the Probability Density Functions (PDFs) indicated by the +parameter 'pdf_weight' (for instance the Best Index estimation obtained +using the 'PDFBest' function). The weight of each member is proportional to +the probability of its index calculated with the PDF "pdf_weight". +} +\examples{ +# Example for the BEI_Weights function +index_weight <- 1 : (10 * 3 * 5 * 1) +dim(index_weight) <- c(sdate = 10, dataset = 3, member = 5, season = 1) +pdf_weight <- 1 : (10 * 3 * 2 * 1) +dim(pdf_weight) <- c(sdate = 10, dataset = 3, statistic = 2, season = 1) +res <- BEI_Weights(index_weight, pdf_weight) +dim(res) +# sdate dataset member season +# 10 3 5 1 + +} +\author{ +Eroteida Sanchez-Garcia - AEMET, \email{esanchezg@aemet.es} +} +\references{ +Regionally improved seasonal forecast of precipitation through Best +estimation of winter NAO, Sanchez-Garcia, E. et al., +Adv. Sci. Res., 16, 165174, 2019, https://doi.org/10.5194/asr-16-165-2019 +} + diff --git a/man/CST_BEI_Weighting.Rd b/man/CST_BEI_Weighting.Rd new file mode 100644 index 0000000000000000000000000000000000000000..6b9a448ae742037887d2e058e679a08edccf71d0 --- /dev/null +++ b/man/CST_BEI_Weighting.Rd @@ -0,0 +1,74 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_BEI_Weighting.R +\name{CST_BEI_Weighting} +\alias{CST_BEI_Weighting} +\title{Weighting SFSs of a CSTools object.} +\usage{ +CST_BEI_Weighting(var_exp, aweights, type = "ensembleMean", + time_dim_name = "time") +} +\arguments{ +\item{var_exp}{An object of the class 's2dv_cube' containing the variable +(e.g. precipitation, temperature, NAO index) array. +The var_exp object is expected to have an element named \code{$data} with +at least a temporal dimension and a dimension named 'member'.} + +\item{aweights}{Normalized weights array with at least dimensions +(time, member), when 'time' is the temporal dimension as default. +When 'aweights' parameter has any other dimensions (as e.g. 'lat') and +'var_exp' parameter has also the same dimension, they must be equals.} + +\item{type}{A character string indicating the type of output. +If 'type' = 'probs', the function returns, in the element data from +'var_exp' parameter, an array with at least two +or four dimensions depending if the variable is spatially aggregated variable +(as e.g. NAO index), dimension (time, tercil) or it is spatial variable +(as e.g. precipitation or temperature), dimension (time, tercile, lat, lon), +containing the terciles probabilities computing with weighted members. +The first tercil is the lower tercile, the second is the normal tercile and +the third is the upper tercile. +If 'type' = 'ensembleMean', the function returns, in the element data from +'var_exp' parameter, an array with at least one or three dimensions +depending if the variable is a spatially aggregated variable +(as e.g. NAO index)(time) or it is spatial variable (as e.g. precipitation +or temperature) (time, lat, lon), containing the ensemble means computing +with weighted members.} + +\item{time_dim_name}{A character string indicating the name of the +temporal dimension, by default 'time'.} +} +\value{ +CST_BEI_Weighting() returns a CSTools object (i.e., of the +class 's2dv_cube'). +This object has at least an element named \code{$data} +with at least a temporal dimension (and dimension 'tercil' when the output +are tercile probabilities), containing the ensemble means computing with +weighted members or probabilities of terciles. +} +\description{ +Function to apply weights to a 's2dv_cube' object. +It could return a weighted ensemble mean (deterministic output) or +the terciles probabilities (probabilistic output) for Seasonal Forecast +Systems (SFSs). +} +\examples{ +var_exp <- 1 : (2 * 4 * 3 * 2) +dim(var_exp) <- c(time = 2, member = 4, lat = 3, lon = 2) +aweights <- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3, 0.1, 0.2, 0.4, 0.4, 0.1, 0.2, 0.4, 0.2) +dim(aweights) <- c(time = 2, member = 4, dataset = 2) +var_exp <- list(data = var_exp) +class(var_exp) <- 's2dv_cube' +res_CST <- CST_BEI_Weighting(var_exp, aweights) +dim(res_CST$data) +# time lat lon dataset +# 2 3 2 2 +} +\author{ +Eroteida Sanchez-Garcia - AEMET, \email{esanchezg@aemet.es} +} +\references{ +Regionally improved seasonal forecast of precipitation through +Best estimation of winter NAO, Sanchez-Garcia, E. et al., +Adv. Sci. Res., 16, 165174, 2019, https://doi.org/10.5194/asr-16-165-2019 +} + diff --git a/man/CST_BiasCorrection.Rd b/man/CST_BiasCorrection.Rd index 75207e2e6e4db8e934dc8a51ba589e90fd64d81d..485199eae252039266185188df10ee21cf1af52b 100644 --- a/man/CST_BiasCorrection.Rd +++ b/man/CST_BiasCorrection.Rd @@ -42,4 +42,3 @@ Verónica Torralba, \email{veronica.torralba@bsc.es} 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) } \encoding{UTF-8} - diff --git a/tests/testthat/test-BEI_PDFBest.R b/tests/testthat/test-BEI_PDFBest.R new file mode 100644 index 0000000000000000000000000000000000000000..f2415f2e036db0a59f843751d35adbc12562ba32 --- /dev/null +++ b/tests/testthat/test-BEI_PDFBest.R @@ -0,0 +1,167 @@ +context("Generic tests") +test_that("basic use case", { + index_obs<- c(2, 3, 0.6) + dim(index_obs) <- c(time = 3) + index_hind1 <- c(-2.4, 3.6, -5, -1, -0.5, -0.2, 3, 0.6, 1, 0, -0.5, 2) + dim(index_hind1) <- c(time = 3, member = 4) + index_hind2 <- c(4, -10, 1.7, 1, -2, 3, -2, -1.5, 2.9) + dim(index_hind2) <- c(time = 3, member = 3) + index_fcst1 <- c(-2.6, -0.6, 1.4, -2.4, 8.1) + dim(index_fcst1) <- c(time = 1, member = 5) + index_fcst2 <- c(0, 4, -10, 2, -6, 3) + dim(index_fcst2) <- c(time = 1, member = 6) + + result <- array(c(2.574303, 0.8630711), + dim=c(time = 1, statistic = 2)) + + expect_equal(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = 'ME'), result, tolerance=1e-4) + + + result <- array(c(0.1282658, 0.4108197), + dim=c(time = 1, statistic = 2)) + + expect_equal(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = c('LMEV', 'ME')), result, + tolerance=1e-4) + + result <- array(c(0.4923944, 1.73306), + dim=c(time = 1, statistic = 2)) + + expect_equal(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = 'none'), result, + tolerance=1e-4) + + index_obs<- c(2, 3, 0.6) + dim(index_obs) <- c(time = 3) + index_hind1 <- c(-2.4, 3.6, -5, -1, -0.5, -0.2, 3, 0.6, 1, 0, -0.5, 2) + dim(index_hind1) <- c(time = 3, member = 4) + index_hind2 <- c(1.194390, 4.298768, -1.222586, 4.419386, 3.972776, 3.998928) + dim(index_hind2) <- c(time = 3, statistic=2) + index_fcst1 <- c(-2.6, -0.6, 1.4, -2.4, 8.1) + dim(index_fcst1) <- c(time = 1, member = 5) + index_fcst2 <- as.array(c(2.94, 4.08)) + dim(index_fcst2) <- c(time = 1, statistic = 2) + + result <- array(c(1.774329, 1.475816), + dim=c(time = 1, statistic = 2)) + + expect_equal(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = c('none',"ME")), result, + tolerance=1e-4) + + result <- array(c(0.2415013, 4.2232416, -0.9925062, 1.5508469, 0.9097119, + 1.5486845), dim=c(time = 3, statistic = 2)) + + expect_equal(BEI_PDFBest(index_obs, index_hind1, index_hind2, + method_BC = c('none',"ME")), result, + tolerance=1e-4) + +}) + +test_that("Sanity checks", { + index_obs<- c(2, 3, 0.6) + dim(index_obs) <- c(time = 3) + index_hind1 <- c(-2.4, 3.6, -5, -1, -0.5, -0.2, 3, 0.6, 1, 0, -0.5, 2) + dim(index_hind1) <- c(time = 3, member = 4) + index_hind2 <- c(1.194390, 4.298768, -1.222586, 4.419386, 3.972776, 3.998928) + dim(index_hind2) <- c(time = 3, statistic=2) + index_fcst1 <- c(-2.6, -0.6, 1.4, -2.4, 8.1) + dim(index_fcst1) <- c(time = 1, member = 5) + index_fcst2 <- as.array(c(2.94, 4.08)) + dim(index_fcst2) <- c(time = 1, statistic = 2) + expect_error(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = c('none',"ME"), na.rm = 2), + "Parameter 'na.rm' must be a logical value.") + + expect_error(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = c('none',"ME"), + time_dim_name = 2), + "Parameter 'time_dim_name' must be a character string indicating ", + "the name of the temporal dimension.") + + expect_error(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = list('none',"ME")), + paste0("Parameter 'method_BC' must be a character vector.")) + + expect_error(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = c('none','ME','ME')), + "Length of parameter 'method_BC' must be 1 or 2.") + + expect_error(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = c('none','HME')), + "Elements of parameter 'method_BC' must be equals to ", + "'none, 'ME' or 'LMEV'") + + expect_error(BEI_PDFBest(index_obs = 2, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = c('none','ME')), + "Parameter 'index_obs' must be an array.") + + expect_error(BEI_PDFBest(index_obs, index_hind1 = 3, index_hind2, index_fcst1, + index_fcst2, method_BC = c('none','ME')), + "Parameter 'index_hind1' must be an array.") + + index_hind1 <- 1:15 + dim(index_hind1) <- c(3, 5) + expect_error(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = c('none','ME')), + "Parameters 'index_obs', 'index_hind1' and 'index_hind2' ", + "should have dimmension names.") + + dim(index_hind1) <- c(pepe = 3, member = 5) + expect_error(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = c('none','ME')), + "Parameter 'index_hind1' must have temporal dimension.") + + index_hind1 <- 1:8 + dim(index_hind1) <- c(time = 2, member = 4) + expect_error(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = c('none','ME')), + "Length of temporal dimensions ", + "of parameter 'index_obs' and 'index_hind1' must be equals.") + + index_hind1 <- 1:12 + dim(index_hind1) <- c(time = 3, season = 4) + expect_error(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = c('none','ME')), + "Parameter 'index_hind1' must have dimension 'member' or 'statistic'") + + dim(index_hind1) <- c(time = 3, member = 4, statistic = 1) + expect_error(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = c('none','ME')), + "Parameter 'index_hind1' must have at least dimension 'member' ", + "or 'statistic', not 'member' and 'statistic' together.") + index_hind1 <- 1:9 + dim(index_hind1) <- c(time = 3, member = 1, season = 3) + expect_error(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = c('none','ME')), + "Length of dimension 'member' ", + "of parameter 'index_hind1' must be greater than 1.") + + + index_hind1 <- 1:12 + dim(index_hind1) <- c(time = 3, member = 4) + index_hind2 <- 1:9 + dim(index_hind2) <- c(time = 3, statistic=3) + expect_error(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = c('none','ME')), + "Length of dimension 'statistic' ", + "of parameter 'index_hind2' must be equal to 2.") + + index_hind2 <- 1:6 + dim(index_hind2) <- c(time = 3, statistic=2) + index_fcst1 <- 1:8 + dim(index_fcst1) <- c(time = 1, member = 8) + expect_error(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1 = 2, + index_fcst2, method_BC = c('none','ME')), + "Parameter 'index_fcst1' must be an array.") + + index_fcst1 <- 1:16 + dim(index_fcst1) <- c(time = 2, member = 8) + + expect_error(BEI_PDFBest(index_obs, index_hind1, index_hind2, index_fcst1, + index_fcst2, method_BC = c('none','ME')), + "Length of temporal dimension ", + "of parameter 'index_fcst1' must be equal to 1.") + +}) diff --git a/tests/testthat/test-BEI_Weights.R b/tests/testthat/test-BEI_Weights.R new file mode 100644 index 0000000000000000000000000000000000000000..fc234744bed4f93efe085ee069a400f5b622d5d5 --- /dev/null +++ b/tests/testthat/test-BEI_Weights.R @@ -0,0 +1,92 @@ +context("Generic tests") +test_that("basic use case", { + index_weight <- 1 : (3 * 3) + dim(index_weight) <- c(sdate = 3, member = 3) + pdf_weight <- 1 : (3 * 2) + dim(pdf_weight) <- c(sdate = 3, statistic = 2) + + result <- array(c(0.4808867, 0.4306591, 0.4017633, + 0.3629923, 0.3597167, 0.3545549, + 0.1561210, 0.2096243, 0.2436818), + dim = c(sdate = 3, member = 3)) + + expect_equal(BEI_Weights(index_weight, pdf_weight), result, tolerance=1e-4) + +}) + +test_that("Sanity checks", { + expect_error( + BEI_Weights(index_weight, pdf_weight, time_dim_name = 1), + "Parameter 'time_dim_name' must be a character string ", + "indicating the name of the temporal dimension.") + + expect_error( + BEI_Weights(index_weight = 1), + "Parameter 'index_weight' must be an array.") + + index_weight <- 1 : (3 * 3) + dim(index_weight) <- c(sdate = 3, member = 3) + expect_error( + BEI_Weights(index_weight, pdf_weight = 1), + "Parameter 'pdf_weight' must be an array.") + + index_weight <- 1 : (3 * 3) + dim(index_weight) <- c(3, 3) + pdf_weight <- 1 : (3 * 2) + dim(pdf_weight) <- c(3, 2) + + expect_error( + BEI_Weights(index_weight, pdf_weight), + "Parameters 'index_weight' and 'pdf_weight'", + " should have dimmension names.") + + index_weight <- 1 : (3 * 3) + dim(index_weight) <- c(time= 3, season = 3) + pdf_weight <- 1 : (3 * 2) + dim(pdf_weight) <- c(time = 3, statistic = 2) + expect_error( + BEI_Weights(index_weight, pdf_weight), + "Parameter 'index_weight' must have dimension 'member'.") + + index_weight <- 1 : (3 * 3) + dim(index_weight) <- c(time= 3, member = 3) + pdf_weight <- 1 : (3 * 2) + dim(pdf_weight) <- c(time = 3, season = 2) + expect_error( + BEI_Weights(index_weight, pdf_weight), + "Parameter 'pdf_weight' must have dimension 'statistic'.") + + index_weight <- 1 : (3 * 3) + dim(index_weight) <- c(time= 3, member = 3) + pdf_weight <- 1 : 2 + dim(pdf_weight) <- c(statistic = 2) + expect_error( + BEI_Weights(index_weight, pdf_weight), + "Parameter 'pdf_weight' must have temporal dimension.") + + index_weight <- 1 : 3 + dim(index_weight) <- c(member = 3) + pdf_weight <- 1 : 4 + dim(pdf_weight) <- c(time = 2, statistic = 2) + expect_error( + BEI_Weights(index_weight, pdf_weight), + "Parameter 'index_weight' must have temporal dimension.") + + index_weight <- 1 : 12 + dim(index_weight) <- c(time= 4, member = 3) + pdf_weight <- 1 : 4 + dim(pdf_weight) <- c(time = 2, statistic = 2) + expect_error( + BEI_Weights(index_weight, pdf_weight), + "Length of temporal dimension of parameters must be equal") + + index_weight <- 1 : (3 * 3) + dim(index_weight) <- c(sdate = 3, member = 3) + pdf_weight <- 1 : (3 * 3) + dim(pdf_weight) <- c(sdate = 3, statistic = 3) + expect_error( + BEI_Weights(index_weight, pdf_weight), + "Length of dimension 'statistic' ", + "of the parameter 'pdf_weight' must be equal to 2.") + +}) diff --git a/tests/testthat/test-CST_BEI_Weighting.R b/tests/testthat/test-CST_BEI_Weighting.R new file mode 100644 index 0000000000000000000000000000000000000000..48d968e5aaee4d9fb478f2d554fa664c7b1f7f61 --- /dev/null +++ b/tests/testthat/test-CST_BEI_Weighting.R @@ -0,0 +1,166 @@ +context("Generic tests") +test_that("basic use case", { + + var_exp <- 1 : (2 * 4 * 3 * 2) + dim(var_exp) <- c(time = 2, member = 4, lat = 3, lon = 2) + aweights <- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3, 0.1, 0.2, 0.4, 0.4, 0.1, 0.2, 0.4, 0.2) + dim(aweights) <- c(time = 2, member = 4, dataset = 2) + var_exp <- list(data = var_exp) + class(var_exp) <- 's2dv_cube' + + result <- array(c(4.4, 5.4, 12.4, 13.4, 20.4, 21.4, + 28.4, 29.4, 36.4, 37.4, 44.4, 45.4, + 4.6, 4.8, 12.6, 12.8, 20.6, 20.8, + 28.6, 28.8, 36.6, 36.8, 44.6, 44.8), + dim = c(time = 2, lat = 3, lon = 2, dataset =2)) + expect_equal(CST_BEI_Weighting(var_exp, aweights, type = 'ensembleMean', + time_dim_name = 'time')$data, result, tolerance=1e-4) + + var_exp <- 1 : (2 * 3 * 1 * 2) + dim(var_exp) <- c(time = 2, member = 3, lat = 1, lon = 2) + aweights <- c(0.2, 0.1, 0.3, 0.4, 0.5, 0.5) + dim(aweights) <- c(time = 2, member = 3) + var_exp <- list(data = var_exp) + class(var_exp) <- 's2dv_cube' + result <- array(c(0.5, 0.1, 0.0, 0.4, 0.5, 0.5, 0.5, 0.1, 0.0, 0.4, 0.5, 0.5), + dim = c(time = 2, tercil = 3, lat = 1, lon = 2)) + expect_equal(CST_BEI_Weighting(var_exp, aweights, type = 'probs', + time_dim_name = 'time')$data, result, tolerance=1e-4) + + var_exp <- 1 : (2 * 3 * 1 * 2) + dim(var_exp) <- c(sdate = 2, member = 3, lat = 1, lon = 2) + aweights <- c(0.2, 0.1, 0.3, 0.4, 0.5, 0.5) + dim(aweights) <- c(sdate = 2, member = 3) + var_exp <- list(data = var_exp) + class(var_exp) <- 's2dv_cube' + + result <- array(c(0.5, 0.1, 0.0, 0.4, 0.5, 0.5, + 0.5, 0.1, 0.0, 0.4, 0.5, 0.5), + dim = c(sdate = 2, tercil = 3, lat = 1, lon = 2)) + + expect_equal(CST_BEI_Weighting(var_exp, aweights, type = 'probs', + time_dim_name = 'sdate')$data, result, tolerance=1e-4) + + var_exp <- 1 : (2 * 3) + dim(var_exp) <- c(sdate = 2, member = 3) + aweights <- c(0.2, 0.1, 0.3, 0.4, 0.5, 0.5) + dim(aweights) <- c(sdate = 2, member = 3) + var_exp <- list(data = var_exp) + class(var_exp) <- 's2dv_cube' + + result <- array(c(0.5, 0.1, 0.0, 0.4, 0.5, 0.5), + dim = c(sdate = 2, tercil = 3)) + + expect_equal(CST_BEI_Weighting(var_exp, aweights, type = 'probs', + time_dim_name = 'sdate')$data, result, tolerance=1e-4) +}) + +test_that("Sanity checks", { + expect_error( + CST_BEI_Weighting(var_exp, aweights, type = 'probs', time_dim_name = 1), + "Parameter 'time_dim_name' must be a character string indicating", + " the name of the temporal dimension.") + + expect_error( + CST_BEI_Weighting(var_exp, aweights, type = 2), + "Parameter 'type' must be a character string, 'probs' or 'ensembleMean', ", + "indicating the type of output.") + + expect_error( + CST_BEI_Weighting(var_exp = 1, aweights), + "Parameter 'var_exp' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + + var_exp <- 1 : (2 * 3) + dim(var_exp) <- c(sdate = 2, member = 3) + var_exp <- list(data = var_exp) + class(var_exp) <- 's2dv_cube' + expect_error( + CST_BEI_Weighting(var_exp, aweights = 2), + "Parameter 'aweights' must be an array.") + + var_exp <- 1 : (2 * 3) + dim(var_exp) <- c(2, 3) + aweights <- c(0.2, 0.1, 0.3, 0.4, 0.5, 0.5) + dim(aweights) <- c(2, 3) + var_exp <- list(data = var_exp) + class(var_exp) <- 's2dv_cube' + + expect_error( + CST_BEI_Weighting(var_exp, aweights), + "Element 'data' from parameter 'var_exp' and parameter 'aweights'", + " should have dimmension names.") + + var_exp <- 1 : (2 * 3) + dim(var_exp) <- c(sdate = 2, member = 3) + aweights <- c(0.2, 0.1, 0.3, 0.4, 0.5, 0.5) + dim(aweights) <- c(sdate = 2, member = 3) + var_exp <- list(data = var_exp) + class(var_exp) <- 's2dv_cube' + + expect_error( + CST_BEI_Weighting(var_exp, aweights), + "Element 'data' from parameter 'var_exp' must have ", + "temporal dimension.") + + var_exp <- 1 : (2 * 3) + dim(var_exp) <- c(time = 2, member = 3) + aweights <- c(0.2, 0.1, 0.3, 0.4, 0.5, 0.5) + dim(aweights) <- c(sdate = 2, member = 3) + var_exp <- list(data = var_exp) + class(var_exp) <- 's2dv_cube' + + expect_error( + CST_BEI_Weighting(var_exp, aweights), + "Parameter 'aweights' must have temporal dimension.") + + var_exp <- 1 : (2 * 3) + dim(var_exp) <- c(time = 2, season = 3) + aweights <- c(0.2, 0.1, 0.3, 0.4, 0.5, 0.5) + dim(aweights) <- c(time = 2, member = 3) + var_exp <- list(data = var_exp) + class(var_exp) <- 's2dv_cube' + + expect_error( + CST_BEI_Weighting(var_exp, aweights), + "Element 'data' from parameter 'var_exp' must have ", + "dimension 'member'.") + + var_exp <- 1 : (2 * 3) + dim(var_exp) <- c(time = 2, member = 3) + aweights <- c(0.2, 0.1, 0.3, 0.4, 0.5, 0.5) + dim(aweights) <- c(time = 2, season = 3) + var_exp <- list(data = var_exp) + class(var_exp) <- 's2dv_cube' + + expect_error( + CST_BEI_Weighting(var_exp, aweights), + "Parameter 'aweights' must have ", + "dimension 'member'.") + + var_exp <- 1 : (3 * 3) + dim(var_exp) <- c(time = 3, member = 3) + aweights <- c(0.2, 0.1, 0.3, 0.4, 0.5, 0.5) + dim(aweights) <- c(time = 2, member = 3) + var_exp <- list(data = var_exp) + class(var_exp) <- 's2dv_cube' + + expect_error( + CST_BEI_Weighting(var_exp, aweights), + "Length of temporal dimensions ", + "of element 'data' from parameter 'var_exp' and parameter ", + "'aweights' must be equals.") + + var_exp <- 1 : (3 * 4) + dim(var_exp) <- c(time = 3, member = 4) + aweights <- c(0.2, 0.1, 0.3, 0.4, 0.5, 0.5) + dim(aweights) <- c(time = 2, member = 3) + var_exp <- list(data = var_exp) + class(var_exp) <- 's2dv_cube' + + expect_error( + CST_BEI_Weighting(var_exp, aweights), + "Length of temporal dimensions of element 'data' from ", + "parameter 'var_exp' and parameter 'aweights' must be equals.") + +}) diff --git a/vignettes/BestEstimateIndex_vignette.Rmd b/vignettes/BestEstimateIndex_vignette.Rmd new file mode 100644 index 0000000000000000000000000000000000000000..e577eeb2cffbfe417333d123d76b48c3f52d8179 --- /dev/null +++ b/vignettes/BestEstimateIndex_vignette.Rmd @@ -0,0 +1,134 @@ +--- +author: "Eroteida Sánchez-García" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteEngine{knitr::knitr} + %\VignetteIndexEntry{Achiving Best Estimate Index} + %\usepackage[utf8]{inputenc} +--- + +Achiving the Precipitation Best prediction giving the NAO index +------------------------------------------------------------------- + +The boreal winter precipitation forecast, acumulated from November to March, can be improved by considering the NAO index. The first step is to find the best estimation of winter NAO, given by two Seasonal Forecast System (SFS). The second step is to employ the enhanced NAO index pdf to produce weights for a SFS (it could be the same or a different SFS from the previous one). The third step is to apply this weights to a precipitation field. The methodology has been proved to improve the skill on precipitation forecast in the iberian peninsula given the relation between the winter precipitation and the NAO index in seasonal time scale (Sánchez-García, E., Voces-Aboy, J., Navascués, N., & Rodríguez-Camino, E. (2019). Regionally improved seasonal forecast of precipitation through Best estimation of winter NAO, Adv. Sci. Res., 16, 165174, ). + +This document is aim to ilustrate a practical use of the functions included in CSTools applying this methodology. + + +## Loading packages and data + +Open an R sesion and load library CSTools: + +``` +library(CSTools) +``` +The required data to applied this methodology are: +- the winter index NAO for two (or three) different SFS in a reference period (hindcast) and in a future simulation (forecast), +- the observed (reconstructed) NAO index in the reference period and +- the acumulated precipitation field from a SFS that aims to be improved (hindcast or hindcast and forecast) + +Given the memory limitations, the following example uses synthetic data. + +The first SFS is a dynamical model containing 25 ensemble members and its output will be save in the object `NAO_hind1` for the 20 years reference period and `NAO_fcst1` for a next season forecast. +The second SFS is a empirical model, so, the `NAO_hind2` and `NAO_fcst2` are characterized by a mean and standard deviation saved in the 'statistics' dimension. + +The synthetic data is created by running the following lines: + + +``` +NAO_obs <- rnorm(3, sd=3) +dim(NAO_obs) <- c(time = 3) +NAO_hind1 <- rnorm(3 * 6, sd=3) +dim(NAO_hind1) <- c(time = 3, member = 6) +NAO_hind2_mean <- rnorm(3, sd=3) +NAO_hind2_sd <- rnorm(3, mean=5, sd=1) +NAO_hind2 <- cbind(NAO_hind2_mean,NAO_hind2_sd) +dim(NAO_hind2) <- c(time=3, statistic=2) + +``` + +The acumulated precipiation field (only march for memory limiations) is loaded by running: + +``` +lonlat_prec$data <- apply(lonlat_prec$data, c(1, 2, 3, 5, 6), sum) +``` + + +## 1- Best Estimate Index NAO + +The function `PDFBest` performes the next steps: + - improves the PDF NAO index for each SFS (it could use bias correction method) and + - does a statistical combination of these improved NAO indexes. +Its output is an array containing the parameters (mean and starndar deviation) of the PDF for the reference period. + + +``` +pdf_hind_best <- BEI_PDFBest(NAO_obs, NAO_hind1, NAO_hind2, index_fcst1 = NULL, + index_fcst2 = NULL, method_BC = 'none', + time_dim_name = 'time', na.rm = FALSE) +``` + + +## 2- Compute weights using the Best Estimation of Index NAO + + +An array of weights is calculated for a SFS. This SFS could be a different SFS than the ones in section 1. +The function WeightIndex computes these weights for each ensemble's member based on the best NAO PDF estimate. + + +``` +weights_hind <- BEI_Weights(NAO_hind1, pdf_hind_best) +names(dim(weights_hind))[1] <- 'sdate' +``` + +The expected dimensions of these weights are 'member' and temporal dimension. + + +## 3- Apply weights to a precipitation field + + +The function `CST_BEI_Weighting` computes the ensemble mean or the terciles probabilities for a climate variable. +The corrected precipitation field is obtained by running: + +``` +em_prec <- CST_BEI_Weighting(lonlat_prec, weights_hind, type = 'ensembleMean', time_dim_name = 'sdate') +prob_prec <- CST_BEI_Weighting(lonlat_prec, weights_hind, type = 'probs', ) +``` + +## Comparison and visualization + +The original model output can be compared against the BEI corrected field. +To do this a equiprobability weights are created and they are applied to the original precipitation field: + + +``` +aweights_raw <- rep(1/6, 3 * 6) +dim(aweights_raw) <- dim(weights_hind) + +em_prec_raw <- CST_BEI_Weighting(lonlat_prec, aweights_raw, type = 'ensembleMean', time_dim_name = 'sdate') +prob_prec_raw <- CST_BEI_Weighting(lonlat_prec, aweights_raw, type = 'probs', time_dim_name = 'sdate') + +``` + +A map with the probability that the total precipitation will be in the lower/normal/upper tercile based on the Best Estimate Index NAO could be obtained using the 'PlotEquiMap' function or 'PlotMostLikelyQuantileMap' function from 'CSToools' package. + +The following figures show the probabilities for lower tercile for precipitation from November to March 2012/13 for ECMWF S5 system applying the methodology exposed or not, obtained using real data: +- NAO indices from the ECMWF-S5 dynamical model and the S-ClimWaRe empirical model from AEMET from 1997 to 2016 for computing the Best Estimation of Index NAO fro this hindcast period. +- The winter precipitation (from November to March) from 1997 to 2016 over Iberia Peninsula from he ECMWF-S5 dynamical model with resolution 0.5º x 0.5º, to weighting with the previous Best Estimation of Index NAO. + + +![](./Figures/BestEstimateIndex_fig1.png){width=70%} + +![](./Figures/BestEstimateIndex_fig3.png){width=70%} + + + +In a similar way we can plot the map with the probability that the total precipitation from November 2012 to +March 2013, for example, will be in the lower tercile from ECMWF Seasonal Forecast System 5 (raw) to compare results running the code: + + + +![](./Figures/BestEstimateIndex_fig2.png){width=70%} + +![](./Figures/BestEstimateIndex_fig4.png){width=70%} \ No newline at end of file diff --git a/vignettes/Figures/BestEstimateIndex_fig1.png b/vignettes/Figures/BestEstimateIndex_fig1.png new file mode 100644 index 0000000000000000000000000000000000000000..bd960203aaded56d21794650520c887339e0824e Binary files /dev/null and b/vignettes/Figures/BestEstimateIndex_fig1.png differ diff --git a/vignettes/Figures/BestEstimateIndex_fig2.png b/vignettes/Figures/BestEstimateIndex_fig2.png new file mode 100644 index 0000000000000000000000000000000000000000..e10333a577d1a2cba98d4efa702af3f6e8c14bab Binary files /dev/null and b/vignettes/Figures/BestEstimateIndex_fig2.png differ diff --git a/vignettes/Figures/BestEstimateIndex_fig3.png b/vignettes/Figures/BestEstimateIndex_fig3.png new file mode 100644 index 0000000000000000000000000000000000000000..ad43826f00abb0f204b549dff8f54bac8b3d94f2 Binary files /dev/null and b/vignettes/Figures/BestEstimateIndex_fig3.png differ diff --git a/vignettes/Figures/BestEstimateIndex_fig4.png b/vignettes/Figures/BestEstimateIndex_fig4.png new file mode 100644 index 0000000000000000000000000000000000000000..632d5e50cff8ff93b580cf6e03ea545216f83411 Binary files /dev/null and b/vignettes/Figures/BestEstimateIndex_fig4.png differ