bfastmonitor.R

version

rm(list = ls())

getwd()

setwd("D:/Learn/rs4ecologist/")
getwd()

library(MODIS) #facilitates MODIS use
library(bfast) #change detection
library(raster) #facilitates raster use
library(rasterVis) #facilitates raster visualization
library(zoo) # time series handling
library(strucchange) # break detection
library(colorspace) #color pallette libraries

mdir = "D:/Learn/rs4ecologist/MODIS/NDVIMOD13Q1/"

#modis.raw are MODIS NDVI rasters taken every 16 days from 01/01/2002 to 01/17/2014.

modis.raw          = brick(paste(mdir,"ndvi_MOD13Q1_ndvi_brick.grd",sep=""))
modis              = brick(paste(mdir,"cvi_MOD13Q1_ndvi_brick.grd",sep=""))
names(modis)       = names(modis.raw) #assigns same layer names to derived data

#summarizing MODIS layers

summary(modis)

#exploratory analysis

indices.2002       = grep("MOD13Q1A2002", names(modis.raw)) # Find images from 2002
modis.2002         = subset(modis, indices.2002) # Create new Raster-Stack with these layers
names(modis.2002)  = names(modis.raw)[indices.2002] # Update layer names

#calculates the median and variance over all the images in 2002 on a pixel (x,y) basis

modis.2002.median  = calc(modis.2002,fun = function(x) median(x, na.rm = T))
modis.2002.var     = calc(modis.2002,fun = function(x) var(x, na.rm = T))

#plots the mean and median for 2002

plot(modis.2002.median/10000)
plot(modis.2002.var/10000)

#subset for 2005 and 2010

indices.2005       = grep("MOD13Q1A2005", names(modis.raw)) # Find images from 2002
modis.2005         = subset(modis, indices.2005) # Create new Raster-Stack with these layers
names(modis.2005)  = names(modis.raw)[indices.2005] # Update layer names

indices.2010       = grep("MOD13Q1A2010", names(modis.raw)) # Find images from 2002
modis.2010         = subset(modis, indices.2010) # Create new Raster-Stack with these layers
names(modis.2010)  = names(modis.raw)[indices.2010] # Update layer names

#calculates the median over all the images in 2005 and 2010 on a pixel (x,y) basis

modis.2005.median  = calc(modis.2005,fun = function(x) median(x, na.rm = T))
modis.2010.median  = calc(modis.2010,fun = function(x) median(x, na.rm = T))

#plot the medians as R-G-B (2002, 2005, 2010) composite

modis.rgb.median   = stack(modis.2002.median/10000,modis.2005.median/10000,modis.2010.median/10000)
plotRGB(modis.rgb.median,stretch="hist") #can also use linear stretch

#the following function is useful for extracting dates from the layer names

doy2date           = function(year, doy) {
  as.Date(doy - 1, origin = paste0(year, "-01-01")) # Convert year + day of year to Date
}

data.dates         = as.character(doy2date(year = 2002,doy = as.numeric(substr(names(modis.2002), 18, 20)))) #applied to the MODIS 2002 time series

modis.stack        = stack() #empty stack for monthly compositing

#this loop computes the monthly composites

for (i in 1:12) {
  bands_i          =  which(sapply(strsplit(data.dates, "-"),function(x) as.numeric(x[2])) == i) #finds all dates whose month == i
  modis.stack      =  stack(modis.stack, calc(subset(modis.2002,bands_i),fun = mean)) #calculates the mean of the corresponding layers and add it to the stack
}

names(modis.stack) = month.abb #assigns three-letter month abberviation to layer names

#plots monthly composites and corresponding boxplots

levelplot(modis.stack / 10000, par.settings = RdBuTheme)#red to blue color palette selected. White = No data, light-blue = non-forest, dark blue = forest, red = water
bwplot(modis.stack)