[Bfast-commits] Trend extraction

[Jan.Verbesselt at csiro.au]

Hi Brain,

See below for comments

> -----Original Message-----
> From: Brian J. Buma [mailto:Brian.Buma at Colorado.EDU]
> Sent: Monday, 22 February 2010 7:16 AM
> To: bfast-commits at r-forge.wu-wien.ac.at
> Subject: Trend extraction
> 
> I have a couple questions regarding bfast trend extraction- hopefully I'm
> not missing something super obvious.
> 
> I'm looking at a system which experienced several disturbances in a short
> time period, and attempting to look at the NDVI response, specifically the
> magnitude of the disturbance impact (the final disturbance, that is) in
> relation to previous disturbances, and also to look at NDVI recovery post-
> disturbance.
> 
> 
> The first question involves the seasonal component.  I have a situation
> where a disturbance dropped the NDVI in 2002 (using MODIS), so I have two
> years or so pre-disturbance.  The recovery period is on-going.  Bfast fits
> the seasonal component over the entire time series, right, so the more
> severe disturbances are fitted with a lower seasonal component because the
> years following the disturbance are lower, thus comparing the magnitude of
> the break around the disturbance itself doesn't work- the more severe
> disturbances have a reduced breakpoint magnitude, because their seasonal
> component is lower overall.  Hopefully that makes sense.  Any ideas for
> making comparisons more equitable would be appreciated- I'm thinking of
> attempting to get bfast to fit the seasonal component on only the pre-
> disturbance years, and then use that for the following recovery as the
> "baseline."  Is that possible?
> 

I'm working on a BFAST update - which is more sensitive to detect seasonal breaks. As such, you would be able to account for the different seasonality before and after the break. Seasonal break detection depends on how significant the change in seasonality is before and after the disturbance. Do you have an example (script/data) time series on which I can run some tests?

> 
> I'm trying to extract recovery rates in NDVI post-disturbance.  I'd really
> like to get the slope of the trend line component in bfast.  Is there a
> way to extract that number specifically?  Since I'm covering a large
> extent (1500+ time series) it needs to be automated.

Try the following code: it is extracts the most relevant info out of a bfast fit.

fit <- bfast()

# Extract BFAST info
      niter <- length(fit$output) # info from the last 
      out <- fit$output[[niter]]

      if (out$Vt.bp[1]>0) {  # only if there are breakpoints
        (nrbps <- length(out$Vt.bp))    # nr of breaks
        (tnrbps <- time(iNDVI)[(out$Vt.bp)])        # Time of the breaks!
        (coefs <- (coef(out$bp.Vt)))      # output coefficients of the model per period
        (int <- coefs[,1])          # intercepts
        (sl <-  coefs[,2])          # slopes
        (per <- names(coefs[,1]))   # periods
        (confs <- confint(out$bp.Vt))
        (lconfs <- confs$confint[,3]-confs$confint[,1])  # length of the confidence intervals
        (startts <- time(iNDVI)[1])
        (endts <- time(iNDVI)[length(iNDVI)])

        for (k in 1:(nrbps)) {
          if (k==1) {
            (dout <- as.data.frame(t(c(i,gids[i,1],nrbps,k,int[k],sl[k],per[k],lconfs[k],startts,tnrbps[k],int[k+1],sl[k+1]))))
          } else {
            (dout <- as.data.frame(t(c(i,gids[i,1],nrbps,k,int[k],sl[k],per[k],lconfs[k],tnrbps[k-1],tnrbps[k],int[k+1],sl[k+1]))))
          }
          names(dout) <- c("nr","gid","nrbps","nrbp","intercept1","slp1","per","lconfint","tsstart","tsend","intercept2","slp2")
          print(dout)
         }
     }

Please let me know if it works for you. Once tested, I'll include it in the help files of next bfast version.

Cheers,
Jan

> 
> Thanks for all your help!
> 
> -brian
> 
> 
> 
> ---
> 
> Brian Buma
> PhD Student
> Dept. of Ecology and Evolutionary Biology
> University of Colorado, Boulder
> Brian.Buma at colorado.edu