Title: Trend Estimating Tools
Description: The traditional linear regression trend, Modified Mann-Kendall (MK) non-parameter trend and bootstrap trend are included in this package. Linear regression trend is rewritten by '.lm.fit'. MK trend is rewritten by 'Rcpp'. Finally, those functions are about 10 times faster than previous version in R. Reference: Hamed, K. H., & Rao, A. R. (1998). A modified Mann-Kendall trend test for autocorrelated data. Journal of hydrology, 204(1-4), 182-196. <doi:10.1016/S0022-1694(97)00125-X>.
Version: 0.1.5
License: MIT + file LICENSE
Encoding: UTF-8
RoxygenNote: 7.2.3
LinkingTo: Rcpp, RcppArmadillo
Imports: Rcpp, fftwtools, boot, magrittr, matrixStats, lubridate, terra, plyr
Suggests: covr, knitr, rmarkdown, testthat (≥ 3.0.0)
URL: https://github.com/rpkgs/rtrend
BugReports: https://github.com/rpkgs/rtrend/issues
Config/testthat/edition: 3
VignetteBuilder: knitr
NeedsCompilation: yes
Packaged: 2024-01-11 01:58:18 UTC; kongdd
Author: Dongdong Kong ORCID iD [aut, cre], Heyang Song ORCID iD [aut]
Maintainer: Dongdong Kong <kongdd.sysu@gmail.com>
Repository: CRAN
Date/Publication: 2024-01-11 03:30:02 UTC

rtrend: Trend Estimating Tools

Description

The traditional linear regression trend, Modified Mann-Kendall (MK) non-parameter trend and bootstrap trend are included in this package. Linear regression trend is rewritten by '.lm.fit'. MK trend is rewritten by 'Rcpp'. Finally, those functions are about 10 times faster than previous version in R. Reference: Hamed, K. H., & Rao, A. R. (1998). A modified Mann-Kendall trend test for autocorrelated data. Journal of hydrology, 204(1-4), 182-196. doi:10.1016/S0022-1694(97)00125-X.

Author(s)

Maintainer: Dongdong Kong kongdd.sysu@gmail.com (ORCID)

Authors:

See Also

Useful links:

faster autocorrelation based on ffw

Description

This function is 4-times faster than stats::acf()

Usage

acf.fft(x, lag.max = NULL)

Arguments

lag.max

maximum lag at which to calculate the acf. Default is 10\log_{10}(N/m) where N is the number of observations and m the number of series. Will be automatically limited to one less than the number of observations in the series.

Value

An array with the same dimensions as x containing the estimated autocorrelation.

References

  1. https://github.com/santiagobarreda/phonTools/blob/main/R/fastacf.R

  2. https://gist.github.com/FHedin/05d4d6d74e67922dfad88038b04f621c

  3. https://gist.github.com/ajkluber/f293eefba2f946f47bfa

  4. http://www.tibonihoo.net/literate_musing/autocorrelations.html#wikispecd

  5. https://lingpipe-blog.com/2012/06/08/autocorrelation-fft-kiss-eigen

Examples

set.seed(1)
x <- rnorm(100)
r_fast <- acf.fft(x)
r <- acf(x, plot=FALSE, lag.max=100)$acf[,,1]

apply function for 3d array

Description

NA values will be removed automatically

Usage

apply_3d(
  array,
  dim = 3,
  FUN = rowMeans2,
  by = NULL,
  scale = 1,
  na.rm = TRUE,
  ...
)

Arguments

array

A 3d array

dim

giving the subscripts to split up data by.

FUN

function, should only be row applied function, e.g. matrixStats::rowMeans2, matrixStats::rowMins, matrixStats::rowRanges. Because 3d array will be convert to matrix first, with the aggregated dim in the last dimension.

by

  • If not provided (NULL), the aggregated dim will be disappear. For example, daily precipitation ⁠[nrow, ncol, 31-days]⁠ aggregate into monthly ⁠[nrow, ncol]⁠.

  • If provided, by should be equal to the aggregated dim. For example, daily precipitation ⁠[nrow, ncol, 365-days]⁠ aggregate into monthly ⁠[nrow, ncol, 12-months]⁠. In that situation, by should be equal to 365, and be format(date, '%Y%m').

scale

in the same length of by, or a const value, value_returned = FUN(x)*scale. This parameter is designed for converting monthly to yearly, meanwhile multiply days in month. Currently, same group should have the same scale factor. Otherwise, only the first is used.

See Also

apply_row matrixStats::rowRanges

Examples

set.seed(1)
size <- c(10, 8, 31)
arr <- array(rnorm(10 * 8 * 31), dim = size)

by <- c(rep(1, 10), rep(2, 21))
r2 <- apply_3d(arr, 3, by = by, FUN = rowMeans)

## Not run: 
arr_yearly <- apply_3d(arr, by = year(dates), scale = days_in_month(dates))

## End(Not run)

apply_col

Description

Usage

apply_col(mat, by, FUN = colMeans2, scale = 1, ...)

apply_row(mat, by, FUN = rowMeans2, scale = 1, ...)

Arguments

mat

matrix, ⁠[nrow, ncol]⁠

by

integer vector, with the dim of ⁠[ntime]⁠

scale

in the same length of by, or a const value, value_returned = FUN(x)*scale. This parameter is designed for converting monthly to yearly, meanwhile multiply days in month. Currently, same group should have the same scale factor. Otherwise, only the first is used.

Details

For example, setting the dimension of mat is ⁠[ngrid, ntime]⁠, if you want to aggregate by time, apply_row should be used here; if you want to aggregate by region (grids), apply_col should be used.

Note

This function also suits for big.matrix object.

Examples

mat <- matrix(rnorm(4 * 6), 4, 6)
mat_bycol <- apply_col(mat, c(1, 1, 2, 2), colMeans)
mat_byrow <- apply_row(mat, c(1, 1, 2, 2, 3, 3), rowMeans)

array_3dTo2d

Description

array_3dTo2d

Usage

array_3dTo2d(array, I_grid = NULL)

array_2dTo3d(array, I_grid = NULL, dim)

Arguments

array

array with the dimension of ⁠[nlon, nlat, ntime]⁠

I_grid

subindex of ⁠[nrow, ncol]⁠

dim

⁠[nrow, ncol]⁠

Value

⁠[nlat*nlon, ntime]⁠

chunk

Description

chunk

Usage

chunk(x, nchunk = 6)

Arguments

x

a vector or list

nchunk

the number of chunks to be splitted

References

https://stackoverflow.com/questions/3318333/split-a-vector-into-chunks-in-r

parallel apply and llply

Description

parallel apply and llply

Usage

llply_par(X, FUN, ..., byrow = TRUE, .combine = c)

parLapply2(X, FUN, ..., byrow = TRUE, .combine = c)

apply_par(X, .margins = 1, FUN, ..., .progress = "text")

Modified Mann Kendall

Description

If valid observations <= 5, NA will be returned.

Usage

mkTrend_r(y, ci = 0.95, IsPlot = FALSE)

mkTrend(y, x = seq_along(y), ci = 0.95, IsPlot = FALSE)

Arguments

y

numeric vector

ci

critical value of autocorrelation

IsPlot

boolean

x

(optional) numeric vector

Details

mkTrend is 4-fold faster with .lm.fit.

Value

Note

slp is significant, if pval < alpha.

Author(s)

Dongdong Kong

References

Hipel, K.W. and McLeod, A.I. (1994), Time Series Modelling of Water Resources and Environmental Systems. New York: Elsevier Science.

Libiseller, C. and Grimvall, A., (2002), Performance of partial Mann-Kendall tests for trend detection in the presence of covariates. Environmetrics 13, 71–84, doi:10.1002/env.507.

See Also

fume::mktrend and trend::mk.test

Examples

x <- c(4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69)
r <- mkTrend(x)
r_cpp <- mkTrend(x, IsPlot = TRUE)

movmean

Description

NA and Inf values in the y will be ignored automatically.

Usage

movmean(y, halfwin = 1L, SG_style = FALSE, w = NULL)

movmean2(y, win_left = 1L, win_right = 0L, w = NULL)

movmean_2d(mat, win_left = 3L, win_right = 0L)

Arguments

y

A numeric vector.

halfwin

Integer, half of moving window size

SG_style

If true, head and tail values will be in the style of SG (more weights on the center point), else traditional moving mean style.

w

Corresponding weights of y, with the same length.

win_left, win_right

windows size in the left and right

mat

numeric matrix

Examples

x <- 1:100
x[50] <- NA; x[80] <- Inf
s1 <- movmean(x, 2, SG_style = TRUE)
s2 <- movmean(x, 2, SG_style = FALSE)
movmean2(c(4, 8, 6, -1, -2, -3, -1), 2, 0)
movmean2(c(4, 8, NA, -1, -2, Inf, -1), 2, 0)

Set dimensions of an Object

Description

Set dimensions of an Object

Usage

set_dim(x, dim)

set_dimnames(x, value)

Arguments

x

an R object, for example a matrix, array or data frame.

dim

integer vector, see also base::dim()

value

For the default method, either NULL or a numeric vector, which is coerced to integer (by truncation).

See Also

base::dim

Examples

x <- 1:12
set_dim(x, c(3, 4))

slope_arr

Description

slope_arr

Usage

slope_arr(
  arr,
  fun = rtrend::slope_mk,
  return.list = FALSE,
  .progress = "text",
  ...
)

Value

t, A 3d array, with the dim of ⁠[nx, ny, 2]⁠.

calculate slope of rast object

Description

calculate slope of rast object

Usage

slope_rast(
  r,
  period = c(2001, 2020),
  outfile = NULL,
  fun = rtrend::slope_mk,
  ...,
  overwrite = FALSE,
  .progress = "text"
)

rast_filter_time(r, period = c(2001, 2020))

Arguments

r

A yearly rast object, which should have time attribute

period

c(year_begin, year_end)

outfile

The path of outputed tiff file. If specified, slope and pvalue will be written into outfile.

fun

the function used to calculate slope, see slope() for details.

...

other parameters ignored

overwrite

logical. If TRUE, outfile is overwritten.

.progress

name of the progress bar to use, see create_progress_bar

Value

A terra rast object, with bands of slope and pvalue.

See Also

terra::rast()

Examples

library(rtrend)
library(terra)

f <- system.file("rast/MOD15A2_LAI_China_G050_2001-2020.tif", package = "rtrend")
r <- rast(f)
r
time(r)

slp <- slope_rast(r,
  period = c(2001, 2020), 
  outfile = "LAI_trend.tif", overwrite = TRUE,
  fun = rtrend::slope_mk, .progress = "none"
)
# if you want to show progress, set `.progress = "text"`
slp
plot(slp)

file.remove("LAI_trend.tif")

slope

Description

Usage

slope_sen(y, x = NULL)

slope(y, x, ...)

slope_p(y, x, fast = TRUE)

slope_sen_r(y, x = seq_along(y), ...)

slope_mk(y, x = NULL, ...)

slope_boot(y, x = NULL, slope_FUN = slope, times = 100, alpha = 0.1, seed, ...)

Arguments

y

vector of observations of length n, or a matrix with n rows.

x

vector of predictor of length n, or a matrix with n rows.

...

ignored.

fast

Boolean. If true, stats::.lm.fit() will be used, which is 10x faster than stats::lm().

slope_FUN

one of slope(), slope_p(), slope_mk()

times

The number of bootstrap replicates.

alpha

significant level, defalt 0.1

seed

a single value, interpreted as an integer, or NULL (see ‘Details’).

Value

For slope_boot, slope is estimated in many times. The lower, mean, upper and standard deviation (sd) are returned.

Examples

y <- c(4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69)
r <- slope(y)
r_p <- slope_p(y)
r_mk <- slope_mk(y)
r_boot <- slope_boot(y)

Weighted Savitzky-Golay

Description

NA and Inf values in the y has been ignored automatically.

Usage

smooth_wSG(y, halfwin = 1L, d = 1L, w = NULL)

smooth_SG(y, halfwin = 1L, d = 1L)

Arguments

y

colvec

halfwin

halfwin of Savitzky-Golay

d

polynomial of degree. When d = 1, it becomes moving average.

w

colvec of weight

Examples

y <- c(1, 3, 2, 5, 6, 8, 10, 1)
w <- seq_along(y)/length(y)

halfwin = 2
d = 2
s1 <- smooth_wSG(y, halfwin, d, w)
s2 <- smooth_SG(y, halfwin, d)

split_data

Description

split_data

Usage

split_data(x, nchunk = 6, byrow = TRUE)

Arguments

x

a vector, list, or matrix (not support 3d array)

nchunk

the number of chunks to be splitted

byrow

If TRUE, split by row, otherwise by column.

summary_lm

Description

summary method for class ".lm.fit".. It's 200 times faster than traditional lm.

Usage

summary_lm(obj, ...)

Arguments

obj

Object returned by .lm.fit.

...

ignored

Value

a p x 4 matrix with columns for the estimated coefficient, its standard error, t-statistic and corresponding (two-sided) p-value. Aliased coefficients are omitted.

Examples

set.seed(129)
n <- 100
p <- 2
X <- matrix(rnorm(n * p), n, p) # no intercept!
y <- rnorm(n)

obj <- .lm.fit (x = cbind(1, X), y = y) 
info <- summary_lm(obj)