exec/Figure2_AbovegroundPheno_revamp.R

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

library(ggplot2)
library(gridExtra)
library(scales)
library(HarmonicRegression)
setwd("D:/R_Workspace/SDEF/SDEF.analysis/inst/extdata/")
# 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_03_2019.csv")
#load sap flow data
setwd("D:/R_Workspace/SDEF/SDEF.analysis/data")
load("Sapflow_normalizedtomax_14dayavg.Rdata")


#diurnal sap flow for points
load("Adenostoma_2016_mphillips_processed_noavg.rda")
ad <- Adenostoma_2016_mphillips_processed_noavg
ad$date <- as.Date(ad$day, format = "%Y-%m-%d")
ad$Native <- ((ad$max_flux)/max(ad$max_flux, na.rm = TRUE))

# subset out the NDVI that are the extra poimf
NDVI <- NDVI[c(1:34), ]
# fix date structure
NDVI$date <- strptime(NDVI$date, format = "%Y%m%d")
str(NDVI$date)
NDVI$date <- as.Date(NDVI$date, format = "%Y-%m-%d")
#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))
grasspheno <- grasspheno[order(grasspheno$date), ]
grasspheno <- grasspheno[-c(1:4), ] # remove first two september dates
grasspheno$measurement <- c(1:30)
grasspheno <- subset(grasspheno, grasspheno$NDVI > 0.1)

grasspheno$normalpheno <- ((grasspheno$NDVI)/max(grasspheno$NDVI))
#fit harmonic regressiont to data
gp <- grasspheno$NDVI
gt <- grasspheno$measurement
grass.fit <- harmonic.regression(inputts = gp, inputtime = gt)
grasspheno$harm <- grass.fit$fit.vals
grasspheno$Invasive <- ((grasspheno$harm)/max(grasspheno$harm))


# plot


colorpheno <-ggplot(aes(x = date, y = Native), data = sap)+ geom_line(aes(colour = "Native"),  size = 1.2) +
  geom_ribbon(aes(ymin=sap$min, ymax=sap$max, colour = "Native"), linetype=2, alpha=0.1, size = 0.8) +
  geom_point(data = ad, aes(x = date, y = max_flux, colour = "Native"), size = 2) +
  geom_line(data = grasspheno, aes(x =  date, y = Invasive, colour = "Invasive"),  size = 1.2)+
  geom_point(data = grasspheno, aes(x =  date, y = normalpheno, colour = "Invasive"),shape = 17,   size = 2)
colorpheno <- colorpheno + xlab("Time (month)") + scale_x_date(labels = date_format("%b-%y"), date_breaks = "2 month") +
  ylab(expression(paste(" Proportion of maximum phenological activity")))
colorpheno <-colorpheno+ theme_bw(base_family = "sans", base_size = 13) +
  scale_colour_manual("", breaks = c("Native", "Invasive"), values=c("black", "grey50"))
colorpheno

#"#E58601", "#E2D200"
#"black", "grey50"