This is to some extent a follow-up to a previous post on measuring oceanic areas in R. I once again in the past months had to measure oceanic areas but this time in the past, based on paleogeographic (or in this case paleotopographic) maps. I found myself using a widely different method from the one detailed in the previous post, so I thought it would be worth showing it here.

For the purpose of this exercise, we will use the paleotopographic maps known as PaleoDEM of Scotese & Wright (2018), as implemented in Adam Kocsis and Nussaibah Raja's nifty package chronosphere to make data acquisition simpler.

library(raster)
library(chronosphere)
paleodems <- fetch(dat="paleomap", var="dem", res=1, ver="20190719")

Yes that's all there is to it: we just extracted a stack of rasters showing paleotopology for each geological stage from 540Ma onwards! Since we are working with rasters here, we will have to use tools adapted to them rather than tools adapted to shapefiles as in last year's post. Shapefiles describe points, lines or polygons while a raster is basically a grid of cells having a value each.

Let's say we want to compute the area of seafloor above the Carbonate Compensation Depth in the late Eocene for reasons (which should hopefully become obvious in a few months).

ccd <- -3900 #Here for the example we use the one computed in Pälike et al. 2012 for the Equatorial Pacific.
w <- which(names(paleodems)==35) #Which raster corresponds to 35Ma?
eocene_topo <- paleodems[[w]]
bathymetry <- getValues(eocene_topo)
area_cells <- getValues(area(eocene_topo))
sum(area_cells[bathymetry<=0 & bathymetry>=ccd]) #Sum the area of cells below 0m and above the CCD
#136842002 km2, i. e. ~136.8 millions km2

Now let's say we want to only look at the Southern Ocean (which here doesn't makes a lot of sense given we are using the equatorial depth of the CCD but just for the exercise:

SO <- extent(-180,180,-90,-60) #Oversimplified as being anything below 60°S
eocene_so  <- crop(eocene_topo, SO)
bath_so  <- getValues(eocene_so)
area_so  <- getValues(area(eocene_so))
sum(area_so[bath_so <=0 & bath_so>=ccd])
#18642768, i. e. ~18.6 millions km2

eocene_contour <- rasterToContour(eocene_topo, levels=c(0,ccd))
plot(eocene_topo)
plot(eocene_contour,add=TRUE)