GOES_Phenology_Paper.Rmd

---
title: "Utilizing New Spectral Bands on GOES-East Satellite to Monitor Phenological Events"
output: html_document
---
target journal: Remote Sensing of Environment
<center><h2>Introduction</h2></center>
Motivation: Phenology is very important for a lot of things ranging from different scales. Climate change is altering phenology, but our models don't predict the timing of phenological events very well. Uncertainties can play a large role in determining if areas are net sinks or sources for carbon and by how much. There are different phenology observing techniques ranging from satellites to human observations. Two of the major metrics to observe phenological events are NDVI and EVI calculated from spectral wavelengths. Currently the primary satellites used to observe phenology include MODIS, Sentinel-2, and Landsat. These satellites, however, have low temporal resolutions because they are sun-synchronous. There are some products and methods to combine different satellite data. The new NASA product, HLS, provides a 5-day repeated temporal resolution for selected sites. Clouds get in the way and phenology is expected to change often on shorter time scales than those provided by most satellites. The newest version of the GOES satellites, GOES-16, has a new vegetation band and thus is able to provide data to calculate NDVI. It is a geostationary satellite and thus can provide spectral data to calculate NDVI from on the minute scale.   ...etc...


<center><h2>Methods</h2></center>
1. Site Selection:



  

Figure 1: site locations? 


2. Fit Maximum likelihood Models to Data Sources

```{r}
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))
  }
}

harvLat <- 42.5378
harvLong <- -72.1715
harvURL <- "http://phenocam.sr.unh.edu/data/archive/harvard/ROI/harvard_DB_0001_1day.csv"

phenoData <- download.phenocam(harvURL)
phenoData2017 <- subset(phenoData,year==2017)
GCC.spring <- phenoData2017$gcc_mean[1:182]
GCC.autumn <- phenoData2017$gcc_mean[182:365]
PC.time = as.Date(phenoData2017$date)

EVIfilePattern <- paste("EVI_Lat",substr(as.character(harvLat),1,5),sep="")
NDVIfilePattern <- paste("NDVI_Lat",substr(as.character(harvLat),1,5),sep="")
MODIS.EVI = read.csv(list.files(pattern=EVIfilePattern)[1],header=FALSE,as.is=TRUE,na.string="-3000")
MODIS.NDVI = read.csv(list.files(pattern=NDVIfilePattern)[1],header=FALSE,as.is=TRUE,na.string="-3000")
EVI = apply(MODIS.EVI[,11:ncol(MODIS.EVI)],1,mean,na.rm=TRUE)*0.0001
NDVI = apply(MODIS.NDVI[,11:ncol(MODIS.NDVI)],1,mean,na.rm=TRUE)*0.0001
MODIS.time = as.Date(substr(MODIS.EVI[,10],1,7),format="%Y%j")

dat <- GCC.autumn
time <- lubridate::yday(PC.time[182:365])
plot(time,dat)
#dat <- NDVI[10:23]

lklog <- function(param,dat,time){
  mu <- param[3]/(1+exp(param[1]+param[2]*time))+param[4] ## process model
  -sum(dnorm(dat,mu,param[5],log=TRUE)) ## data model
}
slope = 0.1
param <- p.ic <- c(-slope*mean(time),slope,diff(range(dat)),min(dat),0.1)
mu <- param[3]/(1+exp(param[1]+param[2]*time))+param[4]
lines(time,mu)
lklog(p.ic,dat,time)

rseq <- seq(182,365)
fit <- nlm(lklog,p.ic,dat=dat,time=time)
fit
param <- fit$estimate

mu.mle <- param[3]/(1+exp(param[1]+param[2]*time))+param[4]
lines(time,mu.mle,col=2)

plot(rseq,-plogis(rseq,location=265,scale=10))


```


```{r} 
#Code to work on parameter model priors
logistic <- function(c,a,b,d,t){
  #print(t)
  den <- 1+exp(a+b*t)
  print(den)
  return (c/den + d)
}

exponential <- function(a,r,t){
  return(a*exp(r*t))
}
#EVI
times <- seq(185,365)
c <-0.3 #0.3
d <- 0.25 #0.25
b <- 0.11 #0.11 in fall; -0.11 in spring
a <- -35 #-30 in fall; 30 in spring
vals <- logistic(c,a,b,d,t=times)
plot(times,vals)
hist(rbeta(1000,2,8))

a2 <- 0.3
r <- -0.013
times2 <- seq(0,140)
vals2 <- exponential(a2,r,times2)
plot(times2,vals2)
hist(rgamma(100,5))
```

```{r}

```


4. Data fusion? (fitting one curve across all the data from the different sources at a site)

<center><h2>Results</h2></center>
Figure 2: Diurnal patterns of NDVI for different months

Figure 3: gap filling for clouds/decreasing uncertainty for GOES data

Figure 4: comparing the phenological patterns/events from the different sources

Figure 5: Something about a basic data assimilation/random walk state space model 

<center><h2>Discussion</h2></center>
Not just phenology, use it to detect storm damages and other disturbances. Maybe look at data from a few sites that got damaged by hurricanes this year

<center><h2>Conclusions</h2></center>