###################################################################################################################
# This R script compares estimates of PET from Hargreaves equation and PRISM data against measured ET at MP and YP flux site. 
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################
rm(list=ls())
require(lattice)
require(dplyr)
setwd('')

Me2.yoi <- 2002:2012
Me6.yoi <- 2010:2013

tower <- read.csv('cascade_tree_traits/data/eddy_tower_monthly_data_Me2Me6.csv', header =T)
head(tower)
tower$site

Me2 <- subset(tower, site == 'Me2')
Me6 <- subset(tower, site == 'Me6')

clim <- read.csv('cascade_tree_traits/data/plot_summaries/climate_by_plot.csv', header = T)
pet <- subset(clim, clim.par == 'eto')
ppt <- subset(clim, clim.par == 'ppt')

pet.Me2 <- subset(pet, plot.id == 'PP.5')
pet.Me2 <- pet.Me2[pet.Me2$year%in%Me2.yoi,]
pet.Me6 <- subset(pet, plot.id == 'PP.4')
pet.Me6 <- pet.Me6[pet.Me6$year%in%Me6.yoi,]

ppt.Me2 <- subset(ppt, plot.id == 'PP.5')
ppt.Me2 <- ppt.Me2[ppt.Me2$year%in%Me2.yoi,]
ppt.Me6 <- subset(ppt, plot.id == 'PP.4')
ppt.Me6 <- ppt.Me6[ppt.Me6$year%in%Me6.yoi,]

dim(ppt.Me2)

clim.mon.Me2 <- data.frame(year=Me2$year, et.tower=Me2$le, pet.prism=pet.Me2$values, ppt.tower=Me2$prec, ppt.prism=ppt.Me2$values)
clim.mon.Me6 <- data.frame(year=Me6$year, et.tower=Me6$le, pet.prism=pet.Me6$values, ppt.tower=Me6$prec, ppt.prism=ppt.Me6$values)
dim(clim.mon.Me2)
#clim.yrly.Me2 <- clim.mon.Me2 %>% group_by(year) %>% summarise_each(funs(mean))

# PET-monthly
plot(pet.prism~et.tower, clim.mon.Me2)
#cor.test(clim.mon.Me2$pet.prism, clim.mon.Me2$et.tower)
summary(lm(pet.prism~et.tower, clim.mon.Me2))

plot(pet.prism~et.tower, clim.mon.Me6)
#cor.test(clim.mon.Me6$pet.prism, clim.mon.Me6$et.tower)
summary(lm(pet.prism~et.tower, clim.mon.Me6))


# PPT-monthly
plot(ppt.prism~ppt.tower, clim.mon.Me2)
#cor.test(clim.mon.Me2$ppt.tower, clim.mon.Me2$ppt.prism)
summary(lm(ppt.prism~ppt.tower, clim.mon.Me2))

plot(ppt.prism~ppt.tower, clim.mon.Me6)
cor.test(clim.mon.Me6$ppt.tower, clim.mon.Me6$ppt.prism)
summary(lm(ppt.prism~ppt.tower, clim.mon.Me6))
