# Simple plotting options for radiocarbon dates in c14_date_lists

#### 2020-01-12

This vignette shows some basic workflows to plot radiocarbon dates in c14_date_lists. This is only a short compilation to get you started.

So let’s begin by loading the main packages for this endeavour: c14bazAAR and ggplot2. And of course magrittr to enable the pipe (%>%) operator. We will use functions from other packages as well, but address them individually with the :: operator.

library(c14bazAAR)
library(ggplot2)
library(magrittr)

The basis for this example code is the adrac data collection by Dirk Seidensticker. So let’s download the data with the relevant c14bazAAR getter function:

adrac <- get_c14data("adrac")
## Trying to download all dates from the requested databases...
##
|
|                                                  |   0%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++|  99%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%

## Temporal plotting of radiocarbon ages

Radiocarbon dating is a method to determine the absolute age of samples. Therefore one of the main aims for plotting naturally is to display temporal information. Let’s select the dates of one individual site – Batalimo in Central Africa – as a subset to reproduce two of the most common types of radiocarbon date plots.

Batalimo <- adrac %>%
dplyr::filter(site == "Batalimo")

If age modelling and date plotting on a local or regional scale is the major aim of your analysis, you might want to take a look at the oxcAAR package. It serves as an R interface to OxCal and provides powerful default plotting methods

### 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
) +
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) +
ylab("uncalibrated age BC/AD")

## Spatial mapping of radiocarbon dates

Most radiocarbon dates that can be accessed with c14bazAAR have coordinate information for the respective sites where the samples were taken. Spatial maps therefore are an important form of data visualization as well.

c14_date_lists can directly be transformed to objects of class sf. sf objects were introduced by the R package sf which provides a tidy interface to work with spatial data in R.

adrac_sf <- adrac %>% as.sf()

This tabular data structure contains the spatial point information for each date in a column geom, but also the initial columns of the input dataset: data.*

## Simple feature collection with 1756 features and 3 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 8.7167 ymin: -15.27044 xmax: 30.68333 ymax: 12.2
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##    data.labnr data.c14age data.c14std                      geom
## 1    AA-78447        2362          39 POINT (16.92445 3.988639)
## 2    AA-78448        2171          37 POINT (16.92445 3.988639)
## 3    AA-78449         834          35 POINT (16.92445 3.988639)
## 4    AA-94529         215          34       POINT (17.47 3.835)
## 5    AA-94530         168          35       POINT (17.47 3.835)
## 6    AA-94531         207          35       POINT (17.47 3.835)
## 7    AA-94532         148          34       POINT (17.47 3.835)
## 8    AA-94533         231          34       POINT (17.47 3.835)
## 9    AA-94534         187          34       POINT (17.47 3.835)
## 10   AA-94535         412          34    POINT (17.5151 3.7628)

It can be manipulated with the powerful dplyr functions. We filter out all dates from one particular publication (Moga 2008), group the dates by site and apply the summarise command to keep only one value per group. As we do not define an operation to fold the other variables in the input table, they are removed. Only the geometry column remains.

Moga_spatial <- adrac_sf %>%
dplyr::filter(grepl("Moga 2008", data.shortref)) %>%
dplyr::group_by(data.site) %>%
dplyr::summarise()

### 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))