Palaeogeography Quick Guide

Palaeogeographic rasters

Rasterised representations of deep time geography, and extending to other properties like vegetation type or climate zone, are increasingly available to palaeobiologists. Some are archived in paper electronic supplements (e.g., the source of the palaeo-Galapagos data included in rTARDIS), but easy access to some of these datasets is provided by a few different R packages. First among these is chronosphere. This vignette gives a short overview of how to quickly access data suitable for building rTARDIS models for palaeogeographic analysis using this package. First we will load the chronosphere package (prior installation required) and examine the raster datasets it offers:

library(chronosphere)
datasets()[which(datasets()$defaultClass == "RasterArray"),4]

#> [1] "PaleoMAP Paleodigital Elevation Models"                                                                                     
#> [2] "Global Mean Surface Temperature"                                                                                            
#> [3] "PALEOMAP PaleoAtlas for GPlates"                                                                                            
#> [4] "PALEOMAP Paleorainfall (Total Annual Precipitation)"                                                                        
#> [5] "S6. Interpolated  temperature  reconstructions  based  on  the  HadCM3L  climate  models  and  the PaleoMAP Paleocoastlines"
#> [6] "Monthly climate layers (Temperature, Precipitation, Surface Albedo, U and V Wind)"

The first of these options is the most relevant to our purposes - digital elevation models of topography and bathymetry at 5-million-year intervals through the Phanerozoic, created by the PALEOMAP project. We can download this dataset as follows:

dem <- fetch("paleomap", "dem")

The data is returned as a RasterArray, a class from the via R package designed for storing temporal arrays of vector and raster data. A bit like a geoglist from rTARDIS, this simply wraps a SpatRaster and a vector of time intervals which we can extract directly. The via package has functions which can conflict with those from terra, so we will immediately unload it to avoid these issues:

tms <- dem@index
dem <- dem@stack
unloadNamespace("via")
tms

#>   0   5  10  15  20  25  30  35  40  45  50  55  60  65  70  75  80  85  90  95 
#>   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
#> 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 
#>  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
#> 200 205 210 215 220 225 230 235 240 245 250 255 260 265 270 275 280 285 290 295 
#>  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 
#> 300 305 310 315 320 325 330 335 340 345 350 355 360 365 370 375 380 385 390 395 
#>  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80 
#> 400 405 410 415 420 425 430 435 440 445 450 455 460 465 470 475 480 485 490 495 
#>  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 
#> 500 505 510 515 520 525 530 535 540 600 630 690 750 
#> 101 102 103 104 105 106 107 108 109 110 111 112 113

If we examine tms, we can see it contains 113 elements, with each element name denoting the age of each raster layer. These ages are also present in the raster layer names themselves (000, 005, 010 ect.). We can select any contiguous subset of these layers to use with rTARDIS. Indeed, the cretaceous data set in the package is just the time layers for 115 and 110 million years ago, taken from chronosphere.

For this vignette, we will extract 4 time slices for the Devonian from the chronosphere data set, mask those rasters to the land surface, and create a geoglist with hexagonal re-gridding to improve representation of geographic space at the poles:

# load the dataset: slices 49 to 46 correspond to 385, 380, 375, and 370 million years ago
dvn <- dem[[78:75]]

# create the land-sea mask
dvn_l <- classify(dvn, matrix(c(-Inf, 0, NA, 0, Inf, 1), ncol = 3, byrow = T), right = F)

# convert to geoglist
hexes <- rast_to_geoglist(dvn, dvn_l, times = c(387.5, 382.5, 377.5, 372.5, 367.5), as.hex = T, hex = 3)

# plot
plot(hexes, 2)

There are many disconnected islands in this landscape. If we want to be able to move between them, we will need to add links between those islands. We can do this very simply by adding single connections between islands until there is some possible route between all the islands:

# link islands
hexes <- link_islands(hexes, klink = 1)

# plot
plot(hexes, 2)

#> Resolving layers [1/4]Resolving layers [2/4]Resolving layers [3/4]Resolving layers [4/4]

This is a simple, conservative linking scheme and others are available. You can take a look at the vignette Building Geographic Spaces for more general guidance.

Plate tectonic drift

Over geological timescales, locations on Earth’s surface have shifted substantially due to plate tectonic drift. As such, cells in a global palaeogeographic landscape may no longer be spatially homologous between layers. rTARDIS can account for this effect by including plate rotation into landscape construction.

Quite a specific data structure is requiree. It must be a list with nlayers - 1 elements relative to the accompanying geoglist. Each list element is a two-column matrix of cell IDs describing pairs of cells to link between successive time layers. These IDs should be spatial, i.e., row-major raster numeric index or hexagonal grid numeric index, rather than full spatiotemporal indexes. This is to permit more intuitive construction of these matrices if desired, with temporal information added automatically during graph construction.

In practice, these lists may be quite complicated to construct by hand, particularly for plate tectonic rotation models. rTARDIS provides a function to automate this procedure called get_rotations(). This takes a geoglist, its accompanying time vector and a plate rotation model name as its principle arguments, then uses the palaeoverse R package function palaeorotate() to calculate rotations between cells in each successive pair of geoglist layers. The user is directed to the documentation for palaeorotate() for more detailed information on what plate rotation models are available. Crucially, however, this should match the palaeogeographic data used to construct the geoglist. This function may take a little while to run for high resolution grids and relies on an API call, which may further affect run times.

We will implement this workflow below. The Devonian dataset contains 4 time layers, based on the PALEOMAP plate rotation model, at 5 million-year intervals. Firstly, we will resample to hexagons, both to better represent global geography and to speed up downstream function runs. We will derive rotations for the geoglist, using its vector of interval times. Note that the function also requires that times are ordered forwards in time:

# get rotations for the PALEOMAP model
rots <- get_rotations(hexes, hexes$tdat, "PALEOMAP")

# interrogate the output
length(rots)
head(rots[[1]])

#> [1] 3
#>       [,1]  [,2]
#> [1,] 19773 19776
#> [2,] 19648 19887
#> [3,] 19770 19751
#> [4,] 19772 19763
#> [5,] 19636 19770
#> [6,] 19649 19772

As expected, rots contains three list elements, for the four layers in hexes. We can also see the effect of the plate movement, with the third entry in the first element changing between layers. Now we can build our tardis object.

# build tardis
htd <- build_tardis(hexes, rotations = rots)
#> Linking geog [1/4]Linking geog [2/4]Linking geog [3/4]Linking geog [4/4]
#> Building graph

We now have a temporally varying palaeogeography model which we can analyse using other functions from rTARDIS. These are covered by other vignettes.