Breakpoints and spatial plots when using bfastlite

[nikosGeography]

I am using the bfastlite function and I found some breaks in my time-series (TS) which I can plot (in a line graph) using the plot() function. I was wondering if there is a way to plot the breaks spatially (i.e., pixel(s)). I think you discussed this in #5 .

Is there any updates or a workaround rgarding #5 ? My code:

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](http://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. The entire code runs in less than 10 seconds.

Windows 11, R 4-3-2, RStudio 2023.12.0 Build 369.

[GreatEmerald]

Yes, you can run it spatially with terra::app(). But the question is what do you want to plot. You can get e.g. the number of breaks, or the timing of the first break, or the timing of the last break, or the magnitude of the first or last break etc. etc.

[nikosGeography]

I think the most important plots for my analysis are the number of break (I mean where are they located spatially, no matter the timing, magnitute etc), as well as the first break. So two plots.

[GreatEmerald]

Hmm, first break as in the earliest break in the time series? Most users are more interested in the latest break. But both are easy to do, I can make some example code and then include it in the examples in the package.

[nikosGeography]

That would be great. I'm interested in the first break because I am using satellite nighttime lights data and I want to correlate the time of the lockdown to the breakpoint in TS. If that makes sense.

Thank you.

[GreatEmerald]

Here's a Jupyter Notebook showing how to do it: https://colab.research.google.com/drive/1M3ICmcomwE29U-4ux66M4nM9Vhczzwpb?usp=sharing

I set it to output the last break, but you can easily get the first break by changing breaks[length(breaks)] to breaks[1]. I'll see how to best incorporate it into the examples.

[nikosGeography]

That's great, many thanks. I will check it this week on my dataset. Looking forward to the updated version in the examples in the package. If that's okay with you, I'd like to keep the post open until I test it on my dataset.

[GreatEmerald]

Yeap, I'll keep it open until I put it into the examples.