diff --git a/DESCRIPTION b/DESCRIPTION index d2444dab773613047218c7002202b443924f499b..44c5d483bda89aad1874e9342fd837dd4d8b71fe 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -22,7 +22,6 @@ Imports: multiApply (>= 2.0.0), PCICt, plyr, - climdex.pcic, stats Suggests: knitr, diff --git a/NAMESPACE b/NAMESPACE index 47ee4dc92d747aebb9591c68cd42617cd21a0648..dd96fe5720228804e4819bd5704fe63ee6a077c0 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -18,6 +18,5 @@ export(WaveDuration) export(WeightedCells) export(WeightedMean) import(PCICt) -import(climdex.pcic) import(multiApply) importFrom(stats,quantile) diff --git a/NEWS.md b/NEWS.md index d72e8eba87f3961e8e032cfe8eb563b16944b2b4..7a8048e546609fdca02f7e2e4e807e524c4bb79f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,6 @@ +# 0.3.2 (Release date: 2023-05-19) +- Remove climdex.pcic dependency + # 0.3.1 (Release date: 2023-03-23) - Subset(): Prioritize the dimension names from names(dim(x)) rather than attribute 'dimensions'; If the input data doesn't have dimension names, the output doesn't have either. diff --git a/R/Climdex.R b/R/Climdex.R index a13b15ba283bb6455729bd2ea0bd6ce69bb6fc46..74da9045edfe2d33ad55a90a576805b6c7d9941c 100644 --- a/R/Climdex.R +++ b/R/Climdex.R @@ -1,23 +1,35 @@ -#'Wrapper for applying the climdex routine ETCCDI climate change indices to n-dimensional arrays. +#'Wrapper for applying the climdex routine ETCCDI climate change indices to +#'n-dimensional arrays. #' -#'@description This function computes the t90p, t10p, cdd or rx5day indices from n-dimensional arrays. +#'@description This function computes the t90p, t10p, cdd or rx5day indices from +#'n-dimensional arrays. #' -#'@param data A numeric n-dimensional array containing daily maximum or minimum temperature, wind speed or precipitation amount. -#'@param metric The metric to be computed, either 't90p', 't10p', 'Wx', 'cdd' or 'rx5day'. -#'@param threshold For the 't90p' and 't10p' metrics, an array of the 90th/10th percentiles must be included. This parameter can be computed with the \code{Threshold} function. -#'@param base.range The years used for the reference period. If NULL (by default), all years are used. -#'@param dates A vector of dates with a calendar attributes. If NULL (by default), the 'time' attributes of parameter 'data' are considered. -#'@param timedim An integer number indicating the position of the time dimension in the parameter \code{data}. If NULL (by default), the dimension called 'time' in parameter \code{data} is considered as temporal dimension. +#'@param data A numeric n-dimensional array containing daily maximum or minimum +#' temperature, wind speed or precipitation amount. +#'@param metric The metric to be computed, either 't90p', 't10p', 'Wx', 'cdd' or +#' 'rx5day'. +#'@param threshold For the 't90p' and 't10p' metrics, an array of the 90th/10th +#' percentiles must be included. This parameter can be computed with the +#' \code{Threshold} function. +#'@param base.range The years used for the reference period. If NULL +#' (by default), all years are used. +#'@param dates A vector of dates with a calendar attributes. If NULL +#' (by default), the 'time' attributes of parameter 'data' are considered. +#'@param timedim An integer number indicating the position of the time dimension +#' in the parameter \code{data}. If NULL (by default), the dimension called +#' 'time' in parameter \code{data} is considered as temporal dimension. #'@param calendar A character indicating the calendar type. #'@param ncores The number of cores to be used when computing the index. #' #'@return A list of length 2: #'\itemize{ -#' \item\code{$result} {An array with the same dimensions as the input array, except for the temporal dimension which is renamed to 'year', moved to the first dimension position and reduce to annual resolution.} -#' \item\code{$years} {A vector of the corresponding years.}} +#' \item\code{$result} {An array with the same dimensions as the input array, +#' except for the temporal dimension which is renamed to 'year', moved to the +#' first dimension position and reduce to annual resolution.} +#' \item\code{$years} {A vector of the corresponding years.} +#'} #' #'@import multiApply -#'@import climdex.pcic #'@import PCICt #'@examples #'##Example synthetic data: @@ -35,8 +47,10 @@ #' seq(ISOdate(1909, 1, 1), ISOdate(1909, 1, 31), "day"), #' seq(ISOdate(1910, 1, 1), ISOdate(1910, 1, 31), "day"), #' seq(ISOdate(1911, 1, 1), ISOdate(1911, 1, 31), "day")) -#'metadata <- list(time = list(standard_name = 'time', long_name = 'time', calendar = 'gregorian', -#' units = 'days since 1970-01-01 00:00:00', prec = 'double', +#'metadata <- list(time = list(standard_name = 'time', long_name = 'time', +#' calendar = 'gregorian', +#' units = 'days since 1970-01-01 00:00:00', +#' prec = 'double', #' dim = list(list(name = 'time', unlim = FALSE)))) #'attr(time, "variables") <- metadata #'attr(data, 'Variables')$dat1$time <- time @@ -45,15 +59,15 @@ #'dim(thres) <- c(jdays = 31, lon = 2, lat = 3, model = 1) #'str(thres) #' -#' #'clim <- Climdex(data, metric = "t90p", threshold = thres) #'str(clim) #'@references David Bronaugh for the Pacific Climate Impacts Consortium (2015). #' climdex.pcic: PCIC Implementation of Climdex Routines. R package #' version 1.1-6. http://CRAN.R-project.org/package=climdex.pcic #'@export -Climdex <- function(data, metric, threshold = NULL, base.range = NULL, dates = NULL, timedim = NULL, - calendar = NULL, ncores = NULL) { +Climdex <- function(data, metric, threshold = NULL, base.range = NULL, + dates = NULL, timedim = NULL, calendar = NULL, + ncores = NULL) { if (is.null(data) | is.null(metric)) { stop("Parameters 'data' and 'metric' cannot be NULL.") } diff --git a/R/ClimdexPcicFunctions.R b/R/ClimdexPcicFunctions.R new file mode 100644 index 0000000000000000000000000000000000000000..a923b5c922415099808688e28a0daccac2766fef --- /dev/null +++ b/R/ClimdexPcicFunctions.R @@ -0,0 +1,392 @@ +#----------------------------------------------------------- +# Auxiliary functions from package climdex.pcic +# (2023-06-01) The package was removed from CRAN. To maintain ClimProjDiags, +# we copy the necessary functions here. +# We will recover the dependency when climdex.pcic is on CRAN again. +#----------------------------------------------------------- +# Package name: climdex.pcic +# Version: 1.1-11 +# Date: 2020-01-21 +# Title: PCIC Implementation of Climdex Routines +# Author: David Bronaugh for the Pacific Climate Impacts +# Consortium +# Maintainer: James Hiebert +# Depends: +# R (>= 2.12.0), +# PCICt (>= 0.5-4) +# Imports: +# methods, +# Rcpp (>= 0.11.4) +# Suggests: +# compiler, +# RUnit +# LinkingTo: Rcpp +# Description: PCIC's implementation of Climdex routines for computation of +# extreme climate indices. Further details on the extreme climate indices +# can be found at +# and in the package manual. +# License: GPL-3 +# URL: https://www.r-project.org +# LazyData: yes +# BugReports: https://github.com/pacificclimate/climdex.pcic/issues/ +# RoxygenNote: 6.0.1 +#----------------------------------------------------------- + +## Helper functions: + +# Lower overhead version of tapply +tapply.fast <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) { + FUN <- if (!is.null(FUN)) + match.fun(FUN) + + if(!is.factor(INDEX)) + stop("INDEX must be a factor.") + + if (length(INDEX) != length(X)) + stop("arguments must have same length") + + if (is.null(FUN)) + return(INDEX) + + namelist <- levels(INDEX) + ans <- lapply(split(X, INDEX), FUN, ...) + + ans <- unlist(ans, recursive = FALSE) + names(ans) <- levels(INDEX) + return(ans) +} + +# Get year +get.years <- function(dates) { + return(as.POSIXlt(dates)$year + 1900) +} + +# Get NA mask given threshold and split factor +get.na.mask <- function(x, f, threshold) { + return(c(1, NA)[1 + as.numeric(tapply.fast(is.na(x), f, function(y) { return(sum(y) > threshold) } ))]) +} + +#----------------------------------------------------------- + +# Functions + +#'Get series length at ends +#' +#'This function takes a series of boolean values and returns a list of +#'integers of the same length corresponding to the lengths at the ends of +#'sequences of TRUE values. +#' +#'It can often be useful to know how long a series of boolean values is. This +#'function provides a method of knowing where and how long such sequences are. +#' +#'@param x Sequence of booleans. +#'@param na.value Value to replace NAs with. +#'@return A vector consisting of the lengths of sequences of TRUE values at +#'the location of the last TRUE value in the sequence, and zeroes elsewhere. +#'@keywords ts climate +#'@examples +#'# Get lengths of sequences of TRUE values in a sequence +#'series.lengths <- get.series.lengths.at.ends(c(TRUE, TRUE, TRUE, FALSE, +#'TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, FALSE)) +#'@noRd +get.series.lengths.at.ends <- function(x, na.value=FALSE) { + stopifnot(is.logical(x) && is.logical(na.value)) + n <- length(x) + if(n == 1) + return(as.numeric(x)) + + res <- rep(0, n) + x[is.na(x)] <- na.value + + ## Compare series to lag-1 and lag+1 series; false added to trigger state transition from TRUE at ends of series + start <- which(x & !(c(FALSE, x[1:(n - 1)]))) + end <- which(x & !(c(x[2:n], FALSE))) + res[end] <- end - start + 1 + return(res) +} + +#'Lengths of strings of TRUE values +#' +#'Computes fraction of days above or below the baseline threshold for each +#'day, and averages them using the date factor passed in. +#' +#'This function computes fractions of days above or below baseline thresholds +#'for each day, then aggregates them using \code{date.factor}. It is used to +#'implement TN/TX 10/90p. +#' +#'@param temp Sequence of temperature values. +#'@param dates Sequence of associated dates. +#'@param jdays Sequence of associated days of year. +#'@param date.factor Factor to aggregate data using. +#'@param threshold.outside.base Sequence of thresholds to be used for data +#' outside the base period. +#'@param base.thresholds Data structure containing sets of thresholds to be +#' used inside the base period; see \link{climdexInput-class}. +#'@param base.range Date range (type PCICt) of the baseline period. +#'@param op Comparison operator to use. +#'@param max.missing.days Maximum number of NA values per time period. +#'@return A vector consisting of the mean fraction of days above or below the +#' supplied set of thresholds. +#'@note If date.factor is omitted, daily series will be returned. +#'@keywords ts climate +#'@noRd +percent.days.op.threshold <- function(temp, dates, jdays, date.factor, + threshold.outside.base, base.thresholds, + base.range, op='<', max.missing.days) { + f <- match.fun(op) + dat <- f(temp, threshold.outside.base[jdays]) + + inset <- dates >= base.range[1] & dates <= base.range[2] + ## Don't use in-base thresholds with data shorter than two years; no years to replace with. + if(sum(inset) > 0 && length(dates) >= 360 * 2) { + jdays.base <- jdays[inset] + years.base <- get.years(dates[inset]) + + ## Get number of base years, subset temp data to base period only. + temp.base <- temp[inset] + years.base.range <- range(years.base) + byrs <- (years.base.range[2] - years.base.range[1] + 1) + + ## Linearize thresholds, then compare them to the temperatures + bdim <- dim(base.thresholds) + dim(base.thresholds) <- c(bdim[1] * bdim[2], bdim[3]) + yday.byr.indices <- jdays.base + (years.base - get.years(base.range)[1]) * bdim[1] + f.result <- f(rep(temp.base, byrs - 1), base.thresholds[yday.byr.indices,]) + dim(f.result) <- c(length(yday.byr.indices), bdim[3]) + + ## Chop up data along the 2nd dim into a list; sum elements of the list + dat[inset] <- rowSums(f.result, na.rm=TRUE) / (byrs - 1) + } + dat[is.nan(dat)] <- NA + if(missing(date.factor)) + return(dat) + na.mask <- get.na.mask(dat, date.factor, max.missing.days) + ## FIXME: Need to monthly-ize the NA mask calculation, which will be ugly. + ret <- tapply.fast(dat, date.factor, mean, na.rm=TRUE) * 100 * na.mask + ret[is.nan(ret)] <- NA + return(ret) +} + +#'Sum of spell lengths exceeding daily threshold +#' +#'This function returns the number of spells of more than \code{min.length} +#'days which exceed or are below the given threshold. +#' +#'@details +#'This routine compares data to the thresholds using the given operator, +#'generating a series of TRUE or FALSE values; these values are then filtered +#'to remove any sequences of less than \code{min.length} days of TRUE values. +#'It then computes the lengths of the remaining sequences of TRUE values +#'(spells) and sums their lengths. +#' +#'The \code{spells.can.span.years} option controls whether spells must always +#'terminate at the end of a period, or whether they may continue until the +#'criteria ceases to be met or the end of the data is reached. The default for +#'fclimdex is FALSE. +#' +#'@param daily.temp Data to compute index on. +#'@param date.factor Date factor to split by. +#'@param jdays Timeseries of days of year. +#'@param thresholds The thresholds to compare to. +#'@param op The operator to use to compare data to threshold. +#'@param min.length The minimum spell length to be considered. +#'@param spells.can.span.years Whether spells can span years. +#'@param max.missing.days Maximum number of NA values per time period. +#'@return A timeseries of maximum spell lengths for each period. +#'@seealso \code{\link{climdex.wsdi}}. +#'@keywords ts climate +#'@examples +#'prec.dat <- c(0.1, 3.0, 4.3, 1.9, 1.3, 6.0, 0, 0, 4.0, 1) +#'phony.date.factor <- factor(rep(1:2, each=5)) +#' +#'## With spells spanning years... +#'alttedi <- threshold.exceedance.duration.index(prec.dat, +#'phony.date.factor, rep(1:5, 2), rep(1, 5), ">=", 2, TRUE, 1) +#' +#'## Without spells spanning years... +#'tedi <- threshold.exceedance.duration.index(prec.dat, phony.date.factor, +#'rep(1:5, 2), rep(1, 5), ">=", 2, FALSE, 1) +#'@noRd +threshold.exceedance.duration.index <- function(daily.temp, date.factor, jdays, + thresholds, op=">", + min.length=6, + spells.can.span.years=TRUE, + max.missing.days) { + stopifnot(is.numeric(c(daily.temp, thresholds, min.length)), is.factor(date.factor), + is.function(match.fun(op)), + min.length > 0) + f <- match.fun(op) + na.mask <- get.na.mask(is.na(daily.temp + thresholds[jdays]), date.factor, max.missing.days) + + if(spells.can.span.years) { + periods <- select.blocks.gt.length(f(daily.temp, thresholds[jdays]), min.length - 1) + return(tapply.fast(periods, date.factor, sum) * na.mask) + } else { + ## fclimdex behaviour... + return(tapply.fast(1:length(daily.temp), date.factor, function(idx) { sum(select.blocks.gt.length(f(daily.temp[idx], thresholds[jdays[idx]]), min.length - 1)) } ) * na.mask) + } +} + +#'Number of days (less than, greater than, etc) a threshold +#' +#'@description +#'Produces sums of values that exceed (or are below) the specified threshold. +#' +#'@details +#'This function takes a data series, the number of days in the running window, +#'a date factor to aggregate by, and an optional modifier parameter +#'(center.mean.on.last.day). It computes the n-day running sum of +#'precipitation and returns the maximum n-day total precipitation per unit +#'time, as defined by \code{date.factor}. +#' +#'@param daily.prec Daily timeseries of precipitation. +#'@param date.factor Factor to aggregate by. +#'@param ndays Number of days in the running window. +#'@param center.mean.on.last.day Whether to center the n-day running mean on +#' the last day of the series, instead of the middle day. +#'@return A vector consisting of the maximum n-day sum of precipitation per +#'time interval. +#'@keywords ts climate +#'@noRd +nday.consec.prec.max <- function(daily.prec, date.factor, ndays, + center.mean.on.last.day=FALSE) { + if (ndays == 1) { + return(suppressWarnings(tapply.fast(daily.prec, date.factor, max, na.rm=TRUE))) + } + ## Ends of the data will be de-emphasized (padded with zero precip data); NAs replaced with 0 + daily.prec[is.na(daily.prec)] <- 0 + prec.runsum <- running.mean(daily.prec, ndays) + prec.runsum[is.na(prec.runsum)] <- 0 + if(center.mean.on.last.day) { + k2 = ndays %/% 2 + prec.runsum <- c(rep(0, k2), prec.runsum[1:(length(prec.runsum) - k2)]) + } + return(tapply.fast(prec.runsum, date.factor, max) * ndays) +} + +#'Maximum spell length +#' +#'@description +#'This function returns the longest string of days which exceed or are below +#'the given threshold. +#' +#'@details +#'This routine compares data to the threshold using the given operator, +#'generating a series of TRUE or FALSE values. It then computes the lengths of +#'sequences of TRUE values (spells) and chooses the longest spell in each +#'period (as defined by date.factor). +#' +#'The \code{spells.can.span.years} option controls whether spells must always +#'terminate at the end of a period, or whether they may continue until the +#'criteria ceases to be met or the end of the data is reached. The default for +#'fclimdex is TRUE. +#' +#'@param daily.prec Data to compute index on. +#'@param date.factor Date factor to split by. +#'@param threshold The threshold to compare to. +#'@param op The operator to use to compare data to threshold. +#'@param spells.can.span.years Whether spells can span years. +#'@return A timeseries of maximum spell lengths for each period. +#'@seealso \code{\link{climdex.cdd}}. +#'@keywords ts climate +#'@examples +#' +#'prec.dat <- c(0.1, 3.0, 4.3, 1.9, 1.3, 6.0, 0, 0, 4.0, 1) +#'phony.date.factor <- factor(rep(1:2, each=5)) +#' +#'## With spells spanning years... +#'cwd <- spell.length.max(prec.dat, phony.date.factor, 1, ">=", TRUE) +#' +#'## Without spells spanning years... +#'altcwd <- spell.length.max(prec.dat, phony.date.factor, 1, ">=", FALSE) +#'@noRd +spell.length.max <- function(daily.prec, date.factor, threshold, op, + spells.can.span.years) { + bools <- match.fun(op)(daily.prec, threshold) + + if(spells.can.span.years) { + all.true <- tapply.fast(bools, date.factor, all) + max.spell <- tapply.fast(get.series.lengths.at.ends(bools), date.factor, max) + + ## Mask out values which are in the middle of a spell with NA + na.mask <- c(1, NA)[as.integer((max.spell == 0) & all.true) + 1] + return(max.spell * na.mask) + } else { + return(tapply.fast(bools, date.factor, function(x) { max(get.series.lengths.at.ends(x)) })) + } +} + +#'Select blocks of TRUE values of sufficient length. +#' +#'Produces a sequence of booleans of the same length as input, with sequences +#'of TRUE values shorter than n replaced with FALSE. +#' +#'This function takes a series of booleans and returns a sequence of booleans +#'of equal length, with all sequences of TRUE of length \code{n} or shorter +#'replaced with sequences of FALSE. NA values are replaced with +#'\code{na.value}. +#' +#'@param d Sequence of booleans. +#'@param n Longest sequence of TRUE to replace with FALSE. +#'@param na.value Values to replace NAs with. +#'@return A vector of booleans, with the length \code{n} or less sequences of +#'TRUE replaced with FALSE. +#'@keywords ts climate +#'@examples +#' +#'## Return only the first sequence of TRUE... second sequence will be FALSE. +#'foo <- select.blocks.gt.length(c(rep(TRUE, 4), FALSE, rep(TRUE, 3)), 3) +#'@noRd +select.blocks.gt.length <- function(d, n, na.value = FALSE) { + stopifnot(is.logical(d), is.numeric(n)) + + if(n < 1) + return(d) + + if(n >= length(d)) + return(rep(FALSE, length(d))) + + d[is.na(d)] <- na.value + + d2 <- Reduce(function(x, y) { return(c(rep(FALSE, y), d[1:(length(d) - y)]) & x) }, 1:n, d) + return(Reduce(function(x, y) { return(c(d2[(y + 1):length(d2)], rep(FALSE, y)) | x) }, 1:n, d2)) +} + +#'Running Mean of a Vector +#' +#'@description Calculates the running means of a vector with a shifting window +#' +#'@details Returns a new vector the same length as vec, where the ith +#'element is the mean of the bin of elements centered at the ith element +#'of the original vector. Means cannot be calculated for elements less +#'than half the width of the bin from the beginning or end of the vector; +#'the result vector has NA in those positions. +#' +#'@param vec A vector +#'@param bin The number of entries to average over for each mean +#' +#'@return a vector containing the running mean of bin elements of vec +#' +#'@examples +#'running.mean(c(1, 2, 3, 4, 5, 6), 2) +#'running.mean(c(5, 5, 5, 5, 5), 4) +#'@noRd +running.mean <- function(vec, bin) { + vec = as.vector(vec) + len = length(vec) + if (bin<=1) { + return (vec) + } + if (bin > len) { + bin = len + } + left.bin = bin%/%2 + + means = double(len) + + right.bin = bin - left.bin - 1 + means = c( sum(vec[1:bin]), diff(vec,bin) ) # find the first sum and the differences from it + means = cumsum(means)/bin # apply precomputed differences + means = c(rep(NA,left.bin), means, rep(NA,right.bin)) # extend to original vector length + return(means) +} diff --git a/R/Extremes.R b/R/Extremes.R index facb6cf4a568356588948eacf246a706ec4dc446..149fd96c9179c3cf0d066d0977a5d4493efa7191 100644 --- a/R/Extremes.R +++ b/R/Extremes.R @@ -1,46 +1,65 @@ #'Sum of spell lengths exceeding daily threshold for n-dimensional arrays #' -#'@description This function returns the number of spells of more than \code{min.length} days which exceed or are below the given \code{threshold} from daily data. +#'@description This function returns the number of spells of more than +#'\code{min.length} days which exceed or are below the given \code{threshold} +#'from daily data. #' #'@param data A n-dimensional array containing daily data. -#'@param threshold A n-dimensional array with the threshold to be/not to be reach, usually given by the a percentile computed with the \code{Threshold} function. +#'@param threshold A n-dimensional array with the threshold to be/not to be +#' reach, usually given by the a percentile computed with the \code{Threshold} +#' function. #'@param op The operator to use to compare data to threshold. #'@param min.length The minimum spell length to be considered. #'@param spells.can.span.years Whether spells can span years. #'@param max.missing.days Maximum number of NA values per time period. -#'@param dates A vector of dates with a calendar attributes. If NULL (by default), the 'time' attributes of parameter 'data' are considered. -#'@param timedim An integer number indicating the position of the time dimension in the parameter \code{data}. If NULL (by default), the dimension called 'time' in parameter \code{data}. +#'@param dates A vector of dates with a calendar attributes. If NULL +#' (by default), the 'time' attributes of parameter 'data' are considered. +#'@param timedim An integer number indicating the position of the time dimension +#' in the parameter \code{data}. If NULL (by default), the dimension called +#' 'time' in parameter \code{data}. #'@param calendar A character indicating the calendar type. #'@param ncores The number of cores to be used when computing the extreme. #' #'@return A list of length 2: #'\itemize{ -#' \item\code{$output1}{An array with the same dimensions as the original \code{data}, except the time dimension which is reduced to annual resolution given a timeseries of maximum spell lengths for each year.} +#' \item\code{$output1}{An array with the same dimensions as the original +#' \code{data}, except the time dimension which is reduced to annual resolution +#' given a timeseries of maximum spell lengths for each year.} #' \item\code{$year}{A vector indicating the corresponding years.} #'} -#'@details This routine compares data to the thresholds using the given operator, generating a series of TRUE or FALSE values; these values are then filtered to remove any sequences of less than \code{min.length} days of TRUE values. It then computes the lengths of the remaining sequences of TRUE values (spells) and sums their lengths. The \code{spells.can.spa .years} option controls whether spells must always terminate at the end of a period, or whether they may continue until the criteria ceases to be met or the end of the data is reached. The default for fclimdex is FALSE. +#'@details This routine compares data to the thresholds using the given +#'operator, generating a series of TRUE or FALSE values; these values are then +#'filtered to remove any sequences of less than \code{min.length} days of TRUE +#'values. It then computes the lengths of the remaining sequences of TRUE values +#'(spells) and sums their lengths. The \code{spells.can.spa .years} option +#'controls whether spells must always terminate at the end of a period, or +#'whether they may continue until the criteria ceases to be met or the end of +#'the data is reached. The default for fclimdex is FALSE. #' #'@import multiApply #'@import PCICt -#'@import climdex.pcic #'@examples #'##Example synthetic data: #'data <- 1:(2 * 3 * 310 * 1) #'dim(data) <- c(time = 310, lon = 2, lat = 3, model = 1) #'time <- as.POSIXct(paste(sort(rep(1902:1911, 31)), 1, 1:31, sep = "-"), tz = "CET") -#'metadata <- list(time = list(standard_name = 'time', long_name = 'time', calendar = 'noleap', +#'metadata <- list(time = list(standard_name = 'time', long_name = 'time', +#' calendar = 'noleap', #' units = 'days since 1970-01-01 00:00:00', prec = 'double', #' dim = list(list(name = 'time', unlim = FALSE)))) #'attr(time, "variables") <- metadata #'attr(data, 'Variables')$dat1$time <- time #'threshold <- Threshold(data, dates = NULL, base.range = NULL, qtiles = 0.9, ncores = NULL) -#'res <- Extremes(data, threshold = threshold, op = ">", min.length = 6, spells.can.span.years = TRUE, -#' max.missing.days = 5, ncores = NULL) +#'res <- Extremes(data, threshold = threshold, op = ">", min.length = 6, +#' spells.can.span.years = TRUE, max.missing.days = 5, +#' ncores = NULL) #'str(res) #' #'@export -Extremes <- function(data, threshold, op = ">", min.length = 6, spells.can.span.years = TRUE, max.missing.days = 5, - dates = NULL, timedim = NULL, calendar = NULL, ncores = NULL) { +Extremes <- function(data, threshold, op = ">", min.length = 6, + spells.can.span.years = TRUE, max.missing.days = 5, + dates = NULL, timedim = NULL, calendar = NULL, + ncores = NULL) { if (is.null(data) | is.null(threshold)) { stop("Parameter 'data' and 'threshold' cannot be NULL.") } diff --git a/R/WaveDuration.R b/R/WaveDuration.R index 23c2d97e68d0753aee49fe13660c651070f6544b..76874f9775641d5c1509fdb465dc90d62a72e5a5 100644 --- a/R/WaveDuration.R +++ b/R/WaveDuration.R @@ -1,40 +1,58 @@ #'Heat and cold waves duration for n-dimensional arrays #' -#'This function computes the duration of a heat/cold wave as the number of consecutive days for which the maximum/minimum temperature is exceeding/below a threshold over a minimum number of days in month or seasonal resolution. +#'This function computes the duration of a heat/cold wave as the number of +#'consecutive days for which the maximum/minimum temperature is exceeding/below +#'a threshold over a minimum number of days in month or seasonal resolution. #' -#'@param data A numeric n-dimensional array containing daily maximum or minimum temperature -#'@param threshold An array with the threshold to be/not to be reach, usually given by the 90th/10th percentiles for heat/cold waves computed with the \code{Threshold} function. -#'@param op A character ">" (by default) or ">=" for heat waves and "<" or "<=" for cold waves indicating the operator must be used to compare data to threshold. -#'@param spell.length A number indicating the number of consecutive days with extreme temperature to be considered heat or cold wave. -#'@param by.seasons If TRUE (by default), the wave duration is computed for each season (DJF/MAM/JJA/SON) separately. If FALSE is specified, the monthly wave duration is computed. -#'@param dates A vector of dates including calendar attributes. If NULL (by default), the 'time' attributes of parameter 'data' is used. +#'@param data A numeric n-dimensional array containing daily maximum or minimum +#' temperature +#'@param threshold An array with the threshold to be/not to be reach, usually +#' given by the 90th/10th percentiles for heat/cold waves computed with the +#' \code{Threshold} function. +#'@param op A character ">" (by default) or ">=" for heat waves and "<" or "<=" +#' for cold waves indicating the operator must be used to compare data to +#' threshold. +#'@param spell.length A number indicating the number of consecutive days with +#' extreme temperature to be considered heat or cold wave. +#'@param by.seasons If TRUE (by default), the wave duration is computed for each +#' season (DJF/MAM/JJA/SON) separately. If FALSE is specified, the monthly wave +#' duration is computed. +#'@param dates A vector of dates including calendar attributes. If NULL (by +#' default), the 'time' attributes of parameter 'data' is used. #'@param calendar A character indicating the calendar type. #'@param ncores The number of cores to be used when computing the wave duration. #' #'@return A list of length 2: #' \itemize{ -#' \item\code{$result}{An array with the same dimensions as the input \code{data}, but with the time dimension reduce from daily to monthly or seasonal resolution depending on the selected resolution in \code{by.season}.} -#' \item\code{$years}{A vector of the years and season/months corresponding to the resolution selected in \code{by.season} and temporal length of the input \code{data}}} +#' \item\code{$result}{An array with the same dimensions as the input +#' \code{data}, but with the time dimension reduce from daily to monthly or +#' seasonal resolution depending on the selected resolution in \code{by.season}.} +#' \item\code{$years}{A vector of the years and season/months corresponding to +#' the resolution selected in \code{by.season} and temporal length of the input +#' \code{data}}} #' #'@import multiApply #'@import PCICt -#'@import climdex.pcic #'@examples #'##Example synthetic data: #'data <- 1:(2 * 3 * 31 * 5) #'dim(data) <- c(lon = 2, lat = 3, time = 31, model = 5) -#'time <- as.POSIXct(paste(paste(1900, 1, 1:31, sep = "-"), paste(12, 0, 0.0, sep = ":")), tz = "CET") -#'metadata <- list(time = list(standard_name = 'time', long_name = 'time', calendar = 'standard', +#'time <- as.POSIXct(paste(paste(1900, 1, 1:31, sep = "-"), paste(12, 0, 0.0, +#' sep = ":")), tz = "CET") +#'metadata <- list(time = list(standard_name = 'time', long_name = 'time', +#' calendar = 'standard', #' units = 'days since 1970-01-01 00:00:00', prec = 'double', #' dim = list(list(name ='time', unlim = FALSE)))) #'attr(time, "variables") <- metadata #'attr(data, 'Variables')$dat1$time <- time #'threshold <- rep(40, 31) #' -#'a <- WaveDuration(data, threshold, op = ">", spell.length = 6, by.seasons = TRUE, ncores = NULL) +#'a <- WaveDuration(data, threshold, op = ">", spell.length = 6, +#' by.seasons = TRUE, ncores = NULL) #'str(a) #'@export -WaveDuration <- function(data, threshold, op = ">", spell.length = 6, by.seasons = TRUE, dates = NULL, +WaveDuration <- function(data, threshold, op = ">", spell.length = 6, + by.seasons = TRUE, dates = NULL, calendar = NULL, ncores = NULL) { if (is.null(data) | is.null(threshold)) { stop("Parameter 'data' and 'threshold' cannot be NULL.") @@ -203,7 +221,8 @@ WaveDuration <- function(data, threshold, op = ">", spell.length = 6, by.seasons op = op, spell.length = spell.length, ncores = ncores) } else { result <- list() - result$output1 <- .WaveDuration(data = data[[1]], threshold = data[[2]], date.factor = date.factor, + result$output1 <- .WaveDuration(data = data[[1]], threshold = data[[2]], + date.factor = date.factor, jdays = jdays, op = op, spell.length = spell.length) } if (length(years) > 1) { @@ -225,7 +244,8 @@ WaveDuration <- function(data, threshold, op = ">", spell.length = 6, by.seasons .WaveDuration <- function(data, threshold, date.factor = date.factor, jdays = jdays, op = op, spell.length = spell.length) { result <- threshold.exceedance.duration.index(daily.temp = data, - date.factor = date.factor, jdays = jdays, + date.factor = date.factor, + jdays = jdays, thresholds = threshold, op = op, min.length = spell.length, spells.can.span.years = TRUE, 1) diff --git a/man/Climdex.Rd b/man/Climdex.Rd index 372712c240998579f23b70b8fbfd998656ce00f7..7afbf0b0085d5e1a54a8d72743b407013574d67a 100644 --- a/man/Climdex.Rd +++ b/man/Climdex.Rd @@ -2,7 +2,8 @@ % Please edit documentation in R/Climdex.R \name{Climdex} \alias{Climdex} -\title{Wrapper for applying the climdex routine ETCCDI climate change indices to n-dimensional arrays.} +\title{Wrapper for applying the climdex routine ETCCDI climate change indices to +n-dimensional arrays.} \usage{ Climdex( data, @@ -16,17 +17,25 @@ Climdex( ) } \arguments{ -\item{data}{A numeric n-dimensional array containing daily maximum or minimum temperature, wind speed or precipitation amount.} +\item{data}{A numeric n-dimensional array containing daily maximum or minimum +temperature, wind speed or precipitation amount.} -\item{metric}{The metric to be computed, either 't90p', 't10p', 'Wx', 'cdd' or 'rx5day'.} +\item{metric}{The metric to be computed, either 't90p', 't10p', 'Wx', 'cdd' or +'rx5day'.} -\item{threshold}{For the 't90p' and 't10p' metrics, an array of the 90th/10th percentiles must be included. This parameter can be computed with the \code{Threshold} function.} +\item{threshold}{For the 't90p' and 't10p' metrics, an array of the 90th/10th +percentiles must be included. This parameter can be computed with the +\code{Threshold} function.} -\item{base.range}{The years used for the reference period. If NULL (by default), all years are used.} +\item{base.range}{The years used for the reference period. If NULL +(by default), all years are used.} -\item{dates}{A vector of dates with a calendar attributes. If NULL (by default), the 'time' attributes of parameter 'data' are considered.} +\item{dates}{A vector of dates with a calendar attributes. If NULL +(by default), the 'time' attributes of parameter 'data' are considered.} -\item{timedim}{An integer number indicating the position of the time dimension in the parameter \code{data}. If NULL (by default), the dimension called 'time' in parameter \code{data} is considered as temporal dimension.} +\item{timedim}{An integer number indicating the position of the time dimension +in the parameter \code{data}. If NULL (by default), the dimension called +'time' in parameter \code{data} is considered as temporal dimension.} \item{calendar}{A character indicating the calendar type.} @@ -35,11 +44,15 @@ Climdex( \value{ A list of length 2: \itemize{ - \item\code{$result} {An array with the same dimensions as the input array, except for the temporal dimension which is renamed to 'year', moved to the first dimension position and reduce to annual resolution.} - \item\code{$years} {A vector of the corresponding years.}} + \item\code{$result} {An array with the same dimensions as the input array, + except for the temporal dimension which is renamed to 'year', moved to the + first dimension position and reduce to annual resolution.} + \item\code{$years} {A vector of the corresponding years.} +} } \description{ -This function computes the t90p, t10p, cdd or rx5day indices from n-dimensional arrays. +This function computes the t90p, t10p, cdd or rx5day indices from +n-dimensional arrays. } \examples{ ##Example synthetic data: @@ -57,8 +70,10 @@ time <- c(seq(ISOdate(1900, 1, 1), ISOdate(1900, 1, 31), "day"), seq(ISOdate(1909, 1, 1), ISOdate(1909, 1, 31), "day"), seq(ISOdate(1910, 1, 1), ISOdate(1910, 1, 31), "day"), seq(ISOdate(1911, 1, 1), ISOdate(1911, 1, 31), "day")) -metadata <- list(time = list(standard_name = 'time', long_name = 'time', calendar = 'gregorian', - units = 'days since 1970-01-01 00:00:00', prec = 'double', +metadata <- list(time = list(standard_name = 'time', long_name = 'time', + calendar = 'gregorian', + units = 'days since 1970-01-01 00:00:00', + prec = 'double', dim = list(list(name = 'time', unlim = FALSE)))) attr(time, "variables") <- metadata attr(data, 'Variables')$dat1$time <- time @@ -67,7 +82,6 @@ thres <- rep(10, 31 * 2 * 3) dim(thres) <- c(jdays = 31, lon = 2, lat = 3, model = 1) str(thres) - clim <- Climdex(data, metric = "t90p", threshold = thres) str(clim) } diff --git a/man/Extremes.Rd b/man/Extremes.Rd index 8042619d0ee0b4aaf0368b227cb3b80d3f73c552..34a4e7119f7152f9d38a14a1bb0054041d3d7277 100644 --- a/man/Extremes.Rd +++ b/man/Extremes.Rd @@ -20,7 +20,9 @@ Extremes( \arguments{ \item{data}{A n-dimensional array containing daily data.} -\item{threshold}{A n-dimensional array with the threshold to be/not to be reach, usually given by the a percentile computed with the \code{Threshold} function.} +\item{threshold}{A n-dimensional array with the threshold to be/not to be +reach, usually given by the a percentile computed with the \code{Threshold} +function.} \item{op}{The operator to use to compare data to threshold.} @@ -30,9 +32,12 @@ Extremes( \item{max.missing.days}{Maximum number of NA values per time period.} -\item{dates}{A vector of dates with a calendar attributes. If NULL (by default), the 'time' attributes of parameter 'data' are considered.} +\item{dates}{A vector of dates with a calendar attributes. If NULL +(by default), the 'time' attributes of parameter 'data' are considered.} -\item{timedim}{An integer number indicating the position of the time dimension in the parameter \code{data}. If NULL (by default), the dimension called 'time' in parameter \code{data}.} +\item{timedim}{An integer number indicating the position of the time dimension +in the parameter \code{data}. If NULL (by default), the dimension called +'time' in parameter \code{data}.} \item{calendar}{A character indicating the calendar type.} @@ -41,29 +46,42 @@ Extremes( \value{ A list of length 2: \itemize{ - \item\code{$output1}{An array with the same dimensions as the original \code{data}, except the time dimension which is reduced to annual resolution given a timeseries of maximum spell lengths for each year.} + \item\code{$output1}{An array with the same dimensions as the original + \code{data}, except the time dimension which is reduced to annual resolution + given a timeseries of maximum spell lengths for each year.} \item\code{$year}{A vector indicating the corresponding years.} } } \description{ -This function returns the number of spells of more than \code{min.length} days which exceed or are below the given \code{threshold} from daily data. +This function returns the number of spells of more than +\code{min.length} days which exceed or are below the given \code{threshold} +from daily data. } \details{ -This routine compares data to the thresholds using the given operator, generating a series of TRUE or FALSE values; these values are then filtered to remove any sequences of less than \code{min.length} days of TRUE values. It then computes the lengths of the remaining sequences of TRUE values (spells) and sums their lengths. The \code{spells.can.spa .years} option controls whether spells must always terminate at the end of a period, or whether they may continue until the criteria ceases to be met or the end of the data is reached. The default for fclimdex is FALSE. +This routine compares data to the thresholds using the given +operator, generating a series of TRUE or FALSE values; these values are then +filtered to remove any sequences of less than \code{min.length} days of TRUE +values. It then computes the lengths of the remaining sequences of TRUE values +(spells) and sums their lengths. The \code{spells.can.spa .years} option +controls whether spells must always terminate at the end of a period, or +whether they may continue until the criteria ceases to be met or the end of +the data is reached. The default for fclimdex is FALSE. } \examples{ ##Example synthetic data: data <- 1:(2 * 3 * 310 * 1) dim(data) <- c(time = 310, lon = 2, lat = 3, model = 1) time <- as.POSIXct(paste(sort(rep(1902:1911, 31)), 1, 1:31, sep = "-"), tz = "CET") -metadata <- list(time = list(standard_name = 'time', long_name = 'time', calendar = 'noleap', +metadata <- list(time = list(standard_name = 'time', long_name = 'time', + calendar = 'noleap', units = 'days since 1970-01-01 00:00:00', prec = 'double', dim = list(list(name = 'time', unlim = FALSE)))) attr(time, "variables") <- metadata attr(data, 'Variables')$dat1$time <- time threshold <- Threshold(data, dates = NULL, base.range = NULL, qtiles = 0.9, ncores = NULL) -res <- Extremes(data, threshold = threshold, op = ">", min.length = 6, spells.can.span.years = TRUE, - max.missing.days = 5, ncores = NULL) +res <- Extremes(data, threshold = threshold, op = ">", min.length = 6, + spells.can.span.years = TRUE, max.missing.days = 5, + ncores = NULL) str(res) } diff --git a/man/WaveDuration.Rd b/man/WaveDuration.Rd index a43c6b9252de6976d55a2db1096a5b971c375eb5..e639d550ed26c5eca3446473d8fa87510e9bcc2a 100644 --- a/man/WaveDuration.Rd +++ b/man/WaveDuration.Rd @@ -16,17 +16,26 @@ WaveDuration( ) } \arguments{ -\item{data}{A numeric n-dimensional array containing daily maximum or minimum temperature} +\item{data}{A numeric n-dimensional array containing daily maximum or minimum +temperature} -\item{threshold}{An array with the threshold to be/not to be reach, usually given by the 90th/10th percentiles for heat/cold waves computed with the \code{Threshold} function.} +\item{threshold}{An array with the threshold to be/not to be reach, usually +given by the 90th/10th percentiles for heat/cold waves computed with the +\code{Threshold} function.} -\item{op}{A character ">" (by default) or ">=" for heat waves and "<" or "<=" for cold waves indicating the operator must be used to compare data to threshold.} +\item{op}{A character ">" (by default) or ">=" for heat waves and "<" or "<=" +for cold waves indicating the operator must be used to compare data to +threshold.} -\item{spell.length}{A number indicating the number of consecutive days with extreme temperature to be considered heat or cold wave.} +\item{spell.length}{A number indicating the number of consecutive days with +extreme temperature to be considered heat or cold wave.} -\item{by.seasons}{If TRUE (by default), the wave duration is computed for each season (DJF/MAM/JJA/SON) separately. If FALSE is specified, the monthly wave duration is computed.} +\item{by.seasons}{If TRUE (by default), the wave duration is computed for each +season (DJF/MAM/JJA/SON) separately. If FALSE is specified, the monthly wave +duration is computed.} -\item{dates}{A vector of dates including calendar attributes. If NULL (by default), the 'time' attributes of parameter 'data' is used.} +\item{dates}{A vector of dates including calendar attributes. If NULL (by +default), the 'time' attributes of parameter 'data' is used.} \item{calendar}{A character indicating the calendar type.} @@ -35,24 +44,33 @@ WaveDuration( \value{ A list of length 2: \itemize{ - \item\code{$result}{An array with the same dimensions as the input \code{data}, but with the time dimension reduce from daily to monthly or seasonal resolution depending on the selected resolution in \code{by.season}.} - \item\code{$years}{A vector of the years and season/months corresponding to the resolution selected in \code{by.season} and temporal length of the input \code{data}}} + \item\code{$result}{An array with the same dimensions as the input + \code{data}, but with the time dimension reduce from daily to monthly or + seasonal resolution depending on the selected resolution in \code{by.season}.} + \item\code{$years}{A vector of the years and season/months corresponding to + the resolution selected in \code{by.season} and temporal length of the input + \code{data}}} } \description{ -This function computes the duration of a heat/cold wave as the number of consecutive days for which the maximum/minimum temperature is exceeding/below a threshold over a minimum number of days in month or seasonal resolution. +This function computes the duration of a heat/cold wave as the number of +consecutive days for which the maximum/minimum temperature is exceeding/below +a threshold over a minimum number of days in month or seasonal resolution. } \examples{ ##Example synthetic data: data <- 1:(2 * 3 * 31 * 5) dim(data) <- c(lon = 2, lat = 3, time = 31, model = 5) -time <- as.POSIXct(paste(paste(1900, 1, 1:31, sep = "-"), paste(12, 0, 0.0, sep = ":")), tz = "CET") -metadata <- list(time = list(standard_name = 'time', long_name = 'time', calendar = 'standard', +time <- as.POSIXct(paste(paste(1900, 1, 1:31, sep = "-"), paste(12, 0, 0.0, + sep = ":")), tz = "CET") +metadata <- list(time = list(standard_name = 'time', long_name = 'time', + calendar = 'standard', units = 'days since 1970-01-01 00:00:00', prec = 'double', dim = list(list(name ='time', unlim = FALSE)))) attr(time, "variables") <- metadata attr(data, 'Variables')$dat1$time <- time threshold <- rep(40, 31) -a <- WaveDuration(data, threshold, op = ">", spell.length = 6, by.seasons = TRUE, ncores = NULL) +a <- WaveDuration(data, threshold, op = ">", spell.length = 6, + by.seasons = TRUE, ncores = NULL) str(a) }