require(graphics)
setwd('C://Users//dweise//Documents//SERDP//RC2640//FBPs')
rm(list=ls(all=TRUE))
library(MASS);
library(openair)
library(ggplot2);library(cowplot)
library(readxl);library(tidyverse)
library(magrittr)
library(DescTools);library(circular);library(clifro)
library(RColorBrewer)
par(family="sans",font=2,cex.lab=2)
options(digits.sec=2)
#
#  Bret Butler FBP observations from raw files
#
FBPTimes<-read.csv("FJFBP_Times.csv")
colnames(FBPTimes)<-c("Plot",colnames(FBPTimes[2:6]))
FBPTimes[,3:6]<-format(FBPTimes[,3:6],format="%H:%M:%OS",tz="UTC")
#
# Plot 24Bs
#
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24Bs"&FBPTimes$logger=="2"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24Bs"&FBPTimes$logger=="2"),5],format="%H:%M:%OS",tz="UTC")
FJ24Bs_2<-read.csv("FJ24Bs_2.csv") # logger 2
colnames(FJ24Bs_2)<-c("Date",colnames(FJ24Bs_2[2:16]))
FJ24Bs_2$UTC.Time<-format(FJ24Bs_2$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ24Bs_2$logger<-"2"
tFJ24Bs_2<- FJ24Bs_2 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24Bs_5<-read.csv("FJ24Bs_5.csv") # logger 5
colnames(FJ24Bs_5)<-c("Date",colnames(FJ24Bs_5[2:16]))
FJ24Bs_5$HorzWind<-NA
FJ24Bs_5$Temp<-NA
FJ24Bs_5$UTC.Time<-format(FJ24Bs_5$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ24Bs_5$logger<-"5"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24Bs"&FBPTimes$logger=="5"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24Bs"&FBPTimes$logger=="5"),5],format="%H:%M:%OS",tz="UTC")
tFJ24Bs_5<- FJ24Bs_5 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24Bs_7<-read.csv("FJ24Bs_7.csv") # logger 7
colnames(FJ24Bs_7)<-c("Date",colnames(FJ24Bs_7[2:16]))
FJ24Bs_7$UTC.Time<-format(FJ24Bs_7$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ24Bs_7$logger<-"7"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24Bs"&FBPTimes$logger=="7"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24Bs"&FBPTimes$logger=="7"),5],format="%H:%M:%OS",tz="UTC")
tFJ24Bs_7<- FJ24Bs_7 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24Bs_12<-read.csv("FJ24Bs_12.csv") # logger 12
colnames(FJ24Bs_12)<-c("Date",colnames(FJ24Bs_12[2:16]))
FJ24Bs_12$UTC.Time<-format(FJ24Bs_12$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ24Bs_12$logger<-"12"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24Bs"&FBPTimes$logger=="12"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24Bs"&FBPTimes$logger=="12"),5],format="%H:%M:%OS",tz="UTC")
tFJ24Bs_12<- FJ24Bs_12 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24Bs_13<-read.csv("FJ24Bs_13.csv") # logger 13
FJ24Bs_13$logger<-"13"
colnames(FJ24Bs_13)<-c("Date",colnames(FJ24Bs_13[2:17]))
FJ24Bs_13$UTC.Time<-format(FJ24Bs_13$UTC.Time,format="%H:%M:%OS",tz="UTC")
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24Bs"&FBPTimes$logger=="13"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24Bs"&FBPTimes$logger=="13"),5],format="%H:%M:%OS",tz="UTC")
tFJ24Bs_13<- FJ24Bs_13 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24Bs_14<-read.csv("FJ24Bs_14.csv") # logger 14
FJ24Bs_14$logger<-"14"
colnames(FJ24Bs_14)<-c("Date",colnames(FJ24Bs_14[2:17]))
FJ24Bs_14$UTC.Time<-format(FJ24Bs_14$UTC.Time,format="%H:%M:%OS",tz="UTC")
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24Bs"&FBPTimes$logger=="14"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24Bs"&FBPTimes$logger=="14"),5],format="%H:%M:%OS",tz="UTC")
tFJ24Bs_14<- FJ24Bs_14 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
#
#  drop logger 5 since wind velocity sensor not working
#
FJ24Bs<-rbind(tFJ24Bs_2,tFJ24Bs_7,tFJ24Bs_12,tFJ24Bs_13,tFJ24Bs_14)
FJ24Bs$plot<-"24Bsquare"
FJ24Bs$angle<-circular(atan2(FJ24Bs$VertWind,FJ24Bs$HorzWind)*360/2/pi,
                       type="angles",units="degrees",zero=0,rotation="counter")
FJ24Bs$VectWind<-sqrt(FJ24Bs$VertWind^2+FJ24Bs$HorzWind^2)
#
# Plot 24Bt
#
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24Bt"&FBPTimes$logger=="6"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24Bt"&FBPTimes$logger=="6"),5],format="%H:%M:%OS",tz="UTC")
FJ24Bt_6<-read.csv("FJ24Bt_6.csv") # logger 6
colnames(FJ24Bt_6)<-c("Date",colnames(FJ24Bt_6[2:16]))
FJ24Bt_6$UTC.Time<-format(FJ24Bt_6$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ24Bt_6$logger<-"6"
tFJ24Bt_6<- FJ24Bt_6 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24Bt"&FBPTimes$logger=="11"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24Bt"&FBPTimes$logger=="11"),5],format="%H:%M:%OS",tz="UTC")
FJ24Bt_11<-read.csv("FJ24Bt_11.csv") # logger 11
colnames(FJ24Bt_11)<-c("Date",colnames(FJ24Bt_11[2:16]))
FJ24Bt_11$UTC.Time<-format(FJ24Bt_11$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ24Bt_11$logger<-"11"
tFJ24Bt_11<- FJ24Bt_11 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24Bt<-rbind(tFJ24Bt_6,tFJ24Bt_11)
FJ24Bt$plot<-"24Btriangle"
FJ24Bt$angle<-circular(atan2(FJ24Bt$VertWind,FJ24Bt$HorzWind)*360/2/pi,
                       type="angles",units="degrees",zero=0,rotation="counter")
FJ24Bt$VectWind<-sqrt(FJ24Bt$VertWind^2+FJ24Bt$HorzWind^2)
#
# Plot 24As
#
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="2"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="2"),5],format="%H:%M:%OS",tz="UTC")
FJ24As_2<-read.csv("FJ24As_2.csv") # logger 2
colnames(FJ24As_2)<-c("Date",colnames(FJ24As_2[2:16]))
FJ24As_2$UTC.Time<-format(FJ24As_2$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ24As_2$logger<-"2"
tFJ24As_2<- FJ24As_2 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24As_4<-read.csv("FJ24As_4.csv") # logger 4
colnames(FJ24As_4)<-c("Date",colnames(FJ24As_4[2:16]))
FJ24As_4$UTC.Time<-format(FJ24As_4$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ24As_4$logger<-"4"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="4"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="4"),5],format="%H:%M:%OS",tz="UTC")
tFJ24As_4<- FJ24As_4 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24As_6<-read.csv("FJ24As_6.csv") # logger 6
colnames(FJ24As_6)<-c("Date",colnames(FJ24As_6[2:16]))
FJ24As_6$UTC.Time<-format(FJ24As_6$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ24As_6$logger<-"6"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="6"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="6"),5],format="%H:%M:%OS",tz="UTC")
tFJ24As_6<- FJ24As_6 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24As_7<-read.csv("FJ24As_7.csv") # logger 7
colnames(FJ24As_7)<-c("Date",colnames(FJ24As_7[2:16]))
FJ24As_7$UTC.Time<-format(FJ24As_7$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ24As_7$logger<-"7"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="7"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="7"),5],format="%H:%M:%OS",tz="UTC")
tFJ24As_7<- FJ24As_7 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24As_11<-read.csv("FJ24As_11.csv") # logger 11
colnames(FJ24As_11)<-c("Date",colnames(FJ24As_11[2:16]))
FJ24As_11$UTC.Time<-format(FJ24As_11$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ24As_11$logger<-"11"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="11"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="11"),5],format="%H:%M:%OS",tz="UTC")
tFJ24As_11<- FJ24As_11 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24As_12<-read.csv("FJ24As_12.csv") # logger 12
colnames(FJ24As_12)<-c("Date",colnames(FJ24As_12[2:16]))
FJ24As_12$UTC.Time<-format(FJ24As_12$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ24As_12$logger<-"12"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="12"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="12"),5],format="%H:%M:%OS",tz="UTC")
tFJ24As_12<- FJ24As_12 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24As_13<-read.csv("FJ24As_13.csv") # logger 13
FJ24As_13$logger<-"13"
colnames(FJ24As_13)<-c("Date",colnames(FJ24As_13[2:17]))
FJ24As_13$UTC.Time<-format(FJ24As_13$UTC.Time,format="%H:%M:%OS",tz="UTC")
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="13"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="13"),5],format="%H:%M:%OS",tz="UTC")
tFJ24As_13<- FJ24As_13 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24As_14<-read.csv("FJ24As_14.csv") # logger 14
FJ24As_14$logger<-"14"
colnames(FJ24As_14)<-c("Date",colnames(FJ24As_14[2:17]))
FJ24As_14$UTC.Time<-format(FJ24As_14$UTC.Time,format="%H:%M:%OS",tz="UTC")
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="14"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24As"&FBPTimes$logger=="14"),5],format="%H:%M:%OS",tz="UTC")
tFJ24As_14<- FJ24As_14 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24As<-rbind(tFJ24As_2,tFJ24As_4,tFJ24As_6,tFJ24As_7,tFJ24As_11,tFJ24As_12,tFJ24As_13,tFJ24As_14)
FJ24As$plot<-"24Asquare"
FJ24As$angle<-circular(atan2(FJ24As$VertWind,FJ24As$HorzWind)*360/2/pi,
                       type="angles",units="degrees",zero=0,rotation="counter")
FJ24As$VectWind<-sqrt(FJ24As$VertWind^2+FJ24As$HorzWind^2)
#
# Plot 24At
#
FJ24At_6<-read.csv("FJ24At_6.csv") # logger 6
colnames(FJ24At_6)<-c("Date",colnames(FJ24At_6[2:16]))
FJ24At_6$UTC.Time<-format(FJ24At_6$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ24At_6$logger<-"6"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24At"&FBPTimes$logger=="6"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24At"&FBPTimes$logger=="6"),5],format="%H:%M:%OS",tz="UTC")
tFJ24At_6<- FJ24At_6 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24At_7<-read.csv("FJ24At_7.csv") # logger 7
colnames(FJ24At_7)<-c("Date",colnames(FJ24At_7[2:16]))
FJ24At_7$UTC.Time<-format(FJ24At_7$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ24At_7$logger<-"7"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24At"&FBPTimes$logger=="7"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24At"&FBPTimes$logger=="7"),5],format="%H:%M:%OS",tz="UTC")
tFJ24At_7<- FJ24At_7 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24At_13<-read.csv("FJ24At_13.csv") # logger 13
FJ24At_13$logger<-"13"
colnames(FJ24At_13)<-c("Date",colnames(FJ24At_13[2:17]))
FJ24At_13$UTC.Time<-format(FJ24At_13$UTC.Time,format="%H:%M:%OS",tz="UTC")
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24At"&FBPTimes$logger=="13"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24At"&FBPTimes$logger=="13"),5],format="%H:%M:%OS",tz="UTC")
tFJ24At_13<- FJ24At_13 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24At_14<-read.csv("FJ24At_14.csv") # logger 14
FJ24At_14$logger<-"14"
colnames(FJ24At_14)<-c("Date",colnames(FJ24At_14[2:17]))
FJ24At_14$UTC.Time<-format(FJ24At_14$UTC.Time,format="%H:%M:%OS",tz="UTC")
IntStart<-format(FBPTimes[(FBPTimes$Plot=="24At"&FBPTimes$logger=="14"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="24At"&FBPTimes$logger=="14"),5],format="%H:%M:%OS",tz="UTC")
tFJ24At_14<- FJ24At_14 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ24At<-rbind(tFJ24At_6,tFJ24At_7,tFJ24At_13,tFJ24At_14)
FJ24At$plot<-"24Atriangle"
FJ24At$angle<-circular(atan2(FJ24At$VertWind,FJ24At$HorzWind)*360/2/pi,
                       type="angles",units="degrees",zero=0,rotation="counter")
FJ24At$VectWind<-sqrt(FJ24At$VertWind^2+FJ24At$HorzWind^2)
#
# Plot 16D1
#
FJ16D1_2<-read.csv("FJ16D1_2.csv") # logger 2
colnames(FJ16D1_2)<-c("Date",colnames(FJ16D1_2[2:16]))
FJ16D1_2$UTC.Time<-format(FJ16D1_2$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ16D1_2$logger<-"2"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D1"&FBPTimes$logger=="2"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D1"&FBPTimes$logger=="2"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D1_2<- FJ16D1_2 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D1_6<-read.csv("FJ16D1_6.csv") # logger 6
colnames(FJ16D1_6)<-c("Date",colnames(FJ16D1_6[2:16]))
FJ16D1_6$UTC.Time<-format(FJ16D1_6$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ16D1_6$logger<-"6"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D1"&FBPTimes$logger=="6"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D1"&FBPTimes$logger=="6"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D1_6<- FJ16D1_6 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D1_12<-read.csv("FJ16D1_12.csv") # logger 12
colnames(FJ16D1_12)<-c("Date",colnames(FJ16D1_12[2:16]))
FJ16D1_12$UTC.Time<-format(FJ16D1_12$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ16D1_12$logger<-"12"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D1"&FBPTimes$logger=="12"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D1"&FBPTimes$logger=="12"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D1_12<- FJ16D1_12 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D1_13<-read.csv("FJ16D1_13.csv") # logger 13
FJ16D1_13$logger<-"13"
colnames(FJ16D1_13)<-c("Date",colnames(FJ16D1_13[2:17]))
FJ16D1_13$UTC.Time<-format(FJ16D1_13$UTC.Time,format="%H:%M:%OS",tz="UTC")
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D1"&FBPTimes$logger=="13"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D1"&FBPTimes$logger=="13"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D1_13<- FJ16D1_13 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D1_14<-read.csv("FJ16D1_14.csv") # logger 14
FJ16D1_14$logger<-"14"
colnames(FJ16D1_14)<-c("Date",colnames(FJ16D1_14[2:17]))
FJ16D1_14$UTC.Time<-format(FJ16D1_14$UTC.Time,format="%H:%M:%OS",tz="UTC")
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D1"&FBPTimes$logger=="14"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D1"&FBPTimes$logger=="14"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D1_14<- FJ16D1_14 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D1<-rbind(tFJ16D1_2,tFJ16D1_6,tFJ16D1_12,tFJ16D1_13,tFJ16D1_14)
FJ16D1$plot<-"16D1"
FJ16D1$angle<-circular(atan2(FJ16D1$VertWind,FJ16D1$HorzWind)*360/2/pi,
                       type="angles",units="degrees",zero=0,rotation="counter")
FJ16D1$VectWind<-sqrt(FJ16D1$VertWind^2+FJ16D1$HorzWind^2)
#
# Plot 16D5
#
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="2"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="2"),5],format="%H:%M:%OS",tz="UTC")
FJ16D5_2<-read.csv("FJ16D5_2.csv") # logger 2
colnames(FJ16D5_2)<-c("Date",colnames(FJ16D5_2[2:16]))
FJ16D5_2$UTC.Time<-format(FJ16D5_2$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ16D5_2$logger<-"2"
tFJ16D5_2<- FJ16D5_2 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D5_4<-read.csv("FJ16D5_4.csv") # logger 4
colnames(FJ16D5_4)<-c("Date",colnames(FJ16D5_4[2:16]))
FJ16D5_4$UTC.Time<-format(FJ16D5_4$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ16D5_4$logger<-"4"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="4"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="4"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D5_4<- FJ16D5_4 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D5_6<-read.csv("FJ16D5_6.csv") # logger 6
colnames(FJ16D5_6)<-c("Date",colnames(FJ16D5_6[2:16]))
FJ16D5_6$UTC.Time<-format(FJ16D5_6$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ16D5_6$logger<-"6"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="6"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="6"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D5_6<- FJ16D5_6 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D5_7<-read.csv("FJ16D5_7.csv") # logger 7
colnames(FJ16D5_7)<-c("Date",colnames(FJ16D5_7[2:16]))
FJ16D5_7$UTC.Time<-format(FJ16D5_7$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ16D5_7$logger<-"7"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="7"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="7"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D5_7<- FJ16D5_7 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D5_11<-read.csv("FJ16D5_11.csv") # logger 12
colnames(FJ16D5_11)<-c("Date",colnames(FJ16D5_11[2:16]))
FJ16D5_11$UTC.Time<-format(FJ16D5_11$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ16D5_11$logger<-"11"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="11"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="11"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D5_11<- FJ16D5_11 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D5_12<-read.csv("FJ16D5_12.csv") # logger 12
colnames(FJ16D5_12)<-c("Date",colnames(FJ16D5_12[2:16]))
FJ16D5_12$UTC.Time<-format(FJ16D5_12$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ16D5_12$logger<-"12"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="12"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="12"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D5_12<- FJ16D5_12 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D5_13<-read.csv("FJ16D5_13.csv") # logger 13
FJ16D5_13$logger<-"13"
colnames(FJ16D5_13)<-c("Date",colnames(FJ16D5_13[2:17]))
FJ16D5_13$UTC.Time<-format(FJ16D5_13$UTC.Time,format="%H:%M:%OS",tz="UTC")
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="13"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="13"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D5_13<- FJ16D5_13 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D5_14<-read.csv("FJ16D5_14.csv") # logger 14
FJ16D5_14$logger<-"14"
colnames(FJ16D5_14)<-c("Date",colnames(FJ16D5_14[2:17]))
FJ16D5_14$UTC.Time<-format(FJ16D5_14$UTC.Time,format="%H:%M:%OS",tz="UTC")
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="14"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D5"&FBPTimes$logger=="14"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D5_14<- FJ16D5_14 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D5<-rbind(tFJ16D5_2,tFJ16D5_4,tFJ16D5_6,tFJ16D5_7,tFJ16D5_11,tFJ16D5_12,tFJ16D5_13,tFJ16D5_14)
FJ16D5$plot<-"16D5"
FJ16D5$angle<-circular(atan2(FJ16D5$VertWind,FJ16D5$HorzWind)*360/2/pi,
                       type="angles",units="degrees",zero=0,rotation="counter")
FJ16D5$VectWind<-sqrt(FJ16D5$VertWind^2+FJ16D5$HorzWind^2)
#
# Plot 16D6
#
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D6"&FBPTimes$logger=="2"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D6"&FBPTimes$logger=="2"),5],format="%H:%M:%OS",tz="UTC")
FJ16D6_2<-read.csv("FJ16D6_2.csv") # logger 2
colnames(FJ16D6_2)<-c("Date",colnames(FJ16D6_2[2:16]))
FJ16D6_2$UTC.Time<-format(FJ16D6_2$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ16D6_2$logger<-"2"
tFJ16D6_2<- FJ16D6_2 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D6_6<-read.csv("FJ16D6_6.csv") # logger 6
colnames(FJ16D6_6)<-c("Date",colnames(FJ16D6_6[2:16]))
FJ16D6_6$UTC.Time<-format(FJ16D6_6$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ16D6_6$logger<-"6"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D6"&FBPTimes$logger=="6"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D6"&FBPTimes$logger=="6"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D6_6<- FJ16D6_6 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D6_12<-read.csv("FJ16D6_12.csv") # logger 12
colnames(FJ16D6_12)<-c("Date",colnames(FJ16D6_12[2:16]))
FJ16D6_12$UTC.Time<-format(FJ16D6_12$UTC.Time,format="%H:%M:%OS",tz="UTC")
FJ16D6_12$logger<-"12"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D6"&FBPTimes$logger=="12"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D6"&FBPTimes$logger=="12"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D6_12<- FJ16D6_12 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D6_13<-read.csv("FJ16D6_13.csv") # logger 13
FJ16D6_13$logger<-"13"
colnames(FJ16D6_13)<-c("Date",colnames(FJ16D6_13[2:17]))
FJ16D6_13$UTC.Time<-format(FJ16D6_13$UTC.Time,format="%H:%M:%OS",tz="UTC")
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D6"&FBPTimes$logger=="13"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D6"&FBPTimes$logger=="13"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D6_13<- FJ16D6_13 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D6_14<-read.csv("FJ16D6_14.csv") # logger 14
FJ16D6_14$logger<-"14"
colnames(FJ16D6_14)<-c("Date",colnames(FJ16D6_14[2:17]))
FJ16D6_14$UTC.Time<-format(FJ16D6_14$UTC.Time,format="%H:%M:%OS",tz="UTC")
IntStart<-format(FBPTimes[(FBPTimes$Plot=="16D6"&FBPTimes$logger=="14"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="16D6"&FBPTimes$logger=="14"),5],format="%H:%M:%OS",tz="UTC")
tFJ16D6_14<- FJ16D6_14 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
FJ16D6<-rbind(tFJ16D6_2,tFJ16D6_6,tFJ16D6_12,tFJ16D6_13,tFJ16D6_14)
FJ16D6$plot<-"16D6"
FJ16D6$angle<-circular(atan2(FJ16D6$VertWind,FJ16D6$HorzWind)*360/2/pi,
                       type="angles",units="degrees",zero=0,rotation="counter")
FJ16D6$VectWind<-sqrt(FJ16D6$VertWind^2+FJ16D6$HorzWind^2)
#
# Tall Timbers Burn 1 4/18/2018
#
IntStart<-format(FBPTimes[(FBPTimes$Plot=="TTB1"&FBPTimes$logger=="1"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="TTB1"&FBPTimes$logger=="1"),5],format="%H:%M:%OS",tz="UTC")
TTB1_1<-read.csv("TTB1_1.csv") # logger 1
colnames(TTB1_1)<-c("Date",colnames(TTB1_1[2:16]))
TTB1_1$UTC.Time<-format(TTB1_1$UTC.Time,format="%H:%M:%OS",tz="UTC")
TTB1_1$logger<-"1"
TTB1<- TTB1_1 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
TTB1$plot<-"TTB1"
TTB1$angle<-circular(atan2(TTB1$VertWind,TTB1$HorzWind)*360/2/pi,
                       type="angles",units="degrees",zero=0,rotation="counter")
TTB1$VectWind<-sqrt(TTB1$VertWind^2+TTB1$HorzWind^2)
#
# TT Burn 2 4/18/2018
#
TTB2_1<-read.csv("TTB2_1.csv") # logger 2
colnames(TTB2_1)<-c("Date",colnames(TTB2_1[2:16]))
TTB2_1$UTC.Time<-format(TTB2_1$UTC.Time,format="%H:%M:%OS",tz="UTC")
TTB2_1$logger<-"1"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="TTB2"&FBPTimes$logger=="1"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="TTB2"&FBPTimes$logger=="1"),5],format="%H:%M:%OS",tz="UTC")
tTTB2_1<- TTB2_1 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
TTB2_4<-read.csv("TTB2_4.csv") # logger 2
colnames(TTB2_4)<-c("Date",colnames(TTB2_4[2:16]))
TTB2_4$UTC.Time<-format(TTB2_4$UTC.Time,format="%H:%M:%OS",tz="UTC")
TTB2_4$logger<-"4"
IntStart<-format(FBPTimes[(FBPTimes$Plot=="TTB2"&FBPTimes$logger=="4"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="TTB2"&FBPTimes$logger=="4"),5],format="%H:%M:%OS",tz="UTC")
tTTB2_4<- TTB2_4 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
TTB2<-rbind(tTTB2_1,tTTB2_4)
TTB2$plot<-"TTB2"
TTB2$angle<-circular(atan2(TTB2$VertWind,TTB2$HorzWind)*360/2/pi,
                       type="angles",units="degrees",zero=0,rotation="counter")
TTB2$VectWind<-sqrt(TTB2$VertWind^2+TTB2$HorzWind^2)
#
# TT Burn 3 4/20/2018
#
IntStart<-format(FBPTimes[(FBPTimes$Plot=="TTB3"&FBPTimes$logger=="1"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="TTB3"&FBPTimes$logger=="1"),5],format="%H:%M:%OS",tz="UTC")
TTB3_1<-read.csv("TTB3_1.csv") # logger 1
colnames(TTB3_1)<-c("Date",colnames(TTB3_1[2:16]))
TTB3_1$UTC.Time<-format(TTB3_1$UTC.Time,format="%H:%M:%OS",tz="UTC")
TTB3_1$logger<-"1"
TTB3<- TTB3_1 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
TTB3$plot<-"TTB3"
TTB3$angle<-circular(atan2(TTB3$VertWind,TTB3$HorzWind)*360/2/pi,
                     type="angles",units="degrees",zero=0,rotation="counter")
TTB3$VectWind<-sqrt(TTB3$VertWind^2+TTB3$HorzWind^2)
#
# TT Burn 4 4/20/2018
#
IntStart<-format(FBPTimes[(FBPTimes$Plot=="TTB4"&FBPTimes$logger=="4"),4],format="%H:%M:%OS",tz="UTC")
IntEnd<-format(FBPTimes[(FBPTimes$Plot=="TTB4"&FBPTimes$logger=="4"),5],format="%H:%M:%OS",tz="UTC")
TTB4_4<-read.csv("TTB4_4.csv") # logger 4
colnames(TTB4_4)<-c("Date",colnames(TTB4_4[2:16]))
TTB4_4$UTC.Time<-format(TTB4_4$UTC.Time,format="%H:%M:%OS",tz="UTC")
TTB4_4$logger<-"4"
TTB4<- TTB4_4 %>%
  filter(UTC.Time >= IntStart & UTC.Time <= IntEnd)
TTB4$plot<-"TTB4"
TTB4$angle<-circular(atan2(TTB4$VertWind,TTB4$HorzWind)*360/2/pi,
                     type="angles",units="degrees",zero=0,rotation="counter")
TTB4$VectWind<-sqrt(TTB4$VertWind^2+TTB4$HorzWind^2)
#
#
FJ<-rbind(FJ24Bs,FJ24Bt,FJ24As,FJ24At,FJ16D1,FJ16D5,FJ16D6)
FJ$Temp<-ifelse(FJ$Temp<0,NA,FJ$Temp)
#
TT<-rbind(TTB1,TTB2,TTB3,TTB4)
TT$Temp<-ifelse(TT$Temp<0,NA,TT$Temp)
#
FJ<-rbind(FJ,TT)
#
#
#
# find max, min observations from time traces for each logger by plot
#
# then calculate harmonic mean of observations and then finally the plot mean and confidence interval
#
#
maxMT.total<-FJ %>% group_by(plot,logger) %>% slice_max(order_by=MT.total)
max.total<-maxMT.total %>% 
  filter(plot != "TTB1" & plot != "TTB3" & plot != "TTB4") %>%
  group_by(plot,logger) %>%
  summarise(total = Hmean(MT.total,method="classic", na.rm=TRUE),
                                                               rad = Hmean(MT.rad,method="classic", na.rm=TRUE),
                                                               conv = Hmean(Conv,method="classic", na.rm=TRUE),
                                                               NAR = Hmean(NAR,method="classic", na.rm=TRUE),
                                                               wind = Hmean(VectWind,method="classic", na.rm=TRUE),
                                                               angle = mean(angle)) 
maxtot<-max.total %>% summarise(maxtotal = Hmean(total,method="classic", conf.level=0.95, sides="two.sided",na.rm=TRUE),
                                                 rad = Hmean(rad,method="classic", na.rm=TRUE),
                                                 conv = Hmean(conv,method="classic", na.rm=TRUE),
                                                 NAR = Hmean(NAR,method="classic", na.rm=TRUE),
                                                 wind = Hmean(wind,method="classic", na.rm=TRUE),
                                                 angle = mean(angle,na.rm=TRUE),.groups="drop")
max.total2<-maxMT.total %>%
  filter(plot == "TTB1" | plot == "TTB3" | plot == "TTB4") %>%
  group_by(plot,logger) %>%
  summarise(maxtotal = mean(MT.total),
            rad = mean(MT.rad),
            conv = mean(Conv),
            NAR = mean(NAR),
            wind = mean(VectWind),
            angle = mean(angle),.groups="drop") %>%
  select(-logger)
maxtot<-rbind(maxtot,max.total2)
write.csv(maxtot,file="maxtot.csv")
maxMT.rad<-FJ %>% group_by(plot,logger) %>% slice_max(order_by=MT.rad)
max.rad<-maxMT.rad %>% 
  filter(plot != "TTB1" & plot != "TTB3" & plot != "TTB4") %>%
  group_by(plot,logger) %>% summarise(total = Hmean(MT.total,method="classic", na.rm=TRUE),
                                                           rad = Hmean(MT.rad,method="classic", na.rm=TRUE),
                                                           conv = Hmean(Conv,method="classic", na.rm=TRUE),
                                                           NAR = Hmean(NAR,method="classic", na.rm=TRUE),
                                                           wind = Hmean(VectWind,method="classic", na.rm=TRUE),
                                                           angle = mean(angle)) 
maxrad<-max.rad %>% summarise(total = Hmean(total,method="classic",na.rm=TRUE),
                              maxrad = Hmean(rad,method="classic", conf.level=0.95, sides="two.sided", na.rm=TRUE),
                              conv = Hmean(conv,method="classic", na.rm=TRUE),
                              NAR = Hmean(NAR,method="classic", na.rm=TRUE),
                              wind = Hmean(wind,method="classic", na.rm=TRUE),
                              angle = mean(angle),.groups="drop")
max.rad2<-maxMT.rad %>%
  filter(plot == "TTB1" | plot == "TTB3" | plot == "TTB4") %>%
  group_by(plot,logger) %>%
  summarise(total = mean(MT.total),
            maxrad = mean(MT.rad),
            conv = mean(Conv),
            NAR = mean(NAR),
            wind = mean(VectWind),
            angle = mean(angle),.groups="drop") %>%
  select(-logger)
maxrad<-rbind(maxrad,max.rad2)
write.csv(maxrad,file="maxrad.csv")

maxConv<-FJ %>% group_by(plot,logger) %>% slice_max(order_by=Conv)
max.Conv<-maxConv %>% group_by(plot,logger) %>%
  filter(plot != "TTB1" & plot != "TTB3" & plot != "TTB4") %>%
  summarise(total = Hmean(MT.total,method="classic", na.rm=TRUE),
                                                          rad = Hmean(MT.rad,method="classic", na.rm=TRUE),
                                                          conv = Hmean(Conv,method="classic", na.rm=TRUE),
                                                          NAR = Hmean(NAR,method="classic", na.rm=TRUE),
                                                          wind = Hmean(VectWind,method="classic", na.rm=TRUE),
                                                           angle = mean(angle)) 
maxCon<-max.Conv %>% summarise(total = Hmean(total,method="classic",na.rm=TRUE),
                              rad = Hmean(rad,method="classic", na.rm=TRUE),
                              maxconv = Hmean(conv,method="classic", conf.level=0.95, sides="two.sided", na.rm=TRUE),
                              NAR = Hmean(NAR,method="classic", na.rm=TRUE),
                              wind = Hmean(wind,method="classic", na.rm=TRUE),
                              angle = mean(angle),.groups="drop")
max.Con2<-maxConv %>%
  filter(plot == "TTB1" | plot == "TTB3" | plot == "TTB4") %>%
  group_by(plot,logger) %>%
  summarise(total = mean(MT.total),
            rad = mean(MT.rad),
            maxconv = mean(Conv),
            NAR = mean(NAR),
            wind = mean(VectWind),
            angle = mean(angle),.groups="drop") %>%
  select(-logger)
maxCon<-rbind(maxCon,max.Con2)
write.csv(maxCon,file="maxconv.csv")

maxNAR<-FJ %>% group_by(plot,logger) %>% slice_max(order_by=NAR)
max.NAR<-maxNAR %>% group_by(plot,logger) %>%
  filter(plot != "TTB1" & plot != "TTB3" & plot != "TTB4") %>%
  summarise(total = Hmean(MT.total,method="classic", na.rm=TRUE),
                                                        rad = Hmean(MT.rad,method="classic", na.rm=TRUE),
                                                        conv = Hmean(Conv,method="classic", na.rm=TRUE),
                                                        NAR = Hmean(NAR,method="classic", na.rm=TRUE),
                                                        wind = Hmean(VectWind,method="classic", na.rm=TRUE),
                                                        angle = mean(angle)) 
maxNarrow<-max.NAR %>% summarise(total = Hmean(total,method="classic",na.rm=TRUE),
                               rad = Hmean(rad,method="classic", na.rm=TRUE),
                               conv = Hmean(conv,method="classic", na.rm=TRUE),
                               maxNAR = Hmean(NAR,method="classic", conf.level=0.95, sides="two.sided", na.rm=TRUE),
                               wind = Hmean(wind,method="classic", na.rm=TRUE),
                               angle = mean(angle),.groups="drop")
max.NAR2<-maxNAR %>%
  filter(plot == "TTB1" | plot == "TTB3" | plot == "TTB4") %>%
  group_by(plot,logger) %>%
  summarise(total = mean(MT.total),
            rad = mean(MT.rad),
            conv = mean(Conv),
            maxNAR = mean(NAR),
            wind = mean(VectWind),
            angle = mean(angle),.groups="drop") %>%
  select(-logger)
maxNarrow<-rbind(maxNarrow,max.NAR2)
write.csv(maxNarrow,file="maxNAR.csv")

minMT.total<-FJ %>% group_by(plot,logger) %>% slice_min(order_by=na_if(MT.total,0))
minMT.rad<-FJ %>% group_by(plot,logger) %>% slice_min(order_by=na_if(MT.rad,0))
minConv<-FJ %>% group_by(plot,logger) %>% slice_min(order_by=na_if(Conv,0))
minNAR<-FJ %>% group_by(plot,logger) %>% slice_min(order_by=na_if(NAR,0))
maxTemp<-FJ %>% group_by(plot,logger) %>% slice_min(order_by=na_if(Temp,0))
dTds<-by(FJ,FJ[,c("plot","logger")],function(x) max(diff(x$Temp,lag=1,differences=1)/diff(x$Elaps,lag=1,differences=1)),simplify=TRUE) 
#
# Calculating harmonic means and confidence intervals for maximum heat fluxes
#
BretFBP<-read.csv("ButlerFBP.csv")
colnames(BretFBP)<-c("Date",colnames(BretFBP[2:30]))
BretCI<-BretFBP %>% 
  filter(Plot != "TTB1" & Plot != "TTB3" & Plot != "TTB4") %>%
  group_by(Plot) %>%
  summarise(PeakTF = Hmean(PeakTF,method="classic",conf.level=0.95, sides="two.sided", na.rm=TRUE),
            PeakRF = Hmean(PeakRF,method="classic",conf.level=0.95, sides="two.sided", na.rm=TRUE),
            PeakCF = Hmean(PeakCF,method="classic",conf.level=0.95, sides="two.sided", na.rm=TRUE),
            PeakNAR =Hmean(PeakNAR,method="classic", conf.level=0.95, sides="two.sided", na.rm=TRUE),
            PeakAirT = MeanCI(PeakAirT,conf.level=0.95,na.rm=TRUE))

BretCI$Stat<-rep(c("Mean","Lower95","Upper95"),8)
factor(BretCI$Plot,levels=c("24Bt","24Bs","24As","24At","16D5","16D6","16D1","TTB2"))
write.csv(BretCI,file="BretCI.csv")

#
# reformatted in Excel and read back in
#
DF<-read.csv("refBretCI.csv")
NoRep<-BretFBP %>%
  filter(Plot == "TTB1" | Plot == "TTB3" | Plot == "TTB4") %>%
  select(c("Plot","PeakTF","PeakRF","PeakCF","PeakNAR","PeakAirT")) %>%
  pivot_longer(cols=2:6,names_to="Type",values_to="Flux")
NoRep$Lower95<-NA
NoRep$Upper95<-NA
DF<-rbind(DF,NoRep)
factor(DF$Type,levels=c("PeakCF","PeakRF","PeakTF","PeakNAR","PeakAirT"),
               labels=c("Convective","Radiant","Total","NARadiant","Air Temperature"))
flux.labs<-c("Convective","Hemispheric Radiant","Hemispheric Total","Narrow Aperture Radiant","Air Temperature")
names(flux.labs)<-c("PeakCF","PeakRF","PeakTF","PeakNAR","PeakAirT")
windows(9,9)
ggplot(DF[(DF[,5]!="PeakAirT"),],aes(x=factor(Plot,levels=c("24Bt","24Bs","24As","24At","16D5","16D6","16D1","TTB1",
                                                            "TTB2","TTB3","TTB4")),
              y=Flux,ymin=Lower95,ymax=Upper95,col=Type))+
  geom_point(aes(shape=18),col="black",size=4)+
  geom_errorbar(stat="identity",size=1.5)+
  scale_shape_identity()+
  theme_classic()+
  xlab("")+
  ylab(expression(paste("Flux ",bgroup("(",over(kW,m^2),")"),sep=" ")))+
  facet_wrap(~Type, scales="free",labeller=labeller(Type = flux.labs))+
  
  scale_color_brewer(
    name="Ft Jackson",
    aesthetics="col",
    type="qual",
    guide="legend",
    palette="Set1") +
  theme(title=element_text(size=20,face="bold",hjust=0.5),
        axis.text.x=element_text(size=20,vjust=0.5,hjust=1,angle=90,face="bold",color="black"),
        axis.text.y=element_text(size=20,vjust=0.5,face="bold",color="black"),
        strip.text.x=element_text(size=20,face="bold",color="blue"),
        legend.title=element_text(size=20,face="bold"),
        legend.text=element_text(size=20,face="bold"),
        legend.position="none")

savePlot(filename = "Mean Flux_FJ_FBP",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)
#
# Replot Bret's original summary data using new style
#
DF3<-read.csv("BretOverallDRW.csv")
colnames(DF3)<-c("Burn",colnames(DF3)[2:6])
factor(DF3$Type,levels=c("Convective","Radiant","Total","NAR","AirT","HorzWind","VertWind"),
       labels=c("Convective","Radiant","Total","NARadiant","Air Temperature","Horizontal Flow",
                "Vertical Flow"))
flux.labs<-c("Convective","Hemispheric Radiant","Hemispheric Total","Narrow Aperture Radiant","Air Temperature",
             "Horizontal Flow","Vertical Flow")
names(flux.labs)<-c("Convective","Radiant","Total","NAR","AirT","HorzWind","VertWind")
ggplot(DF3[1:28,],aes(x=factor(Plot,levels=c("24Bt","24Bs","24As","24At","16D5","16D6","16D1")),
                                     y=Average,group=1,col=Type))+
  geom_point(shape=18,col="black",size=4)+
  facet_wrap(~Type, scales="free",labeller=labeller(Type = flux.labs))+
  geom_line(col="black",linetype=1,size=1)+
  geom_point(aes(Plot,Max,shape=18),col="blue",size=4)+
  geom_line(aes(Plot,Max,group=1),col="blue",linetype=2,size=1)+
  geom_point(aes(Plot,Min,shape=18),col="red",size=4)+
  geom_line(aes(Plot,Min,group=1),col="red",linetype=3,size=1)+
  
  #  geom_errorbar(stat="identity",size=1.5)+
  scale_shape_identity()+
  theme_classic()+
  xlab("Plot")+
  ylab(expression(paste("Flux ",bgroup("(",over(kW,m^2),")"),sep=" ")))+
  
  scale_color_brewer(
    name="Ft Jackson",
    aesthetics="col",
    type="qual",
    guide="legend",
    palette="Set1") +
  theme(title=element_text(size=20,face="bold",hjust=0.5),
        axis.text.x=element_text(size=20,vjust=0.5,hjust=1,angle=90,face="bold",color="black"),
        axis.text.y=element_text(size=20,vjust=0.5,face="bold",color="black"),
        strip.text.x=element_text(size=20,face="bold",color="DarkGreen"),
        legend.title=element_text(size=20,face="bold"),
        legend.text=element_text(size=20,face="bold"),
        legend.position="none")

savePlot(filename = "Bret_Flux_FJ_FBP",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)
windows(9,9)
ggplot(DF3[29:35,],aes(x=factor(Plot,levels=c("24Bt","24Bs","24As","24At","16D5","16D6","16D1")),
                      y=Average,group=1,col=Type))+
  geom_point(shape=18,col="black",size=8)+
  facet_wrap(~Type,scales="fixed",labeller=labeller(Type = flux.labs))+
  geom_line(col="black",linetype=1,size=2)+
  geom_point(aes(Plot,Max,shape=18),col="blue",size=8)+
  geom_line(aes(Plot,Max,group=1),col="blue",linetype=2,size=2)+
  geom_point(aes(Plot,Min,shape=18),col="red",size=8)+
  geom_line(aes(Plot,Min,group=1),col="red",linetype=3,size=2)+
  
  scale_shape_identity()+
  theme_classic()+
  xlab("Plot")+
  ylab(expression(paste(degree,"C")))+
  
  scale_color_brewer(
    name="Ft Jackson",
    aesthetics="col",
    type="qual",
    guide="legend",
    palette="Set1") +
  theme(title=element_text(size=40,face="bold",hjust=0.5),
        axis.text.x=element_text(size=40,vjust=0.5,hjust=1,angle=90,face="bold",color="black"),
        axis.text.y=element_text(size=40,vjust=0.5,face="bold",color="black"),
        strip.text.x=element_text(size=40,face="bold",color="DarkGreen"),
        legend.title=element_text(size=20,face="bold"),
        legend.text=element_text(size=20,face="bold"),
        legend.position="none")

savePlot(filename = "Bret_Air_FJ_FBP",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)

ggplot(DF3[36:42,],aes(x=factor(Plot,levels=c("24Bt","24Bs","24As","24At","16D5","16D6","16D1")),
                       y=Average,group=1,col=Type))+
  geom_point(shape=18,col="black",size=8)+
  facet_wrap(~Type,scales="fixed",labeller=labeller(Type = flux.labs))+
  geom_line(col="black",linetype=1,size=2)+
  geom_point(aes(Plot,Max,shape=18),col="blue",size=8)+
  geom_line(aes(Plot,Max,group=1),col="blue",linetype=2,size=2)+
  geom_point(aes(Plot,Min,shape=18),col="red",size=8)+
  geom_line(aes(Plot,Min,group=1),col="red",linetype=3,size=2)+
  
  #  geom_errorbar(stat="identity",size=1.5)+
  scale_shape_identity()+
  theme_classic()+
  xlab("Plot")+
  ylab("(m/s)")+
  geom_hline(yintercept=0,linetype=4,size=2)+  
  scale_color_brewer(
    name="Ft Jackson",
    aesthetics="col",
    type="qual",
    guide="legend",
    palette="Set1") +
  theme(title=element_text(size=40,face="bold",hjust=0.5),
        axis.text.x=element_text(size=40,vjust=0.5,hjust=1,angle=90,face="bold",color="black"),
        axis.text.y=element_text(size=40,vjust=0.5,face="bold",color="black"),
        strip.text.x=element_text(size=40,face="bold",color="DarkGreen"),
        legend.title=element_text(size=20,face="bold"),
        legend.text=element_text(size=20,face="bold"),
        legend.position="none")

savePlot(filename = "Bret_Vert_FJ_FBP",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)

ggplot(DF3[43:49,],aes(x=factor(Plot,levels=c("24Bt","24Bs","24As","24At","16D5","16D6","16D1")),
                       y=Average,group=1,col=Type))+
  geom_point(shape=18,col="black",size=8)+
  facet_wrap(~Type,scales="fixed",labeller=labeller(Type = flux.labs))+
  geom_line(col="black",linetype=1,size=2)+
  geom_point(aes(Plot,Max,shape=18),col="blue",size=8)+
  geom_line(aes(Plot,Max,group=1),col="blue",linetype=2,size=2)+
  geom_point(aes(Plot,Min,shape=18),col="red",size=8)+
  geom_line(aes(Plot,Min,group=1),col="red",linetype=3,size=2)+
  
  #  geom_errorbar(stat="identity",size=1.5)+
  scale_shape_identity()+
  theme_classic()+
  xlab("Plot")+
  ylab("(m/s)")+
  geom_hline(yintercept=0,linetype=4,size=2)+  
  scale_color_brewer(
    name="Ft Jackson",
    aesthetics="col",
    type="qual",
    guide="legend",
    palette="Set1") +
  theme(title=element_text(size=40,face="bold",hjust=0.5),
        axis.text.x=element_text(size=40,vjust=0.5,hjust=1,angle=90,face="bold",color="black"),
        axis.text.y=element_text(size=40,vjust=0.5,face="bold",color="black"),
        strip.text.x=element_text(size=40,face="bold",color="DarkGreen"),
        legend.title=element_text(size=20,face="bold"),
        legend.text=element_text(size=20,face="bold"),
        legend.position="none")

savePlot(filename = "Bret_Horz_FJ_FBP",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)

#
# Wind direction is a circular variable - positive vertical velocity is upward
# (0) and negative vertical velocity is downward (180)
# horizontal wind - into face (270), away from face (90)
#
DF2<-FJ %>% select(plot,angle,VectWind)
windows(9,9)
spdseq <- seq(0,max(DF2$VectWind),3)
spd.breaks <- c(spdseq,max(DF2$VectWind,na.rm = TRUE))
spd.labels <- seq(1.5,max(DF2$VectWind),3)

spd.binned <- cut(x=DF2$VectWind,breaks=spd.breaks,labels=spd.labels,ordered_result = TRUE)
dirres<-180/8
dir.breaks <- c(seq(dirres/2, 360-dirres/2, by = dirres))
DF2$dir.binned <- cut(DF2$angle,
                  breaks = dir.breaks,
                  ordered_result = TRUE)
DF2$plot<-factor(DF2$plot,levels=c("24Btriangle","24Bsquare","24Asquare","24Atriangle",
                                   "16D5","16D6","16D1","TTB1","TTB2","TTB3","TTB4"),
                          labels=c("24Bt","24Bs","24As","24At","16D5","16D6","16D1",
                                   "TTB1","TTB2","TTB3","TTB4"))

spd.colors<-colorRampPalette(brewer.pal(9,"RdYlBu"))
ggplot(DF2,aes(x=dir.binned,fill=spd.binned))+
  geom_bar()+
  facet_wrap( ~plot)+
  ylab("Number of observations")+
  scale_x_discrete(drop = FALSE,
                   labels = seq(0,360-dirres,dirres)) +
#  coord_polar(start = -9*pi/16, direction=-1) +
  scale_fill_manual(name = "Wind Speed", 
                    values = spd.colors(9),
                    drop = FALSE) +
  theme(axis.title.x = element_blank(),
        panel.grid.major = element_line(colour="black"),
        axis.text.x=element_text(size=10,vjust=0.5,hjust=1,angle=90,face="bold",color="black"),
        axis.text.y=element_text(size=10,vjust=0.5,face="bold",color="black"),
        axis.title=element_text(size=20,face="bold",color="black"),
        strip.text.x=element_text(size=20,face="bold",color="blue"),
        legend.title=element_text(size=20,face="bold"),
        legend.title.align = 0.5,
        legend.text=element_text(size=15,face="bold"),
        legend.text.align = 0.5,
        legend.justification =c(0.5,0.5),
        legend.position=c(0.9,0.12))

savePlot(filename = "Wind_Freq_Bar",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)

ggplot(DF2,aes(x=dir.binned,fill=spd.binned))+
  geom_bar()+
  facet_wrap( ~plot)+
  ylab(expression(paste(log[10],"(Number of observations)")))+
  scale_x_discrete(drop = FALSE,
                 labels = seq(0,360-dirres,dirres)) +
  scale_y_log10()+
  coord_polar(start = -9*pi/16, direction=-1) +
  scale_fill_manual(name = "Wind Speed", 
                    values = spd.colors(9),
                    drop = FALSE) +
  theme(axis.title.x = element_blank(),
        panel.grid.major = element_line(colour="black"),
        axis.text.x=element_text(size=10,vjust=0.5,hjust=1,face="bold",color="black"),
        axis.text.y=element_text(size=10,vjust=0.5,face="bold",color="black"),
        axis.title.y=element_text(size=20,face="bold",color="black"),
        strip.text.x=element_text(size=20,face="bold",color="blue"),
        legend.title=element_text(size=20,face="bold"),
        legend.title.align = 0.5,
        legend.text=element_text(size=20,face="bold"),
        legend.text.align = 0.5,
        legend.justification =c(0.5,0.5),
        legend.position=c(0.9,0.15))

savePlot(filename = "Wind_Freq_Polar",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)

#
# Time plots of vertical and horizontal flow for randomly selected loggers from each fire 
# after ordering in increasing order by logger number and as.integer(runif(1,1,n)) where n
# is number of loggers deployed on burn (excluding logger 5 on 24Bs)
#
# FJ16D1_12, FJ16D5_13,FJ16D6_12,FJ24As_7, FJ24At_6, FJ24Bs_13, FJ24Bt_6
#
DF4<-FJ %>% select(plot,logger,Elaps,Temp,HorzWind,VertWind) %>% 
  subset((plot=="16D1"& logger==12) |(plot=="16D5"& logger==13) |(plot=="16D6"& logger==12) |
         (plot=="24Asquare"& logger==7) | (plot=="24Atriangle"& logger==6) | (plot=="24Bsquare"& logger==13) |
         (plot=="24Btriangle"& logger==6) |
         (plot=="TTB1"& logger==1) |
#         (plot=="TTB2"& logger==1) |
         (plot=="TTB2"& logger==4) |
         (plot=="TTB3"& logger==1) |
         (plot=="TTB4"& logger==4))  %>% group_by(plot,logger) %>% 
         pivot_longer(cols=ends_with("Wind"),names_to = "Flow",values_to = "Velocity")
DF4$Flow<-factor(DF4$Flow,levels=c("HorzWind","VertWind"),labels=c("Horizontal","Vertical"))
DF4$plot<-factor(DF4$plot,levels=c("24Btriangle","24Bsquare","24Asquare","24Atriangle","16D5","16D6","16D1",
                                   "TTB1","TTB2","TTB3","TTB4"),
                          labels=c("24Bt","24Bs","24As","24At","16D5","16D6","16D1","TTB1","TTB2","TTB3","TTB4"))
DF4<-DF4 %>% group_by(plot) %>% 
  mutate(ReLaps = (Elaps-min(Elaps,na.rm=TRUE))/(max(Elaps,na.rm=TRUE)-min(Elaps,na.rm=TRUE)))

windows(9,9)
ggplot(DF4,aes(ReLaps,Velocity,col=Flow))+
  geom_line()+
  theme_classic()+
  facet_wrap( ~plot)+
  xlab("Scaled elapsed time")+
  ylab("Mass Flow (m/s)")+
  scale_color_brewer(
    name="Direction",
    aesthetics="col",
    type="qual",
    guide="legend",
    palette="Set1") +
  geom_hline(yintercept=0,linetype=3,size=1)+  
  theme(axis.text.x=element_text(size=20,vjust=0.5,hjust=1,angle=90,face="bold",color="black"),
        axis.text.y=element_text(size=20,vjust=0.5,face="bold",color="black"),
        axis.title=element_text(size=20,face="bold",color="black"),
        strip.text.x=element_text(size=20,face="bold",color="blue"),
        legend.title=element_text(size=20,face="bold"),
        legend.title.align = 0.5,
        legend.text=element_text(size=20,face="bold"),
        legend.text.align =0,
        legend.justification =c(0.5,0.5),
        legend.position=c(0.9,0.1))

savePlot(filename = "HVWind_Time",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)

ggplot(DF4,aes(ReLaps,Temp))+
  geom_line(color="DarkGreen")+
  theme_classic()+
  facet_wrap( ~plot)+
  xlab("Scaled elapsed time")+
  ylab(expression(paste("Air Temperature (",degree,"C)")))+
  geom_hline(yintercept=0,linetype=3,size=1)+  
  theme(axis.text.x=element_text(size=20,vjust=0.5,hjust=1,angle=90,face="bold",color="black"),
        axis.text.y=element_text(size=20,vjust=0.5,face="bold",color="black"),
        axis.title=element_text(size=20,face="bold",color="black"),
        strip.text.x=element_text(size=20,face="bold",color="blue"),
        legend.title=element_text(size=20,face="bold"),
        legend.title.align = 0.5,
        legend.text=element_text(size=20,face="bold"),
        legend.text.align =0,
        legend.justification =c(0.5,0.5),
        legend.position=c(0.7,0.1))

savePlot(filename = "AirTemp_Time",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)

  #
#  Sara McAllister's new heat flux packages
#
hf16d1<-read.csv("New FBPs//HF16D1.csv")
hf16d1$plot<-"16D1"
hf16d1$RadiantHF<-as.numeric(hf16d1$RadiantHF)
colnames(hf16d1)<-c("UTC",colnames(hf16d1[2:7]))
hf16d1$UTC<-format(hf16d1$UTC,format="%H:%M:%OS",tz="GMT")
start<-min(hf16d1$UTC,na.rm=TRUE)
hf16d1$elaps<-as.POSIXct(hf16d1$UTC,format="%H:%M:%OS")-as.POSIXct(start,format="%H:%M:%OS")

hf16d5<-read.csv("New FBPs//HF16D5.csv")
hf16d5$plot<-"16D5"
hf16d5$RadiantHF<-as.numeric(hf16d5$RadiantHF)
colnames(hf16d5)<-c("UTC",colnames(hf16d5[2:7]))
hf16d5$UTC<-format(hf16d5$UTC,format="%H:%M:%OS",tz="GMT")
start2<-min(hf16d5$UTC,na.rm=TRUE)
hf16d5$elaps<-as.POSIXct(hf16d5$UTC,format="%H:%M:%OS")-as.POSIXct(start2,format="%H:%M:%OS")

hf16d6<-read.csv("New FBPs//HF16D6.csv")
hf16d6$plot<-"16D6"
hf16d6$RadiantHF<-as.numeric(hf16d6$RadiantHF)
colnames(hf16d6)<-c("UTC",colnames(hf16d6[2:7]))
hf16d6$UTC<-format(hf16d6$UTC,format="%H:%M:%OS",tz="GMT")
start3<-min(hf16d6$UTC,na.rm=TRUE)
hf16d6$elaps<-as.POSIXct(hf16d6$UTC,format="%H:%M:%OS")-as.POSIXct(start3,format="%H:%M:%OS")

hf24asq<-read.csv("New FBPs//HF24Asqur.csv")
hf24asq$plot<-"24As"
hf24asq$RadiantHF<-as.numeric(hf24asq$RadiantHF)
colnames(hf24asq)<-c("UTC",colnames(hf24asq[2:7]))
hf24asq$UTC<-format(hf24asq$UTC,format="%H:%M:%OS",tz="GMT")
start4<-min(hf24asq$UTC,na.rm=TRUE)
hf24asq$elaps<-as.POSIXct(hf24asq$UTC,format="%H:%M:%OS")-as.POSIXct(start4,format="%H:%M:%OS")

hf24atr<-read.csv("New FBPs//HF24Atri.csv")
hf24atr$plot<-"24At"
hf24atr$RadiantHF<-as.numeric(hf24atr$RadiantHF)
colnames(hf24atr)<-c("UTC",colnames(hf24atr[2:7]))
hf24atr$UTC<-format(hf24atr$UTC,format="%H:%M:%OS",tz="GMT")
start5<-min(hf24atr$UTC,na.rm=TRUE)
hf24atr$elaps<-as.POSIXct(hf24atr$UTC,format="%H:%M:%OS")-as.POSIXct(start5,format="%H:%M:%OS")

hf24btr<-read.csv("New FBPs//HF24Bdiam.csv")
hf24btr$plot<-"24Bs"
hf24btr$RadiantHF<-as.numeric(hf24btr$RadiantHF)
colnames(hf24btr)<-c("UTC",colnames(hf24btr[2:7]))
hf24btr$UTC<-format(hf24btr$UTC,format="%H:%M:%OS",tz="GMT")
start6<-min(hf24btr$UTC,na.rm=TRUE)
hf24btr$elaps<-as.POSIXct(hf24btr$UTC,format="%H:%M:%OS")-as.POSIXct(start6,format="%H:%M:%OS")

DF<-rbind(hf16d1,hf16d5,hf16d6,hf24asq,hf24atr,hf24btr)
DF$plot<-factor(DF$plot,levels=c("24Bs","24As","24At","16D5","16D6","16D1"))
windows(9,9)
ggplot(DF,aes(as.double(elaps),RadiantHF,col=plot))+
#  geom_point(color="red")+
  geom_line(size=1.5)+
  theme_classic()+
  xlab("Elapsed time (s)")+
  ylab(expression(paste("Radiant flux ",bgroup("(",over(kW,m^2),")"),sep=" ")))+
  facet_wrap(~plot)+
  
    scale_color_brewer(
    name="Ft Jackson",
    aesthetics="col",
    type="qual",
    guide="legend",
    palette="Set1") +
  theme(title=element_text(size=20,face="bold",hjust=0.5),
        axis.text.x=element_text(size=15,vjust=0.5,hjust=1,angle=90,face="bold",color="black"),
        axis.text.y=element_text(size=20,vjust=0.5,face="bold",color="black"),
        strip.text.x=element_text(size=15,face="bold",color="blue"),
        legend.title=element_text(size=20,face="bold"),
        legend.text=element_text(size=20,face="bold"),
        legend.position="none")

savePlot(filename = "RadFlux_FJ_HF",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)
#
# Flame Ionization Detector (FID)
#

fid16d1<-read.csv("New FBPs//FID16D1.csv")
colnames(fid16d1)<-c("Time","Volt")
fid16d1$plot<-"16D1"
fid16d1$Time<-format(fid16d1$Time,format="%H:%M:%OS",tz="GMT")
start<-min(fid16d1$Time,na.rm=TRUE)
fid16d1$elaps<-as.POSIXct(fid16d1$Time,format="%H:%M:%OS")-as.POSIXct(start,format="%H:%M:%OS")
fid16d5<-read.csv("New FBPs//FID16D5.csv")
colnames(fid16d5)<-c("Time","Volt")
fid16d5$plot<-"16D5"
fid16d5$Time<-format(fid16d5$Time,format="%H:%M:%OS",tz="GMT")
start2<-min(fid16d5$Time,na.rm=TRUE)
fid16d5$elaps<-as.POSIXct(fid16d5$Time,format="%H:%M:%OS")-as.POSIXct(start2,format="%H:%M:%OS")
fid16d6<-read.csv("New FBPs//FID16D6.csv")
colnames(fid16d6)<-c("Time","Volt")
fid16d6$plot<-"16D6"
fid16d6$Time<-format(fid16d6$Time,format="%H:%M:%OS",tz="GMT")
start3<-min(fid16d6$Time,na.rm=TRUE)
fid16d6$elaps<-as.POSIXct(fid16d6$Time,format="%H:%M:%OS")-as.POSIXct(start3,format="%H:%M:%OS")
fid24As<-read.csv("New FBPs//FID24As.csv")
colnames(fid24As)<-c("Time","Volt")
fid24As$plot<-"24As"
fid24As$Time<-format(fid24As$Time,format="%H:%M:%OS",tz="GMT")
start4<-min(fid24As$Time,na.rm=TRUE)
fid24As$elaps<-as.POSIXct(fid24As$Time,format="%H:%M:%OS")-as.POSIXct(start4,format="%H:%M:%OS")
fid24At<-read.csv("New FBPs//FID24At.csv")
colnames(fid24At)<-c("Time","Volt")
fid24At$plot<-"24At"
fid24At$Time<-format(fid24At$Time,format="%H:%M:%OS",tz="GMT")
start5<-min(fid24At$Time,na.rm=TRUE)
fid24At$elaps<-as.POSIXct(fid24At$Time,format="%H:%M:%OS")-as.POSIXct(start5,format="%H:%M:%OS")
DF5<-rbind(fid24As,fid24At,fid16d5,fid16d6,fid16d1)
DF5$plot<-factor(DF5$plot,levels=c("24As","24At","16D5","16D6","16D1"))
windows(9,9)
fid<-ggplot(DF5,aes(as.double(elaps),Volt,col=plot))+
  #  geom_point(color="red")+
  geom_point(size=1.5)+
  theme_classic()+
  xlab("Elapsed time (s)")+
  ylab("voltage")+
  facet_wrap(~plot)+
  scale_color_brewer(
    name="Residence time",
    aesthetics="col",
    type="qual",
    guide="legend",
    palette="Set1") +
  theme(title=element_text(size=20,face="bold",hjust=0.5),
        axis.text.x=element_text(size=15,vjust=0.5,hjust=1,angle=90,face="bold",color="black"),
        axis.text.y=element_text(size=20,vjust=0.5,face="bold",color="black"),
        strip.text.x=element_text(size=15,face="bold",color="blue"),
        legend.title=element_text(size=20,face="bold"),
        legend.text=element_text(size=20,face="bold",hjust=0),
        legend.position=c(0.8,0.2))
ggdraw(fid)+
  draw_label("= 335 s",x=0.93,y=0.3,size=20,fontface="bold",hjust=1)+
  draw_label("=  49 s",x=0.925,y=0.27,size=20,fontface="bold",hjust=1)+
  draw_label("=  77 s",x=0.925,y=0.24,size=20,fontface="bold",hjust=1)+
  draw_label("=  31 s",x=0.925,y=0.21,size=20,fontface="bold",hjust=1)+
  draw_label("= 234 s",x=0.93,y=0.185,size=20,fontface="bold",hjust=1)

savePlot(filename = "ResTime_FJ_FID",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)



