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.

Introduction

This document contains the R code used to generate the figures and statistical results that I contributed to

Clarkson, C., Z. Jacobs, B. Marwick, R. Fullagar, L. Wallis, M. Smith, R.G. Roberts, E. Hayes, K. Lowe, X. Carah, S.A. Florin, J. McNeil, D. Cox, L.J. Arnold, Q. Hua, J. Huntley, H.E.A. Brand, T. Manne, A. Fairbairn, J. Shulmeister, L. Lyle, M. Salinas, M. Page, K. Connell, G. Park, K. Norman, T. Murphy and C. Pardoe 2017 Human occupation of northern Australia by 65,000 years ago. Nature 547:306–310.

# prepare data, if interactive, setwd to source file location

# before running library(mjbnaturepaper):

# remotes::install_github("matherion/userfriendlyscience")
# BiocManager::install("IRanges")
# remotes::install_github("benmarwick/mjbnaturepaper")

library(mjbnaturepaper)

library(tidyverse)

# read in from raw data, setwd to source file location
total_station_data_2015 <- read_2015_total_station_data("data/total_station_csv_2015")
total_station_data_2012 <- read_2012_total_station_data("data/total_station_csv_2012")

# combine both years
combined_2012_and_2015_total_station_data <- combine_2012_and_2015_total_station_data(total_station_data_2012,
                                         total_station_data_2015)
# exclude points not in excavation
points_in_main_excavation_area <- 
extract_points_in_main_excavation_area(combined_2012_and_2015_total_station_data)

# rotate points to make the excavation square
rotated_points_in_main_excavation_area <- 
rotate_points_in_main_excavation_area(points_in_main_excavation_area)

# clean out spurious and duplicated points and non-artefacts
cleaned_rotated_points_in_main_excavation_area <- 
clean_points_in_main_excavation_area(rotated_points_in_main_excavation_area)

# extract only stone artefacts and related things
stone_artefacts_only <- 
stone_artefacts(cleaned_rotated_points_in_main_excavation_area)


saveRDS(cleaned_rotated_points_in_main_excavation_area,
        "data/cleaned_rotated_points_in_main_excavation_area.rds")

Main paper figures and tables

This document the code to reproduce the figures in the main paper that I'm involved in.

Extended Data Figure 4

The artefacts occur in three dense bands. Detailed piece plotting revealed three dense bands of artefacts at the site with less dense bands of artefacts in between.

Each band corresponds to a change in raw material use and stone working technology.

# section scatter plot, artefacts in all squares
plot_stone_artefacts(stone_artefacts_only)

# section scatter plot, just artefacts in B and C
stone_artefacts_only_B_C <- stone_artefacts_only[grep("B|C", stone_artefacts_only$square), ]
plot_stone_artefacts(stone_artefacts_only_B_C)

# vertical bar plot of raw materials and technology, B6 only
spit_depths_B6_output <- spit_depths_B6(points_in_main_excavation_area)

B6_raw_materials <- readr::read_csv("data/stone_artefact_data/B6_raw_material_table.csv")

plot_raw_materials_technology <- 
raw_materials_technology_plots(B6_raw_materials, spit_depths_B6_output)

# chi-sq for raw material by phase
chi_sq_raw_material_by_phase_tbl <- chi_sq_raw_material_by_phase(plot_raw_materials_technology)

# this creates the table as originally published, which we now know is wrong, and
# was corrected in https://doi.org/10.1080/03122417.2018.1462884
chisq.test(chi_sq_raw_material_by_phase_tbl)

# These plots are arranged in inkscape to create the final figure EDF-4

Extended Data Figure 8

Associations were found between 19 artefacts (lower = 14, middle = 3) indicating refits and conjoins within the lower two dense artefact layers but not between them. Although artefacts increase in mass slightly with depth overall, change in depth explains only 0.20% of changes in artefact mass. If we look into the upper and lower bands of artefacts separately, we see that ...

# Refits -------------------------------------

# plot refit plan and section views, and rose plot of dip
# check aspect ratio of these plots
refit_output <- refits(stone_artefacts_only)

# median distance between refitting artefact fragments
median_refit_dist <- 
median(refit_output$refit_data_coords_wide$refit_dists_m)
mad_refit_dist <- 
mad(refit_output$refit_data_coords_wide$refit_dists_m)


# Vertical distance between refitted pieces
refit_vertical_dist_median <- 
median(abs(with(refit_output$refit_data_coords_wide, Elevation.x  -  Elevation.y)))

refit_vertical_dist_mad <- 
mad(abs(with(refit_output$refit_data_coords_wide, Elevation.x  -  Elevation.y)))
mad_refit_dist
# Size sorting ---------------------------------------------

# size sorting plot with smoother
size_sorting_plot(spit_depths_B6_output)

Extended Data Figure 9

Each pulse in artifact discard corresponds to a peak in mass-normalized magnetic susceptibility, increasing in magnitude towards the surface, and a peak in the quantity of burnt chert in the lower two bands (chert is absent in the upper band), indicating artifact pulses are probably associated with phases of intense anthropogenic burning in the shelter

prepared_geoarchaeology_data <-  
  prepare_geoarchaeology_data(cleaned_rotated_points_in_main_excavation_area)

# grain size and cluster plot is separate
plot_granulometry_data(prepared_geoarchaeology_data)

plot_geoarchaeology_data(prepared_geoarchaeology_data)

# import PDFs into EDF SVG and resize to fit

Extended Data Figure 10

Micromorphology of the excavated deposits shows aeolian and alluvial origins of predominantly lightly clay-coated rounded quartz grains, with less-weathered disaggregation and slow exfoliation from the rockshelter wall sandstone

Extended Data Figure XX

depth age and age of each artefact band, as well as first occupation

c14_ages <- get_c14_ages(cleaned_rotated_points_in_main_excavation_area)
osl_ages <- get_osl_ages(cleaned_rotated_points_in_main_excavation_area)

plot_ages_and_artefacts(osl_ages, c14_ages, stone_artefacts_only)

depth_age_tbl <- 
  plot_ages_and_artefacts(osl_ages, c14_ages, stone_artefacts_only)

saveRDS(depth_age_tbl,
        "data/depth_age_tbl.rds")

Extended Data Table 1

Raw materials are predominantly quartzite, silcrete, exotic fine-grained quartzite, mudstone and dolerite (Extended Data Table 1)

Data delivery to GAC

cleaned_rotated_points_in_main_excavation_area %>% 
  select(-code) %>% 
readr::write_csv("data/Madjedbebe_2012_and_2015_cleaned_rotated_points_in_main_excavation_area.csv")

# data dictionary
Madjedbebe_2012_and_2015_data_dictionary <- 
  tribble(~"variable",                    ~"explanation",
           "Description"               ,   "Code entered into total station to label the point",  
           "Easting"                   ,  "Local grid easting relative to datum, recorded by total station",
           "Northing"                  ,   "Local grid northing, relative to datum, recorded by total station",    
           "Elevation"                 ,   "Depth below datum, recorded by total station",   
           "year"                      ,   "Excavation year",
           "X"                         ,    "Computed coordinate",   
           "Y"                         ,    "Computed coordinate",
           "Ynew"                       ,   "Computed coordinate, used for plotting data",      
           "Xnew"                       ,   "Computed coordinate",
           "Xnew_flipped"              ,    "Computed coordinate, used for plotting data",
           "type"                       ,    "Type of total station point: PF, EL, OSL, etc., extracted from Description",      
           "square"                    ,     "Excavation square containing the point, extracted from Description",   
             "spit"                     ,   "Excavation spit containing the point, extracted from Description", 
           "cnr"                        ,    "Excavation square corner for end level points, extracted from Description", 
           "find"                      ,     "find type, extracted from Description",
           "findn"                      ,    "find number, extracted from Description",
           "depth_below_ground_surface",     "Computed coordinate, used for plotting data in section view"
  )

readr::write_csv(Madjedbebe_2012_and_2015_data_dictionary,
                 "data/Madjedbebe_2012_and_2015_cleaned_rotated_points_in_main_excavation_area_dictionary.csv")