diff --git a/NAMESPACE b/NAMESPACE index 133942a9fa0a8c645c323c2b3755abdc386efe8d..7b752c7bc4e68351dc651a19f39be0b93cb89375 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -31,5 +31,6 @@ importFrom(s2dv,InsertDim) importFrom(s2dv,Reorder) importFrom(stats,approxfun) importFrom(stats,ecdf) +importFrom(stats,embed) importFrom(stats,quantile) importFrom(utils,read.delim) diff --git a/R/Threshold.R b/R/Threshold.R index c09ebe004fe6424f41edcf1623f058ca24ebb225..3ae72e4c40dfd436eaed176548a284532df667d0 100644 --- a/R/Threshold.R +++ b/R/Threshold.R @@ -5,6 +5,7 @@ #' #'@param data an 's2dv_cube' object as provided function \code{CST_Load} in package CSTools. #'@param threshold a single scalar or vector indicating the relative threshold(s). +#'@param window a integer indicating the number of time steps to use as sample when computing the threshold. By default, one time step is used. #'@param start an optional parameter to defined the initial date of the period to select from the data by providing a list of two elements: the initial date of the period and the initial month of the period. By default it is set to NULL and the indicator is computed using all the data provided in \code{data}. #'@param end an optional parameter to defined the final date of the period to select from the data by providing a list of two elements: the final day of the period and the final month of the period. By default it is set to NULL and the indicator is computed using all the data provided in \code{data}. #'@param time_dim a character string indicating the name of the temporal dimension. By default, it is set to 'ftime'. More than one dimension name matching the dimensions provided in the object \code{data$data} can be specified. This dimension is required to subset the data in a requested period. @@ -31,7 +32,7 @@ #' as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day')) #'exp_probs <- CST_Threshold(exp, threshold, start = list(21, 4), end = list(21, 6)) #'@export -CST_Threshold <- function(data, threshold, start = NULL, end = NULL, +CST_Threshold <- function(data, threshold, window = NULL, start = NULL, end = NULL, time_dim = 'ftime', memb_dim = 'member', sdate_dim = 'sdate', na.rm = FALSE, ncores = NULL) { if (!inherits(data, 's2dv_cube')) { @@ -53,7 +54,8 @@ CST_Threshold <- function(data, threshold, start = NULL, end = NULL, } } } - thres <- Threshold(data$data, threshold, data$Dates[[1]], start, end, + thres <- Threshold(data$data, threshold, window = window, + data$Dates[[1]], start, end, time_dim = time_dim, memb_dim = memb_dim, sdate_dim = sdate_dim, na.rm = na.rm, ncores = ncores) data$data <- thres @@ -71,6 +73,7 @@ CST_Threshold <- function(data, threshold, start = NULL, end = NULL, #' #'@param data a multidimensional array with named dimensions. #'@param threshold a single scalar or vector indicating the relative threshold(s). +#'@param window a integer indicating the number of time steps to use as sample when computing the threshold. By default, one time step is used. #'@param dates a vector of dates or a multidimensional array of dates with named dimensions matching the dimensions on parameter 'data'. By default it is NULL, to select a period this parameter must be provided. #'@param start an optional parameter to defined the initial date of the period to select from the data by providing a list of two elements: the initial date of the period and the initial month of the period. By default it is set to NULL and the indicator is computed using all the data provided in \code{data}. #'@param end an optional parameter to defined the final date of the period to select from the data by providing a list of two elements: the final day of the period and the final month of the period. By default it is set to NULL and the indicator is computed using all the data provided in \code{data}. @@ -83,7 +86,7 @@ CST_Threshold <- function(data, threshold, start = NULL, end = NULL, #'@return A multidimensional array with named dimensions. #' #'@import multiApply -#'@importFrom stats quantile +#'@importFrom stats quantile embed #' #'@examples #'threshold <- 0.9 @@ -93,7 +96,8 @@ CST_Threshold <- function(data, threshold, start = NULL, end = NULL, #'data <- array(rnorm(1 * 3 * 214 * 2), c(member = 1, sdate = 3, time = 214, lon = 2)) #'res <- Threshold(data, threshold) #'@export -Threshold <- function(data, threshold, dates = NULL, start = NULL, end = NULL, +Threshold <- function(data, threshold, window = NULL, + dates = NULL, start = NULL, end = NULL, time_dim = 'time', memb_dim = 'member', sdate_dim = 'sdate', na.rm = FALSE, ncores = NULL) { if (is.null(data)) { @@ -125,12 +129,30 @@ Threshold <- function(data, threshold, dates = NULL, start = NULL, end = NULL, data <- SelectPeriodOnData(data, dates, start, end, time_dim = time_dim, ncores = ncores) } - } + } + if (!is.null(window) && window != 1) { + if (!is.numeric(window)) { + stop("Parameter 'window' must be an integer.") + } + if (time_dim %in% names(dim(data))) { + dims <- dim(data) + positions <- embed(1:dim(data)[time_dim], window) + names(dim(positions)) <- c(time_dim, 'window') + data <- ClimProjDiags::Subset(data, along = time_dim, indices = positions) + postime <- which(names(dims) == time_dim) + dims <- dims[-postime] + dims <- append(dims, dim(positions), postime) + dim(data) <- dims + } + } if (!is.null(memb_dim)) { dimensions <- c(sdate_dim, memb_dim) } else { dimensions <- sdate_dim } + if (!is.null(window) && window != 1) { + dimensions <- c(dimensions, 'window') + } if (length(threshold) == 1) { thres <- Apply(data, target_dims = dimensions, fun = function(x) {quantile(as.vector(x), threshold, na.rm)})$output1 diff --git a/man/CST_Threshold.Rd b/man/CST_Threshold.Rd index bff93fa02f40b70111feb0785de243bb378d5d93..738d1883a8a0469cbd893a0d67383f38dc41e3cc 100644 --- a/man/CST_Threshold.Rd +++ b/man/CST_Threshold.Rd @@ -7,6 +7,7 @@ CST_Threshold( data, threshold, + window = NULL, start = NULL, end = NULL, time_dim = "ftime", @@ -21,6 +22,8 @@ CST_Threshold( \item{threshold}{a single scalar or vector indicating the relative threshold(s).} +\item{window}{a integer indicating the number of time steps to use as sample when computing the threshold. By default, one time step is used.} + \item{start}{an optional parameter to defined the initial date of the period to select from the data by providing a list of two elements: the initial date of the period and the initial month of the period. By default it is set to NULL and the indicator is computed using all the data provided in \code{data}.} \item{end}{an optional parameter to defined the final date of the period to select from the data by providing a list of two elements: the final day of the period and the final month of the period. By default it is set to NULL and the indicator is computed using all the data provided in \code{data}.} diff --git a/man/Threshold.Rd b/man/Threshold.Rd index c2f3872c25d77bd4f7a6a1a1ddd6fd30ac7d229d..12cf88f1e494a13332a440ba3b65c5bdb71937ab 100644 --- a/man/Threshold.Rd +++ b/man/Threshold.Rd @@ -7,6 +7,7 @@ Threshold( data, threshold, + window = NULL, dates = NULL, start = NULL, end = NULL, @@ -22,6 +23,8 @@ Threshold( \item{threshold}{a single scalar or vector indicating the relative threshold(s).} +\item{window}{a integer indicating the number of time steps to use as sample when computing the threshold. By default, one time step is used.} + \item{dates}{a vector of dates or a multidimensional array of dates with named dimensions matching the dimensions on parameter 'data'. By default it is NULL, to select a period this parameter must be provided.} \item{start}{an optional parameter to defined the initial date of the period to select from the data by providing a list of two elements: the initial date of the period and the initial month of the period. By default it is set to NULL and the indicator is computed using all the data provided in \code{data}.} diff --git a/tests/testthat/test-Threshold.R b/tests/testthat/test-Threshold.R index 44f12914d3beee5647760f4af6af16d39c68453a..592a3de423e4921888aefb99cb832476121dadd0 100644 --- a/tests/testthat/test-Threshold.R +++ b/tests/testthat/test-Threshold.R @@ -49,3 +49,28 @@ test_that("Seasonal forecasts", { expect_equal(round(exp1_thresholdP[, 2, 2]), c(281, 279, 276)) }) + +test_that("Temporal window", { + + data <- array(1:(2*3*20), c(member = 2, sdate = 3, time = 20)) + quantile_percentile <- 0.8 + ## No moving window (first day) + expect_equal(as.numeric(Threshold(data = data, threshold = quantile_percentile)[1]), + as.numeric(quantile(data[,,1], quantile_percentile))) + ## Moving window of 3 days (first day) + expect_equal(as.numeric(Threshold(data = data, threshold = quantile_percentile, window = 3)[1]), + as.numeric(quantile(data[,,1:3], quantile_percentile))) + ## Moving window of 5 days (first day) + expect_equal(as.numeric(Threshold(data = data, threshold = quantile_percentile, window = 5)[1]), + as.numeric(quantile(data[,,1:5], quantile_percentile))) + ## No moving window (second day) + expect_equal(as.numeric(Threshold(data = data, threshold = quantile_percentile)[2]), + as.numeric(quantile(data[,,2], quantile_percentile))) + ## Moving window of 3 days (second day) + expect_equal(as.numeric(Threshold(data = data, threshold = quantile_percentile, window = 3)[2]), + as.numeric(quantile(data[,,2:4], quantile_percentile))) + ## Moving window of 5 days (second day) + expect_equal(as.numeric(Threshold(data = data, threshold = quantile_percentile, window = 5)[2]), + as.numeric(quantile(data[,,2:6], quantile_percentile))) + +})