Nothing
#' @title Time Series Segmentation of Residual Trends (MAIN FUNCTION)
#'
#' @importFrom stats coef end frequency lm sd start time ts
#' @importFrom graphics abline arrows legend par plot
#' @importFrom utils tail read.csv
#' @importFrom broom glance
#'
#' @description
#' Time Series Segmented Residual Trend (TSS.RESTREND) methodology.Takes in a complete monthly
#' time series of a VI and its corrosponding precipitation (and temperature). It then looks looks
#' for breakpoints using the BFAST function. The significance of the breakpoin in the residuals
#' and the VPR is assessed using a Chow test, then, the total time series change is calculated.
#'
#' @author Arden Burrell, arden.burrell@unsw.edu.au
#'
#' @param CTSR.VI
#' Complete Monthly Time Series of Vegetation Index values.
#' An object of class \code{'ts'} object without NA's.
#' @param ACP.table
#' A table of every combination of offset period and accumulation period.for precipitation
#' ACP.table can be calculated using the \code{\link{climate.accumulator}}.
#' @note if ACP.table = FALSE, CTSR.RF and acu.RF must be provided as well as
#' rf.b4 and rf.af for \code{'ts'} with a breakpoint in the VPR.
#' @param ACT.table
#' A table of every combination of offset period and accumulation period.for temperature
#' ACP.table can be calculated using the \code{\link{climate.accumulator}}.
#' @param CTSR.RF
#' Complete Time Series of Rainfall. An object of class 'ts' object without NA's
#' and be the same length and cover the same time range as CTSR.VI.
#' If ACP.table is provided, CTSR.RF will be automitaclly calculated using the
#' \code{\link{ACP.calculator}}
#' @param CTSR.TM
#' Complete Time Series of temperature. An object of class 'ts' object without NA's
#' and be the same length and cover the same time range as CTSR.VI. Default (CTSR.TM=NULL).
#' If ACT.table is provided, CTSR.RF will be automitaclly calculated using the
#' \code{\link{ACP.calculator}}
#' @param anu.VI
#' The annual (Growing season) max VI. Must be a object of class \code{'ts'} without NA's.
#' if anu.VI=FALSE, it will be calculated from the CTSR.VI using \code{\link{AnMaxVI}}.
#' @param acu.RF
#' The optimal accumulated rainfall for anu.VI. Must be a object of class \code{'ts'} without
#' NA's and be of equal length and temporal range to anu.VI. if anu.RF=FALSE, it will be
#' calculated from ACP.table usingthe \code{\link{AnnualClim.Cal}}
#' @param acu.TM
#' The optimal accumulated rainfall for anu.TM. Must be a object of class \code{'ts'} without
#' NA's and be of equal length and temporal range to anu.TM. if anu.TM=FALSE, it will be
#' calculated from ACT.table usingthe \code{\link{AnnualClim.Cal}}
#' @param VI.index
#' the index of the CTSR.VI ts that the anu.VI values occur at. Must be the same length
#' as anu.VI. NOTE. R indexs from 1 rather than 0.
#' if VI.index=FALSE, it will be calculated from the CTSR.VI using \code{\link{AnMaxVI}}.
#' @param rf.b4
#' If a breakpoint in the VPR is detected this is the optimial accumulated rainfall before
#' the breakpoint. must be the same length as the anu.VI. If ACP.table is provided it will
#' be generated using \code{\link{AnnualClim.Cal}}
#' @param rf.af
#' If a breakpoint in the VPR is detected this is the optimial accumulated rainfall after
#' the breakpoint. must be the same length as the anu.VI. If ACP.table is provided it will
#' be generated using \code{\link{AnnualClim.Cal}}
#' @param sig
#' Significance of all the functions. defualt sig=0.05
#' @param season
#' See \code{\link[bfast]{bfast}}. This season value only applies to bfast done using the CTS
#' VPR. if a non VPR adjusted BFAST is performed.a harmonic season is used.
#' @param exclude
#' A numberic vector containg months excluded from breakpoint detection. This was included to
#' allow sensor transitions to be masked.
#' @param allow.negative
#' If true, will not preference positive slope in either CTSR or VI calculations. default=FALSE is set
#' because negative associations between rainfall and vegetation in water limited ecosystems is unexpected
#' If temperature data is included then this paramter is forced to TRUE.
#' @param allowneg.retest default=FALSE
#' If temperature data is provided but found to not be significant then a retest is performed.
#' This paramter is to allow negative on re-test.
#' @param h
#' See \code{\link[bfast]{bfast}}, The.minimal segment size between potentially detected breaks in the trend model
#' given as fraction relative to the sample size (i.e. the minimal number of observations in each segment
#' divided by the total length of the timeseries. Default h = 0.15.
#' @param retnonsig
#' Bool. New in v0.3.0. Allows TSSRESTREND to return change estimates of values that filed the sig component in the residual analysis.
#' defualt FALSE will give the same result as eralier versions.
#' @return
#' An object of class \code{'TSSRESTREND'} is a list with the following elements:
#' \describe{
#' \item{\bold{\emph{summary}}}{
#' \describe{
#' \item{Method}{The method used to determine total change. (\emph{RESTREND} see \code{\link{RESTREND}},
#' \emph{segmented.RESTREND} see \code{\link{seg.RESTREND}}, \emph{segmented.VPR} see
#' \code{\link{seg.VPR}})}
#' \item{Total.Change}{The total significant change. Residual.Change + VPR.HeightChange. }
#' \item{Residual.Change}{The change in the VPR Residuals over the time period}
#' \item{VPR.HeightChange}{The change in VI at mean rainfall for a "ts" with a significant
#' breakpoint in the VPR}
#' \item{model.p}{p value of the regression model fitted to the VPR. See \code{\link[stats]{lm}}}
#' \item{residual.p}{p value of the regression model fitted to the VPR Residuals. See \code{\link[stats]{lm}}}
#' \item{VPRbreak.p}{the p value associated with the break height. See \code{\link[stats]{lm}}}
#' \item{bp.year}{The Year of the most significant breakpoint}
#' }}
#' \item{\bold{\emph{ts.data}}}{The Time series used in analysis. See Arguments for description
#' \itemize{
#' \item CTSR.VI
#' \item CTSR.RF
#' \item anu.VI
#' \item VI.index
#' \item acu.RF
#' \item StdVar.RF see \code{\link{seg.VPR}})}
#' }
#' \item{\bold{\emph{ols.summary}}}{
#' \describe{
#' \item{chow.summary}{summary of the most significant breakpoint. }
#' \item{chow.ind}{Summary of every detected breakpoint}
#' \item{OLS.table}{A matrix containing the coefficents for the CTS.fit, VPR.fit, RESTREND.fit and segVPR.fit}
#' }}
#' \item{\bold{\emph{TSSRmodels}}}{
#' models of class "lm" \code{\link[stats]{lm}} and class "bfast" \code{\link[bfast]{bfast}} generated.}
#' }
#'
#' @seealso
#' \itemize{
#' \item \code{\link{plot.TSSRESTREND}}
#' \item \code{\link{print.TSSRESTREND}}
#' }
#' @export
#' @examples
#' \dontrun{
#' #To get the latest version of the package (Still in development)
#' install.packages("devtools")
#' library("devtools")
#' install_github("ArdenB/TSSRESTREND", subdir="TSS.RESTREND")
#' library(TSS.RESTREND)
#' #Find the path of the rabbitRF.csv dataset, read it in and turn it into a time series
#' rf.path<- system.file("extdata", "rabbitRF.csv", package = "TSS.RESTREND", mustWork = TRUE)
#' in.RF <- read.csv(rf.path)
#' rf.data <- ts(in.RF, end=c(2013,12), frequency = 12)
#'
#' #Find the path of the rabbitVI.csv dataset and read it in
#' vi.path <- system.file("extdata", "rabbitVI.csv", package = "TSS.RESTREND", mustWork = TRUE)
#' in.VI <- read.csv(vi.path)
#' CTSR.VI <- ts(in.VI, start=c(1982, 1), end=c(2013,12), frequency = 12)
#'
#' #Define the max accumuulation period
#' max.acp <- 12
#' #Define the max offset period
#' max.osp <- 4
#' #Create a table of every possible precipitation value given the max.acp and max.osp
#' ACP.table <- climate.accumulator(CTSR.VI, rf.data, max.acp, max.osp)
#' results <- TSSRESTREND(CTSR.VI, ACP.table)
#' print(results)
#' plot(results, verbose=TRUE)
#' }
#'
TSSRESTREND <- function(
CTSR.VI, ACP.table = FALSE, ACT.table = NULL, CTSR.RF = FALSE, CTSR.TM = NULL,
anu.VI = FALSE, acu.RF = FALSE, acu.TM = NULL, VI.index = FALSE, rf.b4 = FALSE,
rf.af = FALSE, sig = 0.05, season = "none", exclude = 0, allow.negative = FALSE,
allowneg.retest = FALSE, h = 0.15, retnonsig=FALSE){
# ==============================================================================================
# ========== Sanity check the input data ==========
# Description:
# Each check liiks at a different paramter. If the data fails
# the check will stop, else, it breaks after all the checks
while (TRUE) { #Test the variables for consistenty
# check the class, type and tempora range of provided data
if (class(CTSR.VI) != "ts")
stop("CTSR.VI Not a time series object. Please check the data")
if ((class(ACP.table) == "logical") && (!CTSR.RF || acu.RF))
stop("Insufficent Rainfall data provided. Provide either a complete ACP.table or both the CTSR.RF & acu.RF")
if ((!anu.VI) || (!VI.index)) {
# Get the annual Max VI values
max.df <- AnMaxVI(CTSR.VI)
# Pull the key components from the result
anu.VI <- max.df$Max #the VI values
VI.index <- max.df$index #the indes values
Max.Month <- max.df$Max.Month #month if the year the even occured
}else{
if (class(anu.VI) != "ts") {
stop("anu.VI Not a time series object")
}
}
# Change the allow.negative for multivariate regression with temperature
# (Temperature can have a negative or positive impact on veg in drylands)
# if (!is.null(ACT.table)) {
# allow.negative = TRUE
# }
if (!CTSR.RF) {
# ==============================================================================================
# ===== Calculate the Complete time seties Accumulation using ACP.calculator =====
CTS.Str <- ACP.calculator(
CTSR.VI, ACP.table, ACT.table, allow.negative = allow.negative,
allowneg.retest = allowneg.retest
)
# ==============================================================================================
# ===== Determine if BFAST is applied to the CTS residuals or the raw VI time series =====
# If the allow negative is on, this is ignored, else perform a negative slope check and perform a significance check
# This mod is will impact results comparisons before V0.1.04
if ((!allow.negative && as.numeric(CTS.Str$summary)[1] < 0) || as.numeric(CTS.Str$summary)[4] > sig) {
BFraw = TRUE
} else {BFraw = FALSE}
# +++++ Pull out the relevant paramters for use from the CTS accumulation result +++++
CTSR.RF <- CTS.Str$CTSR.precip #RF values
CTSR.TMraw <- CTS.Str$CTSR.rawtemp
# Check for the presence of temperature data
if (is.null(CTSR.TM)) {
CTSR.TM <- CTS.Str$CTSR.tmp #Temperature values, null if temp is not considered
}
details.CTS.VPR <- CTS.Str$summary #Summay of the FIT between the CTS.VRP
CTSR.osp <- CTS.Str$CTSR.osp #CTS Off set period
CTSR.acp <- CTS.Str$CTSR.acp #CTS Sccumulation period
CTSR.tosp <- CTS.Str$CTSR.tosp #CTS Off set period
CTSR.tacp <- CTS.Str$CTSR.tacp #CTS Sccumulation period
}else{
# ===== Check the times of the datasets =====
if (class(CTSR.RF) != "ts") {
stop("CTSR.RF Not a time series object")
}
# get the time data out
start.ti <- time(CTSR.VI)
freq <- frequency(CTSR.VI)
# check the two ts object cover the same time period
start.ti2 <- time(CTSR.RF)
freq2 <- frequency(CTSR.RF)
#Check the start dates and the frequency are correct
if (!identical(start.ti, start.ti2)) {
stop("ts objects do not have the same time, (CTSR.VI & CTSR.RF)")}
if (!identical(freq, freq2)) {
stop("ts objects do not have the same frequency, (CTSR.VI & CTSR.RF)")}
}
# =================================================================================================================
# ===== Calculate the optimal Accumulated Rainfall and temperature for the annual max VI using AnnualClim.Cal =====
if (!acu.RF) { #if annual accumulated precipitation in not provided
precip.df <- AnnualClim.Cal(anu.VI, VI.index, ACP.table, ACT.table = ACT.table, allow.negative = allow.negative)
# +++++ Pull out and store key values from AnnualClim.Cal result ++++++
osp <- precip.df$osp # offset period
acp <- precip.df$acp # Accumulation period
tosp <- precip.df$tosp # offset period
tacp <- precip.df$tacp # Accumulation period
acu.RF <- precip.df$annual.precip # precip values
acu.TM <- precip.df$annual.temp # precip values
details.VPR <- precip.df$summary # The summary of the lm between rainfall and Vegetation
} else {# If the anmax vi is passed
# Check the passed accumulated rainfall
if (class(acu.RF) != "ts")
stop("acu.RF Not a time series object")
# get the time peramaters
st.ti <- time(anu.VI)
st.f <- frequency(anu.VI)
st.ti2 <- time(acu.RF)
st.f2 <- frequency(acu.RF)
#check the two ts object cover the same time period and frequency
if (!identical(st.ti, st.ti2))
stop("ts object do not have the same time, (acu.RF & anu.VI)")
if (!identical(st.f, st.f2))
stop("ts object do not have the same frequency, (acu.RF & anu.VI)")
}
# ===== Check passed breakpoint rainfall data =====
# if breakpoints are used defined rather than determined, test all the data will work
if (class(rf.b4) != "logical") {
if (length(rf.b4) != (length(rf.af)))
stop("rf.b4 and rf.af are different shapes. They must be the same size and be th same lenths as acu.VI")
}
break
}
# ==============================================================================================
# ===== Perform BFAST to look for potential breakpoints using VPR.BFAST =====
# Pass the infomation about the VI and RF as well as the BFAST method to the VPR.BFAST script
bkp = VPR.BFAST(CTSR.VI, CTSR.RF, CTSR.TM=CTSR.TM, season = season, BFAST.raw = BFraw, h = h)
# Extract the key values from the BFAST result
bp <- bkp$bkps
BFAST.obj <- bkp$BFAST.obj #For the models Bin
CTS.lm <- bkp$CTS.lm #For the Models Bin
bp <- bp[!bp %in% exclude] #remove breakpoints in the exclude list (Sensor transitions)
BFT <- bkp$BFAST.type #Type of BFAST used
# +++++ put all the infomation on the offset periods and accumulation period into a dataframe +++++
acum.df <- data.frame(
CTSR.osp = CTSR.osp, CTSR.acp = CTSR.acp, CTSR.tosp = CTSR.tosp, CTSR.tacp = CTSR.tacp,
osp = osp, acp = acp, tosp = tosp, tacp = tacp, osp.b4 = NaN, acp.b4 = NaN, tosp.b4 = NaN,
tacp.b4 = NaN, osp.af = NaN, acp.af = NaN, tosp.af = NaN, tacp.af = NaN
)
# ===== Check and see if there are breakpoint that need to be tested =====
if (class(bp) == "logical" | length(bp) == 0) {#Should catch both the false and the no breakpoints
# no breakpoints detected by the BFAST
bp <- FALSE
test.Method = "RESTREND" # MEthod set to determine further testing
# Chow summary populated with false
chow.sum <- data.frame(abs.index = FALSE, yr.index = FALSE, reg.sig = FALSE, VPR.bpsig = FALSE)
chow.bpi <- FALSE
}else{# Breakpoints detected by the BFAST
# ===== Perform the chow test on the breakpoints using CHOW function =====
bp <- as.numeric(bkp$bkps)
res.chow <- CHOW(anu.VI, acu.RF, VI.index, bp, acu.TM = acu.TM, sig = sig)
# Pull out the key values from the CHOW
brkp <- as.integer(res.chow$bp.summary["yr.index"]) #this isn't right
chow.sum <- res.chow$bp.summary
chow.bpi <- res.chow$allbp.index
# Use the CHOW results to set the testmethod
test.Method = res.chow$n.Method
}
# ==============================================================================================
# ========== Perform a total change calculation ==========
# Note:
# The method is calculated by the CHOW function
if (test.Method == "RESTREND") {
# ===== No breakpoints, Results calculated using the RESTREND function =====
result <- RESTREND(anu.VI, acu.RF, VI.index, acu.TM=acu.TM, sig = sig, retnonsig=retnonsig)
}else if (test.Method == "seg.RESTREND") {
# ===== breakpoints in the VPR/VCR residuals, Results calculated using the seg.RESTREND function =====
breakpoint = as.integer(res.chow$bp.summary[2])
result <- seg.RESTREND(anu.VI, acu.RF, VI.index, brkp, acu.TM=acu.TM, sig=sig, retnonsig=retnonsig)
}else if (test.Method == "seg.VPR") {
# ===== breakpoints in the VPR/VCR, Results calculated using the seg.VPR function =====
if ((!rf.b4) || (!rf.af)) {
# +++++ Calculate the regression coefficents on either side of the breakpoint using AnnualClim.Cal +++++
VPRbp.df <- AnnualClim.Cal(anu.VI, VI.index, ACP.table, ACT.table=ACT.table, Breakpoint = brkp, allow.negative = allow.negative)
rf.b4 <- VPRbp.df$rf.b4
rf.af <- VPRbp.df$rf.af
tm.b4 <- VPRbp.df$tm.b4
tm.af <- VPRbp.df$tm.af
# Check if temp is insignificant either side of the breakpoint in the VPR,
# if yes, remove temp from segmented VPR
if (is.null(tm.b4) && is.null(tm.af)) {acu.TM = NULL}
#Add the segmented offset periods and accumulation periods to the existing dataframe
acum.df$osp.b4 <- VPRbp.df$osp.b4
acum.df$acp.b4 <- VPRbp.df$acp.b4
acum.df$tosp.b4 <- VPRbp.df$tosp.b4
acum.df$tacp.b4 <- VPRbp.df$tacp.b4
acum.df$osp.af <- VPRbp.df$osp.af
acum.df$acp.af <- VPRbp.df$acp.af
acum.df$tosp.af <- VPRbp.df$tosp.af
acum.df$tacp.af <- VPRbp.df$tacp.af
}
# +++++ Perform segmented VPR/VCR calculation +++++
breakpoint = as.integer(res.chow$bp.summary[2])
print(brkp)
result <- seg.VPR(anu.VI, acu.RF, VI.index, brkp, rf.b4, rf.af, acu.TM, tm.b4, tm.af, sig=sig, retnonsig=retnonsig)
}
# ========== New (in version 0.3.0) Sanity check on Total change values ==========
# +++ Checks to see if value for total change fall within a sane range +++
if (abs(result$summary$Total.Change > (max(CTSR.VI)-min(CTSR.VI)))){
print("Non Valid estimate produced, returning zero")
result$summary$Total.Change = 0
result$summary$Method = "InvalidValueError"
}
# else if (result$summary$Total.Change == 0){
# browser("Failure here somewhere, Take a look and see what the options are")
# # result2 <- RESTREND(anu.VI, acu.RF, VI.index, acu.TM=acu.TM, sig = sig, retnonsig=retnonsig)
# } else if (is.na(result$summary$Total.Change)){
# browser("Failure here somewhere, Take a look and see what the options are")
# }
# ==============================================================================================
# ===== Build the results into a list to be returned to user =====
# +++++ add the common variable to the results list ++++++
# the fitted models
result$TSSRmodels$CTS.fit <- CTS.lm
result$TSSRmodels$BFAST <- BFAST.obj
# Complete Time series values
result$ts.data$CTSR.VI <- CTSR.VI
result$ts.data$CTSR.RF <- CTSR.RF
if (!is.null(ACT.table)) {# Add Temperature if present
result$ts.data$CTSR.TMraw <- ts(ACT.table[1, ], start = c(start(CTSR.VI)[1], start(CTSR.VI)[2]), frequency = 12)
result$ts.data$CTSR.TM <- CTSR.TM
}else{
result$ts.data$CTSR.TM <- CTSR.TM
result$ts.data$CTSR.TMraw <- CTSR.TM
}
# add to the ols summary table
result$ols.summary$chow.sum <- chow.sum
result$ols.summary$chow.ind <- chow.bpi
result$ols.summary$OLS.table["CTS.fit",] <- details.CTS.VPR
# Add the accumulation and offset periods
result$acum.df <- acum.df
#add the bfast method to the results summary
result$summary$BFAST.Method <- BFT
#return the results
return(result)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.