###################################################################################################################
# This R script extracts climate time series (e.g. Tmean, Tmax, VPD, PPT, etc.) for each site. Climate data derived from PRISM. 
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################
rm(list=ls())
require(maptools)
require(raster)
require(lattice)
require(reshape2)
setwd('')
wgs84 <- "+proj=longlat +datum=WGS84"

#-------------------------------------------------------------
# load plot locations and prism climate data
#-------------------------------------------------------------
#xy.cords <- read.csv('C:/Users/Logan/Google Drive/research/cascade_tree_traits/data_backup/plot_details/plot_coordinates_center.csv')
xy.cords <- read.csv('projects/cascade_conifers/data/plot_coordinates_center.csv')
xy.cords$plot.id <- paste(xy.cords$strata, xy.cords$plot, sep = '.')
pts.aea <- SpatialPointsDataFrame(coords =xy.cords[,c(6,7)], xy.cords)
pts.wgs84 <- SpatialPointsDataFrame(coords =xy.cords[,c(4,3)], xy.cords)

tmean.files <- list.files('data/climate/prism/tmean_4km/wus_monthly/', pattern = glob2rx('*.img'), full.names = T)
tmax.files <- list.files('data/climate/prism/tmax_4km/wus_monthly/', pattern = glob2rx('*.img'), full.names = T)
tmin.files <- list.files('data/climate/prism/tmin_4km/wus_monthly/', pattern = glob2rx('*.img'), full.names = T)
ppt.files <- list.files('data/climate/prism/ppt_4km/wus_monthly/', pattern = glob2rx('*.img'), full.names = T)
eto.files <- list.files('data/climate/prism/eto_har_4km/wus_monthly/', pattern = glob2rx('*.img'), full.names = T)
cmi.files <- list.files('data/climate/prism/cmi_har_4km/wus_monthly/', pattern = glob2rx('*.img'), full.names = T)
defecit.files <- list.files('data/climate/prism/defecit_har_4km/wus_monthly/', pattern = glob2rx('*.img'), full.names = T)
vpd.files <- list.files('data/climate/prism/vpd_4km/wus_monthly/', pattern = glob2rx('*.img'), full.names = T)
freeze.files <- list.files('data/climate/prism/freeze_4km/wus_monthly/', pattern = glob2rx('*.img'), full.names = T)
soilmoist.files <- list.files('data/vic_hydrology/wus_wgs84_monthly_avg/soilmoist/', pattern = glob2rx('*.img'), full.names = T)

#-------------------------------------------------------------
# establish time periods and sites
#-------------------------------------------------------------
prism.yrs <- 1895:2013
yoi <- 1955:2013
n.yoi <- length(yoi)
months <- c('01','02','03','04','05','06','07','08','09','10','11','12')
yoi.ts <- paste(sort(rep(yoi, 12)), rep(months, n.yoi), sep = '.')

sites <- pts.aea$plot.id
n.sites <- length(sites)
clim.pars <- c('tmean','tmax','tmin','ppt','eto','cmi','defecit','vpd','freeze','soilmoist')
n.clim.pars <- length(clim.pars)

#-------------------------------------------------------------
# build data frame into which climate data will be extracted
#-------------------------------------------------------------
clim.df <- data.frame(matrix(NA, nrow=0, ncol=7))
colnames(clim.df) <- c('strata','plot','plot.id','clim.par','year','month','value')

for (i in 1:n.clim.pars){
  if (clim.pars[i]=='soilmoist'){
    files <- get(paste(clim.pars[i],'.files',sep=''))
    files.yoi <- files[which(sort(rep(1950:2013,12))%in%yoi)]
    clim.par.i <- stack(files.yoi)
    clim.vals <- as.data.frame(t(raster::extract(clim.par.i, pts.wgs84)))
  } else {
    files <- get(paste(clim.pars[i],'.files',sep=''))
    files.yoi <- files[which(sort(rep(prism.yrs,12))%in%yoi)]
    clim.par.i <- stack(files.yoi)
    clim.vals <- as.data.frame(t(raster::extract(clim.par.i, pts.aea)))
  }  
  
  colnames(clim.vals) <- sites
  clim.vals$timeseries <- yoi.ts
  clim.long <- melt(clim.vals, variable.name="plot.id", value.name=c("values"), id='timeseries', na.rm = F)
  clim.long$clim.par <- rep(clim.pars[i], nrow(clim.long))
  clim.long$year <- as.numeric(substr(clim.long$timeseries, 1,4))
  clim.long$month <- as.numeric(substr(clim.long$timeseries, 6,7)) ## FIX ERROR WITH EXTRACTING THE MONTH!!!!!!!!
  clim.long$strata <- substr(clim.long$plot.id, 1,2)
  clim.long$plot <- as.numeric(substr(clim.long$plot.id, 4,4))
  clim.long <- clim.long[,c(7,8,2,4,5,6,3)] # reorder columns
  clim.df <- rbind(clim.df, clim.long)
  print(c('finished', clim.pars[i]))
}

write.table(clim.df, 'projects/cascade_conifers/output/site_prism_data.csv', sep = ',', row.names = F, col.names = T)
