#Data from tests on 30/10/2019 at met city - including recent crack formation even. setwd("/Users/sarcher/Desktop/R data/1_MOSAiC/") CO2_031119 <- read.csv("~/Desktop/R data/1_MOSAiC/031119_CO2.csv", header = TRUE) library(tidyverse) library(dplyr) library(ggplot2) #simplyfying the df name AX <- CO2_031119 head(AX) #narrow down the axes by altering XR1, XR2 and YR1, YR2 #X and Y -axis range lower, higher XR1 <- 0 XR2 <- 4000 YR1 <- 411.5 YR2 <- 412.5 #Plotting the truncated axes Plot_S1 <- ggplot(AX) + geom_point(aes(x=sec, y = CO2)) + xlim(XR1, XR2) + ylim(YR1,YR2) + labs(title = "031119_CO2") My_Theme = theme( axis.title.x = element_text(size = 20), axis.text.x = element_text(size = 16), axis.title.y = element_text(size = 20), axis.text.y = element_text(size = 16), legend.title = element_text(size =28), title = element_text(size = 20)) Plot_S1 + My_Theme #### Mode transition times: judged from change in voltage of switch Start <- 610 T1 <- 790 T2 <- 1110 T3 <- 1435 T4 <- 1750 T5 <- 2200 T6 <- 2510 T7 <- 2835 T8 <- 3145 T9 <- 3460 T10 <- 3770 T11 <- 4080 T12 <- 4400 T13 <- 4715 Finish <- 5035 #XXXXXXXXXXXXXXXXX #Selecting the sequences (modes) between switching valves, at the moment 60s each side M1 <- AX %>% select(sec,CO2, H2O) %>% filter(sec >Start + 60, sec < T1 - 60) M2 <- AX %>% select(sec,CO2, H2O) %>% filter(sec > T1 + 100, sec < T2 - 60) M3 <- AX %>% select(sec,CO2, H2O) %>% filter(sec >T2 + 90, sec < T3 - 60) M4 <- AX %>% select(sec,CO2, H2O) %>% filter(sec > T3 + 90, sec < T4 - 60) M5 <- AX %>% select(sec,CO2, H2O) %>% filter(sec > T4 + 125, sec < T5 - 30) M6 <- AX %>% select(sec,CO2, H2O) %>% filter(sec > T5 + 60, sec < T6 - 60) M7 <- AX %>% select(sec,CO2, H2O) %>% filter(sec > T6 + 60, sec < T7 - 60) M8 <- AX %>% select(sec,CO2, H2O) %>% filter(sec > T7 + 60, sec < T8 - 60) M9 <- AX %>% select(sec,CO2, H2O) %>% filter(sec > T8 + 60, sec < T9 - 60) M10 <- AX %>% select(sec,CO2, H2O) %>% filter(sec > T9 + 140, sec < T10 - 10) M11 <- AX %>% select(sec,CO2, H2O) %>% filter(sec > T10 + 60, sec < T11 - 30) M12 <- AX %>% select(sec,CO2, H2O) %>% filter(sec > T11 + 90, sec < T12 - 60) M13 <- AX %>% select(sec,CO2, H2O) %>% filter(sec > T12 + 60, sec < T13 - 60) M14 <- AX %>% select(sec,CO2, H2O) %>% filter(sec > T13 + 90, sec < Finish - 30) #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX #Trying to account for the inherent trends during the measurements #adding Chamber modes together: M3, M5, M7, M9 # and Inlet modes together: M2, M4, M6, M8 Chamber <- rbind(M3,M5,M7,M9) head(Chamber) str(Chamber) Inlet <- rbind(M2,M4,M6,M8) head(Inlet) str(Inlet) #Assigning mode to each of the sequences and combining them in one dataframe (CHIN) A <- rep('red',866) B <- rep('blue', 701) mode <- c(A,B) CHIN <- rbind(M3,M5,M7,M9, M2,M4,M6,M8) CHIN <- cbind(CHIN,mode) head(CHIN) #Plotting the combined data #Setting axis limits: RX1 <- 800 RX2 <- 2750 RY1 <- 410.75 RY2 <- 411.1 My_Theme = theme( axis.title.x = element_text(size = 20), axis.text.x = element_text(size = 16), axis.title.y = element_text(size = 20), axis.text.y = element_text(size = 16), legend.title = element_text(size =28), title = element_text(size = 20)) #plotting chamber and inlet separately by faceting Plot_In <- ggplot(data = CHIN)+ geom_point(mapping = aes(x = sec, y = CO2), colour = mode) + geom_smooth(data = Inlet, method = loess, se = TRUE, aes(x = sec, y = CO2), colour = 'blue')+ geom_smooth(data = Chamber, method = loess, se = TRUE, aes(x = sec, y = CO2), colour = 'red')+ labs(title = "301019_CO2 loess smoothed") Plot_In+My_Theme #calculating the trends using loess/polynomial i_loess <- loess(CO2 ~ sec, data = Inlet, degree = 2) c_loess <- loess(CO2 ~ sec, data = Chamber, degree = 2) i_pred <- predict(i_loess, data.frame(sec = seq(850, 3400, 1)), se = TRUE) c_pred <- predict(c_loess, data.frame(sec = seq(850, 3400, 1)), se = TRUE) summary(i_loess) #Inlet prediction: i_pred.df <- data.frame(i_pred) time <- seq(850,3400,1) i_pred.df <- cbind(time, i_pred.df) head(i_pred.df) Predicted <- ggplot(data = i_pred.df) + geom_point(mapping = aes(x = time, y = fit), size = 1, colour = 'dark blue')+ labs(title = "CO2 Inlet prediction") Predicted + My_Theme view(i_pred.df) #Data_Range ##trying to detrend the data - subtracting predicted inlet trend from data #M2 head(M2) DR1 = 891 DR2 = DR1+158 IP2<- filter(i_pred.df, time >= DR1, time <= DR2) IP_M2 <- cbind(M2,IP2) IP_M2 IP_M2 <- mutate(IP_M2, Detrend = CO2 - fit) #M3 head(M3) DR1 = 1201 DR2 = DR1+173 IP3<- filter(i_pred.df, time >= DR1, time <= DR2) IP_M3 <- cbind(M3,IP3) IP_M3 IP_M3 <- mutate(IP_M3, Detrend = CO2 - fit) #M4 head(M4) D4S = 1526 D4E = D4S+163 IP4<- filter(i_pred.df, time >= D4S, time <= D4E) IP_M4 <- cbind(M4,IP4) IP_M4 IP_M4 <- mutate(IP_M4, Detrend = CO2 - fit) #M5 head(M5) s = 1876 e = s-1+294 IP5 <- filter(i_pred.df, time >= s, time <= e) IP_M5 <- cbind(M5,IP5) IP_M5 IP_M5 <- mutate(IP_M5, Detrend = CO2 - fit) #M6 head(M6) s = 2261 e = s-1+189 IP6 <- filter(i_pred.df, time >= s, time <= e) IP_M6 <- cbind(M6,IP6) IP_M6 IP_M6 <- mutate(IP_M6, Detrend = CO2 - fit) #M7 head(M7) head(M7) s = 2571 e = s-1+204 IP7 <- filter(i_pred.df, time >= s, time <= e) IP_M7 <- cbind(M7,IP7) IP_M7 IP_M7 <- mutate(IP_M7, Detrend = CO2 - fit) #M8 head(M8) s = 2896 e = s-1+189 IP8 <- filter(i_pred.df, time >= s, time <= e) IP_M8 <- cbind(M8,IP8) IP_M8 IP_M8 <- mutate(IP_M8, Detrend = CO2 - fit) #M9 head(M9) s = 3206 e = s-1+194 IP9 <- filter(i_pred.df, time >= s, time <= e) IP_M9 <- cbind(M9,IP9) IP_M9 IP_M9 <- mutate(IP_M9, Detrend = CO2 - fit) ### putting the detrended data together... C <- rep('dark red', 174+294+204+194) D <- rep('dark blue', 159+164+189+189) mode <- c(C,D) Detrend_C <- rbind(IP_M3,IP_M5,IP_M7,IP_M9) Detrend_I <- rbind(IP_M2,IP_M4,IP_M6, IP_M8) Detrend_A <- rbind(IP_M3,IP_M5,IP_M7,IP_M9,IP_M2,IP_M4,IP_M6, IP_M8) Detrend_B <- cbind(Detrend_A,mode) Plot_Detrend <- ggplot(data = Detrend_B)+ geom_point(mapping = aes(x = sec, y = Detrend), colour = mode)+ labs(title = "301019: Detrended CO2 data") Plot_Detrend + My_Theme str(Detrend_A) str(Detrend_B) # trying to show as a box plot - differences between sequences: Box_plot_Detrend <- ggplot(data = Detrend_B) + geom_point(mapping = aes(x = sec, y = Detrend), colour = mode)+ geom_boxplot(data = IP_M2,mapping = aes(x = sec, y = Detrend))+ geom_boxplot(data = IP_M3,mapping = aes(x = sec, y = Detrend))+ geom_boxplot(data = IP_M4,mapping = aes(x = sec, y = Detrend))+ geom_boxplot(data = IP_M5,mapping = aes(x = sec, y = Detrend))+ geom_boxplot(data = IP_M6,mapping = aes(x = sec, y = Detrend))+ geom_boxplot(data = IP_M7,mapping = aes(x = sec, y = Detrend))+ geom_boxplot(data = IP_M8,mapping = aes(x = sec, y = Detrend))+ geom_boxplot(data = IP_M9,mapping = aes(x = sec, y = Detrend))+ labs(title = "301019: Detrended CO2 data_B") Box_plot_Detrend + My_Theme IP_M9 head(i_pred) (i_pred.df) # to plot modes separately add: facet_grid(rows = vars(mode)) # to define the axes add: xlim(RX1, RX2) + ylim(RY1,RY2) + #Calculating the regressions: z.lm = lm(CO2 ~ sec, data = Chamber) summary(z.lm) v.lm <- lm(CO2 ~ sec, data = Inlet) summary(v.lm) #ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ # Same again but altering which modes are taken into account #Trying to account for the inherent trends during the measurements #adding Chamber modes together: M3, M5, M7, M9 # and Inlet modes together: M2, M4, M6, M8 ChamberB <- rbind(M3,M5) head(ChamberB) str(ChamberB) InletB <- rbind(M2,M4) head(InletB) str(InletB) #Assigning mode to each of the sequences and combining them in one dataframe (CHIN) A <- rep('red',468) B <- rep('blue',323) mode <- c(A,B) CHIN_B <- rbind(M3,M5,M2,M4) CHIN_B <- cbind(CHINB,mode) head(CHIN_B) #Plotting the combined data #Setting axis limits: RX1 <- 800 RX2 <- 2750 RY1 <- 410.75 RY2 <- 411.1 My_Theme = theme( axis.title.x = element_text(size = 20), axis.text.x = element_text(size = 16), axis.title.y = element_text(size = 20), axis.text.y = element_text(size = 16), legend.title = element_text(size =28), title = element_text(size = 20)) #Plotting Chamber and Inlet on the same plot with separte regressions - reduced number of sequences Plot_B <- ggplot(data = CHIN_B)+ geom_point(mapping = aes(x = sec, y = CO2), colour = mode) + geom_smooth(data = ChamberB, method = lm, se = TRUE, aes(x = sec, y = CO2), colour = 'red') + geom_smooth(data = InletB, method = lm, se = TRUE, aes(x = sec, y = CO2), colour = 'blue') Plot_B #plotting chamber and inlet separately by faceting Plot_Both <- ggplot(data = CHIN)+ geom_point(mapping = aes(x = sec, y = CO2), colour = mode) + geom_smooth(method = lm, se = TRUE, aes(x = sec, y = CO2))+ facet_grid(rows = vars(mode)) Plot_Both #Plotting Chamber and Inlet on the same plot with separte regressions Plot_Both <- ggplot(data = CHIN)+ geom_point(mapping = aes(x = sec, y = CO2), colour = mode) + geom_smooth(data = Chamber, method = lm, se = TRUE, aes(x = sec, y = CO2), colour = 'red') + geom_smooth(data = Inlet, method = lm, se = TRUE, aes(x = sec, y = CO2), colour = 'blue') Plot_Both #Plotting each individual mode: Plot_M14 <- ggplot(M14) + geom_point(aes(x=sec, y = CH4)) + labs(title = "M14_CH4") + abline(lm) Plot_M14 # #M1 - looks like increasing trend in methane 1966.7 to 1966.9 - inlet #M2 - flat trend 1966.7 - chamber #M3 - flat trend 1966.65 - inlet #M4 - flat trend - slightly decreasing - 1966.65 - chamber #M5 - flat trend - slight decrease - 1966.5 - inlet #M6 - flat trend - 1966.4 - chamber #M7 - flat trend - 1966.35 - inlet #M8 - flat trend -decreasing - 1966.35 - chamber #M9 - f/d #M10 - decreasing trend - 1966.35 to 1966.25 #M11 - flat trend 1966.275 #M12 - increasing trend - 1966.8 to 1967.0 #M13 - increasing trend - 1966.9 to 1966.95 #M14 - increasing trend - 1967.0 to 1967.1s #Calcualting the stats for each mode and gas M1S <- pull(M1, CO2) M2S <- pull(M2, CO2) M3S<- pull(M3, CO2) M4S<- pull(M4, CO2) M5S<- pull(M5, CO2) M6S<- pull(M6, CO2) M7S<- pull(M7, CO2) M8S<- pull(M8, CO2) M9S<- pull(M9, CO2) M10S<- pull(M10, CO2) M11S<- pull(M11, CO2) M12S<- pull(M12, CO2) M13S<- pull(M13, CO2) M14S<- pull(M14, CO2) M1S #means of each mode x1 <- mean(M1S) x2<- mean(M2S) x3<- mean(M3S) x4<- mean(M4S) x5<- mean(M5S) x6<- mean(M6S) x7<- mean(M7S) x8<- mean(M8S) x9<- mean(M9S) x10<- mean(M10S) x11<- mean(M11S) x12<- mean(M12S) x13<- mean(M13S) x14<- mean(M14S) x1 #corresponding standard deviation sd1 <- sd(M1S) sd2<- sd(M2S) sd3<- sd(M3S) sd4<- sd(M4S) sd5<- sd(M5S) sd6<- sd(M6S) sd7<- sd(M7S) sd8<- sd(M8S) sd9<- sd(M9S) sd10<- sd(M10S) sd11<- sd(M11S) sd12<- sd(M12S) sd13<- sd(M13S) sd14<- sd(M14S) #corrresponding CV (%) cv1 <- 100*sd1/x1 cv2<- 100*sd2/x2 cv3<- 100*sd3/x3 cv4<- 100*sd4/x4 cv5<- 100*sd5/x5 cv6<- 100*sd6/x6 cv7<- 100*sd7/x7 cv8<- 100*sd8/x8 cv9<- 100*sd9/x9 cv10<- 100*sd10/x10 cv11<- 100*sd11/x11 cv12<- 100*sd12/x12 cv13<- 100*sd13/x13 cv14<- 100*sd14/x14 #Stats for single mode #Tabulating the results Xs <- c(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14) Sds <- c(sd1,sd2,sd3,sd4,sd5,sd6,sd7,sd8,sd9,sd10, sd11,sd12,sd13,sd14) CVs <- c(cv1,cv2,cv3,cv4,cv5,cv6,cv7,cv8,cv9,cv10,cv11,cv12,cv13,cv14) Mins <- c(min1,min2,min3,min4,min5,min6) Max <- c(max1,max2,max3,max4,max5,max6) Xs Mode_concs_301019 <- data.frame('Mode' = c(1:14), 'Means'= Xs,'SD'= Sds,'CV'= CVs) Mode_concs_301019 #Plotting means for each mode Plot_Xs <- ggplot(Mode_concs_301019, aes(x = x, y = Means))+ geom_point() Plot_Xs Plot_Xs <- ggplot(M7, aes(x = sec, y = CH4))+ geom_point() Plot_Xs #Stats comparison of mean from each mode #transition 1 t.test(M1S, M2S,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #transtion 2 t.test(M2S, M3S,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #transtion 3 t.test(M3S, M4S,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #transtion 4 t.test(M4S, M5S,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #transtion 5 t.test(M5S, M6S,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #transtion 6 t.test(M6S, M7S,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #transtion 7 t.test(M7S, M8S,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #transtion 8 t.test(M8S, M9S,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #transtion 9 t.test(M9S, M10S,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #transtion 10 t.test(M10S, M11S,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #transtion 11 t.test(M11S, M12S,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #transtion 12 t.test(M12S, M13S,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #transtion 13 t.test(M13S, M14S,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #Check t.test(M1S, M14S,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #########XXXXXXXXx signif(Xs,6) Xs signif(Sds,4) signif(CVs,3) signif(Mins, 6) signif(Max, 6) A <- matrix(Xs, nrow = 6, ncol = 1) B <- matrix(Sds, nrow = 6, ncol = 1) C <- matrix(CVs, nrow = 6, ncol = 1) D <- matrix(Mins, nrow = 6, ncol = 1) E <- matrix(Max, nrow = 6, ncol = 1) AtoE <- cbind(A,B,C,D,E) colnames(AtoE) <- c("Mean", "SD","CV", "Min","Max") rownames(AtoE) <- c('EQ1', 'I1', 'C1','I2', 'TI1', 'TC1') AtoE ??numbers #this worked last night #converting ppb to ng/l #V = volume per mole V = nRT/P, K is Kelvin, 1 atm K <- 273 R <- 0.082057 n <- 1 P <- 1 V <- n*R*K/P W <- 16 #CH4mass is ug/m3 M1 <- mutate(M,CH4mass = W*CH4/V) head(M1) t.test(I1_CH4, C1_CH4,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) t.test(I1_CH4, C1_CH4,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = TRUE, conf.level = 0.95) t.test(I1_CH4, C1_CH4,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = TRUE, conf.level = 0.95) #Add flow rates here: Q1 <- 5.08 Q2 <- 5.09 Q3 <- 5.08 Q4 <- 5.10 Q5 <- 5.09 Flow <- c(Q2,Q3,Q4,Q5) # F1 removed- prior to actual CH4 measurement #converting to (m3 min-1) Qm2 <- Flow/1000 Q_x <- mean(Qm2) Q_sd <- sd(Qm2) Q_cv <- 100 * Q_sd/Q_x Q_x Q_sd Q_cv #Add chamber diameter: inches D1 <- 26+(5/16) D2 <- 26+(0/16) D3 <- 26+(0/16) D4 <- 25+(14/16) Diam <- c(D1,D2,D3,D4) D_x <- mean(Diam) D_sd <- sd(Diam) D_cv <- 100*D_sd/D_x D_x D_sd D_cv #Chamber area (Converts D_x to meters) A is the area in m2. D_xm <- (D_x * 2.25)/100 D_xm A <- pi*(D_xm/2)*(D_xm/2) A #Calculating flux: #F = Q * (C-I) / A #ug/m2/min Flux1 <- Q_x *(x3 - x1)/A Flux1 Flux2 <- Q_x *(x4 - x2)/A Flux2 Avg <- mean(c(Flux1,Flux2)) Avg # Flux umol / m2/ d Flux1B <- Flux1 * 1440/16 Flux2B <- Flux2 *1440/16 Flux1B Flux2B AvgB <- mean(c(Flux1B,Flux2B)) AvgB Flux1 Parameters <- data.frame('x' = c('Flux1','Flux2','Avg','A','Q_x','Q_sd','Q_cv','x1','cv1','x2','cv2','x3','cv3','x4','cv4')) Parameters OutA <- data.frame('y' = c(Flux1,Flux2,Avg,A,Q_x,Q_sd,Q_cv,x1,cv1,x2,cv2,x3,cv3,x4,cv4)) OutA OutB <- data.frame( = c(Flux1B, Flux2B, AvgB,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)) OutB c Results_1 x1 sd1 min1 max1 cv1 I1A <- tibble(x1,sd1,min1,max1,cv1) I1A write.table(I1A,file=pipe("pbcopy"),sep="\t",col.names=NA) surveys %>% filter(!is.na(hindfoot_length)) %>% group-by(species_id) %>% summarise(mean_hindfoot_length = mean(hindfoot_length), min_hindfoot_length = min(hindfoot_length), max_hindfoot_length = max(hindfoot_length), count = (n)) Plot_All CH4L3_1T <- filter(CH4L3_1, FLUX_TIME <3250) head(CH4L3_1T) Plot_T1 <- ggplot(CH4L3_1T) + geom_point(aes(x=FLUX_TIME, y = CO2)) Plot_T1 #defining the start in seconds S <- 1570597200 F <- 1570600360 Duration <- F-S Duration #adding 2 hours (duration of full test), Total is the total secs of the test. End <- S+120*60 End Total <- End-S Total ChamberT1 <- (S+60*6)-S ChamberT1 ChamberT1F <- (S+60*26)-S ChamberT1F #plotting full sequence: CH4L3_All <- filter(CH4L3_1, FLUX_TIME <7200) head(CH4L3_All) CH4L3_All_B <- filter(CH4L3_All, CO2 < 417.5, CO2 > 410) #Filtering the CH4 data CH4L3_All_C<- filter(CH4L3_All, CH4 < 1950, CH4 > 1940) head(CH4L3_All_C) #not required if already done in Excel CH4L3_1 <- mutate(CH4L3, FLUX_TIME = SECONDS-1570597200) head(CH4L3_1) #plot all time series Plot_All <- ggplot(CH4L3_1T) + geom_point(aes(x=FLUX_TIME, y = CH4)) + xlim(500, 3250) + ylim(1940,1950) Plot_All XXXXXXXXXXx # Water trends?? CH4L3_1 <- mutate(CH4L3, FLUX_TIME = SECONDS-1570597200) head(CH4L3_1) CH4L3_1T <- filter(CH4L3_1, FLUX_TIME <3250) Plot_All <- ggplot(CH4L3_1T) + geom_point(aes(x=FLUX_TIME, y = H2O))+ xlim(500, 3250) + ylim(1250, 1650) Plot_All summary(CH4L3_1) #checking thenumber of rows of data nrow(CH4L3_All_B) #Extra stuff scale_x_continuous(breaks = seq(0, 23, 5)) + scale_y_continuous(breaks = seq(0, 500, 10)) xlim(-0.1,23)+