Newer
Older
#-----------------------------------------------------------
# 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.
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
#-----------------------------------------------------------
# Package name: climdex.pcic
# Version: 1.1-11
# Date: 2020-01-21
# Title: PCIC Implementation of Climdex Routines
# Author: David Bronaugh <bronaugh@uvic.ca> for the Pacific Climate Impacts
# Consortium
# Maintainer: James Hiebert <hiebert@uvic.ca>
# 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 <http://etccdi.pacificclimate.org/list_27_indices.shtml>
# 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)