###########################################################################################
##This code reads in deployment data for chamber flux measurements and calculates CO2 and CH4 fluxes,
##requires environmental data in the read file, and all LGR files in the directory in their native format
#JTC, April 18, 2015
rm(list=ls(all=TRUE))
########################################################################################
#read LGR deployment data 
setwd("C:/Users/Nick/Dropbox/FlamerFlux/run5")
metadata=read.csv("run5_metadata_2015.csv", stringsAsFactors = FALSE)
metadata$LGR.file=as.character(metadata$LGR.file)
#metadata$Date_Time=as.POSIXct(metadata$Date_Time, format="%m/%d/%Y %H:%M:%S")
metadata$start1 <- paste0(metadata$date, " ", metadata$start1)
metadata$start2 <- paste0(metadata$date, " ", metadata$start2)
metadata$start3 <- paste0(metadata$date, " ", metadata$start3)
metadata$stop1 <- paste0(metadata$date, " ", metadata$stop1)
metadata$stop2 <- paste0(metadata$date, " ", metadata$stop2)
metadata$stop3 <- paste0(metadata$date, " ", metadata$stop3)              
metadata$start1=as.POSIXct(metadata$start1, format="%m/%d/%Y %H:%M:%S",tz="America/Chicago")
metadata$stop1=as.POSIXct(metadata$stop1, format="%m/%d/%Y %H:%M:%S",tz="America/Chicago")
metadata$start2=as.POSIXct(metadata$start2, format="%m/%d/%Y %H:%M:%S",tz="America/Chicago")
metadata$stop2=as.POSIXct(metadata$stop2, format="%m/%d/%Y %H:%M:%S",tz="America/Chicago")
metadata$start3=as.POSIXct(metadata$start3, format="%m/%d/%Y %H:%M:%S",tz="America/Chicago")
metadata$stop3=as.POSIXct(metadata$stop3, format="%m/%d/%Y %H:%M:%S",tz="America/Chicago")
metadata$LGR_file=as.character(metadata$LGR.file)
str(metadata)
##############################################################################################################
flux.co2.calc=rep(NA, nrow(metadata))
flux.ch4.calc=rep(NA, nrow(metadata))
mean.r2.co2=rep(NA, nrow(metadata))
mean.r2.ch4=rep(NA, nrow(metadata))
for(i in 1:nrow(metadata)){
  #i=55
  #set environmental parameters
  AirTemp= 273.15+metadata$AirTemp_C[i];#Kelvin
  Pressure= 0.000986923266716*metadata$Pressure_hPa[i]#atmospheres
  GasConstant= 44.617516*Pressure; #corrected gas constant in mol/m3
  ChamberHt= metadata$chamberht;#meters
  AtmosphereCO2=metadata$AtmCO2_ppm[i] #ppmv
  AtmosphereCH4=metadata$AtmCH4_LGR[i] #ppmv
  #read file from LGR
  site.name=metadata$Site[i]
  file.name=metadata$LGR.file[i]
  laser=read.csv(file.name, skip=1) #,comment.char="-"
  laser$Time=as.POSIXct(laser$Time, tz="America/Chicago",  format="%m/%d/%y %H:%M:%S")
  laser=laser[!duplicated(laser[,"Time"]),]#remove duplicated timestamps at the second level
  laser$second=seq(from=1, to=nrow(laser), by=1)# add a string of seconds to the dataset
  ################select flux replicates by timestamp
  start1=metadata$start1[i]
  stop1=metadata$stop1[i]
  time1=subset(laser, Time > as.POSIXct(start1,tz="America/Chicago") & Time <= as.POSIXct(stop1,tz="America/Chicago"))
  ##########
  start2=metadata$start2[i]
  stop2=metadata$stop2[i]
  time2=subset(laser, Time > as.POSIXct(start2,tz="America/Chicago") & Time <= as.POSIXct(stop2,tz="America/Chicago"))
  ##########
  start3=metadata$start3[i]
  stop3=metadata$stop3[i]
  time3=subset(laser, Time > as.POSIXct(start3,tz="America/Chicago") & Time <= as.POSIXct(stop3,tz="America/Chicago"))
  ###########################################################################
  #plot/model the slopes
  graphics.off()
  dev.new()
  pdf(file=paste("Site", site.name, file.name,".pdf", sep="_"), height=5, width=5)
  #dev.new(paste("Site",site.name, file.name,".pdf", sep="_"))
  par(mfrow=c(2,3),mar=c(4,4,1,1) + 0.1,oma=c(4,3,1,1))
  CO2con1<-time1$X.CO2._ppm
  plot(time1$second, CO2con1)
  lm1=lm(CO2con1~time1$second)
  summary(lm1)
  abline(lm1, col="red", lwd=2)
  slope1=NA
  slope1=summary(lm1)[[4]][[2]]
  r1=NA
  r1 = summary(lm1)$r.squared
  ##
  CO2con2<-time2$X.CO2._ppm
  plot(time2$second, CO2con2)
  lm2=lm(CO2con2~time2$second)
  summary(lm2)
  abline(lm2, col="red", lwd=2)
  slope2=NA
  slope2=summary(lm2)[[4]][[2]]
  r2=NA
  r2 = summary(lm2)$r.squared
  ##
  CO2con3<-time3$X.CO2._ppm
  plot(time3$second, CO2con3)
  lm3=lm(CO2con3~time3$second)
  summary(lm3)
  abline(lm3, col="red", lwd=2)
  slope3=NA
  slope3=summary(lm3)[[4]][[2]]
  r3=NA
  r3 = summary(lm3)$r.squared
  r2.co2=c(r1, r2, r3)
  mean.r2.co2[i]=mean(na.omit(r2.co2))
  slopes.co2=c(slope1, slope2, slope3)
  slopemean.co2=mean(na.omit(slopes.co2))#
  ############################################
  #repeat plotting/model of slopes for CH4
  CH4con1<-time1$X.CH4._ppm
  plot(time1$second, CH4con1)
  lm1b=lm(CH4con1~time1$second)
  summary(lm1b)
  abline(lm1b, col="red", lwd=2)
  slope1b=NA
  slope1b=summary(lm1b)[[4]][[2]]
  r1b=NA
  r1b = summary(lm1b)$r.squared
  ##
  CH4con2<-time2$X.CH4._ppm
  plot(time2$second, CH4con2)
  lm2b=lm(CH4con2~time2$second)
  summary(lm2b)
  abline(lm2b, col="red", lwd=2)
  slope2b=NA
  slope2b=summary(lm2b)[[4]][[2]]
  r2b=NA
  r2b = summary(lm2b)$r.squared
  ##
  CH4con3<-time3$X.CH4._ppm
  plot(time3$second, CH4con3)
  lm3b=lm(CH4con3~time3$second)
  summary(lm3b)
  abline(lm3b, col="red", lwd=2)
  slope3b=NA
  slope3b=summary(lm3b)[[4]][[2]]
  r3b=NA
  r3b = summary(lm3b)$r.squared
  slopes.ch4=c(slope1b, slope2b, slope3b)
  r2.ch4=c(r1b, r2b, r3b)
  mean.r2.ch4[i]=mean(na.omit(r2.ch4))
  slopemean.ch4=mean(na.omit(slopes.ch4))
  dev.off()
  #CO2 flux calculation, mol/m2/day
  source("CO2flux.r")
  flux.co2.calc[i]=CO2flux(GasConstant, AirTemp, slopemean.co2, ChamberHt)
  #CH4 flux calculation, mol/m2/day
  source("CH4flux.r")
  flux.ch4.calc[i]=CH4flux(GasConstant, AirTemp, slopemean.ch4, ChamberHt)
  #######
  closeAllConnections()
 }
###############################################################################################
############################################################################
output = data.frame(metadata, flux.co2.calc, flux.ch4.calc, mean.r2.co2, mean.r2.ch4)
str(output)
write.table(output, "run5_CalculatedFluxes_2015.csv", sep=",", row.names=FALSE, col.names=TRUE)
###################################################################