bfastSpatial Peru
Jan Verbesselt and Loïc Dutrieux
02 February, 2016
Getting started!
After bfastSpatial is installed:
- define path for reading and saving files
- set directories
path <- '~/locationofyourscript'
library(bfastSpatial)
tmpDir <- rasterOptions()$tmpdir
# set the path the to the location of your script
inDir <- file.path(path, 'data')
# StepDir is where we store intermediary outputs
stepDir <- file.path(inDir, 'datastep')
# directory for landsat data
landsatDir <- file.path(stepDir, 'landsat')
# it is where individual NDVI layers will be stored before being stacked
# ndviDir is a subdirectory of stepDir
ndviDir <- file.path(stepDir, 'ndvi')
# OutDir is where the output will be stored
outDir <- file.path(inDir, 'out')
Check if all directories exist
## check if dir exist
for (i in c(stepDir, ndviDir, outDir, landsatDir)) {
print(dir.exists(i))
}
## [1] TRUE
## [1] TRUE
## [1] TRUE
## [1] TRUE
Download the archive with the data (~100MB)!
if (!file.exists(file.path(inDir, 'hungaryData.tar.gz'))) {
download.file('http://www.geo-informatie.nl/dutri001/hungaryData.tar.gz',
destfile = file.path(inDir, 'hungaryData.tar.gz'))
# Alternative download: https://googledrive.com/host/0B6NK5CdpGEa7ajFuQ1dBZHUwMW8
# Unpack the data
untar(file.path(inDir, 'hungaryData.tar.gz'), exdir = landsatDir)
}
Preprocessing the data: I
Check files:
# Check that we have the right data in the inDir folder
head(getSceneinfo(list.files(landsatDir)))
## sensor path row date
## LC81850272015108 OLI 185 27 2015-04-18
## LC81850272015124 OLI 185 27 2015-05-04
## LC81860272013125 OLI 186 27 2013-05-05
## LC81860272013141 OLI 186 27 2013-05-21
## LC81860272013157 OLI 186 27 2013-06-06
## LC81860272013173 OLI 186 27 2013-06-22
Note: the lines below take some time!
Preprocessing the data: II
if (!file.exists(file.path(inDir, 'ndvi_stack.grd'))) {
# unzip individual file, use the cloudmask, create ndvi (if this does not exists)
processLandsatBatch(x = landsatDir, outdir = ndviDir, srdir = tmpDir,
delete = TRUE, mask = 'fmask', vi = 'ndvi')
# Make temporal ndvi stack
ndviStack <- timeStack(x = ndviDir, pattern = glob2rx('*.grd'),
filename = file.path(inDir, 'ndvi_stack.grd'),
datatype = 'INT2S')
} else {
ndviStack <- brick(file.path(inDir, 'ndvi_stack.grd'))
}
The function ‘processLandsatBatch’ returns a list of TRUE, nothing to worry about, it means that everything went well
Let’s check output
head(list.files(ndviDir))
## [1] "ndvi.LC81850272015108.grd" "ndvi.LC81850272015108.gri"
## [3] "ndvi.LC81850272015124.grd" "ndvi.LC81850272015124.gri"
## [5] "ndvi.LC81860272013125.grd" "ndvi.LC81860272013125.gri"
Stack info
ndviStack
## class : RasterBrick
## dimensions : 176, 177, 31152, 397 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 569805.5, 575115.5, 5275605, 5280885 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : /Users/verbe039/Git/ChangeMonitor-WUR/Talks/bfastSpatial-2016/data/ndvi_stack.grd
## names : LE71850272000123, LE71850272000155, LE71850272000235, LE71850272000283, LE71850272001141, LE71850272001157, LE71850272001189, LE71850272001237, LE71850272002128, LT51850272002136, LT51850272002168, LE71850272002176, LE71850272002208, LE71850272002240, LE71850272003067, ...
## min values : NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## max values : NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## time : 2000-02-03, 2015-06-12 (min, max)
BfmPixel I
Run BfastMonitor for one pixel time series, by clicking on the map:
# Plot a cloud free recent NDVI layer
names(ndviStack)[394]
## [1] "LE71860272015139"
plot(ndviStack, 394)
BfmPixel II
# Call bfmPixel in interactive mode
# to test settings before running it on the brick
bfmPixel(x = ndviStack, start = 2014,
formula=response~harmon, order=1, history = c(2008, 7),
interactive = TRUE,
plot = TRUE)
## $bfm
##
## BFAST monitoring
##
## 1. History period
## Stable period selected: 2008(24)--2013(341)
## Length (in years): 5.868493
## Model fit:
## (Intercept) harmoncos harmonsin
## 5663.495 -3222.949 -1218.873
## R-squared: 0.875512
##
##
## 2. Monitoring period
## Monitoring period assessed: 2014(1)--2015(163)
## Length (in years): 1.443836
## Break detected at: -- (no break)
##
##
## $cell
## [1] 17934
bfastSpatial I
Produce a change map between 2014 (start=) and end of the time series. The process takes about 20 min (can be monitored with the system.time() function)
formula = is a particularly important parameter to tune, we’ll see what influence it may have on the change results when looking at individual pixels
# Run bfastmonitor (that part takes a long time)
if (!file.exists(fn <- file.path(outDir, 'bfm_2015_harmon_1.grd'))) {
# make sure the filename is correct and the folder exists
bfm <- bfmSpatial(x = ndviStack, pptype = 'irregular',
start = 2014, history = c(2008, 7),
formula = response ~ harmon,
order = 1, mc.cores = 3,
filename = file.path(outDir, 'bfm_2015_harmon_1.grd')
)
} else {
bfm <- brick(fn)
}
bfastSpatial II
What’s the output like?
bfm
## class : RasterBrick
## dimensions : 176, 177, 31152, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 569805.5, 575115.5, 5275605, 5280885 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : /Users/verbe039/Git/ChangeMonitor-WUR/Talks/bfastSpatial-2016/data/out/bfm_2015_harmon_1.grd
## names : layer.1, layer.2, layer.3
## min values : 2014.195, -3836.923, NA
## max values : 2015.444, 3215.836, NA
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.
bfastSpatial output
plot(ndviStack[[204]], col = grey.colors(255), legend = F)
plot(bfm[[1]], add=TRUE)
bfastSpatial: timing of the change
# extract change raster
change <- raster(bfm, 1)
# look at the timing of change
months <- changeMonth(change)
# set up labels and colourmap for months
monthlabs <- c("jan", "feb", "mar", "apr", "may", "jun",
"jul", "aug", "sep", "oct", "nov", "dec")
cols <- rainbow(12)
plot(months, col=cols, breaks=c(1:12), legend=FALSE)
legend("bottomright", legend=monthlabs, cex=0.5, fill=cols, ncol=2)
bfastSpatial: Magnitude of the break
# extract magnitude of the raster and scale values between 0-1.
magn <- raster(bfm, 2) / 10000
# make a version showing only breakpoing pixels
magn_bkp <- magn
magn_bkp[is.na(change)] <- NA
opar <- par(mfrow=c(1, 2))
plot(magn_bkp, main="Magnitude of a breakpoint")
plot(magn, main="Magnitude: all pixels")
par(opar)
bfastSpatial: Set a treshold
# extract and rescale magnitude and apply a threshold
magn_bkp_thresh <- magn_bkp
magn_bkp_thresh[magn_bkp > -0.05] <- NA
# compare
library(colorspace) ## package for nice colors
cols <- heat_hcl(12)
plot(magn_bkp_thresh, main="magnitude < -0.05", col = cols)
bfastSpatial: check pixel time series
# one pixel time
bfmPixel(x = ndviStack, start = 2014, cell=11521,
formula=response~harmon, order=1, history = c(2008, 7),
plot = TRUE)
## $bfm
##
## BFAST monitoring
##
## 1. History period
## Stable period selected: 2008(24)--2013(341)
## Length (in years): 5.868493
## Model fit:
## (Intercept) harmoncos harmonsin
## 6363.5043 -865.9808 -610.0195
## R-squared: 0.479110
##
##
## 2. Monitoring period
## Monitoring period assessed: 2014(1)--2015(155)
## Length (in years): 1.421918
## Break detected at: 2015(51)
##
##
## $cell
## [1] 11521
