diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 866209853288e1301b4b6787c5f70af2055b7249..3017964c5df4c06ad9b391c0bb329e6948075805 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -4,7 +4,7 @@ build: stage: build script: #- module load R/3.4.3-foss-2015a-bare Julia - - module load R + - module load R/3.4.2-foss-2015a-bare #- julia -e 'using Pkg; Pkg.add(["RainFARM"])' - R CMD build --resave-data . - R CMD check --as-cran --no-manual --run-donttest CSTools_*.tar.gz diff --git a/DESCRIPTION b/DESCRIPTION index 0e47736a9c3d0ed4e22f81e732b21162c3c1d622..2c87e1472aeed292a148f0e89adbc82e9fc0453b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: CSTools Title: Assessing Skill of Climate Forecasts on Seasonal-to-Decadal Timescales -Version: 2.0.0 +Version: 3.0.0 Authors@R: c( person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-8568-3071")), person("Louis-Philippe", "Caron", , "louis-philippe.caron@bsc.es", role = "aut", comment = c(ORCID = "0000-0001-5221-0147")), @@ -42,7 +42,7 @@ Description: Exploits dynamical seasonal forecasts in order to provide Van Schaeybroeck et al. (2019) . Yiou et al. (2013) . Depends: - R (>= 3.2.0), + R (>= 3.4.2), maps Imports: s2dverification, @@ -70,4 +70,4 @@ VignetteBuilder: knitr License: Apache License 2.0 Encoding: UTF-8 LazyData: true -RoxygenNote: 5.0.0 +RoxygenNote: 7.0.2 diff --git a/NAMESPACE b/NAMESPACE index e7d7c003b822cbb73b70d26b94fb43d55faab57e..bd5d0f1649b1d8529ffb6d16118bc8c87ec186ab 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -11,6 +11,7 @@ export(CST_Calibration) export(CST_CategoricalEnsCombination) export(CST_EnsClustering) export(CST_Load) +export(CST_MergeDims) export(CST_MultiEOF) export(CST_MultiMetric) export(CST_MultivarRMSE) @@ -21,6 +22,7 @@ export(CST_RainFARM) export(CST_SaveExp) export(CST_SplitDim) export(EnsClustering) +export(MergeDims) export(MultiEOF) export(PlotCombinedMap) export(PlotForecastPDF) diff --git a/NEWS.md b/NEWS.md index 7d1ea1e5bdf9c09a05625ac3a06e97e0b60346c5..3e625ed09e2933e7f4255ceef9fc1ad89e8af161 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,14 @@ +### CSTools 3.0.0 +**Submission date to CRAN: 10-02-2020** + +- New features: + + CST_MergeDims and MergeDims + + Version working with R 3.4.2 + + PlotForecastPDF handles independent terciles, extremes and observations for each panel +- Fixes + + CST_Calibration handles missing values + + BEI functions handle missing values + ### CSTools 2.0.0 **Submission date to CRAN: 25-11-2019** diff --git a/R/CST_BEI_Weighting.R b/R/CST_BEI_Weighting.R index 51c8ba32fdaede17917f610ecd685aa6991ebcda..de7470110c4f9d090bc1ad84b6c786448a5c63c4 100644 --- a/R/CST_BEI_Weighting.R +++ b/R/CST_BEI_Weighting.R @@ -19,6 +19,10 @@ #' (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 terciles A numeric array with at least one dimension 'tercil' equal to +#' 2, the first element is the lower tercil for a hindcast period, and the second +#' element is the upper tercile. By default is NULL, the terciles are computed +#' from var_exp data. #' @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 @@ -36,7 +40,7 @@ #' 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} @@ -56,8 +60,10 @@ #' # time lat lon dataset #' # 2 3 2 2 #' @export -CST_BEI_Weighting <- function(var_exp, aweights, type = 'ensembleMean', - time_dim_name = 'time') { + +CST_BEI_Weighting <- function(var_exp, aweights, terciles = NULL, + 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.") @@ -80,6 +86,27 @@ CST_BEI_Weighting <- function(var_exp, aweights, type = 'ensembleMean', stop("Parameter 'var_exp' must be of the class 's2dv_cube', ", "as output by CSTools::CST_Load.") } + if (!is.null(terciles)){ + if(!is.array(terciles)){ + stop("Parameter 'terciles' must be an array.") + } + if (is.null(names(dim(terciles)))) { + stop("Parameters 'terciles' should have dimmension names.") + } + if(!('tercil' %in% names(dim(terciles)))) { + stop("Parameter 'terciles' must have dimension 'tercil'.") + } + if (dim(terciles)['tercil'] != 2) { + stop("Length of dimension 'tercil' ", + "of parameter 'terciles' must be equal to 2.") + } + if(time_dim_name %in% names(dim(terciles))) { + stop("Parameter 'terciles' must not have temporal dimension.") + } + if('member' %in% names(dim(terciles))) { + stop("Parameter 'terciles' must not have dimension 'member'.") + } + } if (!is.array(aweights)) { stop("Parameter 'aweights' must be an array.") } @@ -115,7 +142,10 @@ CST_BEI_Weighting <- function(var_exp, aweights, 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) + if (is.null(terciles)){ + terciles <- BEI_TercilesWeighting(var_exp$data, aweights, time_dim_name) + } + probs <- BEI_ProbsWeighting(var_exp$data, aweights, terciles, time_dim_name) var_exp$data <- probs } else { stop("Parameter 'type' must be a character string ('probs' or ", @@ -265,6 +295,9 @@ BEI_EMWeighting <- function(var_exp, aweights, time_dim_name = 'time') { #' 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 terciles A numeric array with at least one dimension 'tercil' equal to +#' 2, the first element is the lower tercil for a hindcast period, and the second +#' element is the upper tercile. #' @param time_dim_name A character string indicating the name of the #' temporal dimension, by default 'time'. #' @@ -284,7 +317,9 @@ BEI_EMWeighting <- function(var_exp, aweights, time_dim_name = 'time') { #' 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) +#' terciles <- c(2.5,5) +#' dim(terciles) <- c(tercil = 2) +#' res <- BEI_ProbsWeighting(var_exp, aweights, terciles) #' dim(res) #' # time tercil #' # 2 3 @@ -293,12 +328,15 @@ BEI_EMWeighting <- function(var_exp, aweights, time_dim_name = 'time') { #' 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) +#' terciles <- rep(c(48,50), 2*3) +#' dim(terciles) <- c(tercil = 2, lat = 2, lon = 3) +#' res <- BEI_ProbsWeighting(var_exp, aweights, terciles) #' dim(res) #' # time tercil lat lon #' # 2 3 2 3 #' @noRd -BEI_ProbsWeighting <- function(var_exp, aweights, time_dim_name = 'time') { +BEI_ProbsWeighting <- function(var_exp, aweights, terciles, + time_dim_name = 'time') { if (!is.character(time_dim_name)) { stop("Parameter 'time_dim_name' must be a character string indicating", @@ -309,6 +347,28 @@ BEI_ProbsWeighting <- function(var_exp, aweights, time_dim_name = 'time') { "only the first element will be used.") time_dim_name <- time_dim_name[1] } + if (is.null(terciles)){ + stop("Parameter 'terciles' is null") + } + if(!is.array(terciles)){ + stop("Parameter 'terciles' must be an array.") + } + if (is.null(names(dim(terciles)))) { + stop("Parameters 'terciles' should have dimmension names.") + } + if(!('tercil' %in% names(dim(terciles)))) { + stop("Parameter 'terciles' must have dimension 'tercil'.") + } + if (dim(terciles)['tercil'] != 2) { + stop("Length of dimension 'tercil' ", + "of parameter 'terciles' must be equal to 2.") + } + if(time_dim_name %in% names(dim(terciles))) { + stop("Parameter 'terciles' must not have temporal dimension.") + } + if('member' %in% names(dim(terciles))) { + stop("Parameter 'terciles' must not have dimension 'member'.") + } if (!is.array(var_exp)) { stop("Parameter 'var_exp' must be an array.") } @@ -326,10 +386,10 @@ BEI_ProbsWeighting <- function(var_exp, aweights, time_dim_name = 'time') { stop("Parameter 'aweights' must have temporal dimension.") } if(!('member' %in% names(dim(var_exp)))) { - stop("Parameter 'var_exp' must have temporal dimension.") + stop("Parameter 'var_exp' must have dimension 'member'.") } if(!('member' %in% names(dim(aweights)))) { - stop("Parameter 'aweights' must have temporal dimension.") + stop("Parameter 'aweights' must have dimension 'member'.") } if (dim(var_exp)[time_dim_name] != dim(aweights)[time_dim_name]) { stop("Length of temporal dimensions ", @@ -340,8 +400,11 @@ BEI_ProbsWeighting <- function(var_exp, aweights, time_dim_name = 'time') { "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')), + + res <- Apply(list(var_exp, aweights, terciles), + target_dims = list(c(time_dim_name,'member'), + c(time_dim_name,'member'), + c('tercil')), fun = .BEI_ProbsWeighting, time_dim_name)$output1 return(res) } @@ -352,8 +415,12 @@ BEI_ProbsWeighting <- function(var_exp, aweights, time_dim_name = 'time') { #' by default 'time', and dimension 'member'. #' @param aweights Normalized weights array with a temporal dimension, #' by default 'time', and dimension 'member' +#' @param terciles A numeric array with one dimension 'tercil' equal to 2, +#' the first element is the lower tercil for a hindcast period, and the second +#' element is the upper tercile. #' @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. @@ -363,48 +430,207 @@ BEI_ProbsWeighting <- function(var_exp, aweights, time_dim_name = 'time') { #' # 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) +#' aweights <- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3) #' dim(aweights) <- c(stime = 2, member = 4) -#' res <- .BEI_ProbsWeighting(var_exp, aweights, time_dim_name = 'stime') +#' terciles <- quantile(1:8, probs = c(1/3, 2/3)) +#' dim(terciles) <- c(tercil = 2) +#' res <- .BEI_ProbsWeighting(var_exp, aweights, terciles, 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 +.BEI_ProbsWeighting <- function(var_exp, aweights, terciles, + time_dim_name = 'time') { + if(any(is.na(var_exp)) || any(is.na(aweights))){ + probTercile <- array(NA, dim = c(dim(var_exp)[time_dim_name], tercil = 3)) + } else { + if(any(is.na(terciles))) stop("Terciles are NAs") + terciles_exp <- list(lowerTercile = terciles[1], + upperTercile = terciles[2]) + + 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) +} + +#' Computing the weighted terciles for SFSs. +#' @author Eroteida Sanchez-Garcia - AEMET, \email{esanchezg@aemet.es} +#' +#' @description This function implements the computation to obtain the terciles +#' 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_TercilesWeighting() returns an array with at least one +#' dimension depending if the variable is a spatially aggregated variable +#' (as e.g. NAO index)(tercil) or it is spatial variable (as e.g. +#' precipitation or temperature)(tercil, lat, lon), containing the +#' terciles computing with weighted members. +#' The first tercil is the lower tercile, the second 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_TercilesWeighting(var_exp, aweights) +#' dim(res) +#' # tercil +#' # 2 +#' # 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_TercilesWeighting(var_exp, aweights) +#' dim(res) +#' # tercil lat lon +#' # 2 2 3 +#' @noRd +BEI_TercilesWeighting <- function(var_exp, aweights, time_dim_name = 'time') { - pos <- match(names(dim(aTerciles)), c(time_dim_name,'member')) - aTerciles <- aperm(aTerciles,pos) - names(dim(aTerciles)) <- c(time_dim_name,'member') + 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.") + } - 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 + res <- Apply(list(var_exp, aweights), + target_dims = list(c(time_dim_name,'member'), c(time_dim_name,'member')), + fun = .BEI_TercilesWeighting, time_dim_name)$output1 + return(res) +} + +#' Atomic BEI_TercilesWeighting +#' @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_TercilesWeighting returns a numeric array with dimension tercil +#' equal to 2, the first is the lower tercil and the second the upper tercile, +#' computing with weighted members considering all members and all period. +#' If any member value for any period is NA , the terciles are not computed, and +#' the function return NA value as tercile upper and lower. +#' @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.4, 0.3) +#' dim(aweights) <- c(stime = 2, member = 4) +#' res <- .BEI_TercilesWeighting(var_exp, aweights, time_dim_name = 'stime') +#' dim(res) +#' # tercil +#' # 2 +#' @noRd +.BEI_TercilesWeighting <- function(var_exp, aweights, time_dim_name = 'time') { + if(any(is.na(var_exp)) || any(is.na(aweights))){ + terciles_exp <- array(c(NA, NA), dim = c(tercil = 2)) + } else { + l_terciles_exp <- WeightTerciles(var_exp, aweights, time_dim_name) + terciles_exp <- array(c(l_terciles_exp$lowerTercile, + l_terciles_exp$upperTercile), dim = c(tercil = 2)) } - return(probTercile) + return(terciles_exp) } +# Auxiliar function to compute in which tercile is a data value +Data2Tercil_old <- function(x,lt,ut) { + if(is.na(lt) || is.na(ut)){ + y <- rep(NA, length(x)) + } else { + 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 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") + if(is.na(lt) || is.na(ut)){ + y <- rep(NA, length(x)) + } else { + y <- rep(2,length(x)) + y[x <= lt] <- 1 + y[x >= ut] <- 3 + if (lt == ut) { + y <- rep(NA, length(x)) + } } dim(y) <- c(member = length(x)) + y[which(is.na(x))] <- NA return (y) } - # Auxiliar function to convers weighted terciles to probabilities WeightTercil2Prob <- function(aTerciles, aWeights, idTercil) { return(sum(aWeights[which(aTerciles == idTercil)])) @@ -419,10 +645,12 @@ WeightTerciles <- function(data, aweights, time_dim_name = 'time') { 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 + #lSortData <- sort(vectorData,index.return=TRUE) + indSort <- order(vectorData) # index asociated to weight + # indSort <- lSortData$ix # index asociated to weight # corresponding for this data - dataSort <- lSortData$x + dataSort <- vectorData[indSort] + # 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 diff --git a/R/CST_Calibration.R b/R/CST_Calibration.R index 1cd9bab453be826fd54e49150cd4e136d0ed368b..a7f1924b208663909a774c94e0096f2e3d824428 100644 --- a/R/CST_Calibration.R +++ b/R/CST_Calibration.R @@ -2,27 +2,24 @@ #' #'@author Verónica Torralba, \email{veronica.torralba@bsc.es} #'@author Bert Van Schaeybroeck, \email{bertvs@meteo.be} -#'@description Four types of member-by-member bias correction can be performed. The \code{bias} method corrects the bias only, the \code{evmos} method applies a variance inflation technique to ensure the correction of the bias and the correspondence of variance between forecast and observation (Van Schaeybroeck and Vannitsem, 2011). The ensemble calibration methods \code{"mse_min"} and \code{"crps_min"} correct the bias, the overall forecast variance and the ensemble spread as described in Doblas-Reyes et al. (2005) and Van Schaeybroeck and Vannitsem (2015), respectively. While the \code{"mse_min"} method minimizes a constrained mean-squared error using three parameters, the \code{"crps_min"} method features four parameters and minimizes the Continuous Ranked Probability Score (CRPS). -#'@description Both in-sample or our out-of-sample (leave-one-out cross validation) calibration are possible. -#'@references Doblas-Reyes F.J, Hagedorn R, Palmer T.N. The rationale behind the success of multi-model ensembles in seasonal forecasting-II calibration and combination. Tellus A. 2005;57:234-252. doi:10.1111/j.1600-0870.2005.00104.x -#'@references Van Schaeybroeck, B., & Vannitsem, S. (2011). Post-processing through linear regression. Nonlinear Processes in Geophysics, 18(2), 147. doi:10.5194/npg-18-147-2011 -#'@references Van Schaeybroeck, B., & Vannitsem, S. (2015). Ensemble post-processing using member-by-member approaches: theoretical aspects. Quarterly Journal of the Royal Meteorological Society, 141(688), 807-818. doi:10.1002/qj.2397 +#'@description Equivalent to function \code{Calibration} but for objects of class \code{s2dv_cube}. #' #'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}. #'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed data in the element named \code{$data}. #'@param cal.method is the calibration method used, can be either \code{bias}, \code{evmos}, \code{mse_min} or \code{crps_min}. Default value is \code{mse_min}. #'@param eval.method is the sampling method used, can be either \code{in-sample} or \code{leave-one-out}. Default value is the \code{leave-one-out} cross validation. #'@param multi.model is a boolean that is used only for the \code{mse_min} method. If multi-model ensembles or ensembles of different sizes are used, it must be set to \code{TRUE}. By default it is \code{FALSE}. Differences between the two approaches are generally small but may become large when using small ensemble sizes. Using multi.model when the calibration method is \code{bias}, \code{evmos} or \code{crps_min} will not affect the result. -#'@return an object of class \code{s2dv_cube} containing the calibrated forecasts in the element \code{$data} with the same dimensions of the experimental data. +#'@param na.fill is a boolean that indicates what happens in case calibration is not possible or will yield unreliable results. This happens when three or less forecasts-observation pairs are available to perform the training phase of the calibration. By default \code{na.fill} is set to true such that NA values will be returned. If \code{na.fill} is set to false, the uncorrected data will be returned. +#'@param ncores is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one. +#'@return an object of class \code{s2dv_cube} containing the calibrated forecasts in the element \code{$data} with the same dimensions as the one in the exp object. #' #'@import s2dverification #'@import abind #' #'@seealso \code{\link{CST_Load}} #' -#'# Example -#'# Creation of sample s2dverification objects. These are not complete -#'# s2dverification objects though. The Load function returns complete objects. +#'@examples +#'# Example 1: #'mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7) #'dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7) #'obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7) @@ -37,8 +34,12 @@ #'str(a) #'@export - -CST_Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-one-out", multi.model = F) { +CST_Calibration <- function(exp, obs, + cal.method = "mse_min", + eval.method = "leave-one-out", + multi.model = F, + na.fill = T, + ncores = 1) { if (!inherits(exp, "s2dv_cube") || !inherits(obs, "s2dv_cube")) { stop("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ", "as output by CSTools::CST_Load.") @@ -52,9 +53,12 @@ CST_Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "lea if(!missing(multi.model) & !(cal.method == "mse_min")){ warning(paste0("The multi.model parameter is ignored when using the calibration method ", cal.method)) } - exp$data <- .calibration.wrap(exp = exp$data, obs = obs$data, - cal.method = cal.method, eval.method = eval.method, - multi.model = multi.model) + exp$data <- Calibration(exp = exp$data, obs = obs$data, + cal.method = cal.method, + eval.method = eval.method, + multi.model = multi.model, + na.fill = na.fill, + ncores = ncores) exp$Datasets <- c(exp$Datasets, obs$Datasets) exp$source_files <- c(exp$source_files, obs$source_files) @@ -62,7 +66,39 @@ CST_Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "lea } -.calibration.wrap <- function(exp, obs, cal.method, eval.method, multi.model) { + + +#'Forecast Calibration +#' +#'@author Verónica Torralba, \email{veronica.torralba@bsc.es} +#'@author Bert Van Schaeybroeck, \email{bertvs@meteo.be} +#'@description Four types of member-by-member bias correction can be performed. The \code{bias} method corrects the bias only, the \code{evmos} method applies a variance inflation technique to ensure the correction of the bias and the correspondence of variance between forecast and observation (Van Schaeybroeck and Vannitsem, 2011). The ensemble calibration methods \code{"mse_min"} and \code{"crps_min"} correct the bias, the overall forecast variance and the ensemble spread as described in Doblas-Reyes et al. (2005) and Van Schaeybroeck and Vannitsem (2015), respectively. While the \code{"mse_min"} method minimizes a constrained mean-squared error using three parameters, the \code{"crps_min"} method features four parameters and minimizes the Continuous Ranked Probability Score (CRPS). +#'@description Both in-sample or our out-of-sample (leave-one-out cross validation) calibration are possible. +#'@references Doblas-Reyes F.J, Hagedorn R, Palmer T.N. The rationale behind the success of multi-model ensembles in seasonal forecasting-II calibration and combination. Tellus A. 2005;57:234-252. doi:10.1111/j.1600-0870.2005.00104.x +#'@references Van Schaeybroeck, B., & Vannitsem, S. (2011). Post-processing through linear regression. Nonlinear Processes in Geophysics, 18(2), 147. doi:10.5194/npg-18-147-2011 +#'@references Van Schaeybroeck, B., & Vannitsem, S. (2015). Ensemble post-processing using member-by-member approaches: theoretical aspects. Quarterly Journal of the Royal Meteorological Society, 141(688), 807-818. doi:10.1002/qj.2397 +#' +#'@param exp an array containing the seasonal forecast experiment data. +#'@param obs an array containing the observed data. +#'@param cal.method is the calibration method used, can be either \code{bias}, \code{evmos}, \code{mse_min} or \code{crps_min}. Default value is \code{mse_min}. +#'@param eval.method is the sampling method used, can be either \code{in-sample} or \code{leave-one-out}. Default value is the \code{leave-one-out} cross validation. +#'@param multi.model is a boolean that is used only for the \code{mse_min} method. If multi-model ensembles or ensembles of different sizes are used, it must be set to \code{TRUE}. By default it is \code{FALSE}. Differences between the two approaches are generally small but may become large when using small ensemble sizes. Using multi.model when the calibration method is \code{bias}, \code{evmos} or \code{crps_min} will not affect the result. +#'@param na.fill is a boolean that indicates what happens in case calibration is not possible or will yield unreliable results. This happens when three or less forecasts-observation pairs are available to perform the training phase of the calibration. By default \code{na.fill} is set to true such that NA values will be returned. If \code{na.fill} is set to false, the uncorrected data will be returned. +#'@param ncores is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one. +#'@return an array containing the calibrated forecasts with the same dimensions as the \code{exp} array. +#' +#'@import s2dverification +#'@import abind +#' +#'@seealso \code{\link{CST_Load}} +#' + +Calibration <- function(exp, obs, + cal.method = "mse_min", + eval.method = "leave-one-out", + multi.model = F, + na.fill = T, + ncores = 1) { dim.exp <- dim(exp) amt.dims.exp <- length(dim.exp) @@ -93,14 +129,30 @@ CST_Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "lea if (dim(obs)['member']!=1){ stop("Parameter 'obs' must have a member dimension with length=1") } - + obs <- Subset(obs, along = "member", indices = 1, drop = "selected") + + data.set.sufficiently.large.out <- + Apply(data = list(exp = exp, obs = obs), + target_dims = list(exp = target.dim.names.exp, obs = target.dim.names.obs), + ncores = ncores, + fun = .data.set.sufficiently.large)$output1 - calibrated <- Apply(data = list(mod = exp, obs = obs), + if(!all(data.set.sufficiently.large.out)){ + if(na.fill){ + warning("Some forecast data could not be corrected due to data lack and is replaced with NA values") + } else { + warning("Some forecast data could not be corrected due to data lack and is replaced with uncorrected values") + } + } + + calibrated <- Apply(data = list(exp = exp, obs = obs), cal.method = cal.method, eval.method = eval.method, multi.model = multi.model, - target_dims = list(mod = target.dim.names.exp, obs = target.dim.names.obs), + na.fill = na.fill, + target_dims = list(exp = target.dim.names.exp, obs = target.dim.names.obs), + ncores = ncores, fun = .cal)$output1 dexes <- match(names(dim(exp)), names(dim(calibrated))) @@ -111,6 +163,13 @@ CST_Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "lea return(calibrated) } + +.data.set.sufficiently.large <- function(exp, obs){ + amt.min.samples <- 3 + amt.good.pts <- sum(!is.na(obs) & !apply(exp, c(2), function(x) all(is.na(x)))) + return(amt.good.pts > amt.min.samples) +} + .make.eval.train.dexes <- function(eval.method, amt.points){ if(eval.method == "leave-one-out"){ dexes.lst <- lapply(seq(1, amt.points), function(x) return(list(eval.dexes = x, train.dexes = seq(1, amt.points)[-x]))) @@ -122,14 +181,23 @@ CST_Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "lea return(dexes.lst) } -.cal <- function(mod, obs, cal.method, eval.method, multi.model) { +.cal <- function(exp, obs, cal.method, eval.method, multi.model, na.fill) { - var.obs <- as.vector(obs) - var.fc <- mod - dims.fc <- dim(var.fc) + obs <- as.vector(obs) + dims.fc <- dim(exp) amt.mbr <- dims.fc[1] amt.sdate <- dims.fc[2] - var.cor.fc <- NA * var.fc + var.cor.fc <- NA * exp + names(dim(var.cor.fc)) <- c("member", "sdate") + + if(!.data.set.sufficiently.large(exp = exp, obs = obs)){ + if(na.fill){ + return(var.cor.fc) + } else { + var.cor.fc[] <- exp[] + return(var.cor.fc) + } + } eval.train.dexeses <- .make.eval.train.dexes(eval.method, amt.points = amt.sdate) amt.resamples <- length(eval.train.dexeses) @@ -138,9 +206,9 @@ CST_Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "lea eval.dexes <- eval.train.dexeses[[i.sample]]$eval.dexes train.dexes <- eval.train.dexeses[[i.sample]]$train.dexes - fc.ev <- var.fc[ , eval.dexes, drop = FALSE] - fc.tr <- var.fc[ , train.dexes] - obs.tr <- var.obs[train.dexes , drop = FALSE] + fc.ev <- exp[ , eval.dexes, drop = FALSE] + fc.tr <- exp[ , train.dexes] + obs.tr <- obs[train.dexes , drop = FALSE] if(cal.method == "bias"){ var.cor.fc[ , eval.dexes] <- fc.ev + mean(obs.tr, na.rm = TRUE) - mean(fc.tr, na.rm = TRUE) @@ -165,8 +233,11 @@ CST_Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "lea init.par <- c(.calc.mse.min.par(quant.obs.fc.tr), 0.001) init.par[3] <- sqrt(init.par[3]) #calculate regression parameters on training dataset - optim.tmp <- optim(par = init.par, fn = .calc.crps.opt, - quant.obs.fc = quant.obs.fc.tr, method = "Nelder-Mead") + optim.tmp <- optim(par = init.par, + fn = .calc.crps.opt, + gr = .calc.crps.grad.opt, + quant.obs.fc = quant.obs.fc.tr, + method = "BFGS") mbm.par <- optim.tmp$par #correct evaluation subset @@ -175,7 +246,6 @@ CST_Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "lea stop("unknown calibration method: ",cal.method) } } - names(dim(var.cor.fc)) <- c("member", "sdate") return(var.cor.fc) } @@ -294,6 +364,23 @@ CST_Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "lea ) } +.calc.crps.grad.opt <- function(par, quant.obs.fc){ + sgn1 <- with(quant.obs.fc,sign(obs.per.ens - + (par[1] + par[2] * fc.ens.av.per.ens + + ((par[3])^2 + par[4] / spr.abs.per.ens) * fc.dev))) + sgn2 <- with(quant.obs.fc, sign((par[3])^2 + par[4] / spr.abs.per.ens)) + sgn3 <- with(quant.obs.fc,sign((par[3])^2 * spr.abs + par[4])) + deriv.par1 <- mean(sgn1, na.rm = TRUE) + deriv.par2 <- with(quant.obs.fc, mean(sgn1 * fc.dev, na.rm = TRUE)) + deriv.par3 <- with(quant.obs.fc, + mean(2* par[3] * sgn1 * sgn2 * fc.ens.av.per.ens, na.rm = TRUE) - + mean(spr.abs * sgn3, na.rm = TRUE) / 2.) + deriv.par4 <- with(quant.obs.fc, + mean(sgn1 * sgn2 * fc.ens.av.per.ens / spr.abs.per.ens, na.rm = TRUE) - + mean(sgn3, na.rm = TRUE) / 2.) + return(c(deriv.par1, deriv.par2, deriv.par3, deriv.par4)) +} + #Below are the core or elementary functions to correct the evaluation set based on the regression parameters .correct.evmos.fc <- function(fc, par){ quant.fc.mp <- .calc.fc.quant(fc = fc) diff --git a/R/CST_MergeDims.R b/R/CST_MergeDims.R new file mode 100644 index 0000000000000000000000000000000000000000..a9923c581b95e6c48c37d2efa2e32d7588a947f0 --- /dev/null +++ b/R/CST_MergeDims.R @@ -0,0 +1,127 @@ +#'Function to Merge Dimensions +#' +#'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} +#' +#'@description This function merges two dimensions of the array \code{data} in a 's2dv_cube' object into one. The user can select the dimensions to merge and provide the final name of the dimension. The user can select to remove NA values or keep them. +#' +#'@param data a 's2dv_cube' object +#'@param merge_dims a character vector indicating the names of the dimensions to merge +#'@param rename_dim a character string indicating the name of the output dimension. If left at NULL, the first dimension name provided in parameter \code{merge_dims} will be used. +#'@param na.rm a logical indicating if the NA values should be removed or not. +#' +#'@import abind +#'@import s2dverification +#'@examples +#' +#'data <- 1 : c(2 * 3 * 4 * 5 * 6 * 7) +#'dim(data) <- c(time = 7, lat = 2, lon = 3, monthly = 4, member = 6, +#' dataset = 5, var = 1) +#'data[2,,,,,,] <- NA +#'data[c(3,27)] <- NA +#'data <-list(data = data) +#'class(data) <- 's2dv_cube' +#'new_data <- CST_MergeDims(data, merge_dims = c('time', 'monthly')) +#'dim(new_data$data) +#'new_data <- CST_MergeDims(data, merge_dims = c('lon', 'lat'), rename_dim = 'grid') +#'dim(new_data$data) +#'new_data <- CST_MergeDims(data, merge_dims = c('time', 'monthly'), na.rm = TRUE) +#'dim(new_data$data) +#'@export +CST_MergeDims <- function(data, merge_dims = c('ftime', 'monthly'), rename_dim = NULL, + na.rm = FALSE) { + if (!inherits(data, 's2dv_cube')) { + stop("Parameter 'data' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + data$data <- MergeDims(data$data, merge_dims = merge_dims, + rename_dim = rename_dim, na.rm = na.rm) + return(data) +} +#'Function to Split Dimension +#' +#'@author Nuria Perez-Zanon, \email{nuria.perez@bsc.es} +#' +#'@description This function merges two dimensions of an array into one. The user can select the dimensions to merge and provide the final name of the dimension. The user can select to remove NA values or keep them. +#' +#'@param data an n-dimensional array with named dimensions +#'@param merge_dims a character vector indicating the names of the dimensions to merge +#'@param rename_dim a character string indicating the name of the output dimension. If left at NULL, the first dimension name provided in parameter \code{merge_dims} will be used. +#'@param na.rm a logical indicating if the NA values should be removed or not. +#' +#'@import abind +#'@import s2dverification +#'@examples +#' +#'data <- 1 : 20 +#'dim(data) <- c(time = 10, lat = 2) +#'new_data <- MergeDims(data, merge_dims = c('time', 'lat')) +#'@export +MergeDims <- function(data, merge_dims = c('time', 'monthly'), rename_dim = NULL, + na.rm = FALSE) { + # check data + if (is.null(data)) { + stop("Parameter 'data' cannot be NULL.") + } + if (is.null(dim(data))) { + stop("Parameter 'data' must have dimensions.") + } + if (is.null(names(dim(data)))) { + stop("Parameter 'data' must have dimension names.") + } + dims <- dim(data) + # check merge_dims + if (is.null(merge_dims)) { + stop("Parameter 'merge_dims' cannot be NULL.") + } + if (!is.character(merge_dims)) { + stop("Parameter 'merge_dims' must be a character vector ", + "indicating the names of the dimensions to be merged.") + } + if (length(merge_dims) > 2) { + warning("Only two dimensions can be merge, only the first two ", + "dimension will be used. To merge further dimensions ", + "consider to use this function multiple times.") + merge_dims <- merge_dims[1 : 2] + } else if (length(merge_dims) < 2) { + stop("Parameter 'merge_dims' must be of length two.") + } + if (is.null(rename_dim)) { + rename_dim <- merge_dims[1] + } + if (length(rename_dim) > 1) { + warning("Parameter 'rename_dim' has length greater than 1 ", + "and only the first element will be used.") + rename_dim <- as.character(rename_dim[1]) + } + if (!any(names(dims) %in% merge_dims)) { + stop("Parameter 'merge_dims' must match with dimension ", + "names in parameter 'data'.") + } + pos1 <- which(names(dims) == merge_dims[1]) + pos2 <- which(names(dims) == merge_dims[2]) + if (length(pos1) == 0 | length(pos2) == 0) { + stop("Parameter 'merge_dims' must match with dimension ", + "names in parameter 'data'.") + } + if (pos1 > pos2) { + pos1 <- pos1 - 1 + } + data <- lapply(1 : dims[pos2], function(x) {Subset(data, along = pos2, + indices = x, drop = 'selected')}) + data <- abind(data, along = pos1) + names(dim(data)) <- names(dims)[-pos2] + if (!is.null(rename_dim)) { + names(dim(data))[pos1] <- rename_dim + } + if (na.rm) { + nas <- which(is.na(Subset(data, along = -pos1, indices = 1))) + if (length(nas) != 0) { + nas <- unlist(lapply(nas, function(x) { + if(all(is.na(Subset(data, along = pos1, + indices = x)))) { + return(x)}})) + data <- Subset(data, along = pos1, indices = -nas) + } + } +return(data) +} diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index cf531063b23854f54f3b4107482fd499d32b5e33..c40cc32c3b57211e9d14f2555c67541ec18b523d 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -1,12 +1,12 @@ #'Plot one or multiple ensemble forecast pdfs for the same event #' #'@author Llorenç Lledó \email{llledo@bsc.es} -#'@description This function plots the probability distribution function of several ensemble forecasts for the same event, either initialized at different moments or by different models. Probabilities for tercile categories are computed, plotted in colors and annotated. An asterisk marks the tercile with higher probabilities. Probabilities for extreme categories (above P90 and below P10) can also be included as hatched areas. Individual ensemble members can be plotted as jittered points. The observed value is optionally shown as a diamond. +#'@description This function plots the probability distribution function of several ensemble forecasts. Separate panels are used to plot forecasts valid or initialized at different times or by different models or even at different locations. Probabilities for tercile categories are computed, plotted in colors and annotated. An asterisk marks the tercile with higher probabilities. Probabilities for extreme categories (above P90 and below P10) can also be included as hatched areas. Individual ensemble members can be plotted as jittered points. The observed value is optionally shown as a diamond. #' -#'@param fcst a dataframe or array containing all the ensember members for each frecast. If \code{'fcst'} is an array, it should have two labelled dimensions, and one of them should be \code{'members'}. If \code{'fcsts'} is a data.frame, each column shoul be a separate forecast, with the rows beeing the different ensemble members. -#'@param tercile.limits an array with P33 and P66 values that define the tercile categories. -#'@param extreme.limits (optional) an array with P10 and P90 values that define the extreme categories. (Default: extreme categories are not shown). -#'@param obs (optional) a number with the observed value. (Default: observation is not shown). +#'@param fcst a dataframe or array containing all the ensember members for each forecast. If \code{'fcst'} is an array, it should have two labelled dimensions, and one of them should be \code{'members'}. If \code{'fcsts'} is a data.frame, each column shoul be a separate forecast, with the rows beeing the different ensemble members. +#'@param tercile.limits an array or vector with P33 and P66 values that define the tercile categories for each panel. Use an array of dimensions (nforecasts,2) to define different terciles for each forecast panel, or a vector with two elements to reuse the same tercile limits for all forecast panels. +#'@param extreme.limits (optional) an array or vector with P10 and P90 values that define the extreme categories for each panel. Use an array of (nforecasts,2) to define different extreme limits for each forecast panel, or a vector with two elements to reuse the same tercile limits for all forecast panels. (Default: extreme categories are not shown). +#'@param obs (optional) A vector providing the observed values for each forecast panel or a single value that will be reused for all forecast panels. (Default: observation is not shown). #'@param plotfile (optional) a filename (pdf, png...) where the plot will be saved. (Default: the plot is not saved). #'@param title a string with the plot title. #'@param var.name a string with the variable name and units. @@ -23,6 +23,7 @@ #'@importFrom reshape2 melt #'@importFrom plyr . #'@importFrom plyr dlply +#'@import s2dverification #' #'@examples #'fcsts <- data.frame(fcst1 = rnorm(10), fcst2 = rnorm(10, 0.5, 1.2), @@ -37,7 +38,8 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NULL, plotfile = NULL, title = "Set a title", var.name = "Varname (units)", fcst.names = NULL, add.ensmemb = c("above", "below", "no"), color.set = c("ggplot", "s2s4e", "hydro")) { - value <- init <- extremes <- x <- ymin <- ymax <- tercile <- y <- xend <- yend <- yjitter <- MLT <- NULL + value <- init <- extremes <- x <- ymin <- ymax <- tercile <- NULL + y <- xend <- yend <- yjitter <- MLT <- lab.pos <- NULL ggColorHue <- function(n) { hues <- seq(15, 375, length = n + 1) hcl(h = hues, l = 65, c = 100)[1:n] @@ -71,20 +73,6 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N # Check input arguments #------------------------ add.ensmemb <- match.arg(add.ensmemb) - if (length(tercile.limits) != 2) { - stop("Parameter 'tercile.limits' should be an array with two limits for delimiting tercile categories") - } - if (tercile.limits[1] >= tercile.limits[2]) { - stop("The provided tercile limits are in the wrong order") - } - if (!is.null(extreme.limits)) { - if (length(extreme.limits) != 2) { - stop("Parameter 'extreme.limits' should be an array with two limits for delimiting extreme categories") - } - if (extreme.limits[1] >= tercile.limits[1] | extreme.limits[2] <= tercile.limits[2]) { - stop("The provided extreme limits are not consistent with tercile limits") - } - } #------------------------ # Check fcst type and convert to data.frame if needed #------------------------ @@ -103,6 +91,74 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N } else { stop("Parameter 'fcst' should be an array or a data.frame") } + npanels <- ncol(fcst.df) + #------------------------ + # Check observations + #------------------------ + if (!is.null(obs)) { + if (!is.vector(obs)){ + stop("Observations should be a vector") + } else if (!length(obs) %in% c(1, npanels)) { + stop("The number of observations should equal one or the number of forecasts") + } + } + #------------------------ + # Check tercile limits + #------------------------ + if (is.vector(tercile.limits)) { + if (length(tercile.limits) != 2) { + stop("Provide two tercile limits") + } + tercile.limits <- InsertDim(tercile.limits, 1, npanels) + } else if (is.array(tercile.limits)) { + if (length(dim(tercile.limits)) == 2) { + if (dim(tercile.limits)[2] != 2) { + stop("Provide two tercile limits for each panel") + } + if (dim(tercile.limits)[1] != npanels) { + stop("The number of tercile limits does not match the number of forecasts provided") + } + } else { + stop("Tercile limits should have two dimensions") + } + } else { + stop("Tercile limits should be a vector of length two or an array of dimension (nfcsts,2)") + } + # check consistency of tercile limits + if (any(tercile.limits[, 1] >= tercile.limits[, 2])) { + stop("Inconsistent tercile limits") + } + #------------------------ + # Check extreme limits + #------------------------ + if (!is.null(extreme.limits)) { + if (is.vector(extreme.limits)) { + if (length(extreme.limits) != 2) { + stop("Provide two extreme limits") + } + extreme.limits <- InsertDim(extreme.limits, 1, npanels) + } else if (is.array(extreme.limits)) { + if (length(dim(extreme.limits)) == 2) { + if (dim(extreme.limits)[2] != 2) { + stop("Provide two extreme limits for each panel") + } + if (dim(extreme.limits)[1] != npanels) { + stop("The number of extreme limits does not match the number of forecasts provided") + } + } else { + stop("extreme limits should have two dimensions") + } + } else { + stop("Extreme limits should be a vector of length two or an array of dimensions (nfcsts,2)") + } + # Check that extreme limits are consistent with tercile limits + if (any(extreme.limits[, 1] >= tercile.limits[, 1])) { + stop("Inconsistent lower extreme limits") + } + if (any(extreme.limits[, 2] <= tercile.limits[, 2])) { + stop("Inconsistent higher extreme limits") + } + } #------------------------ # Set proper fcst names #------------------------ @@ -132,8 +188,8 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N } else { stop("Cannot find PANELS in ggp object") } - tmp.df$tercile <- factor(ifelse(tmp.df$x < tercile.limits[1], "Below normal", - ifelse(tmp.df$x < tercile.limits[2], "Normal", "Above normal")), levels = c("Below normal", + tmp.df$tercile <- factor(ifelse(tmp.df$x < tercile.limits[tmp.df$PANEL, 1], "Below normal", + ifelse(tmp.df$x < tercile.limits[tmp.df$PANEL, 2], "Normal", "Above normal")), levels = c("Below normal", "Normal", "Above normal")) #------------------------ # Get the height and width of a panel @@ -145,8 +201,8 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N # Compute hatch coordinates for extremes #------------------------ if (!is.null(extreme.limits)) { - tmp.df$extremes <- factor(ifelse(tmp.df$x < extreme.limits[1], "Below P10", - ifelse(tmp.df$x < extreme.limits[2], "Normal", "Above P90")), levels = c("Below P10", + tmp.df$extremes <- factor(ifelse(tmp.df$x < extreme.limits[tmp.df$PANEL, 1], "Below P10", + ifelse(tmp.df$x < extreme.limits[tmp.df$PANEL, 2], "Normal", "Above P90")), levels = c("Below P10", "Normal", "Above P90")) hatch.ls <- dlply(tmp.df, .(init, extremes), function(x) { # close the polygon @@ -254,7 +310,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N pct <- merge(pct, tot, by = "init") pct$pct <- round(100 * pct$pct/pct$tot, 0) pct$MLT <- pct[, .(MLT = pct == max(pct)), by = init]$MLT - lab.pos <- c(tercile.limits[1], mean(tercile.limits), tercile.limits[2]) + pct$lab.pos <- as.vector(apply(tercile.limits, 1, function(x) {c(min(x), mean(x), max(x))})) #------------------------ # Compute probability for extremes #------------------------ @@ -264,7 +320,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N tot2 <- pct2[, .(tot = sum(pct)), by = init] pct2 <- merge(pct2, tot2, by = "init") pct2$pct <- round(100 * pct2$pct/pct2$tot, 0) - lab.pos2 <- c(extreme.limits[1], NA, extreme.limits[2]) + pct2$lab.pos <- as.vector(apply(extreme.limits, 1, function(x) {c(x[1], NA, x[2])})) pct2 <- merge(pct2, max.df, by = c("init", "extremes")) # include potentially missing groups pct2 <- pct2[CJ(levels(pct2$init), factor(c("Below P10", "Normal", "Above P90"), @@ -280,16 +336,16 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N labpos <- 0 vjust <- -0.5 } - plot <- plot + geom_text(data = pct, aes(x = lab.pos[as.integer(tercile)], y = labpos, + plot <- plot + geom_text(data = pct, aes(x = lab.pos, y = labpos, label = paste0(pct, "%"), hjust = as.integer(tercile) * 1.5 - 2.5), vjust = vjust, - angle = -90, size = 3.2) + geom_text(data = pct[MLT == T, ], aes(x = lab.pos[as.integer(tercile)], + angle = -90, size = 3.2) + geom_text(data = pct[MLT == T, ], aes(x = lab.pos, y = labpos, label = "*", hjust = as.integer(tercile) * 3.5 - 5), vjust = 0.1, angle = -90, size = 7, color = "black") #------------------------ # Add probability labels for extremes #------------------------ if (!is.null(extreme.limits)) { - plot <- plot + geom_text(data = pct2[extremes != "Normal", ], aes(x = lab.pos2[as.integer(extremes)], + plot <- plot + geom_text(data = pct2[extremes != "Normal", ], aes(x = lab.pos, y = 0.9 * y, label = paste0(pct, "%"), hjust = as.integer(extremes) * 1.5 - 2.5), vjust = -0.5, angle = -90, size = 3.2, color = rep(colorLab, dim(fcst.df)[2])) diff --git a/man/Analogs.Rd b/man/Analogs.Rd index ee8a737ece7426357f9b4d7dc675d1de35ea79de..06107c078ff49d5978679daa47625b319bec3926 100644 --- a/man/Analogs.Rd +++ b/man/Analogs.Rd @@ -4,9 +4,19 @@ \alias{Analogs} \title{Analogs based on large scale fields.} \usage{ -Analogs(expL, obsL, time_obsL, expVar = NULL, obsVar = NULL, - criteria = "Large_dist", lonVar = NULL, latVar = NULL, region = NULL, - nAnalogs = NULL, return_list = FALSE) +Analogs( + expL, + obsL, + time_obsL, + expVar = NULL, + obsVar = NULL, + criteria = "Large_dist", + lonVar = NULL, + latVar = NULL, + region = NULL, + nAnalogs = NULL, + return_list = FALSE +) } \arguments{ \item{expL}{an array of N named dimensions containing the experimental field @@ -377,11 +387,6 @@ Local_scalecor <- Analogs(expL=expSLP, str(Local_scalecor) Local_scalecor$AnalogsInfo -} -\author{ -M. Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} - -Nuria Perez-Zanon \email{nuria.perez@bsc.es} } \references{ Yiou, P., T. Salameh, P. Drobinski, L. Menut, R. Vautard, @@ -389,4 +394,8 @@ and M. Vrac, 2013 : Ensemble reconstruction of the atmospheric column from surface pressure using analogues. Clim. Dyn., 41, 1419-1437. \email{pascal.yiou@lsce.ipsl.fr} } +\author{ +M. Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +Nuria Perez-Zanon \email{nuria.perez@bsc.es} +} diff --git a/man/BEI_PDFBest.Rd b/man/BEI_PDFBest.Rd index f836ab722211e359009ca03e7844a00a9b6d2162..0ba24a8487daf4c3f47462b96ac7fbbe829d0700 100644 --- a/man/BEI_PDFBest.Rd +++ b/man/BEI_PDFBest.Rd @@ -4,9 +4,16 @@ \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) +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 @@ -113,12 +120,11 @@ 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 } - +\author{ +Eroteida Sanchez-Garcia - AEMET, \email{esanchezg@aemet.es} +} diff --git a/man/BEI_Weights.Rd b/man/BEI_Weights.Rd index 61db33af3a5edcf6a140f35d7a5f3240320b663b..867a4eb0102a41b88f8ac64ac7707d6fb7b3b37a 100644 --- a/man/BEI_Weights.Rd +++ b/man/BEI_Weights.Rd @@ -43,13 +43,12 @@ 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 } - +\author{ +Eroteida Sanchez-Garcia - AEMET, \email{esanchezg@aemet.es} +} diff --git a/man/CST_Analogs.Rd b/man/CST_Analogs.Rd index 7c9a1e6f125266b510a3f9ccd11c4511d7a927b9..d7dd5e14b8ff65624d6b91807f60f0333768a227 100644 --- a/man/CST_Analogs.Rd +++ b/man/CST_Analogs.Rd @@ -4,8 +4,15 @@ \alias{CST_Analogs} \title{Downscaling using Analogs based on large scale fields.} \usage{ -CST_Analogs(expL, obsL, time_obsL, expVar = NULL, obsVar = NULL, - region = NULL, criteria = "Large_dist") +CST_Analogs( + expL, + obsL, + time_obsL, + expVar = NULL, + obsVar = NULL, + region = NULL, + criteria = "Large_dist" +) } \arguments{ \item{expL}{an 's2dv_cube' object containing the experimental field on the @@ -81,11 +88,6 @@ adapted version of the method of Yiou et al 2013. \examples{ res <- CST_Analogs(expL = lonlat_data$exp, obsL = lonlat_data$obs) -} -\author{ -M. Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} - -Nuria Perez-Zanon \email{nuria.perez@bsc.es} } \references{ Yiou, P., T. Salameh, P. Drobinski, L. Menut, R. Vautard, @@ -97,4 +99,8 @@ from surface pressure using analogues. Clim. Dyn., 41, 1419-1437. code{\link{CST_Load}}, \code{\link[s2dverification]{Load}} and \code{\link[s2dverification]{CDORemap}} } +\author{ +M. Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +Nuria Perez-Zanon \email{nuria.perez@bsc.es} +} diff --git a/man/CST_Anomaly.Rd b/man/CST_Anomaly.Rd index e1c31f0ce5d3f09935f16adfa3c1a1fb08e620ca..07691ea7d7ff5a6105b0e0f4e53a9373d4fdccab 100644 --- a/man/CST_Anomaly.Rd +++ b/man/CST_Anomaly.Rd @@ -4,8 +4,7 @@ \alias{CST_Anomaly} \title{Anomalies relative to a climatology along selected dimension with or without cross-validation} \usage{ -CST_Anomaly(exp = NULL, obs = NULL, cross = FALSE, memb = TRUE, - dim_anom = 3) +CST_Anomaly(exp = NULL, obs = NULL, cross = FALSE, memb = TRUE, dim_anom = 3) } \arguments{ \item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}.} @@ -53,13 +52,12 @@ str(anom3) anom4 <- CST_Anomaly(exp = exp, obs = obs, cross = FALSE, memb = FALSE) str(anom4) +} +\seealso{ +\code{\link[s2dverification]{Ano_CrossValid}}, \code{\link[s2dverification]{Clim}} and \code{\link{CST_Load}} } \author{ Perez-Zanon Nuria, \email{nuria.perez@bsc.es} Pena Jesus, \email{jesus.pena@bsc.es} } -\seealso{ -\code{\link[s2dverification]{Ano_CrossValid}}, \code{\link[s2dverification]{Clim}} and \code{\link{CST_Load}} -} - diff --git a/man/CST_BEI_Weighting.Rd b/man/CST_BEI_Weighting.Rd index 6b9a448ae742037887d2e058e679a08edccf71d0..d6f65bb52123d014377037e435f5c33893756b9f 100644 --- a/man/CST_BEI_Weighting.Rd +++ b/man/CST_BEI_Weighting.Rd @@ -4,8 +4,13 @@ \alias{CST_BEI_Weighting} \title{Weighting SFSs of a CSTools object.} \usage{ -CST_BEI_Weighting(var_exp, aweights, type = "ensembleMean", - time_dim_name = "time") +CST_BEI_Weighting( + var_exp, + aweights, + terciles = NULL, + type = "ensembleMean", + time_dim_name = "time" +) } \arguments{ \item{var_exp}{An object of the class 's2dv_cube' containing the variable @@ -18,6 +23,11 @@ at least a temporal dimension and a dimension named 'member'.} 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{terciles}{A numeric array with at least one dimension 'tercil' equal to +2, the first element is the lower tercil for a hindcast period, and the second +element is the upper tercile. By default is NULL, the terciles are computed +from var_exp data.} + \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 @@ -63,12 +73,11 @@ 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 } - +\author{ +Eroteida Sanchez-Garcia - AEMET, \email{esanchezg@aemet.es} +} diff --git a/man/CST_BiasCorrection.Rd b/man/CST_BiasCorrection.Rd index 485199eae252039266185188df10ee21cf1af52b..a1b415fb525a9f2e3c72171e631c2633e0aefcd6 100644 --- a/man/CST_BiasCorrection.Rd +++ b/man/CST_BiasCorrection.Rd @@ -35,10 +35,9 @@ attr(obs, 'class') <- 's2dv_cube' a <- CST_BiasCorrection(exp = exp, obs = obs) str(a) } -\author{ -Verónica Torralba, \email{veronica.torralba@bsc.es} -} \references{ 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} +\author{ +Verónica Torralba, \email{veronica.torralba@bsc.es} +} diff --git a/man/CST_Calibration.Rd b/man/CST_Calibration.Rd index 36171dbde22451fd681d7dbb899f5f54379a410e..891e2e5fc10fb6751d4e0a6b08cd5b71232a0ed4 100644 --- a/man/CST_Calibration.Rd +++ b/man/CST_Calibration.Rd @@ -4,8 +4,15 @@ \alias{CST_Calibration} \title{Forecast Calibration} \usage{ -CST_Calibration(exp, obs, cal.method = "mse_min", - eval.method = "leave-one-out", multi.model = F) +CST_Calibration( + exp, + obs, + cal.method = "mse_min", + eval.method = "leave-one-out", + multi.model = F, + na.fill = T, + ncores = 1 +) } \arguments{ \item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}.} @@ -17,33 +24,19 @@ CST_Calibration(exp, obs, cal.method = "mse_min", \item{eval.method}{is the sampling method used, can be either \code{in-sample} or \code{leave-one-out}. Default value is the \code{leave-one-out} cross validation.} \item{multi.model}{is a boolean that is used only for the \code{mse_min} method. If multi-model ensembles or ensembles of different sizes are used, it must be set to \code{TRUE}. By default it is \code{FALSE}. Differences between the two approaches are generally small but may become large when using small ensemble sizes. Using multi.model when the calibration method is \code{bias}, \code{evmos} or \code{crps_min} will not affect the result.} + +\item{na.fill}{is a boolean that indicates what happens in case calibration is not possible or will yield unreliable results. This happens when three or less forecasts-observation pairs are available to perform the training phase of the calibration. By default \code{na.fill} is set to true such that NA values will be returned. If \code{na.fill} is set to false, the uncorrected data will be returned.} + +\item{ncores}{is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one.} } \value{ -an object of class \code{s2dv_cube} containing the calibrated forecasts in the element \code{$data} with the same dimensions of the experimental data. +an object of class \code{s2dv_cube} containing the calibrated forecasts in the element \code{$data} with the same dimensions as the one in the exp object. } \description{ -Four types of member-by-member bias correction can be performed. The \code{bias} method corrects the bias only, the \code{evmos} method applies a variance inflation technique to ensure the correction of the bias and the correspondence of variance between forecast and observation (Van Schaeybroeck and Vannitsem, 2011). The ensemble calibration methods \code{"mse_min"} and \code{"crps_min"} correct the bias, the overall forecast variance and the ensemble spread as described in Doblas-Reyes et al. (2005) and Van Schaeybroeck and Vannitsem (2015), respectively. While the \code{"mse_min"} method minimizes a constrained mean-squared error using three parameters, the \code{"crps_min"} method features four parameters and minimizes the Continuous Ranked Probability Score (CRPS). - -Both in-sample or our out-of-sample (leave-one-out cross validation) calibration are possible. -} -\author{ -Verónica Torralba, \email{veronica.torralba@bsc.es} - -Bert Van Schaeybroeck, \email{bertvs@meteo.be} -} -\references{ -Doblas-Reyes F.J, Hagedorn R, Palmer T.N. The rationale behind the success of multi-model ensembles in seasonal forecasting-II calibration and combination. Tellus A. 2005;57:234-252. doi:10.1111/j.1600-0870.2005.00104.x - -Van Schaeybroeck, B., & Vannitsem, S. (2011). Post-processing through linear regression. Nonlinear Processes in Geophysics, 18(2), 147. doi:10.5194/npg-18-147-2011 - -Van Schaeybroeck, B., & Vannitsem, S. (2015). Ensemble post-processing using member-by-member approaches: theoretical aspects. Quarterly Journal of the Royal Meteorological Society, 141(688), 807-818. doi:10.1002/qj.2397 +Equivalent to function \code{Calibration} but for objects of class \code{s2dv_cube}. } -\seealso{ -\code{\link{CST_Load}} - -# Example -# Creation of sample s2dverification objects. These are not complete -# s2dverification objects though. The Load function returns complete objects. +\examples{ +# Example 1: mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7) dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7) obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7) @@ -57,4 +50,11 @@ attr(obs, 'class') <- 's2dv_cube' a <- CST_Calibration(exp = exp, obs = obs, cal.method = "mse_min", eval.method = "in-sample") str(a) } -\encoding{UTF-8} +\seealso{ +\code{\link{CST_Load}} +} +\author{ +Verónica Torralba, \email{veronica.torralba@bsc.es} + +Bert Van Schaeybroeck, \email{bertvs@meteo.be} +} diff --git a/man/CST_CategoricalEnsCombination.Rd b/man/CST_CategoricalEnsCombination.Rd index e551c3ecaf7e9b67cbe212fcffdf7bae905d08fe..c23f8341c201921984cdebb6bf05a43997042752 100644 --- a/man/CST_CategoricalEnsCombination.Rd +++ b/man/CST_CategoricalEnsCombination.Rd @@ -4,8 +4,14 @@ \alias{CST_CategoricalEnsCombination} \title{Make categorical forecast based on a multi-model forecast with potential for calibrate} \usage{ -CST_CategoricalEnsCombination(exp, obs, cat.method = "pool", - eval.method = "leave-one-out", amt.cat = 3, ...) +CST_CategoricalEnsCombination( + exp, + obs, + cat.method = "pool", + eval.method = "leave-one-out", + amt.cat = 3, + ... +) } \arguments{ \item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}. The amount of forecasting models is equal to the size of the \code{dataset} dimension of the data array. The amount of members per model may be different. The size of the \code{member} dimension of the data array is equal to the maximum of the ensemble members among the models. Models with smaller ensemble sizes have residual indices of \code{member} dimension in the data array filled with NA values.} @@ -83,9 +89,6 @@ attr(obs, 'class') <- 's2dv_cube' a <- CST_CategoricalEnsCombination(exp = exp, obs = obs, amt.cat = 3, cat.method = "mmw") } } -\author{ -Bert Van Schaeybroeck, \email{bertvs@meteo.be} -} \references{ Rajagopalan, B., Lall, U., & Zebiak, S. E. (2002). Categorical climate forecasts through regularization and optimal combination of multiple GCM ensembles. Monthly Weather Review, 130(7), 1792-1811. @@ -93,4 +96,6 @@ Robertson, A. W., Lall, U., Zebiak, S. E., & Goddard, L. (2004). Improved combin Van Schaeybroeck, B., & Vannitsem, S. (2019). Postprocessing of Long-Range Forecasts. In Statistical Postprocessing of Ensemble Forecasts (pp. 267-290). } - +\author{ +Bert Van Schaeybroeck, \email{bertvs@meteo.be} +} diff --git a/man/CST_EnsClustering.Rd b/man/CST_EnsClustering.Rd index c13bf20506c21b46cf6ad7daae2e4e06efdd1569..154541d55a3adb839e923ef409fe8cb82cdd59ff 100644 --- a/man/CST_EnsClustering.Rd +++ b/man/CST_EnsClustering.Rd @@ -4,10 +4,18 @@ \alias{CST_EnsClustering} \title{Ensemble clustering} \usage{ -CST_EnsClustering(exp, time_moment = "mean", numclus = NULL, - lon_lim = NULL, lat_lim = NULL, variance_explained = 80, - numpcs = NULL, time_percentile = 90, cluster_dim = "member", - verbose = F) +CST_EnsClustering( + exp, + time_moment = "mean", + numclus = NULL, + lon_lim = NULL, + lat_lim = NULL, + variance_explained = 80, + numpcs = NULL, + time_percentile = 90, + cluster_dim = "member", + verbose = F +) } \arguments{ \item{exp}{An object of the class 's2dv_cube', containing the variables to be analysed. @@ -125,4 +133,3 @@ Paolo Davini - ISAC-CNR, \email{p.davini@isac.cnr.it} Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} } - diff --git a/man/CST_Load.Rd b/man/CST_Load.Rd index 1fee022c23b928c072798a4136d927cbf6c2a287..bf03ba42d98810643a98bfca27f4563f6361162e 100644 --- a/man/CST_Load.Rd +++ b/man/CST_Load.Rd @@ -47,4 +47,3 @@ obs <- CSTools::lonlat_data$obs \author{ Nicolau Manubens, \email{nicolau.manubens@bsc.es} } - diff --git a/man/CST_MergeDims.Rd b/man/CST_MergeDims.Rd new file mode 100644 index 0000000000000000000000000000000000000000..0762e83f9567f082040a8692ebf18c2c90ca45fd --- /dev/null +++ b/man/CST_MergeDims.Rd @@ -0,0 +1,44 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_MergeDims.R +\name{CST_MergeDims} +\alias{CST_MergeDims} +\title{Function to Merge Dimensions} +\usage{ +CST_MergeDims( + data, + merge_dims = c("ftime", "monthly"), + rename_dim = NULL, + na.rm = FALSE +) +} +\arguments{ +\item{data}{a 's2dv_cube' object} + +\item{merge_dims}{a character vector indicating the names of the dimensions to merge} + +\item{rename_dim}{a character string indicating the name of the output dimension. If left at NULL, the first dimension name provided in parameter \code{merge_dims} will be used.} + +\item{na.rm}{a logical indicating if the NA values should be removed or not.} +} +\description{ +This function merges two dimensions of the array \code{data} in a 's2dv_cube' object into one. The user can select the dimensions to merge and provide the final name of the dimension. The user can select to remove NA values or keep them. +} +\examples{ + +data <- 1 : c(2 * 3 * 4 * 5 * 6 * 7) +dim(data) <- c(time = 7, lat = 2, lon = 3, monthly = 4, member = 6, + dataset = 5, var = 1) +data[2,,,,,,] <- NA +data[c(3,27)] <- NA +data <-list(data = data) +class(data) <- 's2dv_cube' +new_data <- CST_MergeDims(data, merge_dims = c('time', 'monthly')) +dim(new_data$data) +new_data <- CST_MergeDims(data, merge_dims = c('lon', 'lat'), rename_dim = 'grid') +dim(new_data$data) +new_data <- CST_MergeDims(data, merge_dims = c('time', 'monthly'), na.rm = TRUE) +dim(new_data$data) +} +\author{ +Nuria Perez-Zanon, \email{nuria.perez@bsc.es} +} diff --git a/man/CST_MultiEOF.Rd b/man/CST_MultiEOF.Rd index fb5847515de2b370481c2b31ec8eb90abef17f92..036a6470665320fea218286579ff33c7e80df387 100644 --- a/man/CST_MultiEOF.Rd +++ b/man/CST_MultiEOF.Rd @@ -4,8 +4,14 @@ \alias{CST_MultiEOF} \title{EOF analysis of multiple variables} \usage{ -CST_MultiEOF(datalist, neof_max = 40, neof_composed = 5, minvar = 0.6, - lon_lim = NULL, lat_lim = NULL) +CST_MultiEOF( + datalist, + neof_max = 40, + neof_composed = 5, + minvar = 0.6, + lon_lim = NULL, + lat_lim = NULL +) } \arguments{ \item{datalist}{A list of objects of the class 's2dv_cube', containing the variables to be analysed. @@ -69,4 +75,3 @@ Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} Paolo Davini - ISAC-CNR, \email{p.davini@isac.cnr.it} } - diff --git a/man/CST_MultiMetric.Rd b/man/CST_MultiMetric.Rd index 079a5588fa13982a8e9d39270b1d71b7cc882040..8e3ce593b7024138e191c4fc818516358503be7d 100644 --- a/man/CST_MultiMetric.Rd +++ b/man/CST_MultiMetric.Rd @@ -37,15 +37,14 @@ c(ano_exp, ano_obs) \%<-\% CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb a <- CST_MultiMetric(exp = ano_exp, obs = ano_obs) str(a) } -\author{ -Mishra Niti, \email{niti.mishra@bsc.es} - -Perez-Zanon Nuria, \email{nuria.perez@bsc.es} -} \references{ Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31.\url{http://link.springer.com/10.1007/s00382-018-4404-z} } \seealso{ \code{\link[s2dverification]{Corr}}, \code{\link[s2dverification]{RMS}}, \code{\link[s2dverification]{RMSSS}} and \code{\link{CST_Load}} } +\author{ +Mishra Niti, \email{niti.mishra@bsc.es} +Perez-Zanon Nuria, \email{nuria.perez@bsc.es} +} diff --git a/man/CST_MultivarRMSE.Rd b/man/CST_MultivarRMSE.Rd index 685eaf7712ccabbcef7d6cdea64a2f6f23556a31..24af608c8360985de73763434a16f1a11fd35ffe 100644 --- a/man/CST_MultivarRMSE.Rd +++ b/man/CST_MultivarRMSE.Rd @@ -56,10 +56,9 @@ weight <- c(1, 2) a <- CST_MultivarRMSE(exp = ano_exp, obs = ano_obs, weight = weight) str(a) } -\author{ -Deborah Verfaillie, \email{deborah.verfaillie@bsc.es} -} \seealso{ \code{\link[s2dverification]{RMS}} and \code{\link{CST_Load}} } - +\author{ +Deborah Verfaillie, \email{deborah.verfaillie@bsc.es} +} diff --git a/man/CST_QuantileMapping.Rd b/man/CST_QuantileMapping.Rd index 1c93843e966ab8e3dce0c46b0f71d44d57fe5481..ad8f4b6ccaf59532d503ea314e29a337fca1b2f4 100644 --- a/man/CST_QuantileMapping.Rd +++ b/man/CST_QuantileMapping.Rd @@ -4,9 +4,16 @@ \alias{CST_QuantileMapping} \title{Quantiles Mapping for seasonal or decadal forecast data} \usage{ -CST_QuantileMapping(exp, obs, exp_cor = NULL, sample_dims = c("sdate", - "ftime", "member"), sample_length = NULL, method = "QUANT", - ncores = NULL, ...) +CST_QuantileMapping( + exp, + obs, + exp_cor = NULL, + sample_dims = c("sdate", "ftime", "member"), + sample_length = NULL, + method = "QUANT", + ncores = NULL, + ... +) } \arguments{ \item{exp}{an object of class \code{s2dv_cube}} @@ -77,10 +84,9 @@ res <- CST_QuantileMapping(exp = exp, obs = obs, sample_dims = 'time', method = 'DIST') } } -\author{ -Nuria Perez-Zanon, \email{nuria.perez@bsc.es} -} \seealso{ \code{\link[qmap]{fitQmap}} and \code{\link[qmap]{doQmap}} } - +\author{ +Nuria Perez-Zanon, \email{nuria.perez@bsc.es} +} diff --git a/man/CST_RFSlope.Rd b/man/CST_RFSlope.Rd index d2b5aec01ec78eb05bb8d3fbf41bf06ebf20dfd4..0c4e16710c82a98b1d1245a2321f4675cd5a578e 100644 --- a/man/CST_RFSlope.Rd +++ b/man/CST_RFSlope.Rd @@ -50,4 +50,3 @@ slopes \author{ Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} } - diff --git a/man/CST_RFWeights.Rd b/man/CST_RFWeights.Rd index 08a7b850b64dfdcaf28313d9da09883e9ad47d54..ef5ebe4d5dd8585143c3fb78d2854e0783dbdd1d 100644 --- a/man/CST_RFWeights.Rd +++ b/man/CST_RFWeights.Rd @@ -47,9 +47,6 @@ nf <- 8 ww <- CST_RFWeights("./worldclim.nc", nf, lon, lat, fsmooth = TRUE) } } -\author{ -Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} -} \references{ Terzago, S., Palazzi, E., & von Hardenberg, J. (2018). Stochastic downscaling of precipitation in complex orography: @@ -57,4 +54,6 @@ A simple method to reproduce a realistic fine-scale climatology. Natural Hazards and Earth System Sciences, 18(11), 2825-2840. http://doi.org/10.5194/nhess-18-2825-2018 . } - +\author{ +Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +} diff --git a/man/CST_RainFARM.Rd b/man/CST_RainFARM.Rd index 4a667f9ad3c4f28b68c821a8497b895df0d41f7c..1c609e084ca509c1343802604dddf9899088da3e 100644 --- a/man/CST_RainFARM.Rd +++ b/man/CST_RainFARM.Rd @@ -4,9 +4,20 @@ \alias{CST_RainFARM} \title{RainFARM stochastic precipitation downscaling of a CSTools object} \usage{ -CST_RainFARM(data, nf, weights = 1, slope = 0, kmin = 1, nens = 1, - fglob = FALSE, fsmooth = TRUE, nprocs = 1, time_dim = NULL, - verbose = FALSE, drop_realization_dim = FALSE) +CST_RainFARM( + data, + nf, + weights = 1, + slope = 0, + kmin = 1, + nens = 1, + fglob = FALSE, + fsmooth = TRUE, + nprocs = 1, + time_dim = NULL, + verbose = FALSE, + drop_realization_dim = FALSE +) } \arguments{ \item{data}{An object of the class 's2dv_cube' as returned by `CST_Load`, @@ -95,13 +106,12 @@ dim(res$data) # dataset member realization sdate ftime lat lon # 1 2 3 3 4 64 64 -} -\author{ -Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} } \references{ Terzago, S. et al. (2018). NHESS 18(11), 2825-2840. http://doi.org/10.5194/nhess-18-2825-2018 ; D'Onofrio et al. (2014), J of Hydrometeorology 15, 830-843; Rebora et. al. (2006), JHM 7, 724. } - +\author{ +Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} +} diff --git a/man/CST_SaveExp.Rd b/man/CST_SaveExp.Rd index 17537205f319893de9874fa55b76e6a59ca33829..0e49c11955488a64772269976f01ed95919f4af3 100644 --- a/man/CST_SaveExp.Rd +++ b/man/CST_SaveExp.Rd @@ -29,11 +29,10 @@ destination <- "./path/" CST_SaveExp(data = data, destination = destination) } -} -\author{ -Perez-Zanon Nuria, \email{nuria.perez@bsc.es} } \seealso{ \code{\link{CST_Load}}, \code{\link{as.s2dv_cube}} and \code{\link{s2dv_cube}} } - +\author{ +Perez-Zanon Nuria, \email{nuria.perez@bsc.es} +} diff --git a/man/CST_SplitDim.Rd b/man/CST_SplitDim.Rd index 2019ea7bb537173efec25d1dc7cb2e9e3fa5d952..ee93aedc2075b38292d392a655a57b83955aecbf 100644 --- a/man/CST_SplitDim.Rd +++ b/man/CST_SplitDim.Rd @@ -43,4 +43,3 @@ dim(new_data$data) \author{ Nuria Perez-Zanon, \email{nuria.perez@bsc.es} } - diff --git a/man/Calibration.Rd b/man/Calibration.Rd new file mode 100644 index 0000000000000000000000000000000000000000..9f884671ebb4bdf6b6412a5f56c2f7a28df3b603 --- /dev/null +++ b/man/Calibration.Rd @@ -0,0 +1,54 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_Calibration.R +\name{Calibration} +\alias{Calibration} +\title{Forecast Calibration} +\usage{ +Calibration( + exp, + obs, + cal.method = "mse_min", + eval.method = "leave-one-out", + multi.model = F, + na.fill = T, + ncores = 1 +) +} +\arguments{ +\item{exp}{an array containing the seasonal forecast experiment data.} + +\item{obs}{an array containing the observed data.} + +\item{cal.method}{is the calibration method used, can be either \code{bias}, \code{evmos}, \code{mse_min} or \code{crps_min}. Default value is \code{mse_min}.} + +\item{eval.method}{is the sampling method used, can be either \code{in-sample} or \code{leave-one-out}. Default value is the \code{leave-one-out} cross validation.} + +\item{multi.model}{is a boolean that is used only for the \code{mse_min} method. If multi-model ensembles or ensembles of different sizes are used, it must be set to \code{TRUE}. By default it is \code{FALSE}. Differences between the two approaches are generally small but may become large when using small ensemble sizes. Using multi.model when the calibration method is \code{bias}, \code{evmos} or \code{crps_min} will not affect the result.} + +\item{na.fill}{is a boolean that indicates what happens in case calibration is not possible or will yield unreliable results. This happens when three or less forecasts-observation pairs are available to perform the training phase of the calibration. By default \code{na.fill} is set to true such that NA values will be returned. If \code{na.fill} is set to false, the uncorrected data will be returned.} + +\item{ncores}{is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one.} +} +\value{ +an array containing the calibrated forecasts with the same dimensions as the \code{exp} array. +} +\description{ +Four types of member-by-member bias correction can be performed. The \code{bias} method corrects the bias only, the \code{evmos} method applies a variance inflation technique to ensure the correction of the bias and the correspondence of variance between forecast and observation (Van Schaeybroeck and Vannitsem, 2011). The ensemble calibration methods \code{"mse_min"} and \code{"crps_min"} correct the bias, the overall forecast variance and the ensemble spread as described in Doblas-Reyes et al. (2005) and Van Schaeybroeck and Vannitsem (2015), respectively. While the \code{"mse_min"} method minimizes a constrained mean-squared error using three parameters, the \code{"crps_min"} method features four parameters and minimizes the Continuous Ranked Probability Score (CRPS). + +Both in-sample or our out-of-sample (leave-one-out cross validation) calibration are possible. +} +\references{ +Doblas-Reyes F.J, Hagedorn R, Palmer T.N. The rationale behind the success of multi-model ensembles in seasonal forecasting-II calibration and combination. Tellus A. 2005;57:234-252. doi:10.1111/j.1600-0870.2005.00104.x + +Van Schaeybroeck, B., & Vannitsem, S. (2011). Post-processing through linear regression. Nonlinear Processes in Geophysics, 18(2), 147. doi:10.5194/npg-18-147-2011 + +Van Schaeybroeck, B., & Vannitsem, S. (2015). Ensemble post-processing using member-by-member approaches: theoretical aspects. Quarterly Journal of the Royal Meteorological Society, 141(688), 807-818. doi:10.1002/qj.2397 +} +\seealso{ +\code{\link{CST_Load}} +} +\author{ +Verónica Torralba, \email{veronica.torralba@bsc.es} + +Bert Van Schaeybroeck, \email{bertvs@meteo.be} +} diff --git a/man/EnsClustering.Rd b/man/EnsClustering.Rd index 27aca45312679c56cb32533c1fde0fae4064fc24..2fd8a3f1036e46a2e7b6c0eedc8372a08205cbd3 100644 --- a/man/EnsClustering.Rd +++ b/man/EnsClustering.Rd @@ -4,10 +4,20 @@ \alias{EnsClustering} \title{Ensemble clustering} \usage{ -EnsClustering(data, lat, lon, time_moment = "mean", numclus = NULL, - lon_lim = NULL, lat_lim = NULL, variance_explained = 80, - numpcs = NULL, time_percentile = 90, cluster_dim = "member", - verbose = T) +EnsClustering( + data, + lat, + lon, + time_moment = "mean", + numclus = NULL, + lon_lim = NULL, + lat_lim = NULL, + variance_explained = 80, + numpcs = NULL, + time_percentile = 90, + cluster_dim = "member", + verbose = T +) } \arguments{ \item{data}{A matrix of dimensions 'dataset member sdate ftime lat lon' containing the variables to be analysed.} @@ -67,4 +77,3 @@ Paolo Davini - ISAC-CNR, \email{p.davini@isac.cnr.it} Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} } - diff --git a/man/MergeDims.Rd b/man/MergeDims.Rd new file mode 100644 index 0000000000000000000000000000000000000000..7539ef6ef08a128c11591c5a98f4654dfb7f1f4e --- /dev/null +++ b/man/MergeDims.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_MergeDims.R +\name{MergeDims} +\alias{MergeDims} +\title{Function to Split Dimension} +\usage{ +MergeDims( + data, + merge_dims = c("time", "monthly"), + rename_dim = NULL, + na.rm = FALSE +) +} +\arguments{ +\item{data}{an n-dimensional array with named dimensions} + +\item{merge_dims}{a character vector indicating the names of the dimensions to merge} + +\item{rename_dim}{a character string indicating the name of the output dimension. If left at NULL, the first dimension name provided in parameter \code{merge_dims} will be used.} + +\item{na.rm}{a logical indicating if the NA values should be removed or not.} +} +\description{ +This function merges two dimensions of an array into one. The user can select the dimensions to merge and provide the final name of the dimension. The user can select to remove NA values or keep them. +} +\examples{ + +data <- 1 : 20 +dim(data) <- c(time = 10, lat = 2) +new_data <- MergeDims(data, merge_dims = c('time', 'lat')) +} +\author{ +Nuria Perez-Zanon, \email{nuria.perez@bsc.es} +} diff --git a/man/MultiEOF.Rd b/man/MultiEOF.Rd index 1e822fc40c27ab72fd67339f9ec1d7b9fa9751ba..dd0fc7fe59e6a1c6606c2de2cbbceca6f0b4b15d 100644 --- a/man/MultiEOF.Rd +++ b/man/MultiEOF.Rd @@ -4,9 +4,19 @@ \alias{MultiEOF} \title{EOF analysis of multiple variables starting from an array (reduced version)} \usage{ -MultiEOF(data, lon, lat, time, lon_dim = "lon", lat_dim = "lat", - neof_max = 40, neof_composed = 5, minvar = 0.6, lon_lim = NULL, - lat_lim = NULL) +MultiEOF( + data, + lon, + lat, + time, + lon_dim = "lon", + lat_dim = "lat", + neof_max = 40, + neof_composed = 5, + minvar = 0.6, + lon_lim = NULL, + lat_lim = NULL +) } \arguments{ \item{data}{A multidimensional array with dimension \code{"var"}, @@ -46,4 +56,3 @@ Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} Paolo Davini - ISAC-CNR, \email{p.davini@isac.cnr.it} } - diff --git a/man/PlotCombinedMap.Rd b/man/PlotCombinedMap.Rd index 6857c64dcd95f0d710d1e560acc5318abe302f38..616b84f98b6951040beb33961f285f571c4523ad 100644 --- a/man/PlotCombinedMap.Rd +++ b/man/PlotCombinedMap.Rd @@ -4,11 +4,27 @@ \alias{PlotCombinedMap} \title{Plot Multiple Lon-Lat Variables In a Single Map According to a Decision Function} \usage{ -PlotCombinedMap(maps, lon, lat, map_select_fun, display_range, - map_dim = "map", brks = NULL, cols = NULL, col_unknown_map = "white", - mask = NULL, col_mask = "grey", bar_titles = NULL, legend_scale = 1, - fileout = NULL, width = 8, height = 5, size_units = "in", res = 100, - ...) +PlotCombinedMap( + maps, + lon, + lat, + map_select_fun, + display_range, + map_dim = "map", + brks = NULL, + cols = NULL, + col_unknown_map = "white", + mask = NULL, + col_mask = "grey", + bar_titles = NULL, + legend_scale = 1, + fileout = NULL, + width = 8, + height = 5, + size_units = "in", + res = 100, + ... +) } \arguments{ \item{maps}{List of matrices to plot, each with (longitude, latitude) dimensions, or 3-dimensional array with the dimensions (longitude, latitude, map). Dimension names are required.} @@ -67,12 +83,11 @@ PlotCombinedMap(list(a, b, c), lons, lats, bar_titles = paste('\% of belonging to', c('a', 'b', 'c')), brks = 20, width = 10, height = 8) } +\seealso{ +\code{PlotCombinedMap} and \code{PlotEquiMap} +} \author{ Nicolau Manubens, \email{nicolau.manubens@bsc.es} Veronica Torralba, \email{veronica.torralba@bsc.es} } -\seealso{ -\code{PlotCombinedMap} and \code{PlotEquiMap} -} - diff --git a/man/PlotForecastPDF.Rd b/man/PlotForecastPDF.Rd index bed0bd3137789bc45216ab058c8d07d9b019f167..c04b43c19f053c9280af505742a02bb8cc1c3aa6 100644 --- a/man/PlotForecastPDF.Rd +++ b/man/PlotForecastPDF.Rd @@ -4,19 +4,27 @@ \alias{PlotForecastPDF} \title{Plot one or multiple ensemble forecast pdfs for the same event} \usage{ -PlotForecastPDF(fcst, tercile.limits, extreme.limits = NULL, obs = NULL, - plotfile = NULL, title = "Set a title", var.name = "Varname (units)", - fcst.names = NULL, add.ensmemb = c("above", "below", "no"), - color.set = c("ggplot", "s2s4e", "hydro")) +PlotForecastPDF( + fcst, + tercile.limits, + extreme.limits = NULL, + obs = NULL, + plotfile = NULL, + title = "Set a title", + var.name = "Varname (units)", + fcst.names = NULL, + add.ensmemb = c("above", "below", "no"), + color.set = c("ggplot", "s2s4e", "hydro") +) } \arguments{ -\item{fcst}{a dataframe or array containing all the ensember members for each frecast. If \code{'fcst'} is an array, it should have two labelled dimensions, and one of them should be \code{'members'}. If \code{'fcsts'} is a data.frame, each column shoul be a separate forecast, with the rows beeing the different ensemble members.} +\item{fcst}{a dataframe or array containing all the ensember members for each forecast. If \code{'fcst'} is an array, it should have two labelled dimensions, and one of them should be \code{'members'}. If \code{'fcsts'} is a data.frame, each column shoul be a separate forecast, with the rows beeing the different ensemble members.} -\item{tercile.limits}{an array with P33 and P66 values that define the tercile categories.} +\item{tercile.limits}{an array or vector with P33 and P66 values that define the tercile categories for each panel. Use an array of dimensions (nforecasts,2) to define different terciles for each forecast panel, or a vector with two elements to reuse the same tercile limits for all forecast panels.} -\item{extreme.limits}{(optional) an array with P10 and P90 values that define the extreme categories. (Default: extreme categories are not shown).} +\item{extreme.limits}{(optional) an array or vector with P10 and P90 values that define the extreme categories for each panel. Use an array of (nforecasts,2) to define different extreme limits for each forecast panel, or a vector with two elements to reuse the same tercile limits for all forecast panels. (Default: extreme categories are not shown).} -\item{obs}{(optional) a number with the observed value. (Default: observation is not shown).} +\item{obs}{(optional) A vector providing the observed values for each forecast panel or a single value that will be reused for all forecast panels. (Default: observation is not shown).} \item{plotfile}{(optional) a filename (pdf, png...) where the plot will be saved. (Default: the plot is not saved).} @@ -34,7 +42,7 @@ PlotForecastPDF(fcst, tercile.limits, extreme.limits = NULL, obs = NULL, a ggplot object containing the plot. } \description{ -This function plots the probability distribution function of several ensemble forecasts for the same event, either initialized at different moments or by different models. Probabilities for tercile categories are computed, plotted in colors and annotated. An asterisk marks the tercile with higher probabilities. Probabilities for extreme categories (above P90 and below P10) can also be included as hatched areas. Individual ensemble members can be plotted as jittered points. The observed value is optionally shown as a diamond. +This function plots the probability distribution function of several ensemble forecasts. Separate panels are used to plot forecasts valid or initialized at different times or by different models or even at different locations. Probabilities for tercile categories are computed, plotted in colors and annotated. An asterisk marks the tercile with higher probabilities. Probabilities for extreme categories (above P90 and below P10) can also be included as hatched areas. Individual ensemble members can be plotted as jittered points. The observed value is optionally shown as a diamond. } \examples{ fcsts <- data.frame(fcst1 = rnorm(10), fcst2 = rnorm(10, 0.5, 1.2), @@ -49,4 +57,3 @@ PlotForecastPDF(fcsts2, c(-0.66, 0.66), extreme.limits = c(-1.2, 1.2), \author{ Llorenç Lledó \email{llledo@bsc.es} } -\encoding{UTF-8} diff --git a/man/PlotMostLikelyQuantileMap.Rd b/man/PlotMostLikelyQuantileMap.Rd index 6c92850ead6ac72be3985fc233a273a7969e4060..4c400b18bc79f2eaf8d3ff7fd50ef7c72d19ebb8 100644 --- a/man/PlotMostLikelyQuantileMap.Rd +++ b/man/PlotMostLikelyQuantileMap.Rd @@ -4,8 +4,15 @@ \alias{PlotMostLikelyQuantileMap} \title{Plot Maps of Most Likely Quantiles} \usage{ -PlotMostLikelyQuantileMap(probs, lon, lat, cat_dim = "bin", - bar_titles = NULL, col_unknown_cat = "white", ...) +PlotMostLikelyQuantileMap( + probs, + lon, + lat, + cat_dim = "bin", + bar_titles = NULL, + col_unknown_cat = "white", + ... +) } \arguments{ \item{probs}{a list of bi-dimensional arrays with the named dimensions 'latitude' (or 'lat') and 'longitude' (or 'lon'), with equal size and in the same order, or a single tri-dimensional array with an additional dimension (e.g. 'bin') for the different categories. The arrays must contain probability values between 0 and 1, and the probabilities for all categories of a grid cell should not exceed 1 when added.} @@ -109,11 +116,10 @@ PlotMostLikelyQuantileMap(bins, lons, lats, mask = 1 - (w1 + w2 / max(c(w1, w2))), brks = 20, width = 10, height = 8) -} -\author{ -Veronica Torralba, \email{veronica.torralba@bsc.es}, Nicolau Manubens, \email{nicolau.manubens@bsc.es} } \seealso{ \code{PlotCombinedMap} and \code{PlotEquiMap} } - +\author{ +Veronica Torralba, \email{veronica.torralba@bsc.es}, Nicolau Manubens, \email{nicolau.manubens@bsc.es} +} diff --git a/man/RFSlope.Rd b/man/RFSlope.Rd index 09a24ff546568ead10fb64f7bb386d50aaedc25b..db3f0e1054ff5a4b93289a453c5d1bd272a4f9e8 100644 --- a/man/RFSlope.Rd +++ b/man/RFSlope.Rd @@ -4,8 +4,7 @@ \alias{RFSlope} \title{RainFARM spectral slopes from an array (reduced version)} \usage{ -RFSlope(data, kmin = 1, time_dim = NULL, lon_dim = "lon", - lat_dim = "lat") +RFSlope(data, kmin = 1, time_dim = NULL, lon_dim = "lon", lat_dim = "lat") } \arguments{ \item{data}{Array containing the spatial precipitation fields to downscale. @@ -60,4 +59,3 @@ slopes \author{ Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} } - diff --git a/man/RainFARM.Rd b/man/RainFARM.Rd index 984dcd42a6b7f374e3067ffe878ed19eff5fd181..0db8467926923fbed0f9da2f9625e67d5fe06a7e 100644 --- a/man/RainFARM.Rd +++ b/man/RainFARM.Rd @@ -4,10 +4,24 @@ \alias{RainFARM} \title{RainFARM stochastic precipitation downscaling (reduced version)} \usage{ -RainFARM(data, lon, lat, nf, weights = 1, nens = 1, slope = 0, kmin = 1, - fglob = FALSE, fsmooth = TRUE, nprocs = 1, time_dim = NULL, - lon_dim = "lon", lat_dim = "lat", drop_realization_dim = FALSE, - verbose = FALSE) +RainFARM( + data, + lon, + lat, + nf, + weights = 1, + nens = 1, + slope = 0, + kmin = 1, + fglob = FALSE, + fsmooth = TRUE, + nprocs = 1, + time_dim = NULL, + lon_dim = "lon", + lat_dim = "lat", + drop_realization_dim = FALSE, + verbose = FALSE +) } \arguments{ \item{data}{Precipitation array to downscale. @@ -117,4 +131,3 @@ dim(res$data) \author{ Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it} } - diff --git a/man/SplitDim.Rd b/man/SplitDim.Rd index e36aa8a5d12ceb1cb360150be81d01f110d57e02..f07e4756bd2fdbdeaafd4d7cfae30010ec29009a 100644 --- a/man/SplitDim.Rd +++ b/man/SplitDim.Rd @@ -35,4 +35,3 @@ new_data <- SplitDim(data, indices = time, freq = 'year') \author{ Nuria Perez-Zanon, \email{nuria.perez@bsc.es} } - diff --git a/man/areave_data.Rd b/man/areave_data.Rd index cc79c85c1647d0b906d7554bef50ca0d78524ebb..a772220a9e3a7a19a02186938bd9b4f2e72fbc5c 100644 --- a/man/areave_data.Rd +++ b/man/areave_data.Rd @@ -41,4 +41,3 @@ areave_data <- Nicolau Manubens \email{nicolau.manubens@bsc.es} } \keyword{data} - diff --git a/man/as.s2dv_cube.Rd b/man/as.s2dv_cube.Rd index 13a2a296d97fd9cf9b9786b9d0190627688aba86..c2b8f3a8b5b164a015008326cb1cf9b3012a1fc9 100644 --- a/man/as.s2dv_cube.Rd +++ b/man/as.s2dv_cube.Rd @@ -40,12 +40,11 @@ data <- as.s2dv_cube(data) class(data) } } +\seealso{ +\code{\link{s2dv_cube}}, \code{\link[s2dverification]{Load}}, \code{\link[startR]{Start}} and \code{\link{CST_Load}} +} \author{ Perez-Zanon Nuria, \email{nuria.perez@bsc.es} Nicolau Manubens, \email{nicolau.manubens@bsc.es} } -\seealso{ -\code{\link{s2dv_cube}}, \code{\link[s2dverification]{Load}}, \code{\link[startR]{Start}} and \code{\link{CST_Load}} -} - diff --git a/man/lonlat_data.Rd b/man/lonlat_data.Rd index eca7abacfc5f0934387540ec2b24c50a70e9a576..0c6ee30fbfa8f8374eb1509fc80b4e08d33e0c77 100644 --- a/man/lonlat_data.Rd +++ b/man/lonlat_data.Rd @@ -41,4 +41,3 @@ lonlat_data <- Nicolau Manubens \email{nicolau.manubens@bsc.es} } \keyword{data} - diff --git a/man/lonlat_prec.Rd b/man/lonlat_prec.Rd index 69cb94e81032916b91ac6e171f1284558582e396..345e3cabfe965dd9813f92a38393d67b06090003 100644 --- a/man/lonlat_prec.Rd +++ b/man/lonlat_prec.Rd @@ -29,4 +29,3 @@ lonlat_prec <- CST_Load('prlr', exp = list(infile), obs = NULL, Jost von Hardenberg \email{j.vonhardenberg@isac.cnr.it} } \keyword{data} - diff --git a/man/s2dv_cube.Rd b/man/s2dv_cube.Rd index 48af7bbb7a04235c90a947dbba2a72919033d392..b0ce896613c52a27a1e5fb0039526dbd757b6b4c 100644 --- a/man/s2dv_cube.Rd +++ b/man/s2dv_cube.Rd @@ -4,8 +4,16 @@ \alias{s2dv_cube} \title{Creation of a 's2dv_cube' object} \usage{ -s2dv_cube(data, lon = NULL, lat = NULL, Variable = NULL, - Datasets = NULL, Dates = NULL, when = NULL, source_files = NULL) +s2dv_cube( + data, + lon = NULL, + lat = NULL, + Variable = NULL, + Datasets = NULL, + Dates = NULL, + when = NULL, + source_files = NULL +) } \arguments{ \item{data}{an array with any number of named dimensions, typically an object output from CST_Load, with the following dimensions: dataset, member, sdate, ftime, lat and lon.} @@ -75,10 +83,9 @@ exp8 <- s2dv_cube(data = exp_original, lon = seq(-10, 10, 5), lat = c(45, 50), end = paste0(rep("31", 10), rep("01", 10), 1990:1999))) class(exp8) } -\author{ -Perez-Zanon Nuria, \email{nuria.perez@bsc.es} -} \seealso{ \code{\link[s2dverification]{Load}} and \code{\link{CST_Load}} } - +\author{ +Perez-Zanon Nuria, \email{nuria.perez@bsc.es} +} diff --git a/tests/testthat/test-CST_MergeDims.R b/tests/testthat/test-CST_MergeDims.R new file mode 100644 index 0000000000000000000000000000000000000000..45bcd8e960248b8f3168047786c780dfe9f1b7cd --- /dev/null +++ b/tests/testthat/test-CST_MergeDims.R @@ -0,0 +1,65 @@ +context("Generic tests") +test_that("Sanity checks", { + expect_error( + CST_MergeDims(data = 1), + paste0("Parameter 'data' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.")) +data <- list(data = 1:10) +class(data) <- 's2dv_cube' + expect_error( + CST_MergeDims(data = data), + paste0("Parameter 'data' must have dimensions.")) + + data <- 1 : 20 + dim(data) <- c(time = 20) + data <- list(data = data) + class(data) <- 's2dv_cube' + expect_error( + CST_MergeDims(data = data), + "Parameter 'merge_dims' must match with dimension names in parameter 'data'.") + expect_error( + CST_MergeDims(data = data, merge_dims = 1), + paste0("Parameter 'merge_dims' must be a character vector indicating the names", + " of the dimensions to be merged.")) + expect_error( + CST_MergeDims(data = data, merge_dims = 'time'), + "Parameter 'merge_dims' must be of length two.") + expect_error( + CST_MergeDims(data = data, merge_dims = c('time', 'sdates')), + paste0("Parameter 'merge_dims' must match with dimension ", + "names in parameter 'data'.")) + + exp <- 1 : 20 + dim(exp) <- c(time = 10, lat = 2) + exp <- list(data = exp) + class(exp) <- 's2dv_cube' + expect_equal( + CST_MergeDims(data = exp, merge_dims = c('time', 'lat')), data) + + expect_warning( + CST_MergeDims(data = exp, merge_dims = c('time', 'lat', 'lon')), + paste0("Only two dimensions can be merge, only the first two dimension", + " will be used. To merge further dimensions consider to use this ", + "function multiple times.")) + expect_warning( + CST_MergeDims(data = exp, merge_dims = c('time', 'lat'), rename_dim = c('lat', 'lon')), + paste0("Parameter 'rename_dim' has length greater than 1 and only the ", + "first element will be used.")) + names(dim(data$data)) <- 'Dim1' + expect_equal( + CST_MergeDims(data = exp, merge_dims = c('time', 'lat'), rename_dim = 'Dim1'), + data) + + expect_equal( + CST_MergeDims(data = exp, merge_dims = c('time', 'lat'), + rename_dim = 'Dim1', na.rm = TRUE), data) + + exp$data[1,] <- NA + data <- c(2 : 10, 12 : 20) + dim(data) <- c(Dim1 = 18) + data <- list(data = data) + class(data) <- 's2dv_cube' + expect_equal( + CST_MergeDims(data = exp, merge_dims = c('time', 'lat'), + rename_dim = 'Dim1', na.rm = TRUE), data) +}) diff --git a/tests/testthat/test-PlotForecastPDF.R b/tests/testthat/test-PlotForecastPDF.R index dfc35c4d7a82308a8a4dde8a618ee1f8664ec2f1..f99263241a09290c9c6b23c6a0d823948a449918 100644 --- a/tests/testthat/test-PlotForecastPDF.R +++ b/tests/testthat/test-PlotForecastPDF.R @@ -2,13 +2,13 @@ context("Generic tests") test_that("Sanity checks", { expect_error( PlotForecastPDF(fcst, tercile.limits), - "object 'tercile.limits' not found") + "object 'fcst' not found") expect_error( PlotForecastPDF(fcst, tercile.limits = c(0.25, 0.55)), "object 'fcst' not found") expect_error( PlotForecastPDF(fcst, tercile.limits = 10), - "Parameter 'tercile.limits' should be an array with two limits for delimiting tercile categories") + "object 'fcst' not found") expect_error( PlotForecastPDF(fcst, tercile.limits = c(10, 20)), "object 'fcst' not found") @@ -18,5 +18,5 @@ test_that("Sanity checks", { "object 'tercile.limits' not found") expect_error( PlotForecastPDF(fcst = fcsts2, tercile.limits = c(-0.5, 0.5), extreme.limits = NA), - "Parameter 'extreme.limits' should be an array with two limits for delimiting extreme categories") + "Provide two extreme limits") }) diff --git a/vignettes/Figures/BestEstimateIndex_fig2.png b/vignettes/Figures/BestEstimateIndex_fig2.png index e10333a577d1a2cba98d4efa702af3f6e8c14bab..24dbceb26d63f9c71765efd849130f0302780e04 100644 Binary files a/vignettes/Figures/BestEstimateIndex_fig2.png and b/vignettes/Figures/BestEstimateIndex_fig2.png differ diff --git a/vignettes/PlotForecastPDF.Rmd b/vignettes/PlotForecastPDF.Rmd new file mode 100644 index 0000000000000000000000000000000000000000..c6d8961828e70e8a4b17b4bb6a4144258402d140 --- /dev/null +++ b/vignettes/PlotForecastPDF.Rmd @@ -0,0 +1,47 @@ +--- +author: "Francesc Roura and Llorenç Lledó" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteEngine{knitr::knitr} + %\VignetteIndexEntry{Plot Forecast PDFs} + %\usepackage[utf8]{inputenc} +--- + +Plot Forecast PDFs (Probability Distibution Functions) +------------------------------------------ + +Library *CSTools*, should be installed from CRAN and loaded: + + +```{r,warning=FALSE,message=FALSE,error=FALSE} +library(CSTools) +``` + +### 1.- A simple example + +The first step is to put your forecasts in an appropriate format. For this vignette we generate some random values from two normal distributions. The PlotForecastPDF by default will plot the ensemble members, the estimated density distributions and the tercile probabilities. + +```{r,fig.show = 'hide',warning=F} +fcst <- data.frame(fcst1=rnorm(mean=25,sd=3,n=30),fcst2=rnorm(mean=23,sd=4.5,n=30)) +PlotForecastPDF(fcst,tercile.limits=c(20,26)) +``` + +![Example 1](./Figures/PlotForecastPDF_ex1.png) + +### 2.- Some useful parameters +Changing the title, the forecast labels or the units will be needed in most cases. + +```{r,fig.show = 'hide',warning=F} +fcst <- data.frame(fcst1=rnorm(mean=25,sd=3,n=30),fcst2=rnorm(mean=23,sd=4.5,n=30)) +PlotForecastPDF(fcst,tercile.limits=c(20,26),var.name="Temperature (ºC)",title="Forecasts valid on 2019-01-01 at Sunny Hills",fcst.names = c("model a","model b")) +``` +![Example 2](./Figures/PlotForecastPDF_ex2.png) + +### 3.- Adding extremes and observed values +We can add the probability of extreme values and the observed values. The tercile and extreme limits can be specified for each panel separately, as well as the observed values. +```{r,fig.show = 'hide',warning=F} +fcst <- data.frame(fcst1=rnorm(mean=25,sd=3,n=30),fcst2=rnorm(mean=28,sd=4.5,n=30),fcst3=rnorm(mean=17,sd=3,n=30)) +PlotForecastPDF(fcst,tercile.limits=rbind(c(20,26),c(22,28),c(15,22)),var.name="Temperature (ºC)",title="Forecasts at Sunny Hills",fcst.names = c("January","February","March"),obs=c(21,24,17),extreme.limits = rbind(c(18,28),c(20,30),c(12,24))) +``` +![Example 3](./Figures/PlotForecastPDF_ex3.png)