# This R script 
rm(list=ls())
require(lattice)
require(latticeExtra)
require(plyr)
#setwd("C:/research/projects/cascade_tree_traits/")
setwd("C:/Users/bernerl/Google Drive/research/cascade_tree_traits/")
#setwd("C:/Users/Logan/Google Drive/research/cascade_tree_traits/")

rw.mstr <- read.csv('data/treerings/treering_master_data.csv', header = T)
rw.mstr$plot.id <- paste(rw.mstr$strata, rw.mstr$plot, sep='.')
drop.sites <- c('wj.4','pp.4','gf.3')
rw.mstr$strata <- factor(rw.mstr$strata, levels = c('wj','pp','gf'))
rw.mstr <- rw.mstr[-which(rw.mstr$plot.id%in%drop.sites),]

# total number of cores
ddply(rw.mstr, .(), plyr::summarize, rw.mm.avg=mean(rw.mm, na.rm = T), n.trees=length(unique(id)))

# total number of cores per measurument method
ddply(rw.mstr, .(method), plyr::summarize, rw.mm.avg=mean(rw.mm, na.rm = T), n.trees=length(unique(id)))

# total number of cores by species 
ddply(rw.mstr, .(species), plyr::summarize, rw.mm.avg=mean(rw.mm, na.rm = T), n.trees=length(unique(id)))

# core count by plot
ddply(rw.mstr, .(strata,plot), plyr::summarize, rw.mm.avg=mean(rw.mm, na.rm = T), n.trees=length(unique(id)))

# core count by plot and species
ddply(rw.mstr, .(strata,plot,species), plyr::summarize, rw.mm.avg=mean(rw.mm, na.rm = T), n.trees=length(unique(id)))

# total number of trees with reliable 50-yr data 
ddply(subset(rw.mstr,reliable<=1963), .(), plyr::summarize, n.trees=length(unique(id)), rw.mm.avg=mean(rw.mm, na.rm = T))

# number of trees by species with reliable 50-yr data 
ddply(subset(rw.mstr,reliable<=1963), .(species), plyr::summarize, n.trees=length(unique(id)))

#----------------------------------------------------------------------------------------------
# growth reduction from 1977 drought
#----------------------------------------------------------------------------------------------
rw <- read.csv('data/treerings/treering_master_data.csv', header = T)
head(rw)
rw <- subset(rw, confidence == 1)
rw <- subset(rw, reliable <= 1964)
rw <- subset(rw, year >= 1964)
rw.197677 <- subset(rw, year == 1976 | year == 1977| year == 1985 |year == 1994)
rw.197677.avg <- ddply(rw.197677, .(strata,year), plyr::summarize, strata=unique(strata), year=unique(year), 
      rw.avg=round(mean(rw.mm, na.rm = T),2), rw.sd=round(sd(rw.mm,na.rm = T),2), n.trees=length(rw.mm))
rw.197677.avg$rw.se <- with(rw.197677.avg, round(rw.sd/sqrt(n.trees),2))
rw.197677.avg
((0.66-0.18)/0.66)*100
((1.65-1.08)/1.65)*100
((1.87- 1.74)/1.87)*100

((1.83-1.13)/1.83)*100
((0.66-0.38)/0.66)*100
((2.07-0.55)/2.07)*100
((2.07-1)/2.07)*100

ss <- subset(rw.mstr, plot.id=='gf.5')
unique(ss$treeid)

#----------------------------------------------------------------------------------------------
# 
#----------------------------------------------------------------------------------------------
rw.var <- ddply(rw.mstr, .(strata,plot,species,treeid), plyr::summarize, rw.var=var(rwi, na.rm = T))
rw.var.sp.avg <- ddply(rw.var, .(strata,plot), plyr::summarize, rw.var.avg=round(mean(rw.var, na.rm = T),3),
                       rw.var.sd=round(sd(rw.var, na.rm = T),3))
rw.var.sp.avg

cor.test(pl)


#----------------------------------------------------------------------------------------------
# comparison of juniper ring-widths during
#----------------------------------------------------------------------------------------------

clim <- read.csv('outputs/climate_cmigy_timeseries_by_strata.csv', header = T)
yoi <- c(1973:1977,1986:1990,2004:2013)

wj.cmi <- subset(clim, strata=='Western juniper')

wj.cmi.yoi <- wj.cmi[wj.cmi$year%in%yoi,]
wj.cmi.yoi$time.period <- c(rep(1,5),rep(2,5),rep(3,10))

t.test(cmi.gy.avg~time.period, subset(wj.cmi.yoi, time.period != 2))
t.test(cmi.gy.avg~time.period, subset(wj.cmi.yoi, time.period != 1))

xyplot(cmi.gy.avg~year.gy, wj.cmi, type='b')
