library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 0.9-1, (SVN revision 518)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.10.1, released 2013/08/26
## Path to GDAL shared files: /usr/share/gdal/1.10
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: (autodetected)
library(bfast)
library(classInt)
library(RColorBrewer)
library(rts)
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(stringr)
library(M3)
## Loading required package: ncdf4
## Loading required package: maps
## Loading required package: mapdata
##############################
### Crop raster to study site extent
##############################
rasterstackcrop = "/home/christopher/Documents/git/MasterThesis/spatialAccuracy/manicore.gri"
rasterStack = stack(rasterstackcrop)
spatialExtent = readOGR("/home/christopher/Documents/Master Thesis/PRODES/Manicore/minisample", layer="extentMinisample")
## OGR data source with driver: ESRI Shapefile
## Source: "/home/christopher/Documents/Master Thesis/PRODES/Manicore/minisample", layer: "extentMinisample"
## with 1 features and 1 fields
## Feature type: wkbPolygon with 2 dimensions
#spatialExtent = spTransform(spatialExtent, crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
#crs(spatialExtent)
#crs(rasterStack)
rasterStack = crop(rasterStack, spatialExtent)
# scale the EVI values
rasterStack = calc(rasterStack, function(x) x/10000)
prodesPolygon = readOGR("/home/christopher/Documents/Master Thesis/PRODES/Manicore/minisample", layer="clipProdesPolygon")
## OGR data source with driver: ESRI Shapefile
## Source: "/home/christopher/Documents/Master Thesis/PRODES/Manicore/minisample", layer: "clipProdesPolygon"
## with 1 features and 1 fields
## Feature type: wkbMultiPolygon with 2 dimensions
plot(rasterStack, 1)
plot(prodesPolygon, add=T)
# writeRaster(rasterStack, filename="minisample.grd", format="raster", overwrite="TRUE", options=c("COMPRESS=NONE"))
##############################
### APPLY BFASTMONITOR ON RASTER STACK
##############################
minisample = "/home/christopher/Documents/Master Thesis/PRODES/Manicore/minisample/minisample.gri"
rasterStack = stack(minisample)
# helper function
xbfastmonitor <- function(x, dates) {
evits <- bfastts(x, dates, type = c("16-day"))
# call of bfastmonitor
bfm <- bfastmonitor(data = evits, start=c(2004, 23), history = c("ROC"))
return(cbind(bfm$breakpoint, bfm$magnitude, bfm$history[1], bfm$history[2], bfm$monitor[1], bfm$monitor[2]))
}
ptm = proc.time()
timeofbreak = calc(rasterStack, fun=function(x){
res = t(apply(x, 1, xbfastmonitor, dates))
return(res)
})
runtimeInMinutes = (proc.time() - ptm) / 60
print(runtimeInMinutes)
## user system elapsed
## 1.255966667 0.004716667 1.548433333
# set meaningful layernames
names(timeofbreak) = c("breakpoint", "magnitude", "historyStart", "historyEnd", "monitorStart", "monitorEnd")
# save the resulting raster object
# writeRaster(timeofbreak, filename="minisampleBfm.grd", format="raster", overwrite="TRUE", options=c("COMPRESS=NONE"))
##############################
### Analyze the bfastmonitor result
##############################
minisampleBfm = "/home/christopher/Documents/Master Thesis/PRODES/Manicore/minisample/minisampleBfm.gri"
minisampleBfmStack = stack(minisampleBfm)
# select the magnitude layer
minisampleBfmRaster = raster(minisampleBfmStack, 2)
plot(minisampleBfmRaster)
plot(prodesPolygon, add=T)
##############################
### Create Raster Time Series
##############################
# get layer names from raster stack
layerNames = names(rasterStack)
# extracting julian date of layer name
dates = lapply(layerNames, function(x) str_extract(x, perl("(?<=X).*?(?=.MOD|$)")))
# transform julian date to date
dates = as.Date(decipher.M3.date(as.numeric(unlist(dates))))
rasterStackTs = rts(rasterStack, dates)
##############################
### Look at the time series inside and outside the PRODES polygons
##############################
# click in the plot get the pixel number to anaylze
# click(minisampleBfmRaster, cell=T)
# Time Series of pixels where change was reported by PRODES 2006
timeSeries = ts(as.numeric(rasterStackTs[422,]), start=c(2000, 4), end=c(2014, 17), frequency=23)
# type = c("OLS-CUSUM", "OLS-MOSUM", "RE", "ME", "fluctuation")
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-CUSUM")
plot(timeSeries, main="Change reported by PRODES 2006")
abline(v=2005.5)
abline(v=2006.5)
plot(bfm)
timeSeries = ts(as.numeric(rasterStackTs[332,]), start=c(2000, 4), end=c(2014, 17), frequency=23)
bfm = bfastmonitor(timeSeries, start=c(2004,11), history=c("BP"))
plot(timeSeries, main="Change reported by PRODES 2006")
abline(v=2005.5)
abline(v=2006.5)
plot(bfm)
# Time Series of pixels where no change was reported by PRODES 2006
timeSeries = ts(as.numeric(rasterStackTs[367,]), start=c(2000, 4), end=c(2014, 17), frequency=23)
plot(timeSeries)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"))
plot(bfm)
timeSeries = ts(as.numeric(rasterStackTs[306,]), start=c(2000, 4), end=c(2014, 17), frequency=23)
plot(timeSeries)
timeSeries = ts(as.numeric(rasterStackTs[187,]), start=c(2000, 4), end=c(2014, 17), frequency=23)
plot(timeSeries)
##############################
### bfastmonitor parameter explanation
##############################
# bfastmonitor(data, start, 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,
# hpc = "none", verbose = FALSE, plot = FALSE)
# data: time series
# start: start of the monitoring period
# formula: formula for a regression model
# order: order of the harmonic term, by default 3
# lag: order of the autoregressive term, by default omitted
# slag: order of the seasonal autoregressive term, by default omitted
# history: specifcation of the start of the stable history period. If character, then selection is
# possible between reverse-ordered CUSUM ("ROC", default), Bai and Perron breakpoint estimation ("BP"),
# or all available observations ("all"). If numeric, the start date can be specified in the same form
# as start. If a function is supplied it is called as history(formula, data) to compute a numeric start date.
# type: character specifying the type of monitoring process. By default, a MOSUM process based on OLS residuals
# is employed. Can be "OLS-CUSUM", "OLS-MOSUM", "ME", "RE"
# h: numeric scalar from interval (0,1) specifying the bandwidth relative to the sample size in MOSUM/ME monitoring
# processes.
# end : numeric. Maximum time (relative to the history period) that will be monitored (in MOSUM/ME processes).
# Default is 10 times the history period.
# h: Significance level of the monitoring (and ROC, if selected) procedure, i.e., probability of type I error.
# hpc: character specifying the high performance computing support. Default is "none", can be set to "foreach".
# verbose: logical. Should information about the monitoring be printed during computation?
# plot: logical. Should the result be plotted?
##############################
### Varying bfastmonitor parameters for pixel INSIDE a prodes polygon
##############################
timeSeries = ts(as.numeric(rasterStackTs[422,]), start=c(2000, 4), end=c(2014, 17), frequency=23)
plot(timeSeries, main="Change reported by PRODES 2006")
abline(v=2005.5)
abline(v=2006.5)
# ROC
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-MOSUM")
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-CUSUM")
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-MOSUM", h=0.5)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-CUSUM", h=0.5)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-MOSUM", h=1)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-CUSUM", h=1)
plot(bfm)
## BP and h
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("BP"), type="OLS-MOSUM")
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("BP"), type="OLS-CUSUM")
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("BP"), type="OLS-MOSUM", h=0.5)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("BP"), type="OLS-CUSUM", h=0.5)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("BP"), type="OLS-MOSUM", h=1)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("BP"), type="OLS-CUSUM", h=1)
plot(bfm)
# ALL
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("all"), type="OLS-MOSUM")
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("all"), type="OLS-MOSUM", level=0.01)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("all"), type="OLS-CUSUM")
plot(bfm)
# "good result" breakpoint detected inside the slope
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("all"), type="OLS-MOSUM", h=0.5)
plot(bfm)
abline(v=2005.66, col="orange") # start of PRODES 2006 data set
abline(v=2006.58, col="orange") # end of PRODES 2006 data set
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("all"), type="OLS-CUSUM", h=0.5)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("all"), type="OLS-MOSUM", h=1)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("all"), type="OLS-CUSUM", h=1)
plot(bfm)
##############################
### Varying bfastmonitor parameters for pixel OUTSIDE a prodes polygon
##############################
timeSeries = ts(as.numeric(rasterStackTs[367,]), start=c(2000, 4), end=c(2014, 17), frequency=23)
plot(timeSeries)
# ROC
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-MOSUM")
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-MOSUM", level=0.01)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-CUSUM")
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-CUSUM", level=0.01)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-MOSUM", h=0.5)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-MOSUM", h=0.5)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-CUSUM", h=0.5)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-MOSUM", h=1)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("ROC"), type="OLS-CUSUM", h=1)
plot(bfm)
## BP and h
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("BP"), type="OLS-MOSUM")
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("BP"), type="OLS-MOSUM", level=0.01)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("BP"), type="OLS-CUSUM")
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("BP"), type="OLS-MOSUM", h=0.5)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("BP"), type="OLS-CUSUM", h=0.5)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("BP"), type="OLS-MOSUM", h=1)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("BP"), type="OLS-CUSUM", h=1)
plot(bfm)
# ALL
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("all"), type="OLS-MOSUM")
plot(bfm)
#
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("all"), type="OLS-MOSUM", level=0.01)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("all"), type="OLS-CUSUM")
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("all"), type="OLS-MOSUM", h=0.5)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("all"), type="OLS-CUSUM", h=0.5)
plot(bfm)
# late breakpoint detected
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("all"), type="OLS-MOSUM", h=1)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("all"), type="OLS-CUSUM", h=1)
plot(bfm)