Using data from the PaleoMAP project with chronosphere 0.1.3
Using data from the PaleoMAP project with chronosphere 0.1.3
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.
## Loading required package: raster## Loading required package: spPlate 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 : 85The 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 : 86The 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 chronosphere.version
## 1 20180811 RasterArray nc TRUE 0.2.0
## 2 20160216 RasterArray nc TRUE 0.2.0
## 3 20160216 RasterArray nc TRUE 0.2.0
## 4 20191122 data.frame rds TRUE 0.2.0
## 5 20170317 data.frame csv comma TRUE 0.2.0
## 6 20020000 data.frame csv semicolon TRUE 0.2.0
## 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 KiesslingCurrently 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/Rtmp09AWUb/1_dem_20190719/1_dem_20190719_0.nc
## names : dem_0
## zvar : zwhich 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)
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.00585Using 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.3names(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/Rtmp09AWUb/1_dem_20190719/1_dem_20190719_300.nc
## names : dem_300
## zvar : zCoordinate 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.69Supplying 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.2The 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/
)
Ádám T. Kocsis
Postdoc
A postdoctoral research fellow, with interests in global scale biogeographical and ecological processes, programming, philosophy and music.