ARCHY 499 Data
Ema Bushnell
March 11, 2015
#read in all the necessary libraries
library(httr)
library(dplyr)##
## Attaching package: 'dplyr'
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(stringr)
library(aqp)## This is aqp 1.8
##
## Attaching package: 'aqp'
##
## The following object is masked from 'package:dplyr':
##
## slicelibrary(reshape2)
library(munsell)
library(rioja)## This is rioja 0.9-3#read in the main data set
all_data <- read.csv("/Users/ema/Desktop/sCOOL/DEPTH.csv")Organic Percentage Mass
#read in LOI data
LOI <- read.csv("/Users/ema/Desktop/sCOOL/LOI.csv")
LOI_sample_ID <- LOI$sampleID
M1.1 <- LOI$crucible.mass.run.1.average
M2.1 <- LOI$X105degree.including.crucible..run.1.average
M3.1 <- LOI$after430run1average
M4.1 <- as.numeric(as.character(LOI$after1000run1average))
M1.2 <- LOI$crucible.mass.run.2.average
M2.2 <- LOI$X105degree.including.crucible..run.2.average
M3.2 <- LOI$after430run2average
M4.2 <- LOI$after1000run2average
M1.3 <- LOI$crucible.mass.run.3.average
M2.3 <- LOI$X105degree.including.crucible..run.3.average
M3.3 <- LOI$after430run3average
M4.3 <- LOI$after1000run3average
org_per_mass1 <- 100 * ((M2.1-M3.1)/(M2.1-M1.1))
org_per_mass2 <- 100 * ((M2.2-M3.2)/(M2.2-M1.2))
org_per_mass3 <- 100 * ((M2.3-M3.3)/(M2.3-M1.3))
#Organize Outliers
# >6 probably wrong, <0 clearly wrong
LOI_stack <- stack(list(run1 = org_per_mass1,
run2 = org_per_mass2,
run3 = org_per_mass3))
LOI_stack$values <- with(LOI_stack, as.numeric(ifelse(values > 6 | values < 0, "NA", values)))
ggplot(LOI_stack, aes(ind, values)) +
ggtitle("Loss on Ignition - Organics")+
geom_violin() +
geom_jitter ()+
xlab("Run")+ ylab("Mass Lost")+
theme_minimal(base_size = 14) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
Stratigraphic Loss on Ignitiion- Organics
org_per_mass_df<- unstack(LOI_stack)
# if a row is all NA, delete it
ind <- apply(org_per_mass_df, 1, function(x) all(is.na(x)))
org_per_mass_df<- org_per_mass_df[ !ind, ]
org_per_mass_av <- rowMeans(org_per_mass_df, na.rm=TRUE)
# get depths for plot labels
org_per_mass_av_df <- data.frame(org_per_mass_av = org_per_mass_av,
depths = as.numeric(all_data$Depth.m.))
library("ggplot2")
ggplot(org_per_mass_av_df, aes(org_per_mass_av, depths)) +
ggtitle("Loss on Ignition- Organics")+
geom_point() +
scale_y_reverse() +
xlab("Mass Lost")+ ylab("Depth")+
theme_minimal(base_size = 14) 
Carbonate Percentage Mass
#%C2=100[(M3-M4)/(M3-M1)]
carb_per_mass1 <- 100 * ((M3.1-M4.1)/(M3.1-M1.1))
carb_per_mass2 <- 100 * ((M3.2-M4.2)/(M3.2-M1.2))
carb_per_mass3 <- 100 * ((M3.3-M4.3)/(M3.3-M1.3))
carboutlierID2 <- LOI[(carb_per_mass2 > 3 | carb_per_mass2 < 0), ]
LOI %>%
filter(carb_per_mass2 > 3 | carb_per_mass2 < 0)## sampleID run1a.total.mass..about.15g. run1b.total.mass..about.15g.
## 1 SWS190 15.4344 15.4314
## 2 SWS310 15.1590 15.1587
## 3 SWS315 14.9937 14.9844
## 4 SWS320 15.0938 15.0936
## 5 SWS325 15.0905 15.0908
## 6 SWS330 14.6958 14.6955
## 7 SWS335 14.9534 14.9541
## run1a.average.total.mass..about.15g. crucible.mass.run.1a
## 1 15.43290 8.0154
## 2 15.15885 8.6511
## 3 14.98905 8.3637
## 4 15.09370 8.3020
## 5 15.09065 7.6852
## 6 14.69565 7.8614
## 7 14.95375 8.0608
## crucible.mass.run.1b crucible.mass.run.1.average crucible.mass.run.2a
## 1 8.0155 8.01545 8.3019
## 2 8.6511 8.65110 8.3836
## 3 8.3634 8.36355 7.8903
## 4 8.3016 8.30180 8.9915
## 5 7.6853 7.68525 8.0322
## 6 7.8611 7.86125 7.8943
## 7 8.0607 8.06075 8.3715
## crucible.mass.run.2b crucible.mass.run.2.average crucible.mass.run.3a
## 1 8.3018 8.30185 8.5117
## 2 8.3847 8.38415 8.1840
## 3 7.8909 7.89060 7.9586
## 4 8.9925 8.99200 8.3599
## 5 8.0333 8.03275 7.4779
## 6 8.8941 8.39420 8.1907
## 7 8.3715 8.37150 7.7409
## crucible.mass.run.3b crucible.mass.run.3.average
## 1 8.5116 8.51165
## 2 8.1840 8.18400
## 3 7.9585 7.95855
## 4 8.3598 8.35985
## 5 7.4777 7.47780
## 6 8.1907 8.19070
## 7 7.7408 7.74085
## X105degree.including.crucible..run.1a
## 1 13.3994
## 2 13.7239
## 3 13.2466
## 4 13.1374
## 5 13.2800
## 6 12.9911
## 7 13.5578
## X105degree.including.crucible..run.1b
## 1 14.3999
## 2 13.7240
## 3 13.2467
## 4 13.1372
## 5 13.2799
## 6 13.9910
## 7 13.5583
## X105degree.including.crucible..run.1.average
## 1 13.89965
## 2 13.72395
## 3 13.24665
## 4 13.13730
## 5 13.27995
## 6 13.49105
## 7 13.55805
## X105degree.including.crucible..run.2a
## 1 13.0262
## 2 13.3780
## 3 13.0373
## 4 14.5908
## 5 12.7997
## 6 13.0438
## 7 13.4876
## X105degree.including.crucible..run.2b
## 1 13.0265
## 2 13.3783
## 3 13.0376
## 4 14.5908
## 5 12.8004
## 6 13.0442
## 7 13.4875
## X105degree.including.crucible..run.2.average
## 1 13.02635
## 2 13.37815
## 3 13.03745
## 4 14.59080
## 5 12.80005
## 6 13.04400
## 7 13.48755
## X105degree.including.crucible..run.3a
## 1 13.7255
## 2 13.2138
## 3 12.8104
## 4 12.9423
## 5 12.1486
## 6 12.5572
## 7 11.9789
## X105degree.including.crucible..run.3b
## 1 13.7256
## 2 13.2136
## 3 12.8103
## 4 12.9424
## 5 12.1488
## 6 12.5574
## 7 11.9789
## X105degree.including.crucible..run.3.average after430run1a after430run1b
## 1 13.72555 13.3754 13.37550
## 2 13.21370 13.7097 13.70950
## 3 12.81035 13.2326 13.23250
## 4 12.94235 13.1268 13.12677
## 5 12.14870 13.2723 13.27220
## 6 12.55730 12.9846 12.98440
## 7 11.97890 13.5516 13.55150
## after430run1average after430run2a after430run2b after430run2average
## 1 13.37545 12.9985 13.9984 13.49845
## 2 13.70960 13.3637 13.3637 13.36370
## 3 13.23255 13.0227 13.0226 13.02265
## 4 13.12678 14.5800 14.5801 14.58005
## 5 13.27225 12.7926 12.7926 12.79260
## 6 12.98450 13.0363 13.0363 13.03630
## 7 13.55155 13.4809 13.4810 13.48095
## after430run3a after430run3b after430run3average after1000run1a
## 1 13.6896 13.6896 13.68960 13.3532
## 2 13.1975 13.1975 13.19750 13.1784
## 3 12.7961 12.7961 12.79610 12.7769
## 4 12.9346 12.9345 12.93455 12.9217
## 5 12.1420 12.1420 12.14200 12.1297
## 6 12.5516 12.5516 12.55160 12.5398
## 7 11.9734 11.9733 11.97335 11.9619
## after1000run1b after1000run1average after1000run2a after1000run2b
## 1 13.3534 13.3533 12.9776 12.9775
## 2 13.1786 13.1785 12.8261 12.8262
## 3 12.7772 12.77705 12.8492 12.8494
## 4 12.9219 12.9218 13.3079 13.3084
## 5 12.1297 12.1297 13.5701 13.5708
## 6 12.5497 12.54475 13.1825 13.1829
## 7 11.9618 11.96185 13.0137 13.0138
## after1000run2average after1000run3a after1000run3b after1000run3average
## 1 12.97755 13.6643 13.6641 13.66420
## 2 12.82615 13.6926 13.6924 13.69250
## 3 12.84930 13.2160 13.2159 13.21595
## 4 13.30815 13.1118 13.1120 13.11190
## 5 13.57045 13.2585 13.2586 13.25855
## 6 13.18270 12.9716 12.9716 12.97160
## 7 13.01375 13.5376 13.5376 13.53760LOI_C__stack <- stack(list(run1 =carb_per_mass1,
run2 =carb_per_mass2,
run3 =carb_per_mass3))
# replace problem values with NA
LOI_C__stack$values <- with(LOI_C__stack, as.numeric(ifelse(values > 5 | values < 0, "NA", values)))
# plot to check again
ggplot(LOI_C__stack, aes(ind, values)) +
ggtitle("Loss on Ignition- Carbonate")+
geom_violin() +
theme_minimal(base_size = 14) +
geom_jitter()+
xlab("Run")+ ylab("Mass Lost")
# looks ok, let's get the averages of the runs...
carb_per_mass_df <- unstack(LOI_C__stack)
carb_per_mass_av <- rowMeans(carb_per_mass_df, na.rm = TRUE)
# get depths for plot labels
carb_per_mass_av_df <- data_frame(carb_per_mass_av = carb_per_mass_av[-length(carb_per_mass_av)],
depths = as.numeric(all_data$Depth.m.))
ggplot(carb_per_mass_av_df, aes(carb_per_mass_av, depths)) +
geom_point() +
ggtitle("Loss on Ignition- Carbonate")+
scale_y_reverse() +
xlab("Mass Lost")+ ylab("Depth") +
theme_minimal(base_size = 14) ## Warning: Removed 4 rows containing missing values (geom_point).
Colour
Munsell: hue: YR etc bit chroma: x value: x-hue y-value/z-chroma, aka 1Y 2/3
colour <- read.csv("/Users/ema/Desktop/sCOOL/test_data.csv")
#split value/chroma
Dry.Value <- sapply(strsplit(as.character(colour$Dry.num), "\\."), "[[", 1)
Dry.Chroma <- sapply(strsplit(as.character(colour$Dry.num), "\\."), "[[", 2)
Wet.Value <- sapply(strsplit(as.character(colour$Wet.num), "\\."), "[[", 1)
Wet.Chroma <- sapply(strsplit(as.character(colour$Wet.num), "\\."), "[[", 2)
Value_df <- data_frame(colour$Sample.ID, Dry.Value, Wet.Value)
Chroma_df <- data_frame(colour$Sample.ID, Dry.Chroma, Wet.Chroma)
Dry_df <- data_frame(colour$Sample.ID, colour$Dry.YR, Dry.Value, Dry.Chroma)
Wet_df <- data_frame(colour$Sample.ID, colour$Wet.YR, Wet.Value, Wet.Chroma)#Read in Munsell (as xYR x/x)
munsell <- read.csv ("/Users/ema/Desktop/sCOOL/munsell.csv")
munsell_df <- data_frame(sampleID= munsell$ID,
value= Dry.Value,
munsell= munsell$Dry..number,
depth= all_data$Depth.m.)
Dry_df <- data.frame(sampleID = all_data$Sample.ID,
dry_yr = colour$Dry.YR,
hue = paste0(colour$Dry.YR, "YR"),
dry_value = Dry.Value,
dry_chroma = Dry.Chroma,
munsell_colour = all_data$Dry..number,
stringsAsFactors = FALSE)
Chroma_df <- data.frame(all_data$Sample.ID, colour$Wet.YR, Wet.Value, Wet.Chroma)
Dry_df$munsell_colour_ <- gsub("10.5", "10", Dry_df$munsell_colour)
our_hex <- munsell2rgb(Dry_df$hue, as.numeric(Dry_df$dry_value), as.numeric(Dry_df$dry_chroma))
our_hex_no_na <- na.omit(our_hex)
#Plots a graph of dry munsell colours, with NA values removed
ggplot(munsell_df, aes(value, sampleID, colour= munsell)) +
scale_colour_manual(values=our_hex_no_na)+
theme_minimal() +
geom_point(size=4)
#pH and Electric Conductivity
ph_and_ec <- read.csv("/Users/ema/Desktop/sCOOL/ph.csv")
#ben's code to fix sheet
n <- nrow(ph_and_ec)
ph <- ph_and_ec[seq(1, n, 2), ]
ec <- ph_and_ec[seq(2, n, 2), ]
#make data frame
ec_df <- data_frame(ph$sample.ID, ec)
#plot these values like organic mass
ggplot(ec, aes(EC.result, all_data$Depth))+
geom_point() +
scale_y_reverse() +
xlab("Electric Conductivity")+
ylab("Depth (m)") +
theme_minimal(base_size = 14) 
ggplot(ph, aes(pH.result, all_data$Depth )) +
geom_point() +
scale_y_reverse() +
xlab("pH")+
ylab("Depth (m)") +
theme_minimal(base_size = 14) 
Magnetic Susceptibility
#read in Magnetic susceptibility data
MS<- read.csv("/Users/ema/Desktop/sCOOL/MS.csv")
#creates a data frame for each Fermion
hfdf<- data.frame(MS$sampleID, MS$HFAverage)
lfdf <- data.frame(MS$sampleID, MS$LFAverage)
#stack the two data frames together to label each value as "lf" or "hf"
ms_stack <- stack(list(hf= MS$HFAverage,
lf= MS$LFAverage))
#stack the IDs to match the stacked values
id_stack <- stack(list(depth= all_data$Depth.m,
depth= all_data$Depth.m))
#creates data frame with aligned depths, lf and hf values
ms_df <- data.frame(ms_stack, id_stack)
#Create a point graph with moth lf and hf magnetic values
ggplot(ms_df, aes(values.1, values, colour=ms_df$ind)) +
geom_point(size = 2) +
ylab("Depth") +
xlab("Magnetic Susceptibility") +
ggtitle("HF and LF: Stratigraphic Plot")+
scale_y_reverse()+
theme_minimal()
#Creates a density histogram
Fermion <- ms_df$ind
ggplot(data=ms_df, aes(x=ms_df$values, fill= Fermion )) +
xlab("Magnetic Susceptbility")+
ylab("Density") +
ggtitle("Frequency Distribution")+
geom_density(alpha=.3)
Stratigraphic Plot
PH<- ph$pH.result
EC <- ec$EC.result
Carbon_Loss <- carb_per_mass_av_df$carb_per_mass_av
Organic_Loss <- org_per_mass_av_df$org_per_mass_av
LF <- lfdf$MS.LFAverage
HF <- hfdf$MS.HFAverage
dry <- colour$Dry.num
wet <- gsub("25.1", "2.5", colour$Wet.num)
D <- all_data$Depth
total <- data.frame(PH,
EC,
Carbon_Loss,
Organic_Loss,
LF,
HF,
dry)
depth <- all_data$Depth
# examine the results and compare to original data
strat.plot(yvar=depth, ylabel= "Depth", cex.ylabel= 1,
cex.xlabel= 1 , srt.xlabel= 360,
title= "Madjebebe Section SWS- Stratigraphic Plot",
y.rev=TRUE,
total) 