Bringing in and processing burn chamber data
library packages
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)From raw data
#### Process for cleaning and working up thermocouple data
#Georgia Harrison, Jan 2021
#### before importing into R
#there may be something with package parsedate to do this in R but I do not know?
#Rename all columns without spaces
#naming system: Ch#_Stat
## i.e. Ch1_Avg
#lets separate time into date and time columns
#Add "Timestamp" title to first column
#also make this column TEXT formatted
#in excel, highlight column -> data Text to columns -> fixed width
#create lines for date . T . time . time zone
#label as "Date" " " "Time" "Tmz"
#format date and time columns as such
#delete the T column in between date and time
#remove no data lines in excel
#sort by any of the thermos, remove blank rows
#resort in order by time
#mass
#make sure trial ID matches temp
#trial ID:
#Species: ALWAYS LIST AG first
#mix - follow order of species
#Need to calculate a residual mass in mass data set
#(pre/post)*100Data Input
#temp<- read.csv("C:/Users/harr4718/OneDrive - University of Idaho/Harrison Grass Combustion/Thermocouple and mass data/11jan2021_temp.csv", header=TRUE)
mass<- read.csv("C:/Users/harr4718/OneDrive - University of Idaho/Harrison Grass Combustion/Thermocouple and mass data/23Feb2021_mass.csv", header=TRUE)
temp<- read.csv("C:/Users/harr4718/OneDrive - University of Idaho/Harrison Grass Combustion/Thermocouple and mass data/feb23_temp.csv", header=TRUE)
temp$Trial_ID <- as.factor(temp$Trial_ID) #convert from character to factor
mass$Trial_ID <- as.factor(mass$Trial_ID) #convert from character to factorCreate new colums for all the outputs
mass["dur_ch1"] <- NA
mass["dur_ch2"] <- NA
mass["dur_ch3"] <- NA
mass["dur_ch4"] <- NA
mass["dur_ch5"] <- NA
mass["dur_ch6"] <- NA
mass["max_ch1"] <- NA
mass["max_ch2"] <- NA
mass["max_ch3"] <- NA
mass["max_ch4"] <- NA
mass["max_ch5"] <- NA
mass["max_ch6"] <- NA
mass["area_ch1"] <- NA
mass["area_ch2"] <- NA
mass["area_ch3"] <- NA
mass["area_ch4"] <- NA
mass["area_ch5"] <- NA
mass["area_ch6"] <- NAseperate by trial ID
summary(mass$Trial_ID)## A B C D E F G H
## 1 1 1 1 1 1 1 1A <- filter(temp, Trial_ID == "A")
B <- filter(temp, Trial_ID == "B")
C <- filter(temp, Trial_ID == "C")
D <- filter(temp, Trial_ID == "D")
E <- filter(temp, Trial_ID == "E")
F <- filter(temp, Trial_ID == "F")
G <- filter(temp, Trial_ID == "G")
H <- filter(temp, Trial_ID == "H")Duration of fire
tempthreshold<-60 #temperature you are considering combustion
timestamp<-1 #internal of each timestamp in seconds
current <- B ## CHANGE THIS EVERY TIME ##
which(mass$Trial_ID=="B")## [1] 2one <- subset(current, Ch1_Avg > tempthreshold)
two <- subset(current, Ch2_Avg > tempthreshold)
three <- subset(current, Ch3_Avg > tempthreshold)
four <- subset(current, Ch4_Avg > tempthreshold)
five <- subset(current, Ch5_Avg > tempthreshold)
six <- subset(current, Ch6_Avg > tempthreshold)
which(colnames(mass)=="dur_ch1")## [1] 10which(colnames(mass)=="dur_ch6")## [1] 15mass[2, 10] <- nrow((one))*timestamp
mass[2, 11] <- nrow((two))*timestamp
mass[2, 12] <- nrow((three))*timestamp
mass[2, 13] <- nrow((four))*timestamp
mass[2, 14] <- nrow((five))*timestamp
mass[2, 15] <- nrow((six))*timestamp
mass## Trial_ID Rep Part Date PSSP BRTE Pre Post
## 1 A 1 2 2/23/2021 50 50 40.6 5.4
## 2 B 2 2 2/23/2021 25 75 34.8 6.0
## 3 C 3 2 2/23/2021 75 25 41.1 7.2
## 4 D 4 2 2/23/2021 100 0 41.9 6.7
## 5 E 5 2 2/23/2021 0 100 27.9 4.5
## 6 F 1 1 2/23/2021 25 75 20.6 3.4
## 7 G 2 1 2/23/2021 0 100 19.9 3.2
## 8 H 3 1 2/23/2021 100 0 39.6 6.0
## notes dur_ch1 dur_ch2 dur_ch3
## 1 weird themocouple reading -started reading mid fire? NA NA NA
## 2 0 0 0
## 3 re-do from earlier, in pt 1 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## 7 NA NA NA
## 8 NA NA NA
## dur_ch4 dur_ch5 dur_ch6 max_ch1 max_ch2 max_ch3 max_ch4 max_ch5 max_ch6
## 1 NA NA NA NA NA NA NA NA NA
## 2 417 291 150 NA NA NA NA NA NA
## 3 NA NA NA NA NA NA NA NA NA
## 4 NA NA NA NA NA NA NA NA NA
## 5 NA NA NA NA NA NA NA NA NA
## 6 NA NA NA NA NA NA NA NA NA
## 7 NA NA NA NA NA NA NA NA NA
## 8 NA NA NA NA NA NA NA NA NA
## area_ch1 area_ch2 area_ch3 area_ch4 area_ch5 area_ch6
## 1 NA NA NA NA NA NA
## 2 NA NA NA NA NA NA
## 3 NA NA NA NA NA NA
## 4 NA NA NA NA NA NA
## 5 NA NA NA NA NA NA
## 6 NA NA NA NA NA NA
## 7 NA NA NA NA NA NA
## 8 NA NA NA NA NA NAMax temp
which(colnames(mass)=="max_ch1") #16## [1] 16which(colnames(mass)=="max_ch6") #21## [1] 21mass[2, 16] <- max(one$Ch1_Avg)## Warning in max(one$Ch1_Avg): no non-missing arguments to max; returning -Infmass[2, 17] <- max(two$Ch2_Avg)## Warning in max(two$Ch2_Avg): no non-missing arguments to max; returning -Infmass[2, 18] <- max(three$Ch3_Avg)## Warning in max(three$Ch3_Avg): no non-missing arguments to max; returning -Infmass[2, 19] <- max(four$Ch4_Avg)
mass[2, 20] <- max(five$Ch5_Avg)
mass[2, 21] <- max(six$Ch6_Avg)
head(mass)## Trial_ID Rep Part Date PSSP BRTE Pre Post
## 1 A 1 2 2/23/2021 50 50 40.6 5.4
## 2 B 2 2 2/23/2021 25 75 34.8 6.0
## 3 C 3 2 2/23/2021 75 25 41.1 7.2
## 4 D 4 2 2/23/2021 100 0 41.9 6.7
## 5 E 5 2 2/23/2021 0 100 27.9 4.5
## 6 F 1 1 2/23/2021 25 75 20.6 3.4
## notes dur_ch1 dur_ch2 dur_ch3
## 1 weird themocouple reading -started reading mid fire? NA NA NA
## 2 0 0 0
## 3 re-do from earlier, in pt 1 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## dur_ch4 dur_ch5 dur_ch6 max_ch1 max_ch2 max_ch3 max_ch4 max_ch5 max_ch6
## 1 NA NA NA NA NA NA NA NA NA
## 2 417 291 150 -Inf -Inf -Inf 784.392 1280.636 1156.841
## 3 NA NA NA NA NA NA NA NA NA
## 4 NA NA NA NA NA NA NA NA NA
## 5 NA NA NA NA NA NA NA NA NA
## 6 NA NA NA NA NA NA NA NA NA
## area_ch1 area_ch2 area_ch3 area_ch4 area_ch5 area_ch6
## 1 NA NA NA NA NA NA
## 2 NA NA NA NA NA NA
## 3 NA NA NA NA NA NA
## 4 NA NA NA NA NA NA
## 5 NA NA NA NA NA NA
## 6 NA NA NA NA NA NAArea under the curve
Since timestamp is only 1 second, just multiply by 1. If timestep was different, multiply the sum term by this (in s).
which(colnames(mass)=="area_ch1") #22## [1] 22which(colnames(mass)=="area_ch6") #27## [1] 27mass[2, 22]<- sum(one$Ch1_Avg)
mass[2, 23]<- sum(two$Ch2_Avg)
mass[2, 24]<- sum(three$Ch3_Avg)
mass[2, 25]<- sum(four$Ch4_Avg)
mass[2, 26]<- sum(five$Ch5_Avg)
mass[2, 27]<- sum(six$Ch6_Avg)summary(mass$Trial_ID)## A B C D E F G H
## 1 1 1 1 1 1 1 1FOR BULK PROCESSING
current <- H ## CHANGE THIS EVERY TIME ##
which(mass$Trial_ID=="H")## [1] 8one <- subset(current, Ch1_Avg > tempthreshold)
two <- subset(current, Ch2_Avg > tempthreshold)
three <- subset(current, Ch3_Avg > tempthreshold)
four <- subset(current, Ch4_Avg > tempthreshold)
five <- subset(current, Ch5_Avg > tempthreshold)
six <- subset(current, Ch6_Avg > tempthreshold)
mass[8, 10] <- nrow((one))*timestamp
mass[8, 11] <- nrow((two))*timestamp
mass[8, 12] <- nrow((three))*timestamp
mass[8, 13] <- nrow((four))*timestamp
mass[8, 14] <- nrow((five))*timestamp
mass[8, 15] <- nrow((six))*timestamp
## Max temp
mass[8, 16] <- max(one$Ch1_Avg)
mass[8, 17] <- max(two$Ch2_Avg)## Warning in max(two$Ch2_Avg): no non-missing arguments to max; returning -Infmass[8, 18] <- max(three$Ch3_Avg)
mass[8, 19] <- max(four$Ch4_Avg)
mass[8, 20] <- max(five$Ch5_Avg)
mass[8, 21] <- max(six$Ch6_Avg)## Warning in max(six$Ch6_Avg): no non-missing arguments to max; returning -Inf## Area
mass[8, 22]<- sum(one$Ch1_Avg)
mass[8, 23]<- sum(two$Ch2_Avg)
mass[8, 24]<- sum(three$Ch3_Avg)
mass[8, 25]<- sum(four$Ch4_Avg)
mass[8, 26]<- sum(five$Ch5_Avg)
mass[8, 27]<- sum(six$Ch6_Avg)Convert data from wide to long - compress channels Make sure spelling and capitalization of col headings are correct!
head(mass)## Trial_ID Rep Part Date PSSP BRTE Pre Post
## 1 A 1 2 2/23/2021 50 50 40.6 5.4
## 2 B 2 2 2/23/2021 25 75 34.8 6.0
## 3 C 3 2 2/23/2021 75 25 41.1 7.2
## 4 D 4 2 2/23/2021 100 0 41.9 6.7
## 5 E 5 2 2/23/2021 0 100 27.9 4.5
## 6 F 1 1 2/23/2021 25 75 20.6 3.4
## notes dur_ch1 dur_ch2 dur_ch3
## 1 weird themocouple reading -started reading mid fire? NA NA NA
## 2 0 0 0
## 3 re-do from earlier, in pt 1 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## dur_ch4 dur_ch5 dur_ch6 max_ch1 max_ch2 max_ch3 max_ch4 max_ch5 max_ch6
## 1 NA NA NA NA NA NA NA NA NA
## 2 417 291 150 -Inf -Inf -Inf 784.392 1280.636 1156.841
## 3 NA NA NA NA NA NA NA NA NA
## 4 NA NA NA NA NA NA NA NA NA
## 5 NA NA NA NA NA NA NA NA NA
## 6 NA NA NA NA NA NA NA NA NA
## area_ch1 area_ch2 area_ch3 area_ch4 area_ch5 area_ch6
## 1 NA NA NA NA NA NA
## 2 0 0 0 114949.2 126310.1 70907.78
## 3 NA NA NA NA NA NA
## 4 NA NA NA NA NA NA
## 5 NA NA NA NA NA NA
## 6 NA NA NA NA NA NAdur <- c("dur_ch1", "dur_ch2", "dur_ch3", "dur_ch4", "dur_ch5", "dur_ch6")
max <- c("max_ch1", "max_ch2", "max_ch3", "max_ch4", "max_ch5", "max_ch6")
area <- c("area_ch1", "area_ch2", "area_ch3", "area_ch4", "area_ch5", "area_ch6")
dur <- pivot_longer(mass, cols = c("dur_ch1", "dur_ch2", "dur_ch3", "dur_ch4", "dur_ch5", "dur_ch6"), names_to = c("D", "channel"), names_sep = "_", values_to="duration")
max <- pivot_longer(mass, cols = c("max_ch1", "max_ch2", "max_ch3", "max_ch4", "max_ch5", "max_ch6"), names_to = c("M", "channel"), names_sep = "_", values_to="max")
area <- pivot_longer(mass, cols = c("area_ch1", "area_ch2", "area_ch3", "area_ch4", "area_ch5", "area_ch6"), names_to = c("A", "channel"), names_sep = "_", values_to="area")create an output data frame
#remmove the columns we just combined
which(colnames(dur)=="notes") #9## [1] 9which(colnames(dur)=="D")#22## [1] 22dur1 <- subset(dur, select = -c(9:22))
max1 <- subset(max, select = -c(9:22))
area1 <- subset(area, select = -c(9:22))
#merge duration, max, and area dataframes
cols <- c("Trial_ID", "Rep", "Part", "Date", "PSSP", "BRTE", "Pre", "Post", "channel")
output <- merge (dur1, max1, by.x=cols, by.y=cols)
output <- merge (output, area1, by.x=cols, by.y=cols)
head(output)## Trial_ID Rep Part Date PSSP BRTE Pre Post channel duration max area
## 1 A 1 2 2/23/2021 50 50 40.6 5.4 ch1 NA NA NA
## 2 A 1 2 2/23/2021 50 50 40.6 5.4 ch2 NA NA NA
## 3 A 1 2 2/23/2021 50 50 40.6 5.4 ch3 NA NA NA
## 4 A 1 2 2/23/2021 50 50 40.6 5.4 ch4 NA NA NA
## 5 A 1 2 2/23/2021 50 50 40.6 5.4 ch5 NA NA NA
## 6 A 1 2 2/23/2021 50 50 40.6 5.4 ch6 NA NA NABasic plot will make better later
plot.ts(five$Ch5_Avg, xlab= "time", ylab="Temperature", main="chanel 5")
gg plots
Stil figuring these out
library(ggplot2)
p <- ggplot(five, aes(x=Time, y=Ch5_Avg))
p <- p + geom_point()+theme_classic()
p <- p + labs(y="Temperature (C)", x="Time")
p #only the time over threshold
p <- ggplot(p5, aes(x=Time, y=Ch2_Avg))
p <- p + geom_point()+theme_classic()
p <- p + labs(y="Temperature (C)", x="Time")
p #all the temp data for that trial
pa <- p + scale_fill_grey(name="Climbing", labels = c("Climbed", "Unclimbed"))
plot(pa)