In benmarwick/mjbnaturepaper: Code and data for analysis of the 2012 and 2015 MJB excavations
knitr::opts_chunk$set(echo = FALSE,
warning = FALSE,
message = FALSE,
quiet = TRUE,
progress=FALSE)
# Submission preparation guidelines:
# http://www.nature.com/nature/authors/gta/
# Nature’s standard figure sizes are 89 mm (single column) and 183 mm (double column) and the full depth of the page is 247 mm.
library(tidyverse)
library(ggrepel)
library(knitr)
invisible(purl("supplementary_information.Rmd", output = "temp", quiet=TRUE))
read_chunk("temp")
unlink("temp")
# only plot one point per artefact (some artefacts have multiple total station points)
library(viridis)
stone_artefacts_only_one <-
stone_artefacts_only %>%
group_by(Description, find) %>%
dplyr::summarise(Xnew_flipped = mean(Xnew_flipped),
depth_below_ground_surface = mean(depth_below_ground_surface))
# determined by plotting row C end levels
row_c <- c(2.4, 1.4, 0.4, -0.6, -1.6, -2.6, -3.6)
row_mids <- row_c/2
size = 0.25
p <- ggplot() +
geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "L", ],
aes(Xnew_flipped,
depth_below_ground_surface),
#colour = plasma(6)[1],
size = size) +
geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "HM", ],
aes(Xnew_flipped,
depth_below_ground_surface),
# colour = plasma(6)[2],
size = size) +
geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "GS", ],
aes(Xnew_flipped,
depth_below_ground_surface),
# colour = plasma(6)[3],
size = size) +
geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "AXE", ],
aes(Xnew_flipped,
depth_below_ground_surface),
# colour = plasma(6)[4],
size = size) +
geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "ART", ],
aes(Xnew_flipped,
depth_below_ground_surface),
# colour = "white",
size = size) +
geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "ART", ],
aes(Xnew_flipped,
depth_below_ground_surface),
# colour = plasma(6)[5],
size = size) +
geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "AF", ], # halo
aes(Xnew_flipped,
depth_below_ground_surface),
# colour = "white",
size = size) +
geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "AF", ],
aes(Xnew_flipped,
depth_below_ground_surface),
# colour = plasma(6)[6],
size = size) +
scale_y_reverse(limits = c(3,0)) +
theme_minimal() +
theme(panel.grid.major.x = element_line(colour = "black")) +
scale_x_continuous(breaks = row_c,
labels = NULL) +
xlab("") +
ylab("Depth below \nground surface (m)") +
scale_color_viridis(discrete=TRUE,
"Artefact\\ntype") +
guides(colour = guide_legend(override.aes = list(size = 5))) +
coord_equal()
row_c = c(2.4, 1.4, 0.4, -0.6, -1.6, -2.6, -3.6)
nums = paste0("B", 6:1)
row_mids <- row_c[-length(row_c)] + diff(row_c)/2
library(grid)
for(i in 1:length(row_mids)){
p = p + annotation_custom(grob = textGrob(nums[i], gp=gpar(fontsize=10)),
xmin = row_mids[i],
xmax = row_mids[i],
ymin = -8.5,
ymax = 2)
}
# Code to override clipping
grid.newpage()
gt <- ggplot_gtable(ggplot_build(p))
gt$layout$clip[gt$layout$name=="panel"] <- "off"
grid.draw(gt)
# output to RStudio plot pane, then save as SVG
library(readxl)
artefacts_per_litre <- read_excel("../../data/stone_artefact_data/Artefacts per litre discard BC4-6 by depth.xlsx")
C2_3 <- read_excel("../../data/stone_artefact_data/Artefacts per litre discard BC4-6 by depth.xlsx", sheet = 2)
names(C2_3) <- C2_3[1, ]
C2_3 <- C2_3[-1, ]
C2_3 <- map_df(C2_3, as.numeric)
C4 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 1:4]
names(C4) <- C4[1, ]
C4 <- C4[-1, ]
C4 <- map_df(C4, as.numeric)
B4 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 6:11]
names(B4) <- B4[1, ]
B4 <- B4[-1, ]
B4 <- map_df(B4, as.numeric)
B5 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 13:16]
names(B5) <- B5[1, ]
B5 <- B5[-1, ]
B5 <- map_df(B5, as.numeric)
B6 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 18:22]
names(B6) <- B6[1, ]
B6 <- B6[-1, ]
B6 <- map_df(B6, as.numeric)
C5 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 24:28]
names(C5) <- C5[1, ]
C5 <- C5[-1, ]
C5 <- map_df(C5, as.numeric)
C6 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 30:34]
names(C6) <- C6[1, ]
C6 <- C6[-1, ]
C6 <- map_df(C6, as.numeric)
# for E2 let's use plotted finds
E2_artefacts <- stone_artefacts_only %>%
filter(find == "L") %>%
filter(grepl("E2", .$Description)) %>%
mutate(Spit = as.numeric(gsub("//d", "", spit))) %>%
group_by(Spit) %>%
tally()
E2_depths <- # 42 rows
cleaned_rotated_points_in_main_excavation_area %>%
filter(grepl("E2", .$Description)) %>%
filter(grepl("EL", .$Description)) %>%
group_by(spit) %>%
summarise(m_d = mean(depth_below_ground_surface, na.rm = TRUE)) %>%
mutate(Spit = as.numeric(gsub("[A-B]", "", spit))) %>%
arrange(Spit) %>%
distinct(Spit, .keep_all = TRUE)
E2_depths$E2_depths_diffs <- abs(c(0, diff(E2_depths$m_d)))
E2_depths$vol <- E2_depths$E2_depths_diffs * 100 * 100 / 1000
E2 <-
E2_artefacts %>%
left_join(E2_depths) %>%
mutate(artefacts_per_cubic_m = n / vol) %>%
select(-spit, -n, -vol, -E2_depths_diffs) %>%
rename(`Artefacts/Litre` = artefacts_per_cubic_m) %>%
rename(depth = m_d)
# or let's use whatever Chris sent, not quite sure how he got that
E2 <- read_excel("../../data/stone_artefact_data/E2 7mm and plotted lithics per litre.xlsx")
# combine these into one df
the_list <- list(
"B4" = B4,
"B5" = B5,
"B6" = B6,
"C2/3" = C2_3,
"C4" = C4,
"C5" = C5,
"C6" = C6,
"E2" = E2
)
artefacts_per_litre_long <-
bind_rows(the_list,
.id = "id")
# Depth in two cols...
artefacts_per_litre_long$depth <-
with(artefacts_per_litre_long, ifelse(is.na(depth),
Depth, depth))
# make facets in specific order
order_we_want <- c("E2", "B4", "B5", "B6" ,
"C2/3", "C4", "C5", "C6")
artefacts_per_litre_long$id_f <-
factor(artefacts_per_litre_long$id,
levels= order_we_want)
# draw vertical lines to show phases
# for back we use 1.90 2.30
# for others we use 2.15-2.6
artefacts_per_litre_long <-
artefacts_per_litre_long %>%
mutate(phase_two_upper =
if_else(id_f %in% c("C4", "C5","C6", "B4", "B5", "B6"),
2.15, 1.9)) %>%
mutate(phase_two_lower =
if_else(id_f %in% c("C4", "C5","C6", "B4", "B5", "B6"),
2.6, 2.3))
ggplot(artefacts_per_litre_long,
aes(depth,
`Artefacts/Litre`)) +
geom_line() +
facet_wrap(~id_f, scales = "free_y", nrow = 2
) +
geom_vline(aes(xintercept = phase_two_upper),
colour = "red") +
geom_vline(aes(xintercept = phase_two_lower),
colour = "red") +
geom_text(aes(x = phase_two_upper,
y = -0.5,
label = phase_two_upper),
size = 2) +
geom_text(aes(x = phase_two_lower,
y = -0.5,
label = phase_two_lower),
size = 2) +
ylim(0, NA) +
theme_minimal() +
xlab("Depth below surface (m)") +
ggtitle("MJB artefact density by depth in select squares")
ggsave("../../figures/Artefact density by depth in squares B-C 4-5-6.png")
Show the dated samples from each square, by depth, so we can see the age constraints on each peak of artefacts'
# get OSL ages from published SI, to be sure we have the right ones
si_p_31 <- tabulizer::extract_tables("../../data/ages/Clarkson_Jacobs_Marwick_2017_SI.pdf", pages = 31, method = "data.frame")[[2]][-c(1:5), ]
si_p_32 <- tabulizer::extract_tables("../../data/ages/Clarkson_Jacobs_Marwick_2017_SI.pdf", pages = 32, method = "data.frame")[[2]][-c(1:3), ]
si_p_34 <- tabulizer::extract_tables("../../data/ages/Clarkson_Jacobs_Marwick_2017_SI.pdf", pages = 34,method = "data.frame")[[2]][-c(1:5), ]
# Want two cols: Sample, osl_age
si_p_31_a <-
si_p_31 %>%
mutate(Sample = gsub("\\s.*", "", X038.nature22968.RESEARCH)) %>%
mutate(depth = as.numeric(stringr::str_match(X038.nature22968.RESEARCH,
"^\\w*\\s(\\d.\\d*)")[,2])) %>%
mutate(osl_age = as.numeric(gsub("\\s.*", "", SUPPLEMENTARY.INFORMATION))) %>%
select(-X038.nature22968.RESEARCH, -SUPPLEMENTARY.INFORMATION) %>%
mutate(Sample = if_else(Sample == "38.3", "SW11A", Sample ))
si_p_31_a$depth[which(is.na(si_p_31_a$depth))] <- 2.28
# ensure that we use 51.7 for SW11A
si_p_31_a$osl_age[which(si_p_31_a$Sample == "SW11A")] <- 51.7
si_p_32_a <-
si_p_32 %>%
mutate(Sample = gsub("\\s.*", "", X038.nature22968)) %>%
mutate(osl_age = as.numeric(gsub("\\s.*", "", X.6))) %>%
mutate(depth = as.numeric(stringr::str_match(X038.nature22968,
"^\\w*\\s(\\d.\\d*)")[,2])) %>%
select(-X038.nature22968,
- RESEARCH.SUPPLEMENTARY.INFORMATION,
-starts_with("X")) %>%
mutate(Sample = ifelse(Sample == "", NA, Sample ))
# patch up the ages that have both CAM and MAM
si_p_32_a$Sample[c(15,17,18,20,24,26,28,30)] <-
c("SW3B", "SW3B",
"SW3A", "SW3A",
"NW3", "NW3",
"NW2", "NW2")
si_p_32_a$depth[which(is.na(si_p_32_a$depth))] <-
c(0.85, 0.85,
0.85, 0.85,
0.54, 0.54,
0.33, 0.33)
si_p_32_a <-
si_p_32_a %>%
filter(!is.na(osl_age))
si_p_34_a <-
si_p_34 %>%
mutate(Sample = gsub("\\s.*", "", X038.nature22968.RESEARCH)) %>%
mutate(osl_age = as.numeric(gsub("\\s.*", "", SUPPLEMENTARY.INFORMATION))) %>%
mutate(depth = -as.numeric(stringr::str_match(X038.nature22968.RESEARCH,
"^\\w*\\s(-\\d.\\d*)")[,2])) %>%
select(-X038.nature22968.RESEARCH, -SUPPLEMENTARY.INFORMATION)
si_osl_ages <-
rbind(si_p_31_a,
si_p_32_a,
si_p_34_a) %>%
arrange(osl_age)
# get depth data from the earlier version of osl_ages
osl_ages2 <-
osl_ages %>%
select(Sample, total_station_depth_below_surf) %>%
left_join(si_osl_ages)
# We need to combine OSL and C14, with a col for square, a col for depth, a col for sample ID, and a col for approx age
c14_ages_rugplot_data <-
c14_ages %>%
select(square,
depth_below_ground_surface,
Bchron_Median,
Lab.ID) %>%
dplyr::rename(grid_square = square,
age = Bchron_Median,
depth = depth_below_ground_surface,
id = Lab.ID) %>%
mutate(age = round(age / 1000, 1))
# For the OSL samples, what square was each sample taken from? We don
# have this in our data so far. We can look EDF 8 from the Nature article.
# NE, NW, SWA, SWB, SWC,
# E2, C4, B4, B5, B5,
osl_ages_rugplot_data <-
si_osl_ages %>%
mutate(grid_square = case_when(
grepl("SW8C", Sample) ~ "B6",
grepl("NE1B", Sample) ~ "C6",
grepl("NE", Sample) ~ "E2",
grepl("NW8B", Sample) ~ "B5",
grepl("NW9B", Sample) ~ "B5",
grepl("NW", Sample) ~ "C5",
grepl("SW.{1,2}A", Sample) ~ "B4",
grepl("SW.{1,2}B", Sample) ~ "B5",
grepl("SW.{1,2}C", Sample) ~ "B5",
grepl("KTL158", Sample) ~ "B4", # auger
grepl("KTL162", Sample) ~ "B4", # DEF30
grepl("KTL164", Sample) ~ "B5", # DEF30
grepl("KTL165", Sample) ~ "B5", # DEF30
is.na(Sample) ~ "Unknown"
)) %>%
select(grid_square,
#total_station_depth_below_surf,
depth,
osl_age,
Sample) %>%
dplyr::rename( # depth = total_station_depth_below_surf,
age = osl_age,
id = Sample) %>%
mutate( # depth = abs(depth),
age = as.numeric(age)) %>%
filter(grid_square != 'Unknown')
# combine the two sets of dates
date_for_plots <-
bind_rows(c14_ages_rugplot_data,
osl_ages_rugplot_data)
date_for_plots <-
date_for_plots %>%
dplyr::rename(id = grid_square,
sample_id = id) %>%
filter(id %in% artefacts_per_litre_long$id )
# find max y-axis value for each sq to plot age points
date_for_plots_max_y <-
artefacts_per_litre_long %>%
group_by(id) %>%
summarise(max_artefacts_litre = max(`Artefacts/Litre`, na.rm = TRUE))
# add max height for red lines
vlines <-
artefacts_per_litre_long %>%
select(id, phase_two_upper, phase_two_lower) %>%
gather(value, xintercept, -id) %>%
select(-value) %>%
distinct(.keep_all = TRUE)
vlines <-
vlines %>%
left_join(date_for_plots_max_y)
date_for_plots <-
date_for_plots %>%
left_join(date_for_plots_max_y) %>%
filter(age > 50) %>%
filter(age < 90) %>%
mutate( label = paste0(sample_id, ": ", age, " ka"))
# plot with ages
vlines$id_f <- factor(vlines$id, levels = order_we_want)
date_for_plots$id_f <- factor(date_for_plots$id, levels = order_we_want)
p <-
ggplot(artefacts_per_litre_long,
aes(depth,
`Artefacts/Litre`)) +
geom_line() +
geom_segment(data = vlines,
aes(x = xintercept,
xend = xintercept,
y = 0,
yend = max_artefacts_litre + 2),
colour = "red") +
geom_text(data = vlines,
aes(x = xintercept,
y = - max_artefacts_litre / 30,
label = xintercept),
size = 2) +
geom_point(data = date_for_plots,
aes(depth,
max_artefacts_litre + 2)) +
geom_text(data = date_for_plots,
aes(depth,
max_artefacts_litre + 2,
label = label),
size = 2,
angle = 45,
hjust = -0.25,
vjust = 0) +
facet_wrap(~id_f,
scales = "free_y",
nrow = 2
) +
theme_minimal() +
xlab("Depth below surface (m)") # +
# ggtitle("MJB artefact density by depth in select squares")
library(grid)
svg("../../figures/Artefact density by depth in select squares.svg",
width = 10, height = 7)
gt = ggplot_gtable(ggplot_build(p))
gt$layout$clip = "off"
grid.draw(gt)
dev.off()
# save, then tidy up in inkscape
benmarwick/mjbnaturepaper documentation built on May 9, 2022, 6:25 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE, quiet = TRUE, progress=FALSE) # Submission preparation guidelines: # http://www.nature.com/nature/authors/gta/ # Nature’s standard figure sizes are 89 mm (single column) and 183 mm (double column) and the full depth of the page is 247 mm. library(tidyverse) library(ggrepel)
library(knitr) invisible(purl("supplementary_information.Rmd", output = "temp", quiet=TRUE)) read_chunk("temp")
unlink("temp")
# only plot one point per artefact (some artefacts have multiple total station points) library(viridis) stone_artefacts_only_one <- stone_artefacts_only %>% group_by(Description, find) %>% dplyr::summarise(Xnew_flipped = mean(Xnew_flipped), depth_below_ground_surface = mean(depth_below_ground_surface)) # determined by plotting row C end levels row_c <- c(2.4, 1.4, 0.4, -0.6, -1.6, -2.6, -3.6) row_mids <- row_c/2 size = 0.25 p <- ggplot() + geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "L", ], aes(Xnew_flipped, depth_below_ground_surface), #colour = plasma(6)[1], size = size) + geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "HM", ], aes(Xnew_flipped, depth_below_ground_surface), # colour = plasma(6)[2], size = size) + geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "GS", ], aes(Xnew_flipped, depth_below_ground_surface), # colour = plasma(6)[3], size = size) + geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "AXE", ], aes(Xnew_flipped, depth_below_ground_surface), # colour = plasma(6)[4], size = size) + geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "ART", ], aes(Xnew_flipped, depth_below_ground_surface), # colour = "white", size = size) + geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "ART", ], aes(Xnew_flipped, depth_below_ground_surface), # colour = plasma(6)[5], size = size) + geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "AF", ], # halo aes(Xnew_flipped, depth_below_ground_surface), # colour = "white", size = size) + geom_point(data = stone_artefacts_only_one[stone_artefacts_only_one$find == "AF", ], aes(Xnew_flipped, depth_below_ground_surface), # colour = plasma(6)[6], size = size) + scale_y_reverse(limits = c(3,0)) + theme_minimal() + theme(panel.grid.major.x = element_line(colour = "black")) + scale_x_continuous(breaks = row_c, labels = NULL) + xlab("") + ylab("Depth below \nground surface (m)") + scale_color_viridis(discrete=TRUE, "Artefact\\ntype") + guides(colour = guide_legend(override.aes = list(size = 5))) + coord_equal() row_c = c(2.4, 1.4, 0.4, -0.6, -1.6, -2.6, -3.6) nums = paste0("B", 6:1) row_mids <- row_c[-length(row_c)] + diff(row_c)/2 library(grid) for(i in 1:length(row_mids)){ p = p + annotation_custom(grob = textGrob(nums[i], gp=gpar(fontsize=10)), xmin = row_mids[i], xmax = row_mids[i], ymin = -8.5, ymax = 2) } # Code to override clipping grid.newpage() gt <- ggplot_gtable(ggplot_build(p)) gt$layout$clip[gt$layout$name=="panel"] <- "off" grid.draw(gt) # output to RStudio plot pane, then save as SVG
library(readxl) artefacts_per_litre <- read_excel("../../data/stone_artefact_data/Artefacts per litre discard BC4-6 by depth.xlsx") C2_3 <- read_excel("../../data/stone_artefact_data/Artefacts per litre discard BC4-6 by depth.xlsx", sheet = 2) names(C2_3) <- C2_3[1, ] C2_3 <- C2_3[-1, ] C2_3 <- map_df(C2_3, as.numeric) C4 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 1:4] names(C4) <- C4[1, ] C4 <- C4[-1, ] C4 <- map_df(C4, as.numeric) B4 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 6:11] names(B4) <- B4[1, ] B4 <- B4[-1, ] B4 <- map_df(B4, as.numeric) B5 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 13:16] names(B5) <- B5[1, ] B5 <- B5[-1, ] B5 <- map_df(B5, as.numeric) B6 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 18:22] names(B6) <- B6[1, ] B6 <- B6[-1, ] B6 <- map_df(B6, as.numeric) C5 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 24:28] names(C5) <- C5[1, ] C5 <- C5[-1, ] C5 <- map_df(C5, as.numeric) C6 <- artefacts_per_litre[1:nrow(artefacts_per_litre) , 30:34] names(C6) <- C6[1, ] C6 <- C6[-1, ] C6 <- map_df(C6, as.numeric) # for E2 let's use plotted finds E2_artefacts <- stone_artefacts_only %>% filter(find == "L") %>% filter(grepl("E2", .$Description)) %>% mutate(Spit = as.numeric(gsub("//d", "", spit))) %>% group_by(Spit) %>% tally() E2_depths <- # 42 rows cleaned_rotated_points_in_main_excavation_area %>% filter(grepl("E2", .$Description)) %>% filter(grepl("EL", .$Description)) %>% group_by(spit) %>% summarise(m_d = mean(depth_below_ground_surface, na.rm = TRUE)) %>% mutate(Spit = as.numeric(gsub("[A-B]", "", spit))) %>% arrange(Spit) %>% distinct(Spit, .keep_all = TRUE) E2_depths$E2_depths_diffs <- abs(c(0, diff(E2_depths$m_d))) E2_depths$vol <- E2_depths$E2_depths_diffs * 100 * 100 / 1000 E2 <- E2_artefacts %>% left_join(E2_depths) %>% mutate(artefacts_per_cubic_m = n / vol) %>% select(-spit, -n, -vol, -E2_depths_diffs) %>% rename(`Artefacts/Litre` = artefacts_per_cubic_m) %>% rename(depth = m_d) # or let's use whatever Chris sent, not quite sure how he got that E2 <- read_excel("../../data/stone_artefact_data/E2 7mm and plotted lithics per litre.xlsx") # combine these into one df the_list <- list( "B4" = B4, "B5" = B5, "B6" = B6, "C2/3" = C2_3, "C4" = C4, "C5" = C5, "C6" = C6, "E2" = E2 ) artefacts_per_litre_long <- bind_rows(the_list, .id = "id") # Depth in two cols... artefacts_per_litre_long$depth <- with(artefacts_per_litre_long, ifelse(is.na(depth), Depth, depth)) # make facets in specific order order_we_want <- c("E2", "B4", "B5", "B6" , "C2/3", "C4", "C5", "C6") artefacts_per_litre_long$id_f <- factor(artefacts_per_litre_long$id, levels= order_we_want) # draw vertical lines to show phases # for back we use 1.90 2.30 # for others we use 2.15-2.6 artefacts_per_litre_long <- artefacts_per_litre_long %>% mutate(phase_two_upper = if_else(id_f %in% c("C4", "C5","C6", "B4", "B5", "B6"), 2.15, 1.9)) %>% mutate(phase_two_lower = if_else(id_f %in% c("C4", "C5","C6", "B4", "B5", "B6"), 2.6, 2.3)) ggplot(artefacts_per_litre_long, aes(depth, `Artefacts/Litre`)) + geom_line() + facet_wrap(~id_f, scales = "free_y", nrow = 2 ) + geom_vline(aes(xintercept = phase_two_upper), colour = "red") + geom_vline(aes(xintercept = phase_two_lower), colour = "red") + geom_text(aes(x = phase_two_upper, y = -0.5, label = phase_two_upper), size = 2) + geom_text(aes(x = phase_two_lower, y = -0.5, label = phase_two_lower), size = 2) + ylim(0, NA) + theme_minimal() + xlab("Depth below surface (m)") + ggtitle("MJB artefact density by depth in select squares") ggsave("../../figures/Artefact density by depth in squares B-C 4-5-6.png")
Show the dated samples from each square, by depth, so we can see the age constraints on each peak of artefacts'
# get OSL ages from published SI, to be sure we have the right ones si_p_31 <- tabulizer::extract_tables("../../data/ages/Clarkson_Jacobs_Marwick_2017_SI.pdf", pages = 31, method = "data.frame")[[2]][-c(1:5), ] si_p_32 <- tabulizer::extract_tables("../../data/ages/Clarkson_Jacobs_Marwick_2017_SI.pdf", pages = 32, method = "data.frame")[[2]][-c(1:3), ] si_p_34 <- tabulizer::extract_tables("../../data/ages/Clarkson_Jacobs_Marwick_2017_SI.pdf", pages = 34,method = "data.frame")[[2]][-c(1:5), ] # Want two cols: Sample, osl_age si_p_31_a <- si_p_31 %>% mutate(Sample = gsub("\\s.*", "", X038.nature22968.RESEARCH)) %>% mutate(depth = as.numeric(stringr::str_match(X038.nature22968.RESEARCH, "^\\w*\\s(\\d.\\d*)")[,2])) %>% mutate(osl_age = as.numeric(gsub("\\s.*", "", SUPPLEMENTARY.INFORMATION))) %>% select(-X038.nature22968.RESEARCH, -SUPPLEMENTARY.INFORMATION) %>% mutate(Sample = if_else(Sample == "38.3", "SW11A", Sample )) si_p_31_a$depth[which(is.na(si_p_31_a$depth))] <- 2.28 # ensure that we use 51.7 for SW11A si_p_31_a$osl_age[which(si_p_31_a$Sample == "SW11A")] <- 51.7 si_p_32_a <- si_p_32 %>% mutate(Sample = gsub("\\s.*", "", X038.nature22968)) %>% mutate(osl_age = as.numeric(gsub("\\s.*", "", X.6))) %>% mutate(depth = as.numeric(stringr::str_match(X038.nature22968, "^\\w*\\s(\\d.\\d*)")[,2])) %>% select(-X038.nature22968, - RESEARCH.SUPPLEMENTARY.INFORMATION, -starts_with("X")) %>% mutate(Sample = ifelse(Sample == "", NA, Sample )) # patch up the ages that have both CAM and MAM si_p_32_a$Sample[c(15,17,18,20,24,26,28,30)] <- c("SW3B", "SW3B", "SW3A", "SW3A", "NW3", "NW3", "NW2", "NW2") si_p_32_a$depth[which(is.na(si_p_32_a$depth))] <- c(0.85, 0.85, 0.85, 0.85, 0.54, 0.54, 0.33, 0.33) si_p_32_a <- si_p_32_a %>% filter(!is.na(osl_age)) si_p_34_a <- si_p_34 %>% mutate(Sample = gsub("\\s.*", "", X038.nature22968.RESEARCH)) %>% mutate(osl_age = as.numeric(gsub("\\s.*", "", SUPPLEMENTARY.INFORMATION))) %>% mutate(depth = -as.numeric(stringr::str_match(X038.nature22968.RESEARCH, "^\\w*\\s(-\\d.\\d*)")[,2])) %>% select(-X038.nature22968.RESEARCH, -SUPPLEMENTARY.INFORMATION) si_osl_ages <- rbind(si_p_31_a, si_p_32_a, si_p_34_a) %>% arrange(osl_age) # get depth data from the earlier version of osl_ages osl_ages2 <- osl_ages %>% select(Sample, total_station_depth_below_surf) %>% left_join(si_osl_ages)
# We need to combine OSL and C14, with a col for square, a col for depth, a col for sample ID, and a col for approx age c14_ages_rugplot_data <- c14_ages %>% select(square, depth_below_ground_surface, Bchron_Median, Lab.ID) %>% dplyr::rename(grid_square = square, age = Bchron_Median, depth = depth_below_ground_surface, id = Lab.ID) %>% mutate(age = round(age / 1000, 1)) # For the OSL samples, what square was each sample taken from? We don # have this in our data so far. We can look EDF 8 from the Nature article. # NE, NW, SWA, SWB, SWC, # E2, C4, B4, B5, B5, osl_ages_rugplot_data <- si_osl_ages %>% mutate(grid_square = case_when( grepl("SW8C", Sample) ~ "B6", grepl("NE1B", Sample) ~ "C6", grepl("NE", Sample) ~ "E2", grepl("NW8B", Sample) ~ "B5", grepl("NW9B", Sample) ~ "B5", grepl("NW", Sample) ~ "C5", grepl("SW.{1,2}A", Sample) ~ "B4", grepl("SW.{1,2}B", Sample) ~ "B5", grepl("SW.{1,2}C", Sample) ~ "B5", grepl("KTL158", Sample) ~ "B4", # auger grepl("KTL162", Sample) ~ "B4", # DEF30 grepl("KTL164", Sample) ~ "B5", # DEF30 grepl("KTL165", Sample) ~ "B5", # DEF30 is.na(Sample) ~ "Unknown" )) %>% select(grid_square, #total_station_depth_below_surf, depth, osl_age, Sample) %>% dplyr::rename( # depth = total_station_depth_below_surf, age = osl_age, id = Sample) %>% mutate( # depth = abs(depth), age = as.numeric(age)) %>% filter(grid_square != 'Unknown') # combine the two sets of dates date_for_plots <- bind_rows(c14_ages_rugplot_data, osl_ages_rugplot_data) date_for_plots <- date_for_plots %>% dplyr::rename(id = grid_square, sample_id = id) %>% filter(id %in% artefacts_per_litre_long$id ) # find max y-axis value for each sq to plot age points date_for_plots_max_y <- artefacts_per_litre_long %>% group_by(id) %>% summarise(max_artefacts_litre = max(`Artefacts/Litre`, na.rm = TRUE)) # add max height for red lines vlines <- artefacts_per_litre_long %>% select(id, phase_two_upper, phase_two_lower) %>% gather(value, xintercept, -id) %>% select(-value) %>% distinct(.keep_all = TRUE) vlines <- vlines %>% left_join(date_for_plots_max_y) date_for_plots <- date_for_plots %>% left_join(date_for_plots_max_y) %>% filter(age > 50) %>% filter(age < 90) %>% mutate( label = paste0(sample_id, ": ", age, " ka")) # plot with ages vlines$id_f <- factor(vlines$id, levels = order_we_want) date_for_plots$id_f <- factor(date_for_plots$id, levels = order_we_want) p <- ggplot(artefacts_per_litre_long, aes(depth, `Artefacts/Litre`)) + geom_line() + geom_segment(data = vlines, aes(x = xintercept, xend = xintercept, y = 0, yend = max_artefacts_litre + 2), colour = "red") + geom_text(data = vlines, aes(x = xintercept, y = - max_artefacts_litre / 30, label = xintercept), size = 2) + geom_point(data = date_for_plots, aes(depth, max_artefacts_litre + 2)) + geom_text(data = date_for_plots, aes(depth, max_artefacts_litre + 2, label = label), size = 2, angle = 45, hjust = -0.25, vjust = 0) + facet_wrap(~id_f, scales = "free_y", nrow = 2 ) + theme_minimal() + xlab("Depth below surface (m)") # + # ggtitle("MJB artefact density by depth in select squares") library(grid) svg("../../figures/Artefact density by depth in select squares.svg", width = 10, height = 7) gt = ggplot_gtable(ggplot_build(p)) gt$layout$clip = "off" grid.draw(gt) dev.off() # save, then tidy up in inkscape
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.