Offline Reconstruction Part II

This is an example of how to use the chronosphere package to reconstruct the paleogeography across the Late Paleozoic Ice Age (LPIA), 360 to 260 Ma. Occurrence data from the Paleobiology Database ( is downloaded and prepared using the divDyn package and then plotted onto these maps.

library(chronosphere) #load chronosphere package
library(divDyn)    #load divDyn package to prepare data

Part Two: Occurrence data

After downloading the occurrence data for all marine mammals during the LPIA from the Paleobiology Database using the fetch() function, the data is prepared following instructions from the divDyn package.

# occurrence data of all marine animals during LPIA 
# downloaded from Paleobiology Database
dat <- fetch("pbdb") 
data(tens) #load 10 million year time scale 
colnames(tens)[colnames(tens)=="X10"] <- "name" #assign column name "name" to table
data(keys) #load keys object which contains categories to group occurrences

#subset for marine animals only
marineNoPlant <- c("",

#get subsetted data
dat <- dat[dat$phylum %in% marineNoPlant,]

#subset to only include occurrences assigned to a genus or species
dat <- dat[dat$accepted_rank %in% c("genus", "species"), ] 

Now categorize entries for stratigraphic binning (details in divDyn handout):

binMin <- categorize(dat[ ,"early_interval"], keys$tenInt)
binMax <- categorize(dat[ ,"late_interval"], keys$tenInt)
binMin <- as.numeric(binMin)
binMax <- as.numeric(binMax)
dat$bin <- rep(NA, nrow(dat))
binCondition <- c(which(binMax==binMin),  which(binMax==-1))
dat$bin[binCondition] <- binMin[binCondition]

Returning to the instructions for the chronosphere package:

dems <- fetch(dat="paleomap", var="dem") #download and install paleomap dataset
#matches two vectors of ages 
demord <- matchtime(dems, tens$mid) 

I then subsetted to include only the required bins (View(tens) to check time bin vs. Ma).

interColl <- dat[dat$bin%in%c(17,18,19,21,23,24),
                 c("collection_no","lng", "lat", "bin", "paleolat" ,"paleolng")]

Remove duplicates, add age of map to data, and reconstruct the paleocoordinates:

interColl <- unique(interColl) #remove duplicates
interColl$mapage <- names(demord)[interColl$bin] #add column to interColl with map age

lpia <- c(360,340,320,300,280,260) #six time intervals during LPIA, 360-260 Ma

#reconstruct paleocoordinates for collections
reconstruct(interColl[1:2, c("lng", "lat")], age=lpia, model=pm)

#reconstruct every row to corresponding age in the ``age`` vector mapage
newCoords <- reconstruct(interColl[, c("lng", "lat")], 
                         age=interColl[, "mapage"], 
                         enumerate=FALSE, verbose=FALSE, model=pm) 

colnames(newCoords) <- c("plng", "plat") #rename columns
allColl <- cbind(interColl,newCoords) #combine vectors

Then I plot all of the occurrences onto the individual paleomaps. First I subset the collections by time slice:

coll360 <- allColl[allColl$bin==17,] # collections in each bin
coll340 <- allColl[allColl$bin==18,]
coll320 <- allColl[allColl$bin==19,]
coll300 <- allColl[allColl$bin==21,]
coll280 <- allColl[allColl$bin==23,]
coll260 <- allColl[allColl$bin==24,]

Use the paleomap dataset:

pr <- fetch(dat="paleomap", var="paleoatlas", res=0.5) #download and install paleomap dataset
## If you use the data in publications, please cite its
## reference, as well as that of the 'chronosphere' package.
## - Scotese, C. R. (2016) Tutorial: PALEOMAP PaleoAtlas for GPlates and the PaleoData Plotter Program. URL:
#accessing the different elements corresponding to each age separately
pa360 <- pr["360", ]  
pa340 <- pr["340", ]
pa320 <- pr["320", ]
pa300 <- pr["300", ]
pa280 <- pr["280", ]
pa260 <- pr["260", ]

And again, I would like to have all maps in one figure:

mapplot(pa360, rgb=TRUE) #Map at 360 Ma
points(coll360$plng, coll360$plat, pch=3, col="red")
mapplot(pa340, rgb=TRUE) #Map at 340 Ma
points(coll340$plng, coll340$plat, pch=3, col="red")
mapplot(pa320, rgb=TRUE) #Map at 320 Ma
points(coll320$plng, coll320$plat, pch=3, col="red")
mapplot(pa300, rgb=TRUE) #Map at 300 Ma
points(coll300$plng, coll300$plat, pch=3, col="red")
mapplot(pa280, rgb=TRUE) #Map at 280 Ma
points(coll280$plng, coll280$plat, pch=3, col="red")
mapplot(pa260, rgb=TRUE) #Map at 260 Ma
points(coll260$plng, coll260$plat, pch=3, col="red")