install.packages(c("chronosphere", "sf", "lwgeom", "via", "divDyn"))Plotting and calculating shelf area through the Phanerozoic with the Paleomap Paleocoastlines
Ádám T. Kocsis 
Preparations
This guide will go through the steps of deriving a time series of flooded shelf area through time using the PaleoMAP Paleocoastlines data product. You have to make sure that the relevant R packages are installed. These are:
| Package | Purpose |
|---|---|
chronosphere |
Accessing original data |
sf |
Vector GIS operations |
lwgeom |
Area calculation (sf dependency) |
via |
Structuring the vector data |
divDyn |
Geological timescale |
You can install the most-up-to date versions of these package with:
Binary distributions of the sf package on Windows and MacOS include the GEOS, GDAL and PROJ libraries. On GNU/Linux operating systems these need to be installed separately so sf can function properly.
1. The Paleomap Paleocoastlines data
This example will use data from the PaleoMAP project a set of reconstructions of the paleocoastlines. These reconstructions are based on the PaleoMAP PaleoDEMs - reconstructed Digital Elevation Models for the Phanerozoic (Scotese and Wright 2018) - as well as on occurrences coming from the Paleobiology Database. The Paleocoastlines meant to express likely maximum transgressive surfaces for the intervals. These are discussed in detail in the paper describing the data (Kocsis and Scotese 2021).
Downloading data using the chronosphere
The Paleocoastlines data product was published on Zenodo, and it was added to the chronosphere for easier access and importing (alternatively, you can download the data directly from its Zenodo repository). This single line of computer code can be used to download and import that data into R:
# Download data from the chronosphere
# set verbose = TRUE to see additional information
pm <- chronosphere::fetch(src="paleomap", ser="paleocoastlines", ver="7", verbose=FALSE)Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUEBesides sf that defines vector data classes, this item requires the via package. If these packages are not available, the chronosphere::fetch() call will throw an error. If they are available, then the packages will be attached automatically, when the item is loaded (hence the message, above).
Exploring the available data
The via package organizes these data in an array-like structure for easier subsetting.
This is an SfcArray-class object, with 2 dimensions that correspond to the their ages and what these vector data represent:
# show method: overiew
pmclass : SfcArray
Element properties:
- class : sfc_MULTIPOLYGON sfc
- geodetic CRS: WGS 84
Array properties:
- dimensions : 81, 2 (nrow, ncol)
- num. layers : 162
- missing : 0
- proxy:
margin coast
0 "0Ma_CM_v7.shx" "0Ma_CS_v7.shx"
5 "5Ma_CM_v7.shx" "5Ma_CS_v7.shx"
10 "10Ma_CM_v7.shx" "10Ma_CS_v7.shx"
15 "15Ma_CM_v7.shx" "15Ma_CS_v7.shx"
20 "20Ma_CM_v7.shx" "20Ma_CS_v7.shx"
25 "25Ma_CM_v7.shx" "25Ma_CS_v7.shx"
30 "30Ma_CM_v7.shx" "30Ma_CS_v7.shx"
35 "35Ma_CM_v7.shx" "35Ma_CS_v7.shx"
40 "40Ma_CM_v7.shx" "40Ma_CS_v7.shx"
45 "45Ma_CM_v7.shx" "45Ma_CS_v7.shx"
50 "50Ma_CM_v7.shx" "50Ma_CS_v7.shx"
60 "60Ma_CM_v7.shx" "60Ma_CS_v7.shx"
65 "65Ma_CM_v7.shx" "65Ma_CS_v7.shx"
70 "70Ma_CM_v7.shx" "70Ma_CS_v7.shx"
80 "80Ma_CM_v7.shx" "80Ma_CS_v7.shx"
85 "85Ma_CM_v7.shx" "85Ma_CS_v7.shx"
90 "90Ma_CM_v7.shx" "90Ma_CS_v7.shx"
95 "95Ma_CM_v7.shx" "95Ma_CS_v7.shx"
105 "105Ma_CM_v7.shx" "105Ma_CS_v7.shx"
120 "120Ma_CM_v7.shx" "120Ma_CS_v7.shx"
125 "125Ma_CM_v7.shx" "125Ma_CS_v7.shx"
130 "130Ma_CM_v7.shx" "130Ma_CS_v7.shx"
135 "135Ma_CM_v7.shx" "135Ma_CS_v7.shx"
140 "140Ma_CM_v7.shx" "140Ma_CS_v7.shx"
150 "150Ma_CM_v7.shx" "150Ma_CS_v7.shx"
155 "155Ma_CM_v7.shx" "155Ma_CS_v7.shx"
160 "160Ma_CM_v7.shx" "160Ma_CS_v7.shx"
165 "165Ma_CM_v7.shx" "165Ma_CS_v7.shx"
170 "170Ma_CM_v7.shx" "170Ma_CS_v7.shx"
180 "180Ma_CM_v7.shx" "180Ma_CS_v7.shx"
185 "185Ma_CM_v7.shx" "185Ma_CS_v7.shx"
195 "195Ma_CM_v7.shx" "195Ma_CS_v7.shx"
200 "200Ma_CM_v7.shx" "200Ma_CS_v7.shx"
205 "205Ma_CM_v7.shx" "205Ma_CS_v7.shx"
220 "220Ma_CM_v7.shx" "220Ma_CS_v7.shx"
230 "230Ma_CM_v7.shx" "230Ma_CS_v7.shx"
240 "240Ma_CM_v7.shx" "240Ma_CS_v7.shx"
245 "245Ma_CM_v7.shx" "245Ma_CS_v7.shx"
250 "250Ma_CM_v7.shx" "250Ma_CS_v7.shx"
255 "255Ma_CM_v7.shx" "255Ma_CS_v7.shx"
260 "260Ma_CM_v7.shx" "260Ma_CS_v7.shx"
265 "265Ma_CM_v7.shx" "265Ma_CS_v7.shx"
270 "270Ma_CM_v7.shx" "270Ma_CS_v7.shx"
275 "275Ma_CM_v7.shx" "275Ma_CS_v7.shx"
285 "285Ma_CM_v7.shx" "285Ma_CS_v7.shx"
295 "295Ma_CM_v7.shx" "295Ma_CS_v7.shx"
300 "300Ma_CM_v7.shx" "300Ma_CS_v7.shx"
305 "305Ma_CM_v7.shx" "305Ma_CS_v7.shx"
310 "310Ma_CM_v7.shx" "310Ma_CS_v7.shx"
320 "320Ma_CM_v7.shx" "320Ma_CS_v7.shx"
325 "325Ma_CM_v7.shx" "325Ma_CS_v7.shx"
340 "340Ma_CM_v7.shx" "340Ma_CS_v7.shx"
355 "355Ma_CM_v7.shx" "355Ma_CS_v7.shx"
365 "365Ma_CM_v7.shx" "365Ma_CS_v7.shx"
375 "375Ma_CM_v7.shx" "375Ma_CS_v7.shx"
385 "385Ma_CM_v7.shx" "385Ma_CS_v7.shx"
390 "390Ma_CM_v7.shx" "390Ma_CS_v7.shx"
400 "400Ma_CM_v7.shx" "400Ma_CS_v7.shx"
410 "410Ma_CM_v7.shx" "410Ma_CS_v7.shx"
415 "415Ma_CM_v7.shx" "415Ma_CS_v7.shx"
420 "420Ma_CM_v7.shx" "420Ma_CS_v7.shx"
425 "425Ma_CM_v7.shx" "425Ma_CS_v7.shx"
430 "430Ma_CM_v7.shx" "430Ma_CS_v7.shx"
435 "435Ma_CM_v7.shx" "435Ma_CS_v7.shx"
440 "440Ma_CM_v7.shx" "440Ma_CS_v7.shx"
445 "445Ma_CM_v7.shx" "445Ma_CS_v7.shx"
450 "450Ma_CM_v7.shx" "450Ma_CS_v7.shx"
455 "455Ma_CM_v7.shx" "455Ma_CS_v7.shx"
465 "465Ma_CM_v7.shx" "465Ma_CS_v7.shx"
470 "470Ma_CM_v7.shx" "470Ma_CS_v7.shx"
475 "475Ma_CM_v7.shx" "475Ma_CS_v7.shx"
480 "480Ma_CM_v7.shx" "480Ma_CS_v7.shx"
485 "485Ma_CM_v7.shx" "485Ma_CS_v7.shx"
490 "490Ma_CM_v7.shx" "490Ma_CS_v7.shx"
495 "495Ma_CM_v7.shx" "495Ma_CS_v7.shx"
500 "500Ma_CM_v7.shx" "500Ma_CS_v7.shx"
505 "505Ma_CM_v7.shx" "505Ma_CS_v7.shx"
510 "510Ma_CM_v7.shx" "510Ma_CS_v7.shx"
515 "515Ma_CM_v7.shx" "515Ma_CS_v7.shx"
525 "525Ma_CM_v7.shx" "525Ma_CS_v7.shx"
535 "535Ma_CM_v7.shx" "535Ma_CS_v7.shx"The available dimensions to this structure can also be accessed with the dimnames() function. The data are available for 81 time intervals through the Phanerozoic, and they represent either the approximate margins of the plates (margin, derived from the 1400m isobath of the PaleoDEMs (Scotese and Wright 2018), and the reconstructed coastline (coast). The ages are in millions of years.
Strategy
In order to make the calculations transparent, it is good practice to:
- Work out the shelf area calculation for one interval, and then
- Repeat for all intervals, to get a time series.
2. Calculating flooded shelf area for a specific interval
Selecting an interval
To prototype our analyses, let’s select an interval first, for instance the Priabonian stage of the Eocene (~35 Ma)!
The margins and the coastlines can be accessed with regular array syntax using character subscripts. Note that the ages have to be submitted as characters, otherwise you subset to the wrong interval!
# rows, then columns (35Ma)
margin <- pm["35", "margin"]
coast <- pm["35", "coast"]The SfcArray objects get reduced to their elements: simple geometry collections from the sf package:
coastGeometry set for 96 features
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 82.82331
Geodetic CRS: WGS 84
First 5 geometries:POLYGON ((180 74.04842, 180 66.41276, 179.9851 ...POLYGON ((180 -84.99437, 180 -90, -180 -90, -18...POLYGON ((164.51 -35.86942, 164.6075 -35.88057,...POLYGON ((162.5601 -34.67409, 162.6511 -34.6827...POLYGON ((159.2816 78.48408, 159.4261 78.49931,...Plotting the data
These can be plotted like any other sf geometry objects. For instance, we can plot these objects simply with:
# use add in the second call to overplot
plot(margin, col="#87cef6", border=NA)
plot(coast, col="#94391cAA", border=NA, add=TRUE)
Hint for plotting multiple intervals
For visualizing all the reconstructions, it is easiest to:
1. Open a pdf() device
2. Write a for loop for all the intervals (e.g. by 1:nrow(pm))
3. Subset and plot the data
4. Close device with dev.off()
Calculating the shelf vector
The shelf area is the geometric difference of the margin (margin) and the coastline object (coast), which can be calculated with the sf::st_difference function. However, to do this calculation, the internals of the vector objects need to be rearranged. Both the coast and the margin objects (also in other intervals!) are composed of many single-part POLYGON items, which sf::st_difference would utilize in an elaborate way - trying to find pairwise differences - which is not what we want.
Therefore it is necessary to first transform the geometry collections so that they are a single, multipart polygon (MULTIPOLYGON) geometry. We can do this with the sf::st_union function. However, since these geometries were hand-drawn, they are not perfect (tiny loops and intersections might be present). The default method to do these calculations (which relies on the S2 spherical library) will not work for all intervals - but they can be done in longitude-latitude planar space. To use this, we first have to switch off the spherical calculations:
sf::sf_use_s2(FALSE)Spherical geometry (s2) switched offNow we can transform the objects with st_union:
marginUnified <- sf::st_union(margin)although coordinates are longitude/latitude, st_union assumes that they are
planarcoastUnified <- sf::st_union(coast)although coordinates are longitude/latitude, st_union assumes that they are
planar
Warning!
The messages here appear because we decided not to use S2. This is a good indication that we have to inspect the results!
The simplest way to check the results is to re-plot the data - now using the unified geometries:
# use add in the second call to overplot
plot(marginUnified, col="#87cef6", border=NA)
plot(coastUnified, col="#94391cAA", border=NA, add=TRUE)
You should not see any artifacts here. This is one more reason to plot the entire series of data (which we will do at a later point).
Now that we have the unified geometries, their difference can be calculated:
shelf <- sf::st_difference(marginUnified, coastUnified)although coordinates are longitude/latitude, st_difference assumes that they
are planarThe message should make more sense now! Again, it is a good idea to plot the result. Let’s draw the shelf outlines and the covered area over our existing map, with a bright, but transparent red color!
plot(shelf, col="#ff000055", border="#ff0000", add=TRUE, lwd=0.5)
If the calculations are correct, you should see a perfect match with the previous blue color.
Calculating the shelf area
The only step that remains is the calcuation of the area. The sf package is very helpful in this regard, the st_area function returns the area of this total vector. If you consult the help page of st_area (?st_area) you can see that depending on which library you want to use, you can either use the default S2 implementation which sometimes does not work at all for geometries that are not perfectly clean. If you don’t want to use S2 then you need to make sure that lwgeom is installed! Note that R will prompt you for this, if this is not installed.
sf::st_area(shelf)4.578451e+13 [m^2]Returns the result of the area calculation in square meters.
Note
You can check the results with S2 (toggling sf_use_s2(TRUE)), which will give you a slighty different result!
3. Repeating the calculation for all intervals
Now that we know how to calculate the shelf area for one time slice, let’s calculate it for all of them!
Set up iterations
First we need to set up a container to store our results. It is good practice to pre-define the structure for this container, and to include for it as many attributes as possible:
# shelf area container
shelfArea<- rep(NA, nrow(pm))
names(shelfArea) <- rownames(pm)It should be clear which value belongs with which time interval:
shelfArea 0 5 10 15 20 25 30 35 40 45 50 60 65 70 80 85 90 95 105 120
NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
125 130 135 140 150 155 160 165 170 180 185 195 200 205 220 230 240 245 250 255
NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
260 265 270 275 285 295 300 305 310 320 325 340 355 365 375 385 390 400 410 415
NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
420 425 430 435 440 445 450 455 465 470 475 480 485 490 495 500 505 510 515 525
NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
535
NA To visually inspect the results of the geometry calculations, it is good idea to plot all the maps. This we can do by setting up a pdf.
pdf(file="allTheMaps.pdf", width=20, height=12)We will populate this pdf all the maps!
Hint
Use the file argument to provide path to the target file. Leaving it as is will create the pdf file in the current working directory (which you can see with getwd()).
Run iteration
Once the preparations are done we just need to take the previous code and iterate it
# iterate for all the intervals
for(i in 1:nrow(pm)){
# use position index to subset, rather than names
margin <- pm[i, "margin"]
coast <- pm[i, "coast"]
# unify geometries
marginUnified <- sf::st_union(margin)
coastUnified <- sf::st_union(coast)
# calculate shelf
shelf <- sf::st_difference(marginUnified, coastUnified)
# plot everything
title <- paste0("Age: ", rownames(pm)[i], "Ma")
plotOrig(marginUnified, col="#87cef6", border=NA, main=title)
plotOrig(coastUnified, col="#94391cAA", border=NA, add=TRUE)
plotOrig(shelf, col="#ff000055", border="#ff0000", add=TRUE, lwd=0.5)
# calculate shelf area - and store
shelfArea[i] <- sf::st_area(shelf)
}If you did make a pdf, don’t forget to switch off the device when the iteration is done!
dev.off()You can inspect the plotted pdf here:
The results of the area calculations are in the shelfArea object:
shelfArea 0 5 10 15 20 25
5.172011e+13 4.415912e+13 4.461499e+13 4.474639e+13 4.927784e+13 3.935787e+13
30 35 40 45 50 60
4.455750e+13 4.578451e+13 5.090562e+13 5.584635e+13 7.363605e+13 6.193659e+13
65 70 80 85 90 95
7.002828e+13 7.559231e+13 9.745919e+13 8.508077e+13 8.393429e+13 7.464471e+13
105 120 125 130 135 140
6.595719e+13 8.045058e+13 6.087005e+13 5.102032e+13 4.826494e+13 5.651167e+13
150 155 160 165 170 180
5.480478e+13 5.758006e+13 7.443396e+13 6.899570e+13 5.997792e+13 4.311000e+13
185 195 200 205 220 230
4.636160e+13 5.196092e+13 6.567870e+13 6.204666e+13 7.150306e+13 6.597844e+13
240 245 250 255 260 265
6.799099e+13 6.972395e+13 6.099348e+13 6.616452e+13 7.852182e+13 7.611595e+13
270 275 285 295 300 305
7.242381e+13 7.376274e+13 7.870736e+13 7.916209e+13 7.783643e+13 7.346216e+13
310 320 325 340 355 365
6.768658e+13 7.836577e+13 8.224693e+13 8.455316e+13 9.161060e+13 9.778632e+13
375 385 390 400 410 415
1.025203e+14 1.037006e+14 9.046260e+13 9.566433e+13 9.094015e+13 9.893962e+13
420 425 430 435 440 445
1.047054e+14 1.001896e+14 1.016129e+14 8.704884e+13 8.503156e+13 8.878136e+13
450 455 465 470 475 480
9.370090e+13 8.867796e+13 9.266572e+13 9.008694e+13 1.058027e+14 9.727434e+13
485 490 495 500 505 510
9.011198e+13 8.728989e+13 8.654702e+13 8.593686e+13 8.544693e+13 8.081970e+13
515 525 535
8.288370e+13 8.335465e+13 7.963508e+13 Plotting the results
You can already visualize the results by doing a scatterplot between the names of shelfArea (cast as numeric) and the calculated values.
plot(as.numeric(names(shelfArea)), shelfArea,
ylab=expression("Shelf Area (m"^2*")"),
xlab="Age (Ma)")
But this plot is a bit sub-par quality. For instance, for instance, the exact values are difficult to evaluate. It makes more sense to evaluate shelf area as a proportion to Earth’s total surface area, which is approximately: 509 600 000 km^2. All we need to do this is to divide our values with this total area:
# the proportion of area
shelfProp <- shelfArea/509600000000000
# should be between 0 and 1
range(shelfProp)[1] 0.07723287 0.20761905Also, To put the ups and downs of the variable into geological context, we need to show the geological timescale, which can be visualized in a number of ways. One way to do it is to use a built-in timescale object of the divDyn package:
data(stages, package="divDyn")
str(stages)'data.frame': 95 obs. of 13 variables:
$ sys : chr "E" "E" "E" "Cm" ...
$ system : chr "Ediacaran" "Ediacaran" "Ediacaran" "Cambrian" ...
$ series : chr "Upper Ediacaran" "Upper Ediacaran" "Upper Ediacaran" "Terreneuvian" ...
$ stage : chr "Avalon Assemblage" "White Sea Assemblage" "Nama Assemblage" "Fortunian" ...
$ short : chr "Ava" "Whs" "Nam" "For" ...
$ bottom : num 580 560 550 539 529 ...
$ mid : num 570 555 544 534 525 ...
$ top : num 560 550 539 529 521 ...
$ dur : num 20 10 11.2 9.8 8 6.5 5.5 4.5 4 3.5 ...
$ stg : int 1 2 3 4 5 6 7 8 9 10 ...
$ systemCol: chr "#fddd88" "#fddd88" "#fddd88" "#94ac72" ...
$ seriesCol: chr "#fddd88" "#fddd88" "#fddd88" "#9db989" ...
$ col : chr "#fcd589" "#fdd587" "#fed583" "#a9be93" ...This object, based on the GTS2020 timescale (Gradstein, Ogg, and Schmitz 2020), can be visualized with the divDyn::tsplot function, that we can use to draw a canvas:
divDyn::tsplot(stages, shading="sys", xlim=4:95, ylim=c(0.05, 0.25),
boxes.col="systemCol", boxes="sys",
ylab="Shelf Area (proportion of Earth's surface)"
)
and then put on the variable with lines:
lines(as.numeric(names(shelfProp)), shelfProp)
There it is, an estimate of shelf area through the Phanerozoic.
Using pre-calculated values
Note that this result is directly available from the chronosphere as:
# download data - delete verbose=FALSE to see additional information
areas <- chronosphere::fetch(src="paleomap", ser="areas", ver="7", verbose=FALSE)
# the total shelf area
areas$shelf_sum [1] 0.10144141 0.08657746 0.08747932 0.08776351 0.09661238 0.07720141
[7] 0.08741413 0.08976380 0.09986567 0.10960176 0.14444422 0.12142108
[13] 0.13727995 0.14829715 0.19111118 0.16684207 0.16459239 0.14634500
[19] 0.12919689 0.15780377 0.11942320 0.10007354 0.09472629 0.11084055
[25] 0.10750778 0.11294771 0.14596706 0.13530968 0.11761877 0.08456972
[31] 0.09095216 0.10195654 0.12881972 0.12166344 0.14022217 0.12948798
[37] 0.13331281 0.13678062 0.11962019 0.12977527 0.15400095 0.14935350
[43] 0.14203516 0.14466921 0.15422069 0.15533068 0.15268436 0.14412440
[49] 0.13283829 0.15377673 0.16138176 0.16584786 0.17968371 0.19180612
[55] 0.20098873 0.20334409 0.17742874 0.18760924 0.17842256 0.19406733
[61] 0.20532905 0.19645994 0.19926132 0.17070369 0.16672680 0.17410159
[67] 0.18370655 0.17384623 0.18169026 0.17657554 0.20743676 0.19069382
[73] 0.17666445 0.17116260 0.16971864 0.16849756 0.16754278 0.15831004
[79] 0.16245405 0.16336793 0.15601866You can compare this to our result by inspecting the absolute difference of the two vectors. Note that there are some differences (calcualtion were done with the old rgeos package, and there is some rounding involved), but these are all extremely small (lower than 0.001):
all(abs((areas$shelf_sum - shelfProp))<0.001)[1] TRUEReferences
Gradstein, Felix M, James George Ogg, and Mark D Schmitz. 2020. The Geologic Time Scale 2020.
Kocsis, Ádám T., and Christopher R. Scotese. 2021. “Mapping Paleocoastlines and Continental Flooding During the Phanerozoic.” Earth-Science Reviews 213 (February): 103463. https://doi.org/10.1016/j.earscirev.2020.103463.
Scotese, Christopher R, and Nicky Wright. 2018. “PALEOMAP Paleodigital Elevation Models (PaleoDEMS) for the Phanerozoic.” https://www.earthbyte.org/paleodem-resource-scotese-and-wright-2018/.