#Data from tests on 30/10/2019 at met city - including recent crack formation event. setwd("/Users/sarcher/Desktop/R data/1_MOSAiC/") AX <- read.csv("~/Desktop/R data/1_MOSAiC/1_Original_Data/1_CH4_System/2019_12_02_CH4.csv", header = TRUE) library(tidyverse) library(dplyr) #simplyfying the df name head(AX) #narrow down the axes by altering XR1, XR2 and YR1, YR2 #X and Y -axis range lower, higher XR1 <- 0 XR2 <- 6500 YR1 <- 1960 YR2 <- 1965 #Plotting the truncated axes 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 =20), legend.text = element_text(size = 14), title = element_text(size = 20)) Plot_S1 <- ggplot(AX) + geom_point(aes(x=sec, y = CH4)) + xlim(XR1, XR2) + ylim(YR1,YR2) + labs(title = "02122019_CH4") Plot_S1 + My_Theme #### Mode transition times: judged from change in voltage of switch #location 1 Start <-675+80 # M1 inlet T1 <- 840+80 # M2 chamber T2 <- 1020+80 # M3 inlet T3 <- 1215+80 # M4 chamber T4 <- 1395+80 # M5 inlet T5 <- 1590+80 # M6 chamber T6 <- 1775+80 # M7 inlet T7 <- 1965+80 # M8 chamber T8 <- 2150+80 # M9 inlet T9 <- 2340+80 # M10 chamber T10 <- 2525+80 # M11 inlet T11 <- 2715+80 # M12 - chamber T12 <- 2900 +80 # M13 - inlet Finish <- 3080 + 80 #Selecting the sequences (modes) between switching valves, at the moment ~ 60s after each transition # and ~ 20 secs before M1 <- AX %>% select(sec, CH4, CO2, H2O) %>% filter(sec > Start +0, sec < T1 - 25)#I M2 <- AX %>% select(sec,CH4,CO2, H2O) %>% filter(sec > T1 + 25, sec < T2 - 25)#C M3 <- AX %>% select(sec,CH4,CO2, H2O) %>% filter(sec >T2 + 25, sec < T3 - 25)#I M4 <- AX %>% select(sec,CH4,CO2, H2O) %>% filter(sec > T3 + 25, sec < T4 - 25)#C M5 <- AX %>% select(sec,CH4,CO2, H2O) %>% filter(sec > T4 + 25, sec < T5 - 25)#I M6 <- AX %>% select(sec,CH4,CO2, H2O) %>% filter(sec > T5 + 25, sec < T6 - 25)#C M7 <- AX %>% select(sec,CH4,CO2, H2O) %>% filter(sec > T6 + 20, sec < T7 - 25)#I M8 <- AX %>% select(sec,CH4,CO2, H2O) %>% filter(sec > T7 + 25, sec < T8 - 25)#C M9 <- AX %>% select(sec,CH4,CO2, H2O) %>% filter(sec > T8 + 25, sec < T9 - 25)#I M10 <- AX %>% select(sec,CH4,CO2, H2O) %>% filter(sec > T9 + 25, sec < T10 - 25)#C M11 <- AX %>% select(sec ,CH4,CO2, H2O) %>% filter(sec > T10 + 25, sec < T11 - 25)#I M12 <- AX %>% select(sec,CH4,CO2, H2O) %>% filter(sec > T11 + 25, sec < T12 - 25)#C M13 <- AX %>% select(sec,CH4,CO2, H2O) %>% filter(sec > T12 + 25, sec < Finish - 25)#I #Plotting each mode - if you want???? Plot_M <- ggplot(M6) + geom_point(aes(x=sec, y = CH4)) + labs(title = "M_CH4") Plot_M #Trying to account for the inherent trends during the measurements #adding Chamber modes together: # and Inlet modes together: Chamber <- rbind(M10, M12) head(Chamber) str(Chamber) Inlet <- rbind(M9,M11,M13) head(Inlet) str(Inlet) #Assigning mode to each of the sequences and combining them in one dataframe (CHIN) A <- rep('red',268) B <- rep('blue',407) mode <- c(A,B) CHIN <- rbind(M10,M12,M9,M11,M13) CHIN <- cbind(CHIN,mode) head(CHIN) #Plotting Chamber and Inlet on the same plot with separte regressions Plot_In <- ggplot(data = CHIN) + geom_point(mapping = aes(x = sec, y = CH4), colour = mode) + geom_smooth(data = Inlet, method = loess, se = TRUE, aes(x = sec, y = CH4), colour = 'blue') + labs(title = '02122019 CH4 Location 1') Plot_In + My_Theme Plot_In <- ggplot(data = CHIN) + geom_point(mapping = aes(x = sec, y = CH4), colour = mode) + geom_smooth(data = Inlet, method = loess, se = TRUE, aes(x = sec, y = CH4),colour = 'blue') + geom_smooth(data = Chamber, method = loess, se = TRUE, aes(x = sec, y = CH4),colour = 'red') + labs(title = '02122019 CH4 Location 1') Plot_In + My_Theme #Plotting theactual curve fits i_loess <- loess(CH4 ~ sec, data = Inlet, degree = 2) c_loess <- loess(CH4 ~ sec, data = Chamber, degree = 2) # not used i_loess i_pred <- predict(i_loess, data.frame(sec = seq(755,3160, 1)), se = TRUE) # generates data frame of prediction i_pred i_pred.df <- data.frame(i_pred) time <- seq(755, 3160,1) #combines time and predicitons i_pred.df <- cbind(time, i_pred.df) i_pred.df Predicted_I <- ggplot(data = i_pred.df) + geom_point(mapping = aes(x = time, y = fit), size = 1, colour = 'dark blue') Predicted_I #Detrending the data? D1 <- filter(i_pred.df, time > Start, time < T1- 25) M1D1 <- cbind(M1,D1) M1DT <- mutate(M1D1, detrend = CH4 - fit) head(M1DT) D2 <- filter(i_pred.df, time > T1+25, time < T2- 25) M2D2 <- cbind(M2,D2) M2DT <- mutate(M2D2, detrend = CH4 - fit) head(M2DT) D3 <- filter(i_pred.df, time > T2+25, time < T3- 25) M3D3 <- cbind(M3,D3) M3DT <- mutate(M3D3, detrend = CH4 - fit) head(M3DT) D4 <- filter(i_pred.df, time > T3+25, time < T4- 25) M4D4 <- cbind(M4,D4) M4DT <- mutate(M4D4, detrend = CH4 - fit) head(M4DT) D5 <- filter(i_pred.df, time > T4+25, time < T5- 25) M5D5 <- cbind(M5,D5) M5DT <- mutate(M5D5, detrend = CH4 - fit) head(M5DT) D6 <- filter(i_pred.df, time > T5 + 25, time < T6 -25) M6D6 <- cbind(M6,D6) M6DT <- mutate(M6D6, detrend = CH4 - fit) head(M6DT) D7 <- filter(i_pred.df, time > T6 + 20, time < T7 -25) M7D7 <- cbind(M7,D7) M7DT <- mutate(M7D7, detrend = CH4 - fit) head(M7DT) D8 <- filter(i_pred.df, time > T7 + 25, time < T8- 25) M8D8 <- cbind(M8,D8) M8DT <- mutate(M8D8, detrend = CH4 - fit) head(M8DT) D9 <- filter(i_pred.df, time > T8 + 25, time < T9- 25) M9D9 <- cbind(M9,D9) M9DT <- mutate(M9D9, detrend = CH4 - fit) head(M9DT) D10 <- filter(i_pred.df, time > T9 + 25, time < T10- 25) M10D10 <- cbind(M10,D10) M10DT <- mutate(M10D10, detrend = CH4 - fit) head(M10DT) D11 <- filter(i_pred.df, time > T10 + 25, time < T11 - 25) M11D11 <- cbind(M11,D11) M11DT <- mutate(M11D11, detrend = CH4 - fit) head(M11DT) D12 <- filter(i_pred.df, time > T11 + 25, time < T12 - 25) M12D12 <- cbind(M12,D12) M12DT <- mutate(M12D12, detrend = CH4 - fit) head(M12DT) D13 <- filter(i_pred.df, time > T12 + 25, time < Finish - 25) M13D13 <- cbind(M13,D13) M13DT <- mutate(M13D13, detrend = CH4 - fit) head(M13DT) ### putting the detrended data together... Detrend_C <- rbind(M10DT, M12DT) Detrend_Cm <- add_column(Detrend_C, valve = 'chamber') Detrend_I <- rbind(M9DT, M11DT, M13DT) Detrend_Im <- add_column(Detrend_I, valve = 'inlet') head(Detrend_Im) Detrend_B <- rbind(Detrend_Cm,Detrend_Im) Plot_DT <- ggplot(Detrend_B) + geom_point(mapping=(aes(x = sec, y = detrend, colour = valve)))+ labs(title = "021219 CH4 detrended_Site 1B") Plot_DT + My_Theme Plot_DT + My_Theme Plot_DT + My_Theme + scale_colour_distiller(palette = 'Spectral') head(Detrend_B) # trying to show as a box chart - differences between sequences: ??brewer Box_plot_Detrend <- ggplot(data = Detrend_B, x = sec, y = detrend)+ geom_boxplot(data = M9DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+ geom_boxplot(data = M10DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+ geom_boxplot(data = M11DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+ geom_boxplot(data = M12DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+ geom_boxplot(data = M13DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+ ylim(-0.375,0.375)+ labs(title = "021219 Detrended CH4 S1B") Box_plot_Detrend + My_Theme c Box_plot_Detrend + My_Theme + scale_colour_manual(values = c("black", "blue")) ?scale_color_manual Cavg <- pull(Detrend_C, detrend) Cx = mean(Cavg) Csd <- sd(Cavg) Iavg <- pull(Detrend_I, detrend) Ix = mean(Iavg) Isd <- sd(Iavg) Cx Csd Ix Isd Cx-Ix #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX #Ambient concentration c CO2ambient <- pull(Detrend_I, CO2) CO2air = mean(CO2ambient) CO2sd <- sd(CO2ambient) CO2air CO2sd H2Oambient <- pull(Detrend_I, H2O) H2Oair = mean(H2Oambient) H2Osd <- sd(H2Oambient) H2Oair H2Osd #ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ #trying to determine whether the Chamber data lies above or below the trend in inlet CO2? #t.test for all chamber versus inlet data: TtestG <- t.test(Cavg, Iavg,alternative = c("greater"), mu = 0, paired = FALSE, var.equal = TRUE, conf.level = 0.95) t.test(Cavg, Iavg,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = TRUE, conf.level = 0.95) t.test(Cavg, Iavg,alternative = c("less"), mu = 0, paired = FALSE, var.equal = TRUE, conf.level = 0.95) #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX #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 V W <- 16 #CO2mass is ug/m3 #for CO2 Cx and Ix uits ppm = umol/mol Cx_mass <- W * Cx #ug/mol Cx_vol <- Cx_mass/V # ug/liter Cx_m3 <- Cx_vol*1000 #ug/m3 Cx_m3 Ix_mass <- W * Ix #ug/mol Ix_vol <- Ix_mass/V # ug/liter Ix_m3 <- Ix_vol*1000 #ug/m3 Ix_m3 #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX $CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC #Flow rate: #Add flow rates here: FR_Avg <- 7.9697 FR_sd <- 0.0357 Q1 <- FR_Avg #SLPM Qsd <- FR_sd Qcv <- 100*Qsd/Q1 #converting to (m3 min-1) Qm3 <- Q1/1000 Qm3 #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.54)/100 D_xm A <- pi*(D_xm/2)*(D_xm/2) A pi #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX #Calculating flux: #F = Q * (C-I) / A #ng/m2/min Flux <- Qm3 *(Cx_m3 - Ix_m3)/A Flux #ng/m2/min Flux_day <- Flux*24*60 Flux_day #ng/m2/d Flux_nmol_day <- Flux_day/16 Flux_nmol_day #nmol/m2/d #Listing the results Date <- '02Dec19 2200 to 3600 s' Location <- 'Site 1B' Gas <- "CH4" Vs <- c(Date, Location, Gas, Flux_day, Flux_nmol_day, Q1, Qsd, CH4air, CH4sd, H2Oair, H2Osd, Cx, Csd, Ix, Isd, T5, FR_Avg, FR_sd, Start, Finish) list(Vs) Parameters <- c('Date', 'Location', 'Gas', 'Flux/day (ng/m2)', 'Flux/day (nmol/m2)', 'Q1 (SLPM)', 'Qsd', 'CH4 (ppb)', 'CH4sd', 'H2Oair (ppm)', 'H2Osd', 'Cx (ppb)', 'Csd', 'Ix (ppb)', 'Isd', 'T5', 'Flow', 'Flow_sd','Start', 'Finish') Results <- data.frame('variables'= Parameters, 'Obs' = Vs) head(Results) View(Results) # extracting the coefficients and calcualting growth rates etc intercept <- round(summary(hyperfit)$coefficients) summary(hyperfit)$coefficients Intercept <- round(summary(hyperfit)$coefficients [1], 6) Intercept Intercept_SE <- round(summary(hyperfit)$coefficients [1,2], 6) Intercept_SE Max_U <- 1/Intercept Max_U Slope <- summary(hyperfit)$coefficient [2,1] Slope Net_U <- 1/(Intercept+Slope) Net_U summary() #ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ #T-tests between detrended CO2 data for the different modes: should maybe do a matrix comparison: #test 1 t.test(MD1, MD2,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 2 t.test(M4D, M6D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 3 t.test(M4D, M7D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 4 t.test(M4D, M8D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 5 t.test(M4D, M9D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 6 t.test(M4D, M10D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 7 t.test(M5D, M6D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 8 t.test(M5D, M7D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 9 t.test(M5D, M8D,alternative =c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 10 t.test(M5D, M9D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 11 t.test(M5D, M10D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 12 t.test(M6D, M7D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 13 t.test(M6D, M8D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 14 t.test(M6D, M9D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 15 t.test(M6D, M10D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 16 t.test(M7D, M8D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 17 t.test(M7D, M9D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 18 t.test(M7D, M10D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 19 t.test(M8D, M9D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 20 t.test(M8D, M10D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) #test 21 t.test(M9D, M10D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)