Visualising the LR04 Benthic Stack

1. Introduction

The aim of this document is to refresh my memory on using R and creating plots. Although I have done this in R scripts, I also wanted to practice creating markdown files as well so I have written it up into here too. I’ve left all of the code on show (echo = TRUE) so that I can see how I have done everything when I look back over it.

The data I am using is the LR04 benthic stack. I am particularly interested in Marine Isotope Stages (MIS) 7 to 9 as this is the timeframe I am investigating in my PhD research into Mediterranean palaeoclimate.

2. Set up

Libraries required:

• ggplot2

Data:

• lr04.txt

library(ggplot2)
lr04 <- read.delim("./data/lr04.txt")

3. Check and view data

str(lr04)

## 'data.frame':    2115 obs. of  3 variables:
##  $ time: num  0 1 2 3 4 5 6 7 8 9 ...
##  $ d18o: num  3.23 3.23 3.18 3.29 3.3 3.26 3.33 3.37 3.42 3.38 ...
##  $ se  : num  0.03 0.04 0.03 0.03 0.03 0.03 0.04 0.04 0.03 0.04 ...

head(lr04)

##   time d18o   se
## 1    0 3.23 0.03
## 2    1 3.23 0.04
## 3    2 3.18 0.03
## 4    3 3.29 0.03
## 5    4 3.30 0.03
## 6    5 3.26 0.03

The data have been loaded correctly. The variables are:

• time – time (ka, thousands of years ago)
• d18O – δ¹⁸O (‰, per mil), the variable I am wanting to plot
• se – standard error (‰, per mil)

4. Plots

4.1 All data (ka)

First of all I just want to plot up the data as it is so I can have a look and see what I want to do with it.

ggplot(data = lr04, aes(x = time, y = d18o)) +        # set up plot
  geom_line(colour = "#00AFBB", size = 0.5) +         # plot line
  labs(title = "LR04 Benthic Stack",                    # plot labels 
       x = "Time (ka)", 
       y = expression(δ^{18}*O~"(‰)"),
       caption = "Data from Lisiecki & Raymo (2005)") + 
  scale_y_reverse() +                                 # reverse scale on y-axis
  theme_classic(base_size = 15)

Some notes on what I’ve done:

Labels

It took me quite a bit of time to figure out how to plot δ¹⁸O as the y-axis label and get it to be displayed correctly. I’ve done it now so I can just copy and reuse the relevant bit of code whenever I need it.

Axes

Lower δ¹⁸O (‰) represents warmer interglacial periods so I reversed the y-axis as it is more intuitive to see warmer temperatures towards the top. I had originally reversed the x-axis using scale_x_reverse() as this makes more sense (going from older on the left to younger on the right) but this doesn’t seem to be the done thing for this type of data. All plots I have seen show it going from younger to older, left to right. Time is plotted in thousands of years but I think it would best easier to understand if it was in millions of years.

Style

I choose the classic theme as I think it is the nicest and the cleanest. There are several other built in themes for ggplot and it is also possible to tweak the themes) too.

4.2 All data (Ma)

It would be easier to understand the plot if the x-axis time scale was in millions rather than thousands of years. This following line of code will do that.

lr04$time <- lr04$time / 1000

Now I can plot it up as I did before.

ggplot(data = lr04, aes(x = time, y = d18o)) +
  geom_line(colour = "#00AFBB", size = 0.5) +
  labs(title = "LR04 Benthic Stack",
       x = "Time (Ma)", 
       y = expression(δ^{18}*O~"(‰)"),                 
       caption = "Data from Lisiecki & Raymo (2005)") + 
  scale_y_reverse() +
  scale_x_continuous(breaks = seq(0, 6, by = 1)) +     # change tick frequency
  theme_classic(base_size = 15) 

As well as converting the time scale to millions of years, I have changed the tick frequency so that there is one everyone million years using scale_x_continuous(breaks = seq(0, 6, by = 1)). They were places similarly to the ka plot with a tick at 2 Ma and 4 Ma before I made this change.
 

It would be nice to plot the line colour as a gradient so cooler conditions are represented by a cold colour and warmer conditions represented by a warm colour.

ggplot(data = lr04, aes(x = time, y = d18o, colour = d18o)) + 
  geom_line(size = 0.5) + 
  scale_colour_gradient(low = "#DB5824", high = "#1A5A95") +
  labs(title = "LR04 Benthic Stack",                   
       x = "Time (Ma)", 
       y = expression(δ^{18}*O~"(‰)"),                
       caption = "Data from Lisiecki & Raymo (2005)") + 
  scale_y_reverse() +
  scale_x_continuous(breaks = seq(0, 6, by = 1)) +
  theme_classic(base_size = 15) +
  theme(legend.position = "none")

I really like this! It clearly reflects how the planet has cooled over the last 5 Ma. It works best for the last million years where parts of the line go blue. It would be good to be able to set the change to blue slighty higher up the plot as some of the earlier glacial periods still appear red. But this is something I can investigate at a later time if I need to. It would be more useful to zoom in on certain parts of the data set next.

4.3 Last 500 ka

Zooming in on part of the data just requires subsetting it to the section of it you are after. Here I have selected the last half a million years. I have also annotated the interglacial periods with their MIS numbers.

# subset last 500 ka
lr04_500 <- subset(lr04, time < 0.5)

# convert time from back from ma to ka
lr04_500$time <- lr04_500$time * 1000

ggplot(data = lr04_500, 
       aes(x = time, y = d18o, colour = d18o)) + 
  geom_line(size = 0.5) + 
  scale_colour_gradient(low = "#DB5824", high = "#1A5A95") +
  labs(title = "LR04 Benthic Stack (Last 500 ka)",
       x = "Age (ka)", 
       y = expression(δ^{18}*O~"(‰)"),                
       caption = "Data from Lisiecki & Raymo (2005)") + 
  scale_y_reverse() +                                 
  annotate("text", 
           x = c(5, 125, 205, 237, 327, 410), 
           y = c(3, 3, 3.3,3.3, 3, 3), 
           label = c("1", "5", "7 a-c", "7e", "9", "11")) +
  theme_classic(base_size = 15) +
  theme(legend.position = "none")

4.4 MIS 7–9

In this example I have focused in on the section of the record that corresponds to the time frame of my PhD research. I’ve left the sections either side of MIS 7–9 in for context. I’ve also labeled the glacial stages on this as well (as best I can quickly). I’ve also chaned the line thickness so it is a bit more clear (geom_line(size = 0.8)).

# subset MIS 7–9 
lr04_7to9 <- subset(lr04_500, time > 130 & time < 374) 

ggplot(data = lr04_7to9, 
       aes(x = time, y = d18o, colour = d18o)) + 
  geom_line(size = 0.8) + 
  scale_colour_gradient(low = "#DB5824", high = "#1A5A95") +
  labs(title = "LR04 Benthic Stack (MIS 7–9)",
       x = "Age (ka)", 
       y = expression(δ^{18}*O~"(‰)"),                
       caption = "Data from Lisiecki & Raymo (2005)") + 
  annotate("text", 
           x = c(160, 200, 225, 237, 260, 283, 297, 310, 320, 327, 360), 
           y = c(4.9, 3.4, 4.5, 3.3, 4.7, 3.75, 4.4, 3.6, 3.95, 3.1, 4.9), 
           label = c("6", "7 a-c", "7d", "7e", "8", "9a", "9b", "9c?", "9d?", 
                     "9e", "10")) +
  scale_y_reverse() +    
  scale_x_continuous(breaks=seq(120, 380, by = 20)) +
  theme_classic(base_size = 15) +
  theme(legend.position = "none")

4.5 Vertical plots

These data are often plotted vertically. All it takes is this small section of code to flip the axes around: scale_y_reverse(position = "right") + scale_x_reverse() + coord_flip()

All data

ggplot(data = lr04, aes(x = time, y = d18o)) +
  geom_line(colour = "#00AFBB", size = 0.5) +
  labs(title = "LR04 Benthic Stack",
       x = "Time (Ma)", 
       y = expression(δ^{18}*O~"(‰)"),                 
       caption = "Data from Lisiecki & Raymo (2005)") + 
  scale_y_reverse(position = "right") +  
  scale_x_reverse(breaks = seq(0, 6, by = 1)) +
  coord_flip() +  
  theme_classic(base_size = 15)

Last 500 ka

ggplot(data = lr04_500, 
       aes(x = time, y = d18o, colour = d18o)) + 
  geom_line(size = 0.5)+ 
  scale_colour_gradient(low = "#DB5824", high = "#1A5A95") +
  labs(title = "LR04 Benthic Stack (Last 500 ka)",
       x = "Age (ka)", 
       y = expression(δ^{18}*O~"(‰)"),                
       caption = "Data from Lisiecki & Raymo (2005)") +
  annotate("text", 
           x = c(5, 125, 205, 237, 327, 410), 
           y = c(3, 3, 3.3,3.3, 3, 3), 
           label = c("1", "5", "7 a-c", "7e", "9", "11")) +
  scale_y_reverse(position = "right") +  
  scale_x_reverse() +
  coord_flip() +  
  theme_classic(base_size = 15) +
  theme(legend.position = "none")

MIS 7-9

ggplot(data = lr04_7to9, aes(x = time, y = d18o, colour = d18o)) + 
  geom_line(size = 0.8) + 
  scale_colour_gradient(low = "#DB5824", high = "#1A5A95") +
  labs(title = "LR04 Benthic Stack (MIS 7–9)",
       x = "Age (ka)", 
       y = expression(δ^{18}*O~"(‰)"),                
       caption = "Data from Lisiecki & Raymo (2005)") + 
  annotate("text", 
           x = c(160, 200, 225, 237, 260, 283, 297, 310, 320, 327, 360), 
           y = c(4.9, 3.4, 4.5, 3.3, 4.7, 3.75, 4.4, 3.6, 3.95, 3.1, 4.9), 
           label = c("6", "7 a-c", "7d", "7e", "8", "9a", "9b", "9c?", "9d?", 
                     "9e", "10")) +
  scale_y_reverse(position = "right") +  
  scale_x_reverse(breaks = seq(120, 380, by = 20)) +
  coord_flip() +                             
  theme_classic(base_size = 15) +
  theme(legend.position = "none")

References

Lisiecki, L. E. & Raymo, M. E. (2005) A Pliocene-Pleistocene stack of 57 globally distributed benthic d18O records. Paleoceanography, 20, PA1003. Available online: http://lorraine-lisiecki.com/stack.html [Downloaded 2019-06-28].