From e78f9c1dd67fda68c93c057ab3a293c49f1d4002 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 10 May 2021 12:04:04 +0200 Subject: [PATCH 1/3] Allow to choose a window of time steps --- NAMESPACE | 1 + R/Threshold.R | 32 +++++++++++++++++++++++++++----- man/CST_Threshold.Rd | 3 +++ man/Threshold.Rd | 3 +++ 4 files changed, 34 insertions(+), 5 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 133942a..7b752c7 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 c09ebe0..3ae72e4 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 bff93fa..738d188 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 c2f3872..12cf88f 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}.} -- GitLab From 500ed4f5b804654fa350e900be538afe8a30ec5c Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Fri, 4 Jun 2021 16:04:41 +0200 Subject: [PATCH 2/3] included test for the temporal window --- tests/testthat/test-Threshold.R | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/tests/testthat/test-Threshold.R b/tests/testthat/test-Threshold.R index 44f1291..46ae8e1 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_perc)[1]), + as.numeric(quantile(data[,,1], quantile_perc))) + ## Moving window of 3 days (first day) + expect_equal(as.numeric(Threshold(data = data, threshold = quantile_perc, window = 3)[1]), + as.numeric(quantile(data[,,1:3], quantile_perc))) + ## Moving window of 5 days (first day) + expect_equal(as.numeric(Threshold(data = data, threshold = quantile_perc, window = 5)[1]), + as.numeric(quantile(data[,,1:5], quantile_perc))) + ## No moving window (second day) + expect_equal(as.numeric(Threshold(data = data, threshold = quantile_perc)[2]), + as.numeric(quantile(data[,,2], quantile_perc))) + ## Moving window of 3 days (second day) + expect_equal(as.numeric(Threshold(data = data, threshold = quantile_perc, window = 3)[2]), + as.numeric(quantile(data[,,2:4], quantile_perc))) + ## Moving window of 5 days (second day) + expect_equal(as.numeric(Threshold(data = data, threshold = quantile_perc, window = 5)[2]), + as.numeric(quantile(data[,,2:6], quantile_perc))) + +}) -- GitLab From dc4111ebb8b6e3ba33395b6f37890487a35e7104 Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Fri, 4 Jun 2021 16:18:27 +0200 Subject: [PATCH 3/3] fixed variable name --- tests/testthat/test-Threshold.R | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/tests/testthat/test-Threshold.R b/tests/testthat/test-Threshold.R index 46ae8e1..592a3de 100644 --- a/tests/testthat/test-Threshold.R +++ b/tests/testthat/test-Threshold.R @@ -55,22 +55,22 @@ 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_perc)[1]), - as.numeric(quantile(data[,,1], quantile_perc))) + 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_perc, window = 3)[1]), - as.numeric(quantile(data[,,1:3], quantile_perc))) + 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_perc, window = 5)[1]), - as.numeric(quantile(data[,,1:5], quantile_perc))) + 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_perc)[2]), - as.numeric(quantile(data[,,2], quantile_perc))) + 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_perc, window = 3)[2]), - as.numeric(quantile(data[,,2:4], quantile_perc))) + 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_perc, window = 5)[2]), - as.numeric(quantile(data[,,2:6], quantile_perc))) + expect_equal(as.numeric(Threshold(data = data, threshold = quantile_percentile, window = 5)[2]), + as.numeric(quantile(data[,,2:6], quantile_percentile))) }) -- GitLab