gistfile1.txt

# This script was prepared for ARCHY 486 AU23. It will draw a plot of particle size distribution 
# of multiple sediment samples on a log scale with major size categories indicated for easy 
# comparison. The data should be formatted as in this sheet:
# https://docs.google.com/spreadsheets/d/11RfkGzjpeAT1MAt1w-L5HFgdqzX_kNcXeX3EymvIBpo/edit#gid=390081773
# with column names exactly as found there. 

# get raw data on mass of sediment in the sieves
sieve_measurements <- 
  lab_data_my_group %>% 
  select(sieve_starting_mass_g, 
         matches("m_mass_g")) %>% 
  mutate(across(everything(), 
                ~ per_starting_mass(., sieve_starting_mass_g))) %>% 
  select(-sieve_starting_mass_g) 

# transpose this table to one row per sieve, samples in columns
sieve_measurements_t <-  
  data.frame(t(sieve_measurements) )

# put the sample IDs as column names
names(sieve_measurements_t) <- 
  paste0(lab_data_my_group$`Bag label`,
         lab_data_my_group$Context)

# get the sieve size units so we can make them all the same unit
seive_size_unit <- 
sieve_measurements %>% 
  names() %>% 
  str_extract(., "mm|um")

# get the sieve size values
seive_size_value <- 
  sieve_measurements %>% 
  names() %>% 
  parse_number()

# get the sieve size values all in the same unit, mm
seive_size <- 
tibble(seive_size_unit = seive_size_unit,
      seive_size_value = seive_size_value ) %>% 
  mutate(seive_size_value_mm = case_when(
    seive_size_unit == "um" ~ seive_size_value / 1000,
    TRUE ~ seive_size_value
  )) %>% 
  select(-seive_size_unit, -seive_size_value) %>% 
  pull(seive_size_value_mm)

# combine sieve sizes with raw data from samples
sieve_measurements_t_sizes_long <- 
bind_cols(sieve_measurements_t) %>% 
  tibble( seive_size = seive_size) %>% 
  pivot_longer(-seive_size,
               names_to = "Sample") %>% 
  mutate(value = ifelse(value == 0, 0.001, value))

# make a table of standard particle size classes so we can put these on the
# plot and make it easier to interpret,  from 
# https://en.wikipedia.org/wiki/Soil_texture, https://en.wikipedia.org/wiki/Grain_size
soil_texture_tbl <- 
tribble(~texture,   ~lower_size_usda, ~upper_size_usda, ~lower_size_wrb, ~upper_size_wrb,
          "Clay",                 0,         0.002,  0,       0.002,
          "Silt",                 0.002,     0.05,   0.002,   0.063,
          "Very fine sand",       0.05,      0.10,   0.063,   0.125,
          "Fine sand",            0.10,      0.25,   0.125,   0.20,
          "Medium sand",          0.25,      0.50,   0.20 ,   0.63,
          "Coarse sand",          0.50,      1.00,   0.63 ,   1.25,
          "Very coarse sand",     1.00,      2.00,   1.25 ,   2.00,
          "Very fine gravel",     2,         4,      2,       4,
          "Fine gravel",          4,         8,      4,       8) %>% 
  mutate(midpoint = lower_size_wrb + (upper_size_wrb - lower_size_wrb) / 2) 

# draw the plot
ggplot(sieve_measurements_t_sizes_long) +
  aes(seive_size,
      value,
      colour = Sample) +
  geom_line() +
  # draw lines for the texture classes, omit clay
  geom_vline(xintercept = soil_texture_tbl$upper_size_wrb[-1],
             colour = "grey80") +
    # draw text labels for the texture classes, omit clay
  annotate("text",
           label = soil_texture_tbl$texture[-1],
           x =soil_texture_tbl$midpoint[-1],
           y = 80,
           size = 3,
           angle = 90) +
  scale_x_continuous(breaks = soil_texture_tbl$upper_size_wrb[-1],
                     name = "Particle size, D (mm)") +
  coord_trans(x = "log10") +
  scale_y_continuous(limits = c(0, 100),
                     name = "Percent by mass") +
  theme_minimal()

Here's an example of the figure that this code produces, using data from https://docs.google.com/spreadsheets/d/11RfkGzjpeAT1MAt1w-L5HFgdqzX_kNcXeX3EymvIBpo/edit#gid=390081773