bfastSpatial Peru

Jan Verbesselt and Loïc Dutrieux

02 February, 2016

Getting started!

After bfastSpatial is installed:

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)

Derive the overal median value

if (!file.exists(fn <- file.path(outDir, 'medianVI.grd')))
  # median values for all layers
  medVI <- summaryBrick(ndviStack, fun=median, na.rm=TRUE, 
              filename = fn) else {
  medVI <- brick(fn)
}
plot(medVI/10000)

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)