Figure of Cenozoic (last 66 million years) and future global temperatures expected with climate change

cc_paleotemp.R

library(ggplot2)
library(cowplot)
library(latex2exp)


#### PALEOTEMPERATURES ####

## Data from Hansen: http://www.columbia.edu/~jeh1/mailings/2012/20120508_ClimateSensitivity.pdf
## Available here: http://www.columbia.edu/~mhs119/Sensitivity+SL+CO2/Table.txt

hansen <- read.table("http://www.columbia.edu/~mhs119/Sensitivity+SL+CO2/Table.txt",
                     header = FALSE, dec = ".", nrows = 17604, skip = 9)
hansen <- hansen[, c(3,6)]
names(hansen) <- c("MyrBP", "Tabs")
hansen$logtime <- log10(hansen$MyrBP)


timebreaks <- c(0.001, 0.01, 0.1, 1, 10, 66)  # in MyrBP
timebreaks.log <- log10(timebreaks)
time.labels <- TeX(c("10^{-3}", "10^{-2}",
                  "10^{-1}", "1", "10", "66"))

temp <- ggplot(hansen, aes(x = logtime, y = Tabs)) +
  ylim(9, 30) +
  labs(x = "Millions of years BP", y = "Temperature (ºC)") +
  theme(axis.text.x = element_text(size = 10)) +
  geom_line(colour = "Dark Red") +
  scale_x_continuous(breaks = timebreaks.log,
                     labels = time.labels,
                     trans = "reverse")


epochs.start <- c(0.0117, 2.58, 5.333, 23.03, 33.9, 56, 66)  # from geoscale

temp.paleo <- temp +
  geom_vline(xintercept = log10(epochs.start), linetype = "dashed", size = 0.2) +
  annotate("text", label = c("P", "Eo", "Ol", "Mi", "Pli", "Ple", "Hol"),
            x = c(1.78, 1.63, 1.44, 1.07, 0.58, -0.7, -2.9),
           y = 30, size = 3)
temp.paleo


#### FUTURE #####

## Climate projections (CMIP5 averages from KNMI Climate Explorer)

rcp26 <- read.table("http://climexp.knmi.nl/CMIP5/Tglobal/global_tas_Amon_modmean_rcp26_000.dat",
                    skip = 5, col.names = c("year", paste("M", 1:12, sep = "")))
rcp85 <- read.table("http://climexp.knmi.nl/CMIP5/Tglobal/global_tas_Amon_modmean_rcp85_000.dat",
                    skip = 5, col.names = c("year", paste("M", 1:12, sep = "")))

rcp26$tm <- apply(rcp26[, -1], 1, function(x) {mean(x) - 273.15}) # convert K to Celsius
rcp85$tm <- apply(rcp85[, -1], 1, function(x) {mean(x) - 273.15}) # convert K to Celsius

rcps <- data.frame(year = rcp26$year, rcp26 = rcp26$tm, rcp85 = rcp85$tm)
rcps <- rcps[rcps$year > 1999, ]


impdates <- c(2000, 2100)

temp.futu <- ggplot(rcps, aes(x = year, y = rcp26)) +
  labs(x = "Year") +
  theme(axis.text.x = element_text(size = 10)) +
  geom_line(colour = "Royal Blue 2", size = 2) +
  scale_x_continuous(breaks = impdates, limits = c(1970, 2130)) +
  scale_y_continuous(position = "right", limits = c(9, 30)) +
  theme(axis.title.y = element_blank()) +
  geom_line(aes(y = rcp85), colour = "Light Salmon 1", size = 2) +
  annotate("text", label = "RCP2.6", x = 2080, y = 14.4, size = 4, colour = "Royal Blue 2") +
  annotate("text", label = "RCP8.5", x = 2080, y = 19, size = 4, colour = "Light Salmon 1")




#### COMBINED FIGURE ####

fig <- plot_grid(temp.paleo, temp.futu, nrow = 1, rel_widths = c(1, 0.2),
                 align = "h")
fig