I am using the bfastlite funtion from the bfast package for time-series
analysis of nighttime light (NTL) imagery. I have monthly data from April
2013 to December 2022 (129 monthly images). Using the code below, the
bfastlite found 2 breaks in the time-series data:
library(bfast)
library(terra)
wd <- "path/"
l <- list.files(paste0(wd), pattern = '.tif$',
all.files = TRUE, full.names = FALSE)
rr <- lapply(paste0(wd, l), rast)
standard <- rr[[1]]
rs <- list(standard)
for (i in 2:length(rr)) {
rr[[i]] <- terra::crop(rr[[i]], standard, extend = TRUE)
}
s <- rast(rr)
m <- terra::as.matrix(s)
m <- m[!rowSums(is.na(m)), ]
# convert the matrix to timeseries matrix
tm <- ts(data = m, start = c(2013, 4), end = c(2022, 12), frequency = 12,
class = "ts")
bf <- bfastlite(
tm,
formula = response ~ harmon,
order = 3,
breaks = "all",
lag = 17,
slag = 7,
na.action = na.omit,
stl = "none",
decomp = c("stl", "stlplus"),
sbins = 1,
h = .25,
level = 0.5,
type = c("OLS-MOSUM", "RE"),
plot = TRUE,
hpc = "foreach"
)
bf
plot(bf)
How can I plot the location (i.e., pixels) of the breakpoints? You can
download the dataset from my GoogleDrive
<https://drive.google.com/drive/folders/1UloTKYHa7I3Wi_2pt3n3xbZ6vFvVz3_Y?usp=sharing>.
The entire code runs in less than 10 seconds.
--
Tziokas Nikolaos
Cartographer
Tel:(+44)07561120302
LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130>
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
