# 2. Reconstructions using modern coordinates (online)

Source:`vignettes/2_online_reconstruction_backwards.Rmd`

`2_online_reconstruction_backwards.Rmd`

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:

#### 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:

#### 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.