BFAST Monitor support

[ceholden]

Create timeseries driver for BFAST monitor.

Plan:

1. [ ] Create proof of concept dataset to validate algorithm implementation against
2. [ ] @ceholden: Implement Python script to run BFAST monitor in R from Python using rpy2
3. [ ] Figure out what parameters, settings, etc. we'll want to be able to change while running BFAST monitor
4. [ ] Use code from #2 to create a timeseries driver for running BFAST monitor
5. [ ] Decide what results beyond the model fits and break points should be plotted and how we should visualize them

[ceholden]

Hi @bendv, @dutri001, @verbe039, @johanez, @bullocke, and @PontusOlofsson

As discussed on the phone, here are the steps I envision to produce a "timeseries driver" for running BFAST within QGIS:

1. [x] Someone from WUR: Create proof of concept dataset to validate algorithm implementation against
   • This can be as simple as linking me to one of your help pages in any of the BFAST packages. I just want to make sure I have the correct method for pre-processing the data prior to ingest into one of the BFAST functions.
   • Ideally this would be a "minimal reproducible example" R script that I can wrap further down the line
2. [ ] @ceholden: Implement Python script to run BFAST monitor in R from Python using rpy2
   • This step simply involves wrapping the R code from #1 in the rpy2 Python package
   • rpy2 handles translation from np.array in Python to the vector R datatype. Further conversion to zoo or whatever can be done in R code if needed
3. [ ] All of us: Figure out what parameters, settings, etc. we'll want to be able to change while running BFAST monitor so we can create GUI widgets
   • This should be pretty easy since @dutri001 already built a Shiny app that does this
4. [ ] @ceholden: Use code from #2 to create a timeseries driver for running BFAST monitor
   • This just involves me combining the R script that runs BFAST to the other components of the "timeseries driver" API
5. [ ] All of us: Decide what results beyond the model fits and break points should be plotted and how we should visualize them

Please comment and make suggestions or changes if you like.

[bendv]

Thanks for this @ceholden. I guess you meant to mention @johanez as well?

[ceholden]

Whoops! Sorry for the spam, Johannes!

And sorry for confusing your username, @johanez!

[loicdtx]

Hi Chris, Regarding #1 I believe there is some material here. The code brings you from the raw landsat SR data as they are delivered by ESPA (these should still be available here) to a bfast friendly multi-layer object. Also conclusions of that training are well in line with the idea of the TStool.

Once you have a vector of values and a vector of dates, running bfastmonitor in R can be as simple as:

runPixel <- function(x, dates, ...) {
    ts <- bfastts(x, dates=dates, type='irregular')
    bfm <- bfastmonitor(ts, ...)
    return(bfm)
}

In the shiny app, it looks like I had chosen for start=, formula=, order =, and history= to be passed to the ellipsis. I feel that's really the bare minimum though.

[janverbesselt]

Hi Chris, thanks! You can also check examples include within the R BFAST package when you go to the help for the bfastmonitor function (see ?bfastmonitor()). As for visualisation see e.g. the standard/basic plotting of the bfastmonitor() result in the example or see also Ben's nice ggplot visualisation (https://github.com/bendv/bfastPlot) with some examples.

[ceholden]

Great! Thanks Loïc and Jan. I'll start up on the second and fourth tasks and tag commits with this issue.

How do you all store your data once it's been downloaded and pre-processed? Do you store each observation in separate files or some other way? I guess the most important bit to know would be if it's "stacked" (same bands, same raster size, same extent, and same projection / geo-transform). It's obviously not hard to deal with a non-stacked scenario, but I haven't had the need yet to write code for that situation.

We've been storing everything in more or less the same method since forever. 8 band raster images (7 Landsat bands + Fmask) that are the same size/extent/projection/datatype/etc, each observation stored in a folder named according to the Landsat ID. There's an example dataset on this repo: https://github.com/ceholden/landsat_stack