Simple plotting options for radiocarbon dates in c14_date_lists
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
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")
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.
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
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: 7 ## sites: 1 ## countries: 1 ## uncalBP: 2000 ― 200 ## ## # A data frame: 7 × 17 ## sourcedb sourcedb_ve…¹ method labnr c14age c14std calpr…² c13val site feature ## <chr> <date> <chr> <chr> <int> <int> <list> <dbl> <chr> <chr> ## 1 adrac 2023-03-05 NA Bdy-… 1890 130 <df> NA Bata… Z III ## 2 adrac 2023-03-05 NA Bdy-… 1730 120 <df> NA Bata… Z IV ## 3 adrac 2023-03-05 NA Bdy-… 1990 210 <df> NA Bata… Z III ## 4 adrac 2023-03-05 NA Bdy-… 240 75 <df> NA Bata… Z VI ## 5 adrac 2023-03-05 NA Bdy-… 1270 125 <df> NA Bata… Z VI ## 6 adrac 2023-03-05 NA Bdy-… 1798 101 <df> NA Bata… Z VI ## 7 adrac 2023-03-05 NA Gif-… 1590 90 <df> NA Bata… ZII ## # … with 7 more variables: period <chr>, material <chr>, country <chr>, ## # lat <dbl>, lon <dbl>, shortref <chr>, comment <chr>, and abbreviated ## # variable names ¹sourcedb_version, ²calprobdistr
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.
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")
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
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.
## [] ## dens from to ## 1 91.9 1521 2148 ## 2 3.3 2275 2287 ## ## [] ## dens from to ## 1 93.7 1363 1889 ## 2 1.5 1913 1917 ## ## [] ## dens from to ## 1 94.1 1407 2430 ## 2 1.2 2453 2460
The resulting table can also be unnested to make the list column content available in the main table.
Now we can plot the calibration curve and – on top – error bars with the
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")
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 objects were introduced by the R package sf which provides a tidy interface to work with spatial data in R.
## Warning in remove_dates_without_coordinates(., quiet): Dates without ## coordinates were removed.
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 1913 features and 3 fields ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: 8.717 ymin: -15.27 xmax: 30.683 ymax: 12.2 ## Geodetic CRS: WGS 84 ## First 10 features: ## data.labnr data.c14age data.c14std geom ## 1 AA-33225 10015 55 POINT (28.6 2.17) ## 2 AA-33226 6025 70 POINT (28.6 2.17) ## 3 AA-78447 2362 39 POINT (16.924 3.989) ## 4 AA-78448 2171 37 POINT (16.924 3.989) ## 5 AA-78449 834 35 POINT (16.924 3.989) ## 6 AA-94529 215 34 POINT (17.47 3.835) ## 7 AA-94530 168 35 POINT (17.47 3.835) ## 8 AA-94531 207 35 POINT (17.47 3.835) ## 9 AA-94532 148 34 POINT (17.47 3.835) ## 10 AA-94533 231 34 POINT (17.47 3.835)
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(.groups = "drop")
Interactive map view
sf object can be plotted interactively with the mapview package.
# Moga_spatial %>% mapview::mapview()
Static map plot
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.