GE509_ClassProject.Rmd

log(---
title: "GE509_ClassProject Preliminary Analysis"
author: "Kathryn Wheeler"
output: html_document
---

```{r}
## Settings
library(rjags)
library(coda)
library("MODISTools")
library("numDeriv")

```

```{r}
##' Download Phenocam data
##' 
##' @param URL  web address where data is located
download.phenocam <- function(URL) {
  ## check that we've been passed a URL
  if (length(URL) == 1 & is.character(URL) & substr(URL,1,4)=="http") {
    
    ## read data
    dat <- read.csv(URL,skip = 22)
    
    ## convert date
    dat$date <- as.Date(as.character(dat$date))
    
    return(dat)
  } else {
    print(paste("download.phenocam: Input URL not provided correctly",URL))
  }
}



##' Create a Bayes Model for a deciduous broadleaf site
##' 
##' @param Lat  latitude of desired site in decimals
##' @param Long longitude of desired site in decimals
##' @param data.source data source (GOES.NDVI, MODIS.NDVI, PC.GCC)
##' @param site.name Site Name
##' @param URL PhenoCam network URL
createBayesModel.DB <- function(Lat=0, Long=0, data.source,site.name="",URL="",niter=30000) {
  nchain = 5
  inits <- list()
  if(data.source=="PC.GCC"){
    data = PC_data(URL)
    print(data$x)
    for(i in 1:nchain){
      inits[[i]] <- list(a=-30,b=rnorm(1,0.10,0.015),c=rnorm(1,0.05,0.01),d=rnorm(1,0.33,0.03))
    }
  }
  else if(data.source == "MODIS.NDVI"){
    data = MODIS_data(Lat,Long,data.source,site.name=site.name)
    for(i in 1:(nchain)){
      inits[[i]] <- list(a=rnorm(1,-30,3),b=rnorm(1,0.11,0.05),c=0.2,d=0.7)
    }
  }
  else if(data.source=="GOES.NDVI"){
    data = GOES_data(site.name)
    for(i in 1:(nchain)){
      inits[[i]] <- list(a=rnorm(1,-29.9,0.3),b=rnorm(1,0.25,0.1),c=rnorm(1,0.4,0.1),d=rnorm(1,0.4,0.1))
      #inits[[nchain]] <- list(a=-29.9,b=0.25,c=0.40,d=0.40)
    }
  }
  data$s1 <- 0.001#0.5
  data$s2 <- 0.00001#0.2
  data$v.a <- 1#0.00001
  data$v.b <- .001#0.00001
  data$mean.a <- -30
  data$mean.b <- 0.11
  
    #data$s1 <- 0.5
  #data$s2 <- 0.2
  #data$v.a <- 3
  #data$v.b <- 0.001
    #data$mean.a <- -30
  #data$mean.b <- 0.11
  
  DB_model <- "
  model{
    ##priors
    a ~ dnorm(mean.a,v.a)
    b ~ dnorm(mean.b,v.b)
    d ~ dbeta(alpha.d,beta.d)
    c ~ dbeta(alpha.c,beta.c)
    prec ~ dgamma(s1,s2)
  
    for(i in 1:n){
    mu[i] <- c/(1 + exp(a+b*x[i]))+d   	## process model
    y[i]  ~ dnorm(mu[i],prec)		## data model (will need to change to beta eventually)
    }
  }
  "
   j.model   <- jags.model(file = textConnection(DB_model),
                           data = data,
                           inits=inits,
                           n.chains = nchain)
   var.out   <- coda.samples (model = j.model,
                               variable.names = c("a","b","c","d","prec"),
                               n.iter = niter)
  output <- list()
  output$var.out <- var.out
  output$x <- data$x
  output$y <- data$y
  return(output)
}

##' For GOES NDVI data, construct the data object for input into MCMC
##' 
##' @param site.name Site Name
GOES_data <- function(site.name) {
  ##Data
  NDVI.fileName1 <- paste("GOES_NDVI_",site.name,"2017_kappaDQF2.csv",sep="")
  GOES1 <- read.csv(NDVI.fileName1,header=FALSE) 
  NDVI.fileName2 <- paste("GOES_NDVI_",site.name,"2018_kappaDQF.csv",sep="")
  GOES2 <- read.csv(NDVI.fileName2,header=FALSE)
  GOES2[1,] <- GOES2[1,]+365
  GOES <- cbind(GOES1,GOES2)
  
  GOES_Days <- as.numeric(GOES[1,])
  GOES_NDVI <- as.numeric(GOES[2,])
  
  sep.val <- min(which(GOES_Days>182)) 
  y <- GOES_NDVI[sep.val:length(GOES_Days)]
  x <- GOES_Days[sep.val:length(GOES_Days)]
  
  data <- list(x=x,y=y,n=length(y))
  ##Specify Priors
  data$alpha.c <- 4
  data$beta.c <- 3
  data$alpha.d <- 2
  data$beta.d <- 3
  data$s1 <- 0.001#0.5
  data$s2 <- 0.00001#0.2
  data$v.a <- 1#0.00001
  data$v.b <- 1#0.00001
  data$mean.a <- -30
  data$mean.b <- 0.11
  
  #data$s1 <- 0.5
  #data$s2 <- 0.2
  #data$v.a <- 3
  #data$v.b <- 2
  #data$mean.a <- -30
  #data$mean.b <- 0.4

  return(data)
}

##' For PhenoCam data, construct the data object for input into MCMC
##' 
##' @param URL PhenoCam network URL
PC_data <- function(URL) {
  ##Data
  PC.data <- subset(download.phenocam(URL),year%in%c(2017,2018))
  PC.data <- PC.data[1:425,]
  PC.time = as.Date(PC.data$date)
  y <- PC.data$gcc_mean[185:425]
  x <- lubridate::yday(PC.time[185:425])
  for(i in 1:length(x)){
    if(x[i]<100){
      x[i] <- x[i]+365
    }
  }
  data <- list(x=x,y=y,n=length(y))
  ##Specify Priors
  data$beta.c <- 5
  data$alpha.c <- 1
  data$alpha.d <- 3
  data$beta.d <- 7
  data$s1 <- 0.5
  data$s2 <- 0.2
  data$v.a <- 3
  data$v.b <- 0.001
  data$mean.a <- -30
  data$mean.b <- 0.11
  return(data)
}

##' For MODIS EVI data, construct the data object for input into MCMC
##' 
##' @param Lat  latitude of desired site in decimals
##' @param Long longitude of desired site in decimals
##' @param data.source data source (GOES.NDVI, MODIS.NDVI, MODIS.EVI, PC.GCC)
MODIS_data <- function(Lat,Long,data.source,site.name) {
  ##Data
  #site.name <- "HarvardForest"
  fileName <- paste(site.name,"_MODIS_NDVI2.csv",sep="")
  #options(scipen=999)
  #options(scipen=0)
  MODIS = read.csv(fileName,header=FALSE)
  y <- MODIS[,7]
  x <- as.integer(MODIS[,5])
  for(i in 1:length(x)){
    if(x[i]<100){
      x[i] <- as.numeric(x[i]) + 365
    }
  }

  data <- list(x=x,y=y,n=length(y))
  ##Specify Priors
  data$alpha.c <- 1
  data$beta.c <- 5
  data$alpha.d <- 3.5
  data$beta.d <- 5
  data$s1 <- 0.001#0.5
  data$s2 <- 0.00001#0.2
  data$v.a <- 0.1#0.00001
  data$v.b <- 1#0.00001
  data$mean.a <- -30
  data$mean.b <- 0.11

  return(data)
}

pheno.logistic <- function(a,b,c,d,xseq){
  return(c/(1 + exp(a+b*xseq))+d)
}

K.prime <- function(a,b,c,d,t){
  z <- exp(a+b*t)
  num1 <- 3*z*(1-z)*(1+z)**3*(2*(1+z)**3+b**2*c**2*z)
  den1 <- ((1+z)**4+(b*c*z)**2)**(5/2)
  num2 <- (1+z)**2*(1+2*z-5*z**2)
  den2 <- ((1+z)**4+(b*c*z)**2)**(3/2)
  return(b**3*c*z*((num1/den1)-num2/den2))
}

```

```{r}
siteData <- read.csv("GE509_Project_Sites.csv",header=FALSE)
results.table <- data.frame(matrix(nrow=nrow(siteData)*3,ncol=7))
colnames(results.table) <- c("Burnin","Eff. Sample Size","a+-SE","b+-SE","c+-SE","d+-SE","precision+-SE")
#rownames(results.table) <- c("HarvForest GOES","HarvForest MODIS","HarvForest PhenoCam","Bartlett GOES","Bartlett MODIS","Bartlett PhenoCam","missouriozarks GOES","missouriozarks MODIS","missouriozarks PhenoCam","willowCreek GOES","willowCreek MODIS","willowCreek PhenoCam")
PC.trans.dates <- matrix(nrow=nrow(siteData),ncol=3)
colnames(PC.trans.dates) <- c("Trans1","Trans2","Trans3")
MODIS.trans.dates <- matrix(nrow=nrow(siteData),ncol=3)
colnames(MODIS.trans.dates) <- c("Trans1","Trans2","Trans3")
GOES.trans.dates <- matrix(nrow=nrow(siteData),ncol=3)
colnames(GOES.trans.dates) <- c("Trans1","Trans2","Trans3")
nrows <- 1
iseq <- c(1,3,4,7)
for (i in iseq){
  i <- 1
  siteName <- as.character(siteData[i,1])
  Lat <- as.character(siteData[i,2])
  Long <-as.character(siteData[i,3])
  URL <- as.character(siteData[i,4])
  
#GOES
  # if(i==2){
  #   out.GOES <- createBayesModel.DB(data.source="GOES.NDVI",site.name =  siteName,niter=400000)
  # }
  # else{
    out.GOES <- createBayesModel.DB(data.source="GOES.NDVI",site.name =  siteName)
  #}
#sum.GOES <- summary(out.GOES$var.out)
#sum.GOES
#plot(out.GOES$var.out)
gelman.diag(out.GOES$var.out)
GBR <- gelman.plot(out.GOES$var.out)
burnin <- GBR$last.iter[tail(which(apply(GBR$shrink[,,2]>1.05,1,any)),1)+1]
results.table[nrows,1] <-  burnin
results.table[nrows,2] <-  3*(20000-burnin)
if(length(burnin) == 0) burnin = 1
var.burn <- window(((out.GOES$var.out)),start=burnin)
burn.sum <- summary(var.burn)

#howland.GOES.varBurn <- var.burn
fileName <- paste(siteName,"_GOES_varBurn.RData",sep="")
save(howland.GOES.varBurn,file=fileName)


var.GOES.a <- burn.sum$statistics[1,1]
var.GOES.b <- burn.sum$statistics[2,1]
var.GOES.c <- burn.sum$statistics[3,1]
var.GOES.d <- burn.sum$statistics[4,1]

results.table[nrows,3] <-paste(as.character(burn.sum$statistics[1,1]),"+-",as.character(burn.sum$statistics[1,2]))
results.table[nrows,4] <- paste(as.character(burn.sum$statistics[2,1]),"+-",as.character(burn.sum$statistics[2,2]))
results.table[nrows,5] <- paste(as.character(burn.sum$statistics[3,1]),"+-",as.character(burn.sum$statistics[3,2]))
results.table[nrows,6] <-paste(as.character(burn.sum$statistics[4,1]),"+-",as.character(burn.sum$statistics[4,2]))
results.table[nrows,7] <-paste(as.character(burn.sum$statistics[5,1]),"+-",as.character(burn.sum$statistics[5,2]))

x.GOES <- out.GOES$x
var.GOES.pred <- pheno.logistic(var.GOES.a,var.GOES.b,var.GOES.c,var.GOES.d,xseq=x.GOES)
var.mat<-as.matrix(var.burn)

a<-var.mat[,1]
b<-var.mat[,2]
c <- var.mat[,3]
d <- var.mat[,4]
xpred <- seq(182,425, by =1)
ycred <- matrix(0,nrow=10000,ncol=244)

for(g in 1:10000){
  Ey <- (c[g]/(1 + exp(a[g]+b[g]*xpred))+d[g])
  ycred[g,] <- Ey
}
ci <- apply(ycred,2,quantile,c(0.025,0.5, 0.975), na.rm= TRUE)
xpred.GOES <- xpred
ci.GOES <- ci

xseq <- seq(182,365)
K.prime.vals.mean <- K.prime(var.GOES.a,var.GOES.b,var.GOES.c,var.GOES.d,xseq)

max.index <- which.max(K.prime.vals.mean)
GOES.trans.dates[i,2] <- xseq[max.index]
GOES.trans.dates[i,1] <- xseq[which.min(K.prime.vals.mean[1:max.index])]
GOES.trans.dates[i,3] <- xseq[which.min(K.prime.vals.mean)]
nrows <- nrows + 1
#MODIS:
if(i==10){
  out.MODIS <- createBayesModel.DB(Lat = Lat, Long=Long,data.source = "MODIS.NDVI",site.name=siteName)
}
else{
  out.MODIS <- createBayesModel.DB(Lat = Lat, Long=Long,data.source = "MODIS.NDVI",site.name=siteName,niter=300000)
}

#plot(out.MODIS$var.out)
#summary(out.MODIS$var.out)

gelman.diag(out.MODIS$var.out)
GBR <- gelman.plot(out.MODIS$var.out)
burnin <- GBR$last.iter[tail(which(apply(GBR$shrink[,,2]>1.05,1,any)),1)+1]
results.table[nrows,1] <-  burnin
results.table[nrows,2] <-  3*(20000-burnin)
if(length(burnin) == 0) burnin = 1
var.burn <- window(((out.MODIS$var.out)),start=burnin)
burn.sum <- summary(var.burn)
#burn.sum

howland.MODIS.varBurn <- var.burn
fileName <- paste(siteName,"_MODIS_varBurn.RData",sep="")
save(howland.MODIS.varBurn,file=fileName)

var.a <- burn.sum$statistics[1,1]
var.b <- burn.sum$statistics[2,1]
var.c <- burn.sum$statistics[3,1]
var.d <- burn.sum$statistics[4,1]

results.table[nrows,3] <-paste(as.character(burn.sum$statistics[1,1]),"+-",as.character(burn.sum$statistics[1,2]))
results.table[nrows,4] <- paste(as.character(burn.sum$statistics[2,1]),"+-",as.character(burn.sum$statistics[2,2]))
results.table[nrows,5] <- paste(as.character(burn.sum$statistics[3,1]),"+-",as.character(burn.sum$statistics[3,2]))
results.table[nrows,6] <-paste(as.character(burn.sum$statistics[4,1]),"+-",as.character(burn.sum$statistics[4,2]))
results.table[nrows,7] <-paste(as.character(burn.sum$statistics[5,1]),"+-",as.character(burn.sum$statistics[5,2]))

x.MODIS <- out.MODIS$x
var.MODIS.pred <- pheno.logistic(var.a,var.b,var.c,var.d,xseq=as.numeric(x.MODIS))

var.mat<-as.matrix(var.burn)

a<-var.mat[,1]
b<-var.mat[,2]
c <- var.mat[,3]
d <- var.mat[,4]
xpred <- seq(x.MODIS[1],x.MODIS[length(x.MODIS)], by =1)
ycred <- matrix(0,nrow=10000,ncol=length(xpred))

for(g in 1:10000){
  Ey <- (c[g]/(1 + exp(a[g]+b[g]*xpred))+d[g])
  ycred[g,] <- Ey
}
ci.MODIS <- apply(ycred,2,quantile,c(0.025,0.5, 0.975), na.rm= TRUE)
xpred.MODIS <- xpred

xseq <- seq(182,365)
K.prime.vals.mean <- K.prime(var.a,var.b,var.c,var.d,xseq)

max.index <- which.max(K.prime.vals.mean)
MODIS.trans.dates[i,2] <- xseq[max.index]
MODIS.trans.dates[i,1] <- xseq[which.min(K.prime.vals.mean[1:max.index])]
MODIS.trans.dates[i,3] <- xseq[which.min(K.prime.vals.mean)]

nrows <- nrows+1
#PC:
out.PC <- createBayesModel.DB(data.source="PC.GCC",URL=URL)#,niter=400000)
#plot(out.PC$var.out)
#summary(out.PC$var.out)

gelman.diag(out.PC$var.out)
GBR <- gelman.plot(out.PC$var.out)
burnin <- GBR$last.iter[tail(which(apply(GBR$shrink[,,2]>1.05,1,any)),1)+1]
burnin
results.table[nrows,1] <-  burnin
results.table[nrows,2] <-  3*(20000-burnin)

if(length(burnin) == 0) burnin = 1
var.burn <- window(((out.PC$var.out)),start=burnin)
burn.sum <- summary(var.burn)

Coweeta.PC.varBurn <- var.burn
fileName <- paste(siteName,"_PC_varBurn.RData",sep="")
save(Coweeta.PC.varBurn,file=fileName)

var.a <- burn.sum$statistics[1,1]
var.b <- burn.sum$statistics[2,1]
var.c <- burn.sum$statistics[3,1]
var.d <- burn.sum$statistics[4,1]

results.table[nrows,3] <-paste(as.character(burn.sum$statistics[1,1]),"+-",as.character(burn.sum$statistics[1,2]))
results.table[nrows,4] <- paste(as.character(burn.sum$statistics[2,1]),"+-",as.character(burn.sum$statistics[2,2]))
results.table[nrows,5] <- paste(as.character(burn.sum$statistics[3,1]),"+-",as.character(burn.sum$statistics[3,2]))
results.table[nrows,6] <-paste(as.character(burn.sum$statistics[4,1]),"+-",as.character(burn.sum$statistics[4,2]))
results.table[nrows,7] <-paste(as.character(burn.sum$statistics[5,1]),"+-",as.character(burn.sum$statistics[5,2]))

x.PC <- out.PC$x

var.PC.pred <- pheno.logistic(var.a,var.b,var.c,var.d,xseq=x.PC)
var.mat<-as.matrix(var.burn)

a<-var.mat[,1]
b<-var.mat[,2]
c <- var.mat[,3]
d <- var.mat[,4]
xpred <- seq(x.PC[1],x.PC[length(x.PC)], by =1)
ycred <- matrix(0,nrow=10000,ncol=length(xpred))

for(g in 1:10000){
  Ey <- (c[g]/(1 + exp(a[g]+b[g]*xpred))+d[g])
  ycred[g,] <- Ey
}
ci.PC <- apply(ycred,2,quantile,c(0.025,0.5, 0.975), na.rm= TRUE)
xpred.PC <- xpred

xseq <- seq(182,365)
K.prime.vals.mean <- K.prime(var.a,var.b,var.c,var.d,xseq)

max.index <- which.max(K.prime.vals.mean)
PC.trans.dates[i,2] <- xseq[max.index]
PC.trans.dates[i,1] <- xseq[which.min(K.prime.vals.mean[1:max.index])]
PC.trans.dates[i,3] <- xseq[which.min(K.prime.vals.mean)]

plot(x=list(),y=list(),xlim=c(160,420),ylim=c(0,1),ylab="Value",xlab="Day of Year",main=siteName)
lines(x.GOES,var.GOES.pred,col="black")
lines(xpred.GOES,ci.GOES[1,], col="black", lty = 2)
lines(xpred.GOES,ci.GOES[3,],col="black",lty=2)
points(out.GOES$x,out.GOES$y,col="Black")
lines(x.PC,var.PC.pred,col="green")
lines(xpred,ci.PC[1,], col="green", lty = 2)
lines(xpred,ci.PC[3,],col="green",lty=2)
points(out.PC$x,out.PC$y,col="green")
lines(x.MODIS,var.MODIS.pred,col="red")
lines(xpred.MODIS,ci.MODIS[1,], col="red", lty = 2)
lines(xpred.MODIS,ci.MODIS[3,],col="red",lty=2)
points(out.MODIS$x,out.MODIS$y,col="red")
nrows <- nrows+1
}

```


Fig. 1-4: Plots showing logistic curve fits (solid line) of data points (open circles) with 95% credible intervals (dashed line). Green indicates PhenoCam GCC, red indicates MODIS NDVI, and black indicates GOES NDVI. 



```{r}
results.table
```


```{r}
plot(PC.trans.dates[,1],GOES.trans.dates[,1],xlim=c(200,350),ylim=c(200,350),col="Black",ylab="GOES Date",xlab="PhenoCam Date")
abline(a=0,b=1)
points(PC.trans.dates[,2],GOES.trans.dates[,2],col="Green")
points(PC.trans.dates[,3],GOES.trans.dates[,3],col="Red")
```

Fig 5. Comparison between the phenological transition dates from GOES data vs PhenoCam data for the four different sites. Black circles indicate the start of senescence, green indicates the senescence inflection point, and red indicates the end of senescence. The black line denotes a 1:1 relationship. 


```{r}
plot(PC.trans.dates[,1],MODIS.trans.dates[,1],xlim=c(200,350),ylim=c(200,350),col="black",ylab="MODIS Date",xlab="PhenoCam Date")
abline(a=0,b=1)
points(PC.trans.dates[,2],MODIS.trans.dates[,2],col="Green")
points(PC.trans.dates[,3],MODIS.trans.dates[,3],col="Red")
```

Fig 6. Comparison between the phenological transition dates from MODIS data vs PhenoCam data for the four different sites. Black circles indicate the start of senescence, green indicates the senescence inflection point, and red indicates the end of senescence. The black line denotes a 1:1 relationship. 


```{r}
plot(GOES.trans.dates[,1],MODIS.trans.dates[,1],xlim=c(200,350),ylim=c(200,350),col="black",ylab="MODIS Date",xlab="GOES Date")
abline(a=0,b=1)
points(GOES.trans.dates[,2],MODIS.trans.dates[,2],col="green")
points(GOES.trans.dates[,3],MODIS.trans.dates[,3],col="red")
```

Fig 7. Comparison between the phenological transition dates from MODIS data vs GOES data for the four different sites. Black circles indicate the start of senescence, green indicates the senescence inflection point, and red indicates the end of senescence. The black line denotes a 1:1 relationship. 

_Results:_

Logistic curves were fitted for autumn phenology changes in four broadleaf deciduous forests based on three remote sensing sources: MODIS NDVI, GOES NDVI, and PhenoCam GCC. Burnin values, effective sample size, and mean parameter values for all sites and data sources can be found in Table 1. For all sites, the credible interval for PhenoCam was narrowest followed by GOES, with MODIS resulting in the widest credible interval (Figures 1-4). Reasonable logistic curves were fit for both PhenoCam and GOES, but the MODIS NDVI fit does not appear reasonable for all sites. GOES predicted a later date for the start of senescence and the inflection point of senescence than PhenoCam for all sites and predicted an earlier end of senescence for all sites except for one (Figure 5). For the majority of transition dates, PhenoCam predicted later dates than MODIS, except for one senescence inflection date and one end of senescence date (Figure 6).  Likewise, GOES predicted later dates for all transitions, except for one senescence inflection date and one start senescence date (Figure 7).