Normalized Difference Vegetation Index (NDVI) is a satellite-derived global vegetation indicator based on vegetation reflectance that proivdes information on vegetation productivity and phenology (Hamel et al. 2009). NDVI data comes in 15-day composite images, values range from 0.0 to 1.0, and data manipulation is required to identify the parameter of interest (e.g., peak timing of high quality vegetation). For our purpose, we will determine maximum increase between successive NDVI sampling periods and the sum of bimonthly values for each month they area available (i.e., May, July, September) similar to previous research (Hamel et al. 2009).
- Exercise 9.4 - Download and extract zip folder into your preferred location
- Set working directory to the extracted folder in R under File - Change dir...
-
First we need to load the packages needed for the exercise
We will start by importing individual rasters of our study site for each time period NDVI was available. We will then combine all layers into a raster "stack" for ease of manipulation.
library(sp)
library(lattice)
library(rgdal)#readOGR
library(rgeos)#gIntersection
library(raster)#to use "raster" function
library(adehabitatHR)
library(maptools)#readAsciiGrid
library(zoo)
ndviMay09 <- raster("may9_final")
ndviMay25 <- raster("may25_final")
ndviJune10 <- raster("june10_final")
ndviJuly12 <- raster("july12_final")
ndviJuly28 <- raster("july28_final")
ndviAug29 <- raster("aug29_final")
ndviSep14 <- raster("sept14_final")
ndviSep30 <- raster("sept30_final")
proj4string(ndviMay09)
#Create a Stack of all Rasters
r <- stack(list(ndviMay09=ndviMay09, ndviMay25=ndviMay25, ndviJune10=ndviJune10, ndviJuly12=ndviJuly12,ndviJuly28=ndviJuly28,ndviAug29=ndviAug29,
ndviSep14=ndviSep14, ndviSep30=ndviSep30))
nlayers(r)
plot(r)
names(r) -
Now open the script "SnowyNDVIprep.R" and run code directly from the script
-
We will start by importing individual rasters of our study site for each time period NDVI was available. We will then combine all layers into a raster "stack" for ease of manipulation.
ndviMay09 <- raster("snowyMay09.tif")
ndviMay25 <- raster("snowyMay25.tif")
ndviJune10 <- raster("snowyJun10.tif")
ndviJuly12 <- raster("snowyJuly12.tif")
ndviJuly28 <- raster("snowyJuly28.tif")
ndviAug29 <- raster("snowyAug29.tif")
ndviSep14 <- raster("snowySep14.tif")
ndviSep30 <- raster("snowySep30.tif")
proj4string(ndviMay09)
#Create a Stack of all Rasters
r <- stack(list(ndviMay09=ndviMay09, ndviMay25=ndviMay25, ndviJune10=ndviJune10,
ndviJuly12=ndviJuly12,ndviJuly28=ndviJuly28,ndviAug29=ndviAug29,
ndviSep14=ndviSep14, ndviSep30=ndviSep30))
nlayers(r)
plot(r)
names(r) - Then we can get a mean for each NDVI layer for each 13 km2 cell of our sampling grid by extracting all rasters within each sampling grid cell. We will start first by creating our own functions using the raster package.
#We will start by creating a function to determine maximum increase between successive NDVI periods
#NOTE: x[1] refers to the first layer in your raster stack
maxmay <- function(x){x[1]-x[2]}
maxjune <- function(x){x[2]-x[3]}
maxjuly1 <- function(x){x[3]-x[4]}
maxjuly2 <- function(x){x[4]-x[5]}
maxaug <- function(x){x[5]-x[6]}
maxsept1 <- function(x){x[6]-x[7]}
maxsept2 <- function(x){x[7]-x[8]}
#Now perform the function on each raster in the stack
may <- calc(r,maxmay)
june <- calc(r,maxjune)
july1 <- calc(r,maxjuly1)
july2 <- calc(r,maxjuly2)
aug <- calc(r,maxaug)
sept1 <- calc(r,maxsept1)
sept2 <- calc(r,maxsept2) - We can plot out each layer in 4 x 3 dimensions if we want to look over what the function created and determine values that resulted from performing the functions on each raster group.
#Set up the figure layout
par(mfcol=c(3,3),mar=c(2,3.5,2,2),oma=c(3,3,3,3)) #Bottom,Left,Top,Right.
plot(may)
# Create a title with a red, bold, italic font
title(main="May", col.main="black", font.main=4)
plot(june)
# Create a title with a red, bold, italic font
title(main="June", col.main="black", font.main=4)
plot(july1)
# Create a title with a red, bold, italic font
title(main="July 1", col.main="black", font.main=4)
plot(july2)
# Create a title with a red, bold, italic font
title(main="July 2", col.main="black", font.main=4)
plot(aug)
# Create a title with a red, bold, italic font
title(main="Aug", col.main="black", font.main=4)
plot(sept1)
# Create a title with a red, bold, italic font
title(main="Sept 1", col.main="black", font.main=4)
plot(sept2)
# Create a title with a red, bold, italic font
title(main="Sept 2", col.main="black", font.main=4) - Now we can use the raster package to create functions to determine sum of the bimonthly values for May, July, and September.
#Start by creating a function to sum the bimonthly values for May, July, and September
summay <- function(x){x[1]+x[2]}
sumjuly <- function(x){x[4]+x[5]}
sumsept <- function(x){x[7]+x[8]}
#Now perform the function on each month that we have 2 layers of NDVI
maysum <- calc(r,summay)
julysum <- calc(r,sumjuly)
septsum <- calc(r,sumsept) - Plot out each layer to look over what the function created and determine values that resulted from performing the functions for each month.
windows()#opens a new graphic window if needed.
par(mfcol=c(2,2))
plot(maysum)
# Create a title with a red, bold, italic font
title(main="May", col.main="black", font.main=4)
plot(julysum)
# Create a title with a red, bold, italic font
title(main="July", col.main="black", font.main=4)
plot(septsum)
# Create a title with a red, bold, italic font
title(main="Sept", col.main="black", font.main=4) - Now we need to get a mean for each covariate for each 13 km2 cell of our sampling grid similar to our DEM layer means.
#Means for maximum increase
extmay <- extract(may, grid, weights=TRUE, fun = mean)
extjune <- extract(june, grid, weights=TRUE, fun = mean)
extjuly1 <- extract(july1, grid, weights=TRUE, fun = mean)
extjuly2 <- extract(july2, grid, weights=TRUE, fun = mean)
extaug <- extract(aug, grid, weights=TRUE, fun = mean)
extsept1 <- extract(sept1, grid, weights=TRUE, fun = mean)
extsept2 <- extract(sept2, grid, weights=TRUE, fun = mean)
#Means for bimontly sums
extmaysum <- extract(maysum, grid, weights=TRUE, fun = mean)
extjulysum <- extract(julysum, grid, weights=TRUE, fun = mean)
extseptsum <- extract(septsum, grid, weights=TRUE, fun = mean) - Now convert each matrix to a data frame so it is easier to manipulate then combine into a single dataset.
mayNDVImax <- as.data.frame(extmay)
juneNDVImax <- as.data.frame(extjune)
jul1NDVImax <- as.data.frame(extjuly1)
jul2NDVImax <- as.data.frame(extjuly2)
augNDVImax <- as.data.frame(extaug)
sep1NDVImax <- as.data.frame(extsept1)
sep2NDVImax <- as.data.frame(extsept2)
#Means for bimontly sums
mayNDVIsum <- as.data.frame(extmaysum)
julyNDVImax <- as.data.frame(extjulysum)
sepNDVImax <- as.data.frame(extseptsum)
#Bind all NDVI layers created above
Snowy_NDVI <- cbind(maySnowymax,juneSnowymax,jul1Snowymax, jul2Snowymax, augSnowymax,sep1Snowymax,sep2Snowymax,maySnowysum,julySnowysum,sepSnowysum)
colnames(Snowy_NDVI)<-c("maySnowymax","juneSnowymax","jul1Snowymax",
"jul2Snowymax","augSnowymax","sep1Snowymax","sep2Snowymax","maySnowysum","julySnowysum","sepSnowysum")
head(Snowy_NDVI)
#Now assign actual sampling grid IDs to the covariate summaries in "Snowy_NDVI"
Snowy_NDVI$id <- paste(snowygrid@data$id)
head(Snowy_NDVI)
write.table(Snowy_NDVI,"SnowyNDVI.txt", sep="\t")
#Save workspace so all data is available
save.image("Snowy_NDVI.RData") -
Now let's combine Habitat and DEM data from previous exercise to NDVI data just created.
#First let's import the previous NLCD and DEM data we created
snowy2 <- read.csv("SnowyData.txt", sep="\t")
snowy2[1:5,]
#Now we need to "join" the appropriate covariate data to each moose sampled based on the sampling grid cell it occurred in (e.g., g953)
Snowy_NDVI[1:5,]
SnowyFinal <- merge(snowy2, Snowy_NDVI, by=c("id"))
SnowyFinal[1:10,]
write.table(SnowyFinal,"SnowyFinal.txt", sep="\t")
#Copy file into R2WinBUGS exercise 9.5
