###################################################################################################################
# This R script injests monthly climtae timeseries for each site and derives growing-year time series, as well as 
# monthly and growing-year climatologiles.
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################
rm(list=ls())
require(plyr) # for ddply 
require(raster)
require(dplyr)
require(lattice)
setwd("")

strata <- c('WJ','PP','GF')
plots <- 1:5
subplots <- 1:4
yoi <- 1964:2013
n.yoi <- length(yoi)

clim <- read.csv('data/plot_summaries/climate_by_plot.csv', header=T)
clim <- subset(clim, year>=1963)
clim$month <- rep(1:12, nrow(clim)/12)
clim$strata <- factor(clim$strata, levels=strata)

#-------------------------------------------------------------------------------------------
# create vector of indices that need to be droped before calculating growing-year values
#-------------------------------------------------------------------------------------------
index.drop.onepar <- (c(1:9,(((n.yoi+1)*12)-2):((n.yoi+1)*12)))
index.drop.full <- c()
for (i in seq(0,nrow(clim)-612,612)){
  index.drop.full <- c(index.drop.full, index.drop.onepar+i)
}
index.drop.full
clim.gy <- clim[-index.drop.full,]
clim.gy$year.gy <- rep(sort(rep(yoi,12)),15)
head(clim.gy)

#-------------------------------------------------------------------------------------------
# Monthly climatologies by site
#-------------------------------------------------------------------------------------------
site.clim.monthly.avgs <- ddply(clim.gy, .(clim.par,strata,plot,month), plyr::summarize, plot.id=unique(plot.id),
                           clim.avg = mean(values), clim.sd = sd(values))
write.table(site.clim.monthly.avgs, 'outputs/climate_monthly_avgs_by_site.csv', sep = ',', col.names = T, row.names = F)

#-------------------------------------------------------------------------------------------
# Monthly climatologies by strata
#-------------------------------------------------------------------------------------------
strata.clim.monthly.avgs <- ddply(clim.gy, .(clim.par,strata,month), plyr::summarize, plot.id=unique(strata),
                                clim.avg = mean(values), clim.sd = sd(values))
write.table(strata.clim.monthly.avgs, 'outputs/climate_monthly_avgs_by_strata.csv', sep = ',', col.names = T, row.names = F)

#-------------------------------------------------------------------------------------------
# Annual climatologies by site
#-------------------------------------------------------------------------------------------
site.clim.sum.annual.avgs <- ddply(clim.gy, .(clim.par,strata,plot,year.gy), plyr::summarize, clim.sum = sum(values))
site.clim.sum.mean.annual <- ddply(site.clim.sum.annual.avgs, .(clim.par,strata,plot), plyr::summarize, clim.avg=round(mean(clim.sum),2), clim.sd=round(sd(clim.sum),2))         
write.table(site.clim.sum.mean.annual, 'outputs/climate_sum_annual_avgs_by_site.csv', sep = ',', col.names = T, row.names = F)

site.clim.mean.annual.avgs <- ddply(clim.gy, .(clim.par,plot.id,year), plyr::summarize, plot.id=unique(plot.id),strata=unique(strata), tmean.avg = mean(values))
site.clim.mean.mean.annual <- ddply(site.clim.mean.annual.avgs, .(clim.par,plot.id), plyr::summarize, plot.id=unique(plot.id), strata=unique(strata),
                                     clim.avg=round(mean(tmean.avg),2), clim.sd=round(sd(tmean.avg, na.rm = T),2))         
write.table(site.clim.mean.mean.annual, 'outputs/climate_mean_annual_avgs_by_site.csv', sep = ',', col.names = T, row.names = F)

kruskal.test(clim.avg~strata, site.clim.tmean.mean.annual)
kruskal.test(clim.avg~strata, subset(site.clim.sum.mean.annual, clim.par=='ppt'))
kruskal.test(clim.avg~strata, subset(site.clim.sum.mean.annual, clim.par=='cmi'))
kruskal.test(clim.avg~strata, subset(site.clim.sum.mean.annual, clim.par=='eto'))

#-------------------------------------------------------------------------------------------
# Annual climatologies by strata
#-------------------------------------------------------------------------------------------
strata.clim.sum.annual.avgs <- ddply(site.clim.sum.annual.avgs, .(clim.par,strata,year.gy), plyr::summarize, clim.avg = mean(clim.sum), clim.sd=sd(clim.sum))
strata.clim.sum.mean.annual <- ddply(strata.clim.sum.annual.avgs, .(clim.par,strata), plyr::summarize, gy.avg=round(mean(clim.avg)), gy.sd=round(sd(clim.avg)))         
write.table(strata.clim.sum.mean.annual, 'outputs/climate_water_annual_avgs_by_strata.csv', sep = ',', col.names = T, row.names = F)

strata.clim.mean.annual.avgs <- ddply(clim.gy, .(clim.par,strata,year), plyr::summarize, tmean.avg = mean(values))
strata.clim.mean.mean.annual <- ddply(strata.clim.mean.annual.avgs, .(clim.par,strata), plyr::summarize, clim.avg=round(mean(tmean.avg),2), clim.sd=round(sd(tmean.avg, na.rm = T),2))         
write.table(strata.clim.mean.mean.annual, 'outputs/climate_mean_annual_avgs_by_strata.csv', sep = ',', col.names = T, row.names = F)

#-------------------------------------------------------------------------------------------
# growing-year values and anomalies by site/strata
#-------------------------------------------------------------------------------------------
vars.sum <- c('cmi','deficit','eto','freeze','ppt')
vars.avg <- c('tmax','tmean','tmin','vpd')

site.water.mons <- clim.gy[clim.gy$clim.par%in%vars.sum,]
site.water.gy <- ddply(site.water.mons, .(strata,plot,clim.par,year.gy), plyr::summarize, clim.sum = sum(values))

# compute anomalies
site.water.gy <- site.water.gy %>%
  group_by(strata,plot,clim.par) %>%
  arrange(year.gy) %>%
  mutate(clim.50yr.avg = mean(clim.sum)) %>%
  mutate(clim.anom = clim.sum-clim.50yr.avg) %>%
  mutate(clim.anom.pcnt = (clim.sum-clim.50yr.avg)/clim.50yr.avg*100)

strata.water.gy <- summarise(group_by(site.water.gy,strata,clim.par,year.gy), clim.avg=mean(clim.sum),clim.sd=sd(clim.sum),
                             clim.anom.avg=mean(clim.anom), clim.anom.sd=sd(clim.anom),
                             clim.anom.pcnt.avg=mean(clim.anom.pcnt), clim.anom.pcnt.sd=sd(clim.anom.pcnt))

strata.water.gy <- strata.water.gy %>%
  mutate(clim.anom.upper=clim.anom.avg+clim.anom.sd) %>% 
  mutate(clim.anom.lower=clim.anom.avg-clim.anom.sd) %>% 
  mutate(clim.anom.pcnt.upper=clim.anom.pcnt.avg+clim.anom.pcnt.sd) %>%
  mutate(clim.anom.pcnt.lower=clim.anom.pcnt.avg-clim.anom.pcnt.sd)

write.table(strata.water.gy, 'outputs/climate_water_growingyear_timeseries_by_strata.csv', sep = ',', col.names = T, row.names = F)

# heat related variables
site.heat.mons <-  clim.gy[clim.gy$clim.par%in%vars.avg,]
site.heat.gy <- ddply(site.heat.mons, .(strata,plot,clim.par,year.gy), plyr::summarize, clim.avg = mean(values), clim.sd = sd(values))
site.heat.gy <- site.heat.gy %>%
  group_by(strata,plot,clim.par) %>%
  arrange(year.gy) %>%
  mutate(clim.50yr.avg = mean(clim.avg)) %>%
  mutate(clim.anom = clim.avg-clim.50yr.avg) %>%
  mutate(clim.anom.pcnt = (clim.avg-clim.50yr.avg)/clim.50yr.avg*100)


strata.heat.gy <- summarise(group_by(site.heat.gy,strata,clim.par,year.gy), clim.avg=mean(clim.avg),clim.sd=sd(clim.avg),
                            clim.anom.avg=mean(clim.anom), clim.anom.sd=sd(clim.anom),
                            clim.anom.pcnt.avg=mean(clim.anom.pcnt), clim.anom.pcnt.sd=sd(clim.anom.pcnt))
strata.heat.gy <- strata.heat.gy %>% 
  mutate(clim.anom.upper=clim.anom.avg+clim.anom.sd) %>% 
  mutate(clim.anom.lower=clim.anom.avg-clim.anom.sd) %>%
  mutate(clim.anom.pcnt.upper=clim.anom.pcnt.avg+clim.anom.pcnt.sd) %>%
  mutate(clim.anom.pcnt.lower=clim.anom.pcnt.avg-clim.anom.pcnt.sd)

write.table(strata.heat.gy, 'outputs/climate_heat_growingyear_timeseries_by_strata.csv', sep = ',', col.names = T, row.names = F)

#-------------------------------------------------------------------------------------------
# growing-year moving averages/sd
#-------------------------------------------------------------------------------------------
# compute moving averages/sd
site.water.gy.mavg = site.water.gy %>%
  group_by(strata,plot,clim.par) %>%
  arrange(year.gy) %>%
  mutate(clim.5yr.avg = movingFun(clim.sum,3,mean, type='to'))

site.water.gy.msd = site.water.gy %>%
  group_by(strata,plot,clim.par) %>%
  arrange(year.gy) %>%
  mutate(clim.5yr.sd = movingFun(clim.sum,3,sd, type='to'))

site.water.gy.mavg <- data.frame(site.water.gy.mavg)
site.water.gy.mavg$clim.5yr.sd <- site.water.gy.msd$clim.5yr.sd

# plot moving averages/sd
xyplot(clim.5yr.sd~year.gy|strata, subset(site.water.gy.mavg, clim.par == 'cmi'), groups = plot, layout=c(1,3),
       panel = function(x, y, ...){
         panel.xyplot(x, y, type=c('l'), cex=0.6, lty=1,lwd=2,col='black', ...)
         panel.abline(v=c(1985,1994), lty = 3, lwd=0.7)}
)

xyplot(clim.5yr.avg~year.gy|strata, subset(site.water.gy.mavg, clim.par == 'cmi'), groups = plot, layout=c(1,3),
       panel = function(x, y, ...){
         panel.xyplot(x, y, type=c('l'), cex=0.6, lty=1,lwd=2,col='black', ...)
         panel.abline(v=c(1985,1994), lty = 3, lwd=0.7)}
)

#-------------------------------------------------------------------------------------------
# growing-year moving sums
#-------------------------------------------------------------------------------------------
# compute moving averages/sd
site.water.gy.msum = site.water.gy %>%
  group_by(strata,plot,clim.par) %>%
  arrange(year.gy) %>%
  mutate(clim.5yr.sum = movingFun(clim.sum,4,sum, type='to'))

xyplot(clim.5yr.sum~year.gy|strata, subset(site.water.gy.msum, clim.par == 'cmi'), groups = plot, layout=c(1,3),
       panel = function(x, y, ...){
         panel.xyplot(x, y, type=c('l'), cex=0.6, lty=1,lwd=2,col='black', ...)
         panel.abline(v=c(1985,1994), lty = 3, lwd=0.7)}
)
