R package: Extract Remote Sensing Vegetation Phenology by TIMESAT Fortran library.
Note that TSF_process.exe and TSF_fit2time should be added to system environment path.
You can install the released version of rTIMESAT from GitHub with:
# install.packages("rTIMESAT")
devtools::install_github("kongdd/rTIMESAT")This is a basic example which shows you how to use rTIMESAT:
TSF_main <- function(d, nptperyear = 23, cache = T){
## 1. Prepare inputs
sitename <- d$site[1]
nyear <- floor(nrow(d)/nptperyear)
npt <- nyear * nptperyear
d <- d[1:npt, ]
# pls make sure it's complete year in input
file_y <- sprintf("TSM_%s_y.txt", sitename)
file_w <- sprintf("TSM_%s_w.txt", sitename)
file_set <- sprintf("TSM_%s.set", sitename)
write_input(d$y , file_y, nptperyear)
write_input(d$SummaryQA %>% as.numeric(), file_w, nptperyear)
## 2. Update options
options <- list(
job_name = "",
file_y = file_y, # Data file list/name
file_w = file_w, # Mask file list/name
nyear_and_nptperear = c(nyear, nptperyear), # No. years and no. points per year
ylu = c(0, 9999), # Valid data range (lower upper)
qc_1 = c(0, 0, 1), # Quality range 1 and weight
qc_2 = c(1, 1, 0.5), # Quality range 2 and weight
qc_3 = c(2, 3, 0.2), # Quality range 3 and weight
A = 0.1, # Amplitude cutoff value
output_type = c(1, 1, 0), # Output files (1/0 1/0 1/0), 1: seasonality data; 2: smoothed time-series; 3: original time-series
seasonpar = 1.0, # Seasonality parameter (0-1)
iters = 2, # No. of envelope iterations (3/2/1)
FUN = 2, # Fitting method (1/2/3): (SG/AG/DL)
half_win = 7, # half Window size for Sav-Gol.
meth_pheno = 1, # (1: seasonal amplitude, 2: absolute value, 3: relative amplitude, 4: STL trend)
trs = c(0.5, 0.5) # Season start / end values
)
options$job_name <- sitename
opt <- update_setting(options)
write_setting(opt, file_set)
TSF_process(file_set) # call TSF_process.exe
file_tts <- sprintf("%s_fit.tts", opt$job_name)
file_tpa <- sprintf("%s_TS.tpa", opt$job_name)
# note: only suit for ascii
tidy_tts <- function(d_tts){
sites <- d_tts$row %>% paste0("v", .)
npt <- ncol(d_tts) - 2
d <- d_tts %>% {.[, 3:ncol(.)]} %>% as.matrix() %>% t() %>% data.frame() %>%
set_colnames(sites) %>% cbind(t = 1:npt, .)
d
}
d_tts <- read_tts(file_tts) %>% tidy_tts()
d_tpa <- read_tpa(file_tpa)
if (!cache){
status1 <- file.remove(c(file_tts, file_tpa, file_y, file_w, file_set))
status2 <- dir(".", "*.ndx", full.names = T) %>% file.remove()
}
list(fit = d_tts, pheno = d_tpa)
}library(rTIMESAT)
library(magrittr)
data("MOD13A1")
sitename <- "US-KS2"
# sitename <- "CA-NS6"
df <- subset(MOD13A1$dt, date >= as.Date("2004-01-01") & date <= as.Date("2010-12-31"))
d <- subset(df, site == sitename)
d$y <- d$EVI/1e4
r <- TSF_main(d, cache = F)
print(str(r))
#> List of 2
#> $ fit :'data.frame': 161 obs. of 2 variables:
#> ..$ t : int [1:161] 1 2 3 4 5 6 7 8 9 10 ...
#> ..$ v1: num [1:161] 0.371 0.362 0.359 0.359 0.362 ...
#> $ pheno:'data.frame': 6 obs. of 16 variables:
#> ..$ row : num [1:6] 1 1 1 1 1 1
#> ..$ col : num [1:6] 1 1 1 1 1 1
#> ..$ season : num [1:6] 1 2 3 4 5 6
#> ..$ time_start : num [1:6] 7.18 29.79 53.33 75.63 98.19 ...
#> ..$ time_end : num [1:6] 18.4 43.4 61.9 88.7 102.7 ...
#> ..$ time_peak : num [1:6] 12.8 35.8 57.6 82 100.2 ...
#> ..$ len : num [1:6] 11.23 13.57 8.6 13.03 4.47 ...
#> ..$ val_start : num [1:6] 0.392 0.462 0.458 0.444 0.452 ...
#> ..$ val_end : num [1:6] 0.413 0.459 0.466 0.435 0.458 ...
#> ..$ val_peak : num [1:6] 0.425 0.522 0.52 0.476 0.511 ...
#> ..$ val_base : num [1:6] 0.38 0.399 0.404 0.403 0.4 ...
#> ..$ ampl : num [1:6] 0.0453 0.1236 0.1157 0.0734 0.1111 ...
#> ..$ der_l : num [1:6] 0.018 0.0305 0.0322 0.0133 0.0492 ...
#> ..$ der_r : num [1:6] 0.00617 0.01563 0.02708 0.01506 0.02284 ...
#> ..$ integ_large: num [1:6] 5.44 7.96 4.98 6.94 2.88 ...
#> ..$ integ_small: num [1:6] 0.499 1.583 0.942 0.898 0.484 ...
#> NULL