Using data from the PaleoMAP project

true

Goals and introduction

The purpose of this vignette is to demonstrate the newly implemented elements of the chronosphere. The primary feature of the package is to load published datasets of deep time paleogeography and implement functions that ease interacting with them. Paleogeographic data are available through the GPlates Web Service and through the repositories of the chronosphere package. For instructions on installing the package, please click here.

Plate reconstructions

The ‘static plate’ polygons (plates rotated with time) can be accessed directly from the GPlates Web Service using the reconstruct() function. For the demo we have chosen two ages that reflect a good match between PaleoDB occurrences and landmass distribution.

maps <- reconstruct("plates", age=c(300,295))

The function call create a list class object with as many elements as many ages are given. Have a look at the object:

maps
## $`300`
## class       : SpatialPolygonsDataFrame 
## features    : 87 
## extent      : -180, 180, -90, 69.6  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 1
## names       : FID 
## min values  :   0 
## max values  :  86 
## 
## $`295`
## class       : SpatialPolygonsDataFrame 
## features    : 86 
## extent      : -180, 180, -90, 71.61  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 1
## names       : FID 
## min values  :   0 
## max values  :  85

The different elements corresponding to each age can be accessed separately with:

# accessing the first polygon object (300 Ma)
ma300 <- maps[[1]]
ma300
## class       : SpatialPolygonsDataFrame 
## features    : 87 
## extent      : -180, 180, -90, 69.6  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 1
## names       : FID 
## min values  :   0 
## max values  :  86

The new ma300 object behaves like a regular SpatialPolygonsDataFrame object from the sp package. You can use its default methods to plot and manipulate the object

#plot
plot(ma300, col="gray")

The different elements can also be accessed by referring to their ages with character subscripts, i.e. quotation marks are required.

ma295 <- maps[["295"]]

plot(ma295, col="gray")

Coastlines

Coastlines can also be reconstructed via the web service. These can be accessed with.

coast <- reconstruct("coastlines", age=c(300,295))

This object has the same structure as the maps object before, and also has named elements, that you can use for subsetting:

ma295 <- maps[["295"]]
coast295 <- coast[["295"]]

Plotting works with the methods implemented in the sp package:

plot(ma295, col="gray", border=NA)
lines(coast295)

The reconstruct() function allows reconstructions with multiple plate tectonic models. The default setting is PALEOMAP that refers to the PaleoMAP project of C. Scotese (v3, 2018), but other models are available such as that of Seton et al. (2012).

Environmental variables

The primary goal of the package is to make raster-type environmental variables available directly from the R environment. These are hosted on a remote server, and can be downloaded whenever required. You can get the list of the currently available maps with the function dataindex().

ind <- dataindex()
head(ind)
##                                                  name        dat
## 1                    PALEOMAP Digital Elevation model   paleomap
## 2                    PALEOMAP Paleoatlas Paleorasters   paleomap
## 3                    PALEOMAP Paleoatlas Paleorasters   paleomap
## 4                               Paleobiology Database       pbdb
## 5 Continent Fragmentation index of Zaffos et al. 2017 timeseries
## 6             J.J. Sepkoski's genus level compendium    sepkoski
##          var res               ver rotation     licence proj4string
## 1        dem 1.0          20190719       NA CC-BY-NC-ND          NA
## 2 paleoatlas 0.1        20160216v3       NA CC-BY-NC-ND          NA
## 3 paleoatlas 0.5        20160216v3       NA CC-BY-NC-ND          NA
## 4             NA          20191122       NA       CC-BY          NA
## 5   fragment  NA zaffos_et_al_2017       NA  not formal          NA
## 6     genera  NA    2002_kiessling       NA  not formal          NA
##       date        type format separator public
## 1 20180811 RasterArray     nc             TRUE
## 2 20160216 RasterArray     nc             TRUE
## 3 20160216 RasterArray     nc             TRUE
## 4 20191122  data.frame    rds             TRUE
## 5 20170317  data.frame    csv     comma   TRUE
## 6 20020000  data.frame    csv semicolon   TRUE
##                                                                    url
## 1 https://www.earthbyte.org/paleodem-resource-scotese-and-wright-2018/
## 2           https://www.earthbyte.org/paleomap-paleoatlas-for-gplates/
## 3           https://www.earthbyte.org/paleomap-paleoatlas-for-gplates/
## 4                                              https://paleobiodb.org/
## 5                      https://github.com/UW-Macrostrat/PNAS_201702297
## 6                                 http://strata.geology.wisc.edu/jack/
##                                                                                                                                                                                                                         citation
## 1                                            Scotese, C. R. Wright, N. (2018). PALEOMAP Paleodigital Elevation Models (PaleoDEMS) for the Phanerozoic. URL: https://www.earthbyte.org/paleodem-resource-scotese-and-wright-2018/
## 2                                                             Scotese, C. R. (2016) Tutorial: PALEOMAP PaleoAtlas for GPlates and the PaleoData Plotter Program. URL: https://www.earthbyte.org/paleomap-paleoatlas-for-gplates/
## 3                                                             Scotese, C. R. (2016) Tutorial: PALEOMAP PaleoAtlas for GPlates and the PaleoData Plotter Program. URL: https://www.earthbyte.org/paleomap-paleoatlas-for-gplates/
## 4                                                                                                                          Please remember to acknowledge the Paleobiology Database (http://paleobiodb.org) in your publication.
## 5 Zaffos, A., Finnegan, S., & Peters, S. E. (2017). Plate tectonic regulation of global marine animal diversity. Proceedings of the National Academy of Sciences, 114(22), 5653\x965658. https://doi.org/10.1073/pnas.1702297114
## 6                                                                                         Sepkoski, J. J., Jablonski, D. & Foote, M. (2002). Sepkoski\x92s Genus Compendium\x97Editors Introduction & Timescale. no.363 (2002). 
##                                                     comment
## 1                                                          
## 2   \\Caradocian map omitted -  same date as Darriwillian\\
## 3                                        Resampled from 0.1
## 4                                                          
## 5                                                          
## 6 Numeric stratigraphic entries added by Wolfgang Kiessling

Currently the available variables are limited to the PaleoDEM (Scotese and Wright, 2018) and PaleoRaster series. The registry object includes the the name of the dataset, the name of the variable in question, spatial resolution, if applicable and the variables version. These are input parameters to the fetch() function, that downloads and attached the selected data. The PaleoDEMS can be downloaded using the following command. Note: the argument dat represents the downloaded dataset and var represents the desired variable. All available ages will be downloaded with the above function call.

dems <- fetch(dat="paleomap", var="dem")
## If you use the data in publications, please cite its
## reference, as well as that of the 'chronosphere' package.
## - Scotese, C. R. Wright, N. (2018). PALEOMAP Paleodigital Elevation Models (PaleoDEMS) for the Phanerozoic. URL: https://www.earthbyte.org/paleodem-resource-scotese-and-wright-2018/
dems
## class         : RasterArray 
## RasterLayer properties: 
## - dimensions  : 181, 361  (nrow, ncol)
## - resolution  : 1, 1  (x, y)
## - extent      : -180.5, 180.5, -90.5, 90.5  (xmin, xmax, ymin, ymax)
## - coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## Array properties: 
## - dimensions   : 109  (vector)
## - num. layers    : 109
## - proxy:
##          0         5        10        15        20        25        30 
##   "dem_0"   "dem_5"  "dem_10"  "dem_15"  "dem_20"  "dem_25"  "dem_30" 
##        35        40        45        50        55        60        65 
##  "dem_35"  "dem_40"  "dem_45"  "dem_50"  "dem_55"  "dem_60"  "dem_65" 
##        70        75        80        85        90        95       100 
##  "dem_70"  "dem_75"  "dem_80"  "dem_85"  "dem_90"  "dem_95" "dem_100" 
##       105       110       115       120       125       130       135 
## "dem_105" "dem_110" "dem_115" "dem_120" "dem_125" "dem_130" "dem_135" 
##       140       145       150       155       160       165       170 
## "dem_140" "dem_145" "dem_150" "dem_155" "dem_160" "dem_165" "dem_170" 
##       175       180       185       190       195       200       205 
## "dem_175" "dem_180" "dem_185" "dem_190" "dem_195" "dem_200" "dem_205" 
##       210       215       220       225       230       235       240 
## "dem_210" "dem_215" "dem_220" "dem_225" "dem_230" "dem_235" "dem_240" 
##       245       250       255       260       265       270       275 
## "dem_245" "dem_250" "dem_255" "dem_260" "dem_265" "dem_270" "dem_275" 
##       280       285       290       295       300       305       310 
## "dem_280" "dem_285" "dem_290" "dem_295" "dem_300" "dem_305" "dem_310" 
##       315       320       325       330       335       340       345 
## "dem_315" "dem_320" "dem_325" "dem_330" "dem_335" "dem_340" "dem_345" 
##       350       355       360       365       370       375       380 
## "dem_350" "dem_355" "dem_360" "dem_365" "dem_370" "dem_375" "dem_380" 
##       385       390       395       400       405       410       415 
## "dem_385" "dem_390" "dem_395" "dem_400" "dem_405" "dem_410" "dem_415" 
##       420       425       430       435       440       445       450 
## "dem_420" "dem_425" "dem_430" "dem_435" "dem_440" "dem_445" "dem_450" 
##       455       460       465       470       475       480       485 
## "dem_455" "dem_460" "dem_465" "dem_470" "dem_475" "dem_480" "dem_485" 
##       490       495       500       505       510       515       520 
## "dem_490" "dem_495" "dem_500" "dem_505" "dem_510" "dem_515" "dem_520" 
##       525       530       535       540 
## "dem_525" "dem_530" "dem_535" "dem_540"

The function returns an object of the RasterArray class. This S4 class (still under intense development) is built on RasterStacks and simulates base R arrays with RasterLayers as elements. In other words, you can organize and subset individual RasterLayers as if they were simple logical/numeric/character values that can be combined to build vectors, matrices, or multidimensional arrays. Typing the object name to the console will display the proxy object that symbolizes the structure of the RasterArray representing every RasterLayer with its name.

The dems object is a one-dimensional array (practically a vector), where every element is a time slice-specific RasterLayer (DEM). Subsetting rules are implemented to access these elements the same way as you would do it with base R objects. For instance, if you want to get the first element, you can use:

dems[1]
## class      : RasterLayer 
## dimensions : 181, 361, 65341  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : -180.5, 180.5, -90.5, 90.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## source     : /tmp/RtmpFEqfyV/1_dem_20190719/1_dem_20190719_0.nc 
## names      : dem_0 
## zvar       : z

which immediately reduces to a simple RasterLayer, so you can treat it as such:

plot(dems[1])

The elements have the ages of the DEMs as names, and can also be subsetted as regular R vectors.

plot(dems["300"])

To access additional variables, you only have to vary the dat, var and sometimes the res arguments of the package, depending on the resolutions are available.

pr <- fetch(dat="paleomap", var="paleoatlas", res=0.5)
## 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: https://www.earthbyte.org/paleomap-paleoatlas-for-gplates/
pr
## class         : RasterArray 
## RasterLayer properties: 
## - dimensions  : 360, 720  (nrow, ncol)
## - resolution  : 0.5, 0.5  (x, y)
## - extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## - coord. ref. : +proj=longlat 
## Array properties: 
## - dimensions  : 89, 3  (nrow, ncol)
## - num. layers    : 267
## - proxy:
##      paleoatlas_R       paleoatlas_G       paleoatlas_B      
## 0   "paleoatlas_0_R"   "paleoatlas_0_G"   "paleoatlas_0_B"  
## 1   "paleoatlas_1_R"   "paleoatlas_1_G"   "paleoatlas_1_B"  
## 4   "paleoatlas_4_R"   "paleoatlas_4_G"   "paleoatlas_4_B"  
## 6   "paleoatlas_6_R"   "paleoatlas_6_G"   "paleoatlas_6_B"  
## 10  "paleoatlas_10_R"  "paleoatlas_10_G"  "paleoatlas_10_B" 
## 15  "paleoatlas_15_R"  "paleoatlas_15_G"  "paleoatlas_15_B" 
## 20  "paleoatlas_20_R"  "paleoatlas_20_G"  "paleoatlas_20_B" 
## 25  "paleoatlas_25_R"  "paleoatlas_25_G"  "paleoatlas_25_B" 
## 30  "paleoatlas_30_R"  "paleoatlas_30_G"  "paleoatlas_30_B" 
## 35  "paleoatlas_35_R"  "paleoatlas_35_G"  "paleoatlas_35_B" 
## 40  "paleoatlas_40_R"  "paleoatlas_40_G"  "paleoatlas_40_B" 
## 45  "paleoatlas_45_R"  "paleoatlas_45_G"  "paleoatlas_45_B" 
## 50  "paleoatlas_50_R"  "paleoatlas_50_G"  "paleoatlas_50_B" 
## 55  "paleoatlas_55_R"  "paleoatlas_55_G"  "paleoatlas_55_B" 
## 60  "paleoatlas_60_R"  "paleoatlas_60_G"  "paleoatlas_60_B" 
## 66  "paleoatlas_66_R"  "paleoatlas_66_G"  "paleoatlas_66_B" 
## 70  "paleoatlas_70_R"  "paleoatlas_70_G"  "paleoatlas_70_B" 
## 75  "paleoatlas_75_R"  "paleoatlas_75_G"  "paleoatlas_75_B" 
## 80  "paleoatlas_80_R"  "paleoatlas_80_G"  "paleoatlas_80_B" 
## 90  "paleoatlas_90_R"  "paleoatlas_90_G"  "paleoatlas_90_B" 
## 95  "paleoatlas_95_R"  "paleoatlas_95_G"  "paleoatlas_95_B" 
## 100 "paleoatlas_100_R" "paleoatlas_100_G" "paleoatlas_100_B"
## 105 "paleoatlas_105_R" "paleoatlas_105_G" "paleoatlas_105_B"
## 110 "paleoatlas_110_R" "paleoatlas_110_G" "paleoatlas_110_B"
## 115 "paleoatlas_115_R" "paleoatlas_115_G" "paleoatlas_115_B"
## 120 "paleoatlas_120_R" "paleoatlas_120_G" "paleoatlas_120_B"
## 125 "paleoatlas_125_R" "paleoatlas_125_G" "paleoatlas_125_B"
## 130 "paleoatlas_130_R" "paleoatlas_130_G" "paleoatlas_130_B"
## 135 "paleoatlas_135_R" "paleoatlas_135_G" "paleoatlas_135_B"
## 140 "paleoatlas_140_R" "paleoatlas_140_G" "paleoatlas_140_B"
## 145 "paleoatlas_145_R" "paleoatlas_145_G" "paleoatlas_145_B"
## 150 "paleoatlas_150_R" "paleoatlas_150_G" "paleoatlas_150_B"
## 155 "paleoatlas_155_R" "paleoatlas_155_G" "paleoatlas_155_B"
## 160 "paleoatlas_160_R" "paleoatlas_160_G" "paleoatlas_160_B"
## 165 "paleoatlas_165_R" "paleoatlas_165_G" "paleoatlas_165_B"
## 170 "paleoatlas_170_R" "paleoatlas_170_G" "paleoatlas_170_B"
## 175 "paleoatlas_175_R" "paleoatlas_175_G" "paleoatlas_175_B"
## 180 "paleoatlas_180_R" "paleoatlas_180_G" "paleoatlas_180_B"
## 185 "paleoatlas_185_R" "paleoatlas_185_G" "paleoatlas_185_B"
## 190 "paleoatlas_190_R" "paleoatlas_190_G" "paleoatlas_190_B"
## 195 "paleoatlas_195_R" "paleoatlas_195_G" "paleoatlas_195_B"
## 200 "paleoatlas_200_R" "paleoatlas_200_G" "paleoatlas_200_B"
## 210 "paleoatlas_210_R" "paleoatlas_210_G" "paleoatlas_210_B"
## 220 "paleoatlas_220_R" "paleoatlas_220_G" "paleoatlas_220_B"
## 230 "paleoatlas_230_R" "paleoatlas_230_G" "paleoatlas_230_B"
## 240 "paleoatlas_240_R" "paleoatlas_240_G" "paleoatlas_240_B"
## 245 "paleoatlas_245_R" "paleoatlas_245_G" "paleoatlas_245_B"
## 250 "paleoatlas_250_R" "paleoatlas_250_G" "paleoatlas_250_B"
## 255 "paleoatlas_255_R" "paleoatlas_255_G" "paleoatlas_255_B"
## 260 "paleoatlas_260_R" "paleoatlas_260_G" "paleoatlas_260_B"
## 270 "paleoatlas_270_R" "paleoatlas_270_G" "paleoatlas_270_B"
## 275 "paleoatlas_275_R" "paleoatlas_275_G" "paleoatlas_275_B"
## 280 "paleoatlas_280_R" "paleoatlas_280_G" "paleoatlas_280_B"
## 290 "paleoatlas_290_R" "paleoatlas_290_G" "paleoatlas_290_B"
## 295 "paleoatlas_295_R" "paleoatlas_295_G" "paleoatlas_295_B"
## 300 "paleoatlas_300_R" "paleoatlas_300_G" "paleoatlas_300_B"
## 305 "paleoatlas_305_R" "paleoatlas_305_G" "paleoatlas_305_B"
## 310 "paleoatlas_310_R" "paleoatlas_310_G" "paleoatlas_310_B"
## 315 "paleoatlas_315_R" "paleoatlas_315_G" "paleoatlas_315_B"
## 320 "paleoatlas_320_R" "paleoatlas_320_G" "paleoatlas_320_B"
## 330 "paleoatlas_330_R" "paleoatlas_330_G" "paleoatlas_330_B"
## 340 "paleoatlas_340_R" "paleoatlas_340_G" "paleoatlas_340_B"
## 350 "paleoatlas_350_R" "paleoatlas_350_G" "paleoatlas_350_B"
## 360 "paleoatlas_360_R" "paleoatlas_360_G" "paleoatlas_360_B"
## 370 "paleoatlas_370_R" "paleoatlas_370_G" "paleoatlas_370_B"
## 380 "paleoatlas_380_R" "paleoatlas_380_G" "paleoatlas_380_B"
## 390 "paleoatlas_390_R" "paleoatlas_390_G" "paleoatlas_390_B"
## 395 "paleoatlas_395_R" "paleoatlas_395_G" "paleoatlas_395_B"
## 400 "paleoatlas_400_R" "paleoatlas_400_G" "paleoatlas_400_B"
## 410 "paleoatlas_410_R" "paleoatlas_410_G" "paleoatlas_410_B"
## 415 "paleoatlas_415_R" "paleoatlas_415_G" "paleoatlas_415_B"
## 420 "paleoatlas_420_R" "paleoatlas_420_G" "paleoatlas_420_B"
## 425 "paleoatlas_425_R" "paleoatlas_425_G" "paleoatlas_425_B"
## 430 "paleoatlas_430_R" "paleoatlas_430_G" "paleoatlas_430_B"
## 440 "paleoatlas_440_R" "paleoatlas_440_G" "paleoatlas_440_B"
## 445 "paleoatlas_445_R" "paleoatlas_445_G" "paleoatlas_445_B"
## 450 "paleoatlas_450_R" "paleoatlas_450_G" "paleoatlas_450_B"
## 460 "paleoatlas_460_R" "paleoatlas_460_G" "paleoatlas_460_B"
## 470 "paleoatlas_470_R" "paleoatlas_470_G" "paleoatlas_470_B"
## 480 "paleoatlas_480_R" "paleoatlas_480_G" "paleoatlas_480_B"
## 490 "paleoatlas_490_R" "paleoatlas_490_G" "paleoatlas_490_B"
## 500 "paleoatlas_500_R" "paleoatlas_500_G" "paleoatlas_500_B"
## 510 "paleoatlas_510_R" "paleoatlas_510_G" "paleoatlas_510_B"
## 520 "paleoatlas_520_R" "paleoatlas_520_G" "paleoatlas_520_B"
## 530 "paleoatlas_530_R" "paleoatlas_530_G" "paleoatlas_530_B"
## 540 "paleoatlas_540_R" "paleoatlas_540_G" "paleoatlas_540_B"
## 600 "paleoatlas_600_R" "paleoatlas_600_G" "paleoatlas_600_B"
## 690 "paleoatlas_690_R" "paleoatlas_690_G" "paleoatlas_690_B"
## 750 "paleoatlas_750_R" "paleoatlas_750_G" "paleoatlas_750_B"

This function call downloads a RasterArray of the PaleoMap Paleoatlas raster series with the resolution of 0.5 degrees. This RasterArray has 3 columns, each corresponding to one of the RGB channels.

You can refer to the 300Ma map by referring to the row of this RasterArray and then you can plot it with the function mapplot()

pa300 <- pr["300", ]
mapplot(pa300, rgb=TRUE)

Note: in the future the SpArray class will also be implemented, building R arrays from SpatialPolygonsDataFrame objects.

Merging with PaleoDB Data

Loading the data

The Paleobiology Database data for this demo will be processed in the same way as in the divDyn (Kocsis et al., 2019) R package.

## v0.8.0. See latest changes with 'news(package="divDyn")'. Default: 
## Time flows forward with 'bin', and backward with 'age' and 'slice()'.

These data are available from the ddPhanero GitHub repository, and you can import them with the following command. This imports a copy downloaded on 2019-05-31 and assigns it to the name dat.

load(url(
"https://github.com/divDyn/ddPhanero/raw/master/data/PaleoDB/2019-05-31_paleoDB.RData"
))

These data are prepared the same way as in (Kocsis et al., 2019) with a single function call:

# load the function
source("https://github.com/divDyn/assets/raw/master/examples/misc/preparation.R")

# run it
dat <- prepare(dat)

The data are now filtered for marine environment and are assigned to stages.

Reconstructing the coordinates to match the plots

In order to use the provided maps, the data has to be rotated back to their paleopositions, with the same reconstruction model.

For the sake of this demonstration, we will focus on time slices 42 and 43 (Gzhelian and Asselian stages). But how do we know which reconstructions to use? The mean ages of the slices are not the same as the reconstruction dates for maps/environmental variables, so we have to approximate and use the reconstruction which is closest to the stage mid age.

The chronosphere package provides a function matchtime() that does exactly that: it takes a vector of ages or a set of features with ages given as names (such as dems and maps above) and reorders it to best match another given vector of ages. For instance, the dems object is organized backwards in time:

names(dems)
##   [1] "0"   "5"   "10"  "15"  "20"  "25"  "30"  "35"  "40"  "45"  "50" 
##  [12] "55"  "60"  "65"  "70"  "75"  "80"  "85"  "90"  "95"  "100" "105"
##  [23] "110" "115" "120" "125" "130" "135" "140" "145" "150" "155" "160"
##  [34] "165" "170" "175" "180" "185" "190" "195" "200" "205" "210" "215"
##  [45] "220" "225" "230" "235" "240" "245" "250" "255" "260" "265" "270"
##  [56] "275" "280" "285" "290" "295" "300" "305" "310" "315" "320" "325"
##  [67] "330" "335" "340" "345" "350" "355" "360" "365" "370" "375" "380"
##  [78] "385" "390" "395" "400" "405" "410" "415" "420" "425" "430" "435"
##  [89] "440" "445" "450" "455" "460" "465" "470" "475" "480" "485" "490"
## [100] "495" "500" "505" "510" "515" "520" "525" "530" "535" "540"

The time scale we use for the analyses (stages) on the other hand, has ages in the opposite order, and has totally different values:

stages$mid
##  [1] 567.50000 554.50000 545.00000 535.00000 525.00000 517.50000 511.50000
##  [8] 506.75000 502.50000 498.75000 495.50000 491.75000 487.45000 481.55000
## [15] 473.85000 468.65000 462.85000 455.70000 449.10000 444.30000 442.10000
## [22] 439.65000 435.95000 431.95000 428.95000 426.50000 424.30000 421.10000
## [29] 415.00000 409.20000 400.45000 390.50000 385.20000 377.45000 365.55000
## [36] 352.80000 338.80000 327.05000 319.20000 311.10000 305.35000 301.30000
## [43] 297.20000 292.80000 284.70000 275.80000 270.55000 266.95000 262.50000
## [50] 257.03500 253.17000 251.68500 249.20000 244.60000 239.50000 232.50000
## [57] 218.25000 204.90000 200.30000 195.05000 186.75000 178.40000 172.20000
## [64] 169.30000 167.20000 164.80000 160.40000 154.70000 148.55000 142.40000
## [71] 136.35000 131.15000 127.20000 119.00000 106.75000  97.20000  91.85000
## [78]  88.05000  84.95000  77.85000  69.05000  63.80000  58.80000  51.90000
## [85]  44.55000  39.65000  35.95000  31.00000  25.56500  19.50000  13.79500
## [92]   8.47650   3.96050   1.29985   0.00585

Using the dems object as the first, and stages$mid as the second argument of matchtime() will reorder the elements of dems to match stages$mid.

demord <- matchtime(dems, stages$mid)

The resulting object (demord)has the same class as the original one (dems), but the order of elements are changed:

names(demord)
##  [1] "540" "540" "540" "535" "525" "515" "510" "505" "500" "500" "495"
## [12] "490" "485" "480" "475" "470" "465" "455" "450" "445" "440" "440"
## [23] "435" "430" "430" "425" "425" "420" "415" "410" "400" "390" "385"
## [34] "375" "365" "355" "340" "325" "320" "310" "305" "300" "295" "295"
## [45] "285" "275" "270" "265" "260" "255" "255" "250" "250" "245" "240"
## [56] "230" "220" "205" "200" "195" "185" "180" "170" "170" "165" "165"
## [67] "160" "155" "150" "140" "135" "130" "125" "120" "105" "95"  "90" 
## [78] "90"  "85"  "80"  "70"  "65"  "60"  "50"  "45"  "40"  "35"  "30" 
## [89] "25"  "20"  "15"  "10"  "5"   "0"   "0"

The utility of the reorganization is in the referencing of subsets. Now we can use the slice ID as index for subsetting a single map. For instance, slice number 42 is not only the 42th row in the stages object, but the relevant map for this stage is also the 42th element of the map container:

stages$mid[42]
## [1] 301.3
names(demord)[42]
## [1] "300"
demord[42]
## class      : RasterLayer 
## dimensions : 181, 361, 65341  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : -180.5, 180.5, -90.5, 90.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## source     : /tmp/RtmpFEqfyV/1_dem_20190719/1_dem_20190719_300.nc 
## names      : dem_300 
## zvar       : z

Coordinate rotation

To reconstruct the coordinates, it is faster to use collection-level information only. First, let’s select the collection id, longitude and latitude columns of the time slices 42 and 43, and remove duplicate entries.

# subset
interColl <- dat[dat$stg%in%c(42,43),c("collection_no","lng", "lat", "stg", "paleolat" ,"paleolng")]
# duplicates
interColl <- unique(interColl)

To keep things organized, it is a good idea to record which reconstruction is used for which collection.

interColl$mapage <- names(demord)[interColl$stg]

The coordinate columns of the data.frame interColl can be directly fed to the reconstruct() function that sends the coordinates to the GPlates Web Service, which will rotate them back to the given ages. All data in this tutorial uses C. Scotese’s model, therefore the model argument of the reconstruct() function can be left to the default PALEOMAP value. For a single collection/coordinate, we would use the following function call:

reconstruct(interColl[1, c("lng", "lat"), drop=FALSE], age=300)
##      paleolong paleolat
## [1,]     18.84    25.69

Supplying multiple ages to the age argument will reconstruct all the provided collections to all given ages.

reconstruct(interColl[1:2, c("lng", "lat")], age=c(300, 295))
## $`300`
##      paleolong paleolat
## [1,]     18.84    25.69
## [2,]     18.84    25.69
## 
## $`295`
##      paleolong paleolat
## [1,]     18.88     27.2
## [2,]     18.88     27.2

The results are organized in a list, where every element corresponds to an age given in the age argument - the default behaviour of the function. As we have a specific age to which we want to reconstruct every coordinate, we can use the option enumerate=FALSE. This will reconstruct every row to the corresponding age in the age vector, which we have in mapage.

newCoords <- reconstruct(interColl[, c("lng", "lat")], age=interColl[, "mapage"], enumerate=FALSE, verbose=FALSE)

# rename the coordinates
colnames(newCoords) <- c("plng", "plat")

# and join data with the rest
allColl <- cbind(interColl,newCoords)

These paleocoordinates are almost the same as the ones provided by the PBDB, but they use an updated reconstruction model:

plot(allColl$plat, allColl$paleolat, 
  xlab="reconstruct() function", ylab="downloaded with the PBDB")

We can now finally plot the PBDB collections on the reconstructed maps.

mapplot(pr["300",], rgb=TRUE)

# collections in this bin
coll42 <- allColl[allColl$stg==42,]

# the points
points(coll42$plng, coll42$plat, pch=3, col="red")

Extracting values

The coordinates can be used to extract values from the RasterLayers such as an element of demord.

# the 300ma dem
dem42 <- demord[42]

# extracting values with extact()
coll42$elev <- extract(dem42, coll42[, c("plng", "plat")])

You can iterate this procedure with the extract() method implemented for the RasterArray class. This will extract values from multiple RasterLayers in a RasterArray, given an index vector is provided for every coordinate pair.

# extracting values with extact()
allColl$elev <- extract(demord, allColl, by="mapage")

The collection table with the newly added elevation estimates can now be merged with the occurrence data. This is the workflow that will be used with other raster-type environmental variables (GCM reconstructions).

References

Kocsis, Á. T., Reddin, C. J., Alroy, J., & Kiessling, W. (2019). The R package divDyn for quantifying diversity dynamics using fossil sampling data. Methods in Ecology and Evolution, 10(5), 735–743. https://doi.org/10.1111/2041-210X.13161

Seton, M., Mueller, R. D., Zahirovic, S., Gaina, C., Torsvik, T. H., Shephard, G., … Chandler, M. (2012). Global continental and ocean basin reconstructions since 200 Ma. Earth-Science Reviews, 113(3–4), 212–270. https://doi.org/10.1016/j.earscirev.2012.03.002

Scotese, C. and Wright, N. (2018) PaleoDEM Resource – https://www.earthbyte.org/paleodem-resource-scotese-and-wright-2018/

Previous
Next