3

I have a stack of 300+ landsat NDVI images. I am using the BFAST package in R to identify breakpoints. The breaks are often hugely obvious as you can see from this image:
enter image description here

Note how there is a huge dip in NDVI around 1988 followed by a gradual increase. BFAST ignores the obvious break, and instead places a breakpoint around 1994 in the middle of the gradual increase.

I used the following R-Code to run BFAST:

bfast(ndvi.ts, h=.3, season="harmonic", max.iter=1, breaks=1)

(Tweaking the h parameter does not seem to improve the situation).

A few questions for everyone:

  1. Is there a parameter besides h that I can tweak to improve the results?
  2. If not, is there a way to massage the data to get better results?
  3. If not, is there another breakpoint analysis package in R that might yield better results?
rvp
  • 51
  • 5
  • 2
    The h parameter says that your minimal segment size in this case is 30% of the data. But your break occurs after about 15% of the data. Thus, you won't be able to catch it exactly. Hence, the recommendation would be to use a bandwidth h that is small enough to catch relevant breaks but still large enough to get reasonable estimates from the regression model in each segment. – Achim Zeileis Apr 30 '15 at 16:38
  • Also beside the h parameter, it is import to select the right model (e.g. seasonal + trend model, or only a trend model) for this time series. Is a season+trend model appropriate here? Can you make this example reproducible? (see here for how to do this: https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example – Janvb Apr 30 '15 at 17:19
  • @rvp would you be welling to discuss this elsewhere? I'm very interested in this approach with BFAST. – Paulo E. Cardoso May 01 '15 at 14:55
  • Thanks for the help everyone. I misunderstood what h represents. Once I set h = .15, the breakpoint is properly identified as long as I restrict the number of breaks to 1. However, if I don't restrict the number of breaks, BFAST identifies a huge number of spurious breakpoints. Is there any way to specify a minimum magnitude for breaks? – rvp May 06 '15 at 14:04

1 Answers1

1

I have wrestled with the same analyses and ended up using the ecp package http://cran.r-project.org/web/packages/ecp/vignettes/ecp.pdf

In my workflow I extract my index per image date for a number of sites and put them in a data frame. Then from there I have a loop that calculates the changepoint and creates a ggplot of the timeseries with the changepoints indicated. I won't show the whole loop here however the relevant code for identifying the changepoints is something like:

library("ecp")
df <- "your data frame"
df2 <- df[ ,1]#assuming your data values are in first column
ecp.mat <- matrix(df2, ncol = 1)
ecp.out <- e.divisive(ecp.mat, R = 499, sig.lvl = sig, alpha = 1 )
ecp.est <- ecp.out$estimates[c(-1, -length(ecp.out$estimates))]#drop first and last records

The last object above contains the location of your changepoints and I then plot this. Hope that helps.

Bart
  • 136
  • 1
  • 8