Simple plotting options for radiocarbon dates in c14_date_lists

Ridgeplots of density distributions

One way to plot radiocarbon ages is to show the probability density distribution of individual calibrated dates as ridgeplots. To produce a plot like that, we first of all need the age-probability information for each date. We can calculate that with the function c14bazAAR::calibrate().

Batalimo_calibrated <- Batalimo %>%
  calibrate(choices = "calprobdistr")

This adds a list column calprobdistr to the input c14_date_list. The list column contains a nested data.frame for each date with its probability distribution.

##  Radiocarbon date list
##  dates: 8
##  sites: 1
##  countries: 1
##  uncalBP: 2000 ― 200 
## 
## # A tibble: 8 x 16
##   sourcedb sourcedb_version labnr c14age c14std calprobdistr c13val site 
##   <chr>    <date>           <chr>  <int>  <int> <list>        <dbl> <chr>
## 1 adrac    2020-01-12       Bdy-…   1890    130 <df[,2] [1,…      0 Bata…
## 2 adrac    2020-01-12       Bdy-…   1730    120 <df[,2] [1,…      0 Bata…
## 3 adrac    2020-01-12       Bdy-…   1990    210 <df[,2] [1,…      0 Bata…
## 4 adrac    2020-01-12       Bdy-…    240     75 <df[,2] [54…      0 Bata…
## 5 adrac    2020-01-12       Bdy-…   1270    125 <df[,2] [1,…      0 Bata…
## 6 adrac    2020-01-12       Bdy-…   1798    101 <df[,2] [95…      0 Bata…
## 7 adrac    2020-01-12       Gif-…   1590     90 <df[,2] [73…      0 Bata…
## 8 adrac    2020-01-12       OxTL…   1570    220 <df[,2] [1,…      0 Bata…
## # … with 8 more variables: feature <chr>, period <chr>, culture <chr>,
## #   material <chr>, country <chr>, lat <dbl>, lon <dbl>, shortref <chr>

With tidyr::unnest() the list column can be dissolved (“unnested”) and integrated into the initial c14_date_list. Of course the latter looses its original structure and meaning with this step. Each row now represents the probability for one date and year.

Batalimo_cal_dens <- Batalimo_calibrated %>% tidyr::unnest()

## Warning: `cols` is now required.
## Please use `cols = c(calprobdistr)`

A table like that can be used for plotting a ridgeplot.

Batalimo_cal_dens %>%
  ggplot() +
  # a special geom for ridgeplots is provided by the ggridges package
  ggridges::geom_ridgeline(
    # the relevant variables that have to be mapped for this geom are 
    # x (the time -- here the calibrated age transformed to calBC), 
    # y (the individual lab number of the dates) and
    # height (the probability for each year and date) 
    aes(x = -calage + 1950, y = labnr, height = density),
    # ridgeplots lack a scientifically clear y axis for each 
    # distribution plot and we can adjust the scaling to our needs
    scale = 300
  ) +
  xlab("age calBC/calAD") +
  ylab("dates")

Calcurve plot

Another way to plot radiocarbon dates is to project them onto a calibration curve. The Bchron R package contains a data.frame with the intcal13 calibration curve data. We can load the intcal13 table directly from Bchron.

load(system.file('data/intcal13.rda', package = 'Bchron'))

For this kind of plot it is more convenient to work with the simplified calrange output of c14bazAAR::calibrate().

Batalimo_calibrated <- Batalimo %>%
  calibrate(choices = "calrange")

Like the calprobdistr option this also adds a list column to the input c14_date_list, but a much smaller one. For each date only the age ranges that make up the 2-sigma significance interval of the probability distribution are stored.

Batalimo_calibrated$calrange[1:3]

## [[1]]
##   dens from   to
## 1 95.2 1530 2148
## 
## [[2]]
##   dens from   to
## 1 93.1 1377 1899
## 2  2.0 1913 1919
## 
## [[3]]
##   dens from   to
## 1 94.3 1413 2439
## 2  1.0 2450 2455

The resulting table can also be unnested to make the list column content available in the main table.

Batalimo_cal_range <- Batalimo_calibrated %>% tidyr::unnest()

## Warning: `cols` is now required.
## Please use `cols = c(calrange)`

Now we can plot the calibration curve and – on top – error bars with the calrange sequences.

ggplot() +
  # line plot of the intcal curve
  geom_line(
    data = intcal13,
    # again we transform the age information from BP to BC
    mapping = aes(x = -V1 + 1950, y = -V2 + 1950)
  ) +
  # the errorbars are plotted on top of the curve
  geom_errorbarh(
    data = Batalimo_cal_range,
    mapping = aes(y = -c14age + 1950, xmin = -to + 1950, xmax = -from + 1950)
  ) +
  # we define the age range manually -- typically the calcurve
  # is arranged to go from the top left to the bottom right corner
  xlim(-1000, 2000) +
  ylim(2000, -1000) +
  xlab("age calBC/calAD") +
  ylab("uncalibrated age BC/AD")

Interactive map view

The resulting sf object can be plotted interactively with the mapview package.

Moga_spatial %>% mapview::mapview()

Static map plot

The sf object can also be used for a static plot – which is useful for publications. We download some simple country border base map vector data with the rnaturalearth R package and transform it to sf as well.

countries <- rnaturalearth::ne_countries() %>% sf::st_as_sf()

Now we can combine the base layer and our point data to create the prototype of a static map plot.

ggplot() +
  # geom_sf is a special geom to handle spatial data in the sf format
  geom_sf(data = countries) +
  # the explicit mapping of variables is not necessary here, as geom_sf 
  # automatically finds the *geom* column in the input table
  geom_sf_text(data = countries, mapping = aes(label = formal_en), size = 2) +
  geom_sf(data = Moga_spatial) +
  # with geom_sf comes coord_sf to manage the underlying coordinate grid
  coord_sf(xlim = c(10, 30), ylim = c(0, 15))

Please feel free to open an issue here if you have questions about plotting radiocarbon dates.