Plotting the geological timescale in R

Posted on by

225px-Calycoceras_(Newboldiceras)_asiaticum_JIMBO_spinosum_KOSSMAT_microconchI’ve been figuring out nice ways to use R to plot the entire geological timescale. This being R, there is a package (called “geoscale“) which does this already, but I want one where the age is converted linearly to greyscale, and which can fit entirely onto one page, with the short periods of the Holocene and Pleistocene crammed in somehow. I also want to be able to update it easily when new values are published.

Below is my effort, with code: it also puts in the Tertiary and the American names for the subperiods of the Carboniferous. I’ve only gone to the base of the Archaean, and the names need some tidying to fit in some places (e.g. the Lower Triassic). But it’s nicer in some respects than most I’ve seen online (click for a pdf).
Geoscale
Here’s the code. It uses the font family STIXGeneral as that’s the only one on my system that has the unicode symbol for the Cambrian (Ꞓ), so you might need to change the “fam” variable to something elseaNB: I’m unsure whether to use the US symbols for Triassic and Palaeogene, or simply the letters “T” and “E” respectively. There seems to me to be no reason to reserve “T” for the tertiary any more..

geoscale.r (code in the public domain)

fam = "STIXGeneral-Regular"
Phanerozoic=list(
  CZ=list(Q=c(Holocene=0.0117,Pleistocene=2.58),
          Tertiary=list(N=c(Pliocene=5.333,Miocene=23.03), 
          PG=c(Oligocene=33.9,Eocene=56,Paleocene=66))),
  MZ=list(K=c(Upper=100.5,Lower=145),
          J=c(Upper=163.5,Middle=174.1,Lower=201.3),
          TR=c(Upper=237,Middle=247.2,Lower=252.17)),
          PZ=list(P=c(Lopingian=259.8,Guadalupian=272.3,Cisuralian=298.9),      
         C=list(Penn.=c(Upper=307,Middle=315.2,Lower=323.2),Miss.=c(Upper=330.9,Middle=346.7,Lower=358.9)), 
          D=c(Upper=382.7,Middle=393.3,Lower=419.2),   
          S=c(Pridoli=423,Ludlow=427.4,Wenlock=433.4,Llandovery=443.8),
          O=c(Upper=458.4,Middle=470,Lower=485.4), 
          Cb=c(Furongian=497,'Series 3'=509,'Series 2'=521,Terreneuvian=541))
)
 
Precambrian = list(PROTEROZOIC=c(NP=1000,MP=1600,PP=2500),
                   ARCHEAN=c(NEOARCHEAN=2800,MESOARCHEAN=3200,PALAEOARCHEAN=3600,EOARCHEAN=4000))
 
abbrev=c(NP='NEOPROTEROZOIC',MP='MESOPROTEROZOIC',PP='PALAEOPROTEROZOIC',
         CZ='CENOZOIC',MZ='MESOZOIC',PZ='PALEOZOIC',
         Q='Quaternary',N='Neogene',PG='Palaeogene',T='Tertiary',
         K='Cretaceous',J='Jurassic',TR='Triassic',
         P='Permian',C='Carboniferous',D='Devonian',S='Silurian',O='Ordovician',Cb='Cambrian')
 
pos=c(-2,7,13,31.5,48,55)
diagonal_frac=0.2 #for diagonal boxes in the holocene etc
textheight = 10
 
frac = function(mn,mx) mean(c(mn,mx))/4500
 
greyscale <- function(mn,mx) do.call(rgb,as.list(c(rep(1-frac(mn,mx),3),maxColorValue = 1)))
 
rmatch <- function(x, name) {
  pos <- match(name, names(x))
  if (!is.na(pos)) return(x[[pos]])
  for (el in x) if (length(names(el))) {
      out <- rmatch(el, name)
      if (!is.null(out)) return(out)
  }
}
 
draw_levels <- function(name, children, x, mindate, depth) {
    maxdate=mindate
    truemindate=ifelse(mindate<0,0,mindate)
    y.add = 0.5
    x.add = 0.4
    #exceptions here    
    switch(name,
           Tertiary =   { #Plot Tertiary (as in USGS) although not in the official stratigraphy
                            x <- c(x[1],weighted.mean(x[1:2],c(4,1)),x[-1])
                            depth=depth-1
                         },
           C =          { #Plot American names for Carboniferous divisions
                            x <- c(x[1],weighted.mean(x[1:2],c(1,10)),weighted.mean(x[2:3],c(5,2)),x[-1:-2])
                        },
           S =          { #Omit Silurian names for lack of space
                            names(children) <- rep("", length(names(children)))
                        }
           )
 
 
    if (length(unlist(children)) > 1) {
        for (n in 1:length(names(children))) {
            maxdate=draw_levels(names(children)[n], children[[n]], x[-1], maxdate, depth+1)
        }
    } else {
        maxdate = children[[1]]
    }
 
    if (!is.na(name)) {
    pos2 = sum(c(diagonal_frac, 1-diagonal_frac)*x[1:2])
    switch(name,
           Holocene =   {
                            polygon(c(x[1],pos2,x[2],x[2],pos2,x[1]),c(mindate+textheight,mindate+textheight,maxdate,truemindate,mindate,mindate),col= greyscale(truemindate, maxdate))
                        },
           Pleistocene ={
                            mindate <- maxdate - textheight
                            polygon(c(x[1],x[2],x[2],pos2,x[1]),c(maxdate,maxdate,truemindate,mindate,mindate),col=greyscale(truemindate, maxdate))
                        },
           Pliocene    = {
                            polygon(c(x[1],pos2,x[2],x[2],x[1]),c(textheight+mindate,textheight+ mindate,maxdate,mindate,mindate),col=greyscale(truemindate, maxdate))
                            maxdate <- textheight + mindate                          
                        },
           Miocene    = {
                            truemindate = rmatch(Phanerozoic,'Pliocene')
                            polygon(c(x[1],x[2],x[2],pos2,x[1]),c(maxdate,maxdate,truemindate,mindate,mindate),col=greyscale(truemindate, maxdate))
                        },
           rect(x[1],mindate,x[2],maxdate,col=greyscale(truemindate, maxdate))
    )
    if(name %in% names(abbrev)) {
        if (name=="Cb")
            fullname=list("Cambrian","(\uA792)")
        else if(nchar(name)>1)
            fullname=list(abbrev[name], bquote("("*.(substring(name,1,1))*scriptscriptstyle(bold(.(substring(name,2))))*")"))
        else
            fullname=list(abbrev[name], bquote("("*.(name)*")"))
    } else {
            fullname=name
    }
    cl = ifelse(frac(truemindate, maxdate)>0.55,'white','black')
    sz = ifelse(name=="Precambrian",1.5,1)
    if (name == "PG" || name == "C") {
        text(x[1]+x.add, mindate+y.add, fullname[[1]], adj=c(0,1.2), family=fam, col=cl,cex=sz)
        text(x[1]+x.add, mindate+textheight+y.add, fullname[[2]], adj=c(0,1.2), family=fam, col=cl,cex=sz)
    } else {
        if (length(fullname)>1) fullname = bquote(.(fullname[[1]])~.(fullname[[2]]))
        if(depth > 1)
            text(x[1]+x.add, mindate+y.add, fullname, adj=c(0,1.2), family=fam, col=cl,cex=sz)
        else
            text(mean(x[1:2]), mean(c(mindate,maxdate))+y.add, fullname, srt=90, family=fam, col=cl,cex=sz)
    }
    if (depth==2) {
        rect(rev(x)[2],truemindate,rev(x)[1],maxdate)
        text(rev(x)[2]+x.add,maxdate + y.add,maxdate,adj=c(0,-0.3),cex=0.8)
    }
    }
    return(maxdate)
}
 
 
par(mar=c(0.5,1,2,1), xpd=NA)
layout(t(1:2))
plot(1,1,type="n", xlim=range(pos), ylim=rev(range(unlist(Phanerozoic))), bty="n", xlab="",ylab="", yaxs="i", axes=FALSE)
draw_levels("PHANEROZOIC", Phanerozoic, pos, rmatch(data,'Pleistocene')-textheight*2, 0)
plot(1,1,type="n", xlim=range(pos), ylim=rev(c(0,4500)), bty="n", xlab="",ylab="", yaxs="i", axes=FALSE)
draw_levels("Precambrian", Precambrian, pos[-4], max(unlist(Phanerozoic)), 0)
shortPhan <- list(PHANEROZOIC=sapply(Phanerozoic,function(x) max(unlist(x))))
names(shortPhan[[1]]) <- rep("", length(shortPhan[[1]]))
draw_levels(NA, shortPhan, pos[c(1,1,3,5,6)], 0, 0)
lines(c(-12,-2),c(0,0), lty=3, xpd=NA)
lines(c(-12,-2),c(4500,max(unlist(Phanerozoic))), lty=3, xpd=NA)

Notes   [ + ]

a. ↑ NB: I’m unsure whether to use the US symbols for Triassic and Palaeogene, or simply the letters “T” and “E” respectively. There seems to me to be no reason to reserve “T” for the tertiary any more.

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.