exec/Figure2_AbovegroundPheno.R

In bmcnellis/SDEF.analysis: Data Analysis for 2016 San Dimas Experimental Forest Project

library(ggplot2)
library(gridExtra)
library(scales)
setwd("D:/R_Workspace/SDEF/SDEF.analysis/inst/raw_data/")
# there is an extra point in the dataset cause i was lazy about fucking with my GEE script
# plot is in path_row_ 40_36 and 41_36
NDVI <- read.csv("SDEF_NDVI.csv")

# import sap flow data
setwd("D:/R_Workspace/SDEF/SDEF.analysis/data")
# this data is on a 0-1 scale... relative to maximum flux
load("Adenostoma_2016_mphillips_processed.rda")
sap <- Adenostoma_2016_mphillips_processed

# subset out the NDVI that are the extra poimf
NDVI <- NDVI[c(1:21), ]
# fix date structure
NDVI$date <- strptime(NDVI$date, format = "%Y%m%d")
str(NDVI$date)
NDVI$date <- as.POSIXct(NDVI$date)
#keep only parts to calulate NDVI
grasspheno <- NDVI[, c(5,6,9)]
# NDVI = (IR-R)/(IR+R) Band is IR and Band 4 is R for Landsat8
grasspheno$NDVI <- ((NDVI$B5 - NDVI$B4)/ (NDVI$B5 + NDVI$B4))

# this is to relativize the NDVI to 0-1-- max NDVI is 1
grasspheno$normalpheno <- ((grasspheno$NDVI)/max(grasspheno$NDVI))
# fix date stucture
sap$day <- as.Date(sap$day)
grasspheno$date <- as.Date(grasspheno$date)
grasspheno <- grasspheno[order(grasspheno$date),]
grasspheno$mesurement <- c(1:16)
sap$mesurement <- c(1:191)
# remove first five rows from grass pheno to match timing up with sap flow
grasspheno <- grasspheno[-c(1:5), ]


pheno<-ggplot(aes(x = day, y = max_flux  ), data = sap)+ geom_smooth(size = 1, color = "black") +
  geom_smooth(data = grasspheno, aes(x =  date, y = normalpheno ),  size = 1, color = "gray")
pheno <- pheno + xlab("Time (month)") + scale_x_date(labels = date_format("%b-%y"), date_breaks = "month") +
  ylab(expression(paste(" Proportion of maximum phenological activity")))
pheno<-pheno+ theme_bw(base_family = "sans", base_size = 13)
pheno