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)
##
## Attaching package: 'bfast'
##
## The following object is masked _by_ '.GlobalEnv':
##
## dates
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, 11), history=c("all"), type="OLS-MOSUM", level=0.01)
bfm <- bfastmonitor(data = evits, start=c(2004, 11), formula = response ~ harmon, level=0.01, history=c("all"))
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
## 0.789816667 0.001616667 0.793866667
# set meaningful layernames
names(timeofbreak) = c("breakpoint", "magnitude", "historyStart", "historyEnd", "monitorStart", "monitorEnd")
# save the resulting raster object
# writeRaster(timeofbreak, filename="minisampleBfm2.grd", format="raster", overwrite=TRUE, options=c("COMPRESS=NONE"))
##############################
### Analyze the bfastmonitor result: Plotting the magnitude
##############################
### minisampleBfm.gri
### bfm <- bfastmonitor(data = evits, start=c(2004, 11), history=c("all"), type="OLS-MOSUM", level=0.01)
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)
options(scipen=999)
magVals = na.omit(getValues(minisampleBfmRaster))
plot(minisampleBfmRaster)
plot(minisampleBfmRaster, col = heat.colors(4), zlim=c(-0.2, 0))
plot(prodesPolygon, add=TRUE)
plot(minisampleBfmRaster, col = rainbow(4), zlim=c(-0.2, 0), breaks=c(0, -0.05, -0.1, -0.15, -0.2))
plot(prodesPolygon, add=TRUE)
# negative magnitudes
redPurples = c("#f1eef6", "#d7b5d8", "#df65b0", "#dd1c77", "#980043")
fj5NegMag = classIntervals(magVals[magVals < 0], n=5, style="fisher", dataPrecision=2)
plot(minisampleBfmRaster, col = rev(redPurples), zlim=c(round(fj5NegMag$brks[1], 2), 0),
breaks=round(fj5NegMag$brks, 3),
main="Change magnitudes of bfastmonitor", legend=FALSE)
# positive magnitudes
yellowGreens = c("#addd8e", "#78c679", "#41ab5d", "#238443", "#005a32")
fj5PosMag = classIntervals(magVals[magVals > 0], n=5, style="fisher", dataPrecision=2)
plot(minisampleBfmRaster, col = (yellowGreens), zlim=c(0, round(fj5PosMag$brks[5], 2)),
breaks=round(fj5PosMag$brks, 3), add=TRUE, legend=FALSE)
plot(prodesPolygon, add=TRUE)
plot(minisampleBfmRaster, legend.only=TRUE, legend.shrink=1, legend.width=1,
zlim=c(round(fj5NegMag$brks[1], 3), round(fj5PosMag$brks[5], 3)),
breaks=c(round(fj5NegMag$brks, 3), round(fj5PosMag$brks, 3)[-1]),
col=c(rev(redPurples), yellowGreens),
axis.args=list(cex.axis=0.6))
### minisampleBfm2.gri
### bfm <- bfastmonitor(data = evits, start=c(2004, 11), formula = response ~ harmon, level=0.01, history=c("all"))
minisampleBfm = "/home/christopher/Documents/Master Thesis/PRODES/Manicore/minisample/minisampleBfm2.gri"
minisampleBfmStack = stack(minisampleBfm)
# select the magnitude layer
minisampleBfmRaster = raster(minisampleBfmStack, 2)
plot(minisampleBfmRaster)
plot(prodesPolygon, add=T)
options(scipen=999)
magVals = na.omit(getValues(minisampleBfmRaster))
plot(minisampleBfmRaster)
plot(minisampleBfmRaster, col = heat.colors(4), zlim=c(-0.2, 0))
plot(prodesPolygon, add=TRUE)
plot(minisampleBfmRaster, col = rainbow(4), zlim=c(-0.2, 0), breaks=c(0, -0.05, -0.1, -0.15, -0.2))
plot(prodesPolygon, add=TRUE)
# negative magnitudes
redPurples = c("#f1eef6", "#d7b5d8", "#df65b0", "#dd1c77", "#980043")
fj5NegMag = classIntervals(magVals[magVals < 0], n=5, style="fisher", dataPrecision=2)
plot(minisampleBfmRaster, col = rev(redPurples), zlim=c(round(fj5NegMag$brks[1], 2), 0),
breaks=round(fj5NegMag$brks, 3),
main="Change magnitudes of bfastmonitor", legend=FALSE)
# positive magnitudes
yellowGreens = c("#addd8e", "#78c679", "#41ab5d", "#238443", "#005a32")
fj5PosMag = classIntervals(magVals[magVals > 0], n=5, style="fisher", dataPrecision=2)
plot(minisampleBfmRaster, col = (yellowGreens), zlim=c(0, round(fj5PosMag$brks[5], 2)),
breaks=round(fj5PosMag$brks, 3), add=TRUE, legend=FALSE)
plot(prodesPolygon, add=TRUE)
plot(minisampleBfmRaster, legend.only=TRUE, legend.shrink=1, legend.width=1,
zlim=c(round(fj5NegMag$brks[1], 3), round(fj5PosMag$brks[5], 3)),
breaks=c(round(fj5NegMag$brks, 3), round(fj5PosMag$brks, 3)[-1]),
col=c(rev(redPurples), yellowGreens),
axis.args=list(cex.axis=0.6))
plot(minisampleBfmRaster, col=rev(redPurples), zlim=c(-0.112, -0.014), breaks=c(round(fj5NegMag$brks, 3))[-6])
plot(prodesPolygon, add=TRUE)
plot(minisampleBfmRaster, col=rev(redPurples), zlim=c(-0.112, -0.024), breaks=c(round(fj5NegMag$brks, 3))[-5:-6])
plot(prodesPolygon, add=TRUE)
plot(minisampleBfmRaster, col=rev(redPurples), zlim=c(-0.112, -0.039), breaks=c(round(fj5NegMag$brks, 3))[-4:-6])
plot(prodesPolygon, add=TRUE)
plot(minisampleBfmRaster, col=rev(redPurples), zlim=c(-0.112, -0.077), breaks=c(round(fj5NegMag$brks, 3))[-3:-6])
plot(prodesPolygon, add=TRUE)
##############################
### Plotting the time of change point
##############################
minisampleBfmBrkpt = raster(minisampleBfmStack, 1)
plot(minisampleBfmBrkpt, col = heat.colors(8))
plot(prodesPolygon, add=TRUE)
plot(minisampleBfmBrkpt, col = heat.colors(5), zlim=c(2004.4347, 2007))
plot(prodesPolygon, add=TRUE)
##############################
### 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, afterwards press ESC
# 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)
##############################
### 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)
##
# dark green: 1022
# middle green: 385
# light green: 474
# dark purple inside polygon: 239 type 5
# middle purple inside polygon: 283 type 6
# middle purple inside polygon: 600 type 6
# light purple outside polygon: 348 no type
# dark purple inside polygon: 903 type 5
timeSeries = ts(as.numeric(rasterStackTs[238,]), start=c(2000, 4), end=c(2014, 17), frequency=23)
plot(timeSeries)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c("all"), type="OLS-MOSUM", level=0.01, end=10)
plot(bfm$model)
abline(v=2005.66, col="red") # start of PRODES 2006 data set
abline(v=2006.58, col="red") # end of PRODES 2006 data set
abline(v=2)
bf01 = bfast01(timeSeries)
plot(bf01)
bf01cl = bfast01classify(bf01)
bf01cl
## flag_type flag_significance p_segment1 p_segment2 pct_segment1
## 1 4 1 0.0000001005356 0.130925 -1.843347
## pct_segment2 flag_pct_stable
## 1 -0.3216709 NA
###
bfm = bfastmonitor(timeSeries, start=c(2004, 11), formula = response ~ harmon + trend, order=3, level=0.01)
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), formula = response ~ harmon, level=0.01, history=c("all"))
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), formula = response ~ harmon, order=4, level=0.01, history=c("all"))
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), formula = response ~ harmon, order=3, level=0.01, history=c("ROC"))
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), level=0.01, history=c("BP"))
plot(bfm)
bfm = bfastmonitor(timeSeries, start=c(2004, 11), history=c(2000,2))
plot(bfm)
