2. Reconstructions using modern coordinates (online)

Out of the box, the package relies on the GPlates Web Service (GWS), an online service that executes paleogeographic rotations using data provided in URLs. rgplates (1) sends the data to the GWS, which (2) calculates paleocoordinates, and then (3) rgplates reads in the returned result. For this online method, you have to be connected to the internet.

To start, you have to attach the package:

library(rgplates)
#> Loading required package: sf
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE

When attached, rgplates automatically loads the Simple Features for R (sf) package, a standard R package used for processing vector spatial data.

Reconstructing plates

All paleocoordinate rotations are executed with the reconstruct() function. The default model used by this is the Merdith et al. 2021, described here. (If you look at the reference of reconstruct(), you can see that by default model="MERDITH2021").

Every tectonic model relies on plates, which are rotated on the surface of Earth. Plate positions) can be reconstructed to any age that the model covers. The present-day positions of the plates can be queried with the reconstruct() function, with a string passed (static_plates or plate_polygons, depending on the model) as the first argument, and the target age set to 0 (in Ma).

pl0 <- reconstruct("static_polygons", age=0)
pl0
#> Simple feature collection with 2391 features and 0 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -179.99 ymin: -89.99 xmax: 179.99 ymax: 89.99
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>                          geometry
#> 1  POLYGON ((112.0797 39.2348,...
#> 2  POLYGON ((119.2562 36.6774,...
#> 3  POLYGON ((130.47 41.99, 130...
#> 4  POLYGON ((104.7254 40.4621,...
#> 5  POLYGON ((102.6257 37.6989,...
#> 6  POLYGON ((104.6209 38.9799,...
#> 7  POLYGON ((94.0783 73.2609, ...
#> 8  POLYGON ((122.0576 66.9608,...
#> 9  POLYGON ((119.9099 61.8243,...
#> 10 POLYGON ((105.7441 72.3461,...

The rgplates package relies on sf to handle vector spatial data - like the plates. You can do anything with this that you normally can with and sf object: manipulate it, change its map projection, export or use it for calculations. These data are in standard equirectangular projection, registered with longitude and latitude data in the WGS84 Coordinate Reference System.

You can plot the distribution of plates with the plot() function. To focus on the spatial data and not the attributes of the features (which are technical in nature), you can plot the geometries of this object.

plot(pl0$geometry)

Setting the age argument allows you to access the state of the plate tectonic configuration at past points in time. The age argument accepts dates in millions of years. To reconstruct the position of the plates at around the Early/Late Cretaceous boundary (approx. 100Ma), you have to set the age argument to 100:

pl100 <- reconstruct("static_polygons", age=100)
pl100
#> Simple feature collection with 1042 features and 0 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -179.99 ymin: -89.99 xmax: 179.99 ymax: 89.99
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>                          geometry
#> 1  POLYGON ((123.8607 43.0607,...
#> 2  POLYGON ((131.9549 41.4295,...
#> 3  POLYGON ((142.9576 47.937, ...
#> 4  POLYGON ((115.9434 43.2923,...
#> 5  POLYGON ((114.4442 40.2853,...
#> 6  POLYGON ((116.192 41.8191, ...
#> 7  POLYGON ((84.8005 73.1575, ...
#> 8  POLYGON ((121.8271 71.3283,...
#> 9  POLYGON ((123.7043 66.1533,...
#> 10 POLYGON ((97.0729 73.9899, ...

You can plot the results in a similar way:

plot(pl100$geometry)

Reconstructing modern `coastlines`

Note that the result includes oceanic plates as well, which can be quite distracting when you want to visualize past distribution of the continents. To help with this, the paleogeographic positions of the present-day coastlines (and plate margins) can be quite helpful, which can be accessed, by providing "coastlines" as the first argument of the reconstruct() function:

c0 <- reconstruct("coastlines", age=0)
c100 <- reconstruct("coastlines", age=100)
c100
#> Simple feature collection with 2451 features and 0 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -179.99 ymin: -89.3928 xmax: 179.99 ymax: 83.1004
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>                          geometry
#> 1  POLYGON ((132.3511 47.0744,...
#> 2  POLYGON ((141.963 46.5718, ...
#> 3  POLYGON ((134.1119 44.0571,...
#> 4  POLYGON ((134.1819 42.4726,...
#> 5  POLYGON ((133.5996 44.4895,...
#> 6  POLYGON ((142.7117 48.1375,...
#> 7  POLYGON ((133.7077 44.7113,...
#> 8  POLYGON ((133.5384 44.3543,...
#> 9  POLYGON ((132.7962 42.5404,...
#> 10 POLYGON ((132.486 42.1309, ...

You can visalize these in the same way:

plot(c100$geometry)

Again, the pl0 and pl100 and c100 objects are sf-class objects. You can customize their plotting as you would for any other sf object - for instance setting a fill color for the polygons and not plotting their boundaries.

plot(pl100$geometry, col="gray", border=NA)

You can easily plot them on top of each other:

plot(pl100$geometry, col="gray", border=NA)
plot(c100$geometry, col="darkgray", border=NA, add=TRUE)

Individual locations

This is nice, but plotting the plates on their own has only so much use. The true use of reconstruct() is the ability to calculate the paleocoordinates of present-day locations for a given age.

Single present-day point

Let’s consider the location of London (here it is on Google maps)! Either in the URL, or with the user interface of Google maps, you can find that the coordinates of the city center are around 51.52°N (latitude) and 0.38°W (longitude). To figure out where the city was at the Triassic/Jurassic boundary you need to: (1) register these coordinates in R; and (2) provide them to the reconstruct() function.

Note that we are dealing with global-scale, approximate coordinates here. For more precise results, always make sure that the CRS of the points is matching the CRS of the maps (including ellipses)!

Because with the usually used equirectangular projection the x-axis of plotting become longitude, and the y axis becomes latitude, we usually register the coordinates in the order of longitude first, and latitude second (with easting longitudes and northing latitudes registered as positive values).

To make this structure absolutely clear, it is best to register coordinates as 2-column matrices, with longitude being the first, and latitude the second column:

# the coordinates
london <- c(-0.38, 51.52)

# make it a matrix
londonMat <- matrix(london, ncol=2, byrow=TRUE)

# add column names (optional)
colnames(londonMat) <- c("long", "lat")
londonMat
#>       long   lat
#> [1,] -0.38 51.52

Since the coordinate reference system (CRS) of the maps is longitude-latitude, you can use these coordinates directly to indicate the positionx of the city on the present-day map using points() - in this case with red plus signs:

plot(c0$geometry)
points(londonMat, col="red", pch=3)

Paleocoordinates of a single locality

Finding the paleocoordinates of such localities is as easy as calculating those of the plates. You have to use the matrix that you defined earlier (londonMat) as the first argument, of reconstruct() (where "static_polygons" and "coastlines" were given earlier), and provide a target age in million years:

londonMat100 <- reconstruct(londonMat, age=100)
londonMat100
#>      paleolong paleolat
#> [1,]    16.501  43.4247

The result of this calculation is a similar matrix: now including the paleolongitude and paleolatitude columns. If you provide coordinates as a plain matrix, coordinates are inferred to be longitude and latitude.

You can visualize this the same way, as you visualized the present-day position of the location:

plot(pl100$geometry, col="gray", border=NA)
points(londonMat100, col="red", pch=3)

Paleocoordinates of a multiple localities

If you want to reconstruct multiple locations, all you need to provide is more rows in the matrix. For instance, if you also want to calculate the positions of Sydney, AU (33.85°S, 151.11°E) and Montréal (CA) (45.52°N, 73.61°W), you need to add these in a similar fashion.

# coordinates of the two other cities
sydney<- c(151.17, -33.85)
montreal<- c(-73.61, 45.52)

# all cities in a single matrix
cities<- rbind(london, sydney, montreal)

#optional column names
colnames(cities) <- c("long", "lat")
cities
#>            long    lat
#> london    -0.38  51.52
#> sydney   151.17 -33.85
#> montreal -73.61  45.52

Now that we have a matrix of longitudes and latitudes, all we need to do is use this as the first argument of the reconstruct() function:

cities100 <- reconstruct(cities, age=100)
cities100
#>          paleolong paleolat
#> london     16.5010  43.4247
#> sydney    168.1989 -64.7057
#> montreal  -24.3468  43.0799

Note that the order of entities remains the same. We can plot these the same way similar to a single city.

# the background map
plot(c100$geometry, col="gray", border=NA)

# the reconstructed cities 
points(cities100, col="red", pch=3)

Other reconstruction models

If you look into the Details of the reference of reconstruct(), you will see that the model argument can be set to other character strings besides "MERDITH2021". These indicate other models that are accessible via the GPLates Web Service. For instance, if you want to execute the same calculations to reconstruct the position position of the plates at 100Ma, with the Seton et al. 2012 model, all you have to do is set model="SETON2012":

c100seton <- reconstruct("coastlines", age=100, model="SETON2012")
c100seton
#> Simple feature collection with 1893 features and 0 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -179.99 ymin: -89.99 xmax: 179.99 ymax: 85.6335
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>                          geometry
#> 1  POLYGON ((-10.7862 27.3367,...
#> 2  POLYGON ((-141.3289 -23.446...
#> 3  POLYGON ((-60.5945 5.9602, ...
#> 4  POLYGON ((-59.9848 4.8269, ...
#> 5  POLYGON ((-36.5397 16.789, ...
#> 6  POLYGON ((-36.7068 16.8855,...
#> 7  POLYGON ((-36.494 16.7124, ...
#> 8  POLYGON ((-36.3418 16.5887,...
#> 9  POLYGON ((-12.5495 22.9016,...
#> 10 POLYGON ((-10.4138 22.3064,...

You can compare this with the Merdith et al. model, by plotting this result on top of that with some transparency (e.g. a semi-transparent red, in HTML RGBA: "#FF000077"):

plot(c100$geometry, col="gray", border=NA)
plot(c100seton$geometry, col="#FF000077", add=TRUE, border=NA)

You can see that the the two reconstructions differ quite a bit, which becomes more apparent as we go back in time. Also, with the transparency, you can see how the plates overlap in convergent zones.