Implements bfastmonitor function, from the bfast package on any kind of rasterBrick object. Time information is provided as an extra object and the time series can be regular or irregular.
bfmSpatial(x, dates = NULL, pptype = "irregular", start, monend = NULL,
formula = response ~ trend + harmon, order = 3, lag = NULL,
slag = NULL, history = c("ROC", "BP", "all"), type = "OLS-MOSUM",
h = 0.25, end = 10, level = 0.05, mc.cores = 1,
returnLayers = c("breakpoint", "magnitude", "error"), sensor = NULL, ...)rasterBrick or rasterStack object, or file name to a multilayer raster object stored on disk.
A date vector. The number of dates must match the number of layers of x.
Character. Type of preprocessing to be applied to individual time series vectors. The two options are 'irregular' and '16-days'. See bfastts for more details.
See bfastmonitor
Numeric. Optional: end of the monitoring period in the format c(year, julian day). All raster data after this time will be removed before running bfastmonitor
See bfastmonitor
See bfastmonitor
See bfastmonitor
See bfastmonitor
See bfastmonitor
See bfastmonitor
See bfastmonitor
See bfastmonitor
Numeric. Number of cores to be used for the job.
Character. Result layers to be returned. Can be any combination of c("breakpoint", "magnitude", "error", "history", "r.squared", "adj.r.squared", "coefficients"). By default, breakpoint, magnitude and error are returned by the function. See details for more information.
Character. Optional: Limit analysis to one or more particular sensors. Can be any combintation of c("ETM+", "ETM+ SLC-on", "ETM+ SLC-off", "TM", or "OLI")
Arguments to be passed to mc.calc
A rasterBrick with layers depending on what has been supplied to returnLayers. By default, 3 layers are returned: (1) breakpoint: timing of breakpoints detected for each pixel; (2) magnitude: the median of the residuals within the monitoring period; (3) error: a value of 1 for pixels where an error was encountered by the algorithm (see try), and NA where the method was successfully run. See bfastmonitor for more information on the other possible layers.
bfmSpatial applies bfastmonitor over a raster time series. For large raster datasets, processing times can be long. Given the number of parameters that can be set, it is recommended to first run bfmPixel over some test pixels or bfmSpatial over a small test area to gain familiarity with the time series being analyzed and to test several parameters.
Note that there is a difference between the monend argument included here and the end argument passed to bfastmonitor. Supplying a date in the format c(year, Julian day) to monend will result in the time series being trimmed before running bfastmonitor. While this may seem identical to trimming the resulting bfastmonitor object per pixel, trimming the time series before running bfastmonitor will have an impact on the change magnitude layer, which is calculated as the median residual withint the entire monitoring period, whether or not a breakpoint is detected.
While bfmSpatial can be applied over any raster time series with a time dimension (implicit or externally supplied), an additional feature relating to the type of Landsat sensor is also included here. This feature allows the user to specify data from a particular sensor, excluding all others. The sensor argument accepts any combination of the following characters (also see getSceneinfo): "all" - all layers; "TM" - Landsat 5 Thematic Mapper; "ETM+" - Landsat 7 Enhanced Thematic Mapper Plus (all); "ETM+ SLC-on" - ETM+ data before failure of Scan Line Corrector; ETM+ data after failure of the Scan Line Corrector.
Note that for the sensor argument to work, names(x) must correspond to Landsat sceneID's (see getSceneinfo), otherwise any values passed to sensor will be ignored with a warning.
returnLayers can be used to specify which bfasmonitor results to return. Regardless of which parameters are assigned, the output layers will always follow the order: c("breakpoint", "magnitude", "error", "history", "r.squared", "adj.r.squared", "coefficients"). This is important if mc.cores is set to be greater than 1, since this causes the layer names in the output brick to be lost, so it is important to know which layers have been requested and in which order they will be exported. Note that if "coefficients" is included, the output will include the following: "(Intercept)" and any trend and/or harmonic coefficients depending on the values of formula and order.
# NOT RUN {
# load tura dataset
data(tura)
# run BFM over entire time series with a monitoring period starting at the beginning of 2009
t1 <- system.time(bfm <- bfmSpatial(tura, start=c(2009, 1)))
# }
# NOT RUN {
# with multi-core support (see ?mc.calc)
t2 <- system.time(bfm <- bfmSpatial(tura, start=c(2009, 1), mc.cores=4))
# difference processing time
t1 - t2
# }
# NOT RUN {
# plot the result
plot(bfm)
# }
Run the code above in your browser using DataLab