Using pre-assigned plate IDs to rotate coordinates

Q&A
1

Hello!

I’ve been trying to follow the trajectory of palaeo-coordinates (i.e. Scotese’s PaleoDEM coordinates) through time, e.g. from 540 Ma to 535, with Gplately/pygplates. This works without a hitch when Gplately/pygplates assigns the plate IDs, but something goes wrong when using pre-assigned plate IDs (edit: they end up in the wrong location as illustrated in the code/figures). I would like to use pre-assigned IDs as some coordinates do not fall on the plate polygons. If I let gplately. Points assign an ID, the coordinates that fall outside of the plates are either removed or incorrectly rotated.

I’ve managed to work out a partial fix that works if all coordinates fall on the plate polygons but ignores all coordinates that do not. This was done by manipulating the ‘features’ part of a gplately in the same way as if there would not have been a plate ID. The ‘features’ is probably where the problem is, and it might have something to do with coordinates incorrectly being specified to t = 0 Ma instead of t = 540 Ma. However, those are all big questions marks, and I’ve hit a wall.

Any help would be greatly appreciated!

I’ve deposited the code in a Dropbox folder along with the dataset with only continental coordinates, and a csv with all paleoDEM coordinates pre-assigned to a plate ID ( https://www.dropbox.com/scl/fo/e31nc0lx0uam58hggrtdr/h?rlkey=ysxuxm1owkns874u8wm87lx7i&dl=0 )

Cheers

Bouwe

### R code to run to start the interactive python session in Rstudio
# library(reticulate)
# library(tidyverse)
# use_condaenv("pygplates_py310")
# reticulate::repl_python()

# versions:
# python version 3.10.12
# pygplates 0.39
# gplately 1.0

# import packages
import pandas as pd
import numpy as np
import gplately
from gplately import pygplates
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

# load the plate model
gdownload = gplately.download.DataServer("Scotese2016")
rotation_model, topology_features, static_polygons = gdownload.get_plate_reconstruction_files()
coastlines, continents, COBs = gdownload.get_topology_geometries()

model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)
gplot = gplately.PlotTopologies(model, coastlines=coastlines, continents=continents, COBs=COBs)

paleodem_A = pd.read_csv('Continental_coords_pID_assigned.csv')

# first use the pre-assigned plate IDs (I calculated these by with gplately.Points previously and saved the output)
gp = gplately.Points(plate_reconstruction=model, lons=paleodem_A['Longitude'], lats=paleodem_A['Latitude'], time=540, plate_id=paleodem_A['Plate_ID'])
# calculate lons and lats
rlons, rlats = gp.reconstruct(535, return_array=True)

# Set up a GeoAxis plot
fig = plt.figure(figsize=(32,24), dpi=600)
ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude = 0))
ax.set_global()
ax.gridlines(color='0.7',linestyle='--', xlocs=np.arange(-180,180,15), ylocs=np.arange(-90,90,15))
ax.set_title("Projected - plate ID specified")

# Plot shapefile features, subduction zones and MOR boundaries at 0 Ma
gplot.time = 540 # Ma
gplot.plot_coastlines(ax, color='DarkGrey')
sc = ax.scatter(paleodem_A['Longitude'], paleodem_A['Latitude'], color='orange', 
                transform=ccrs.PlateCarree(), label='A', s = 10)
sc = ax.scatter(rlons, rlats, color='blue',
                transform=ccrs.PlateCarree(), label='Projected', s= 10)
ax.legend(frameon=False)
plt.show()

# now calculate the plate_ID with gplately.points
gp2 = gplately.Points(plate_reconstruction=model, lons=paleodem_A['Longitude'], lats=paleodem_A['Latitude'], time=540)
# all plate IDs are identical!
all(gp.plate_id == gp2.plate_id)
# all coords are identical!
all(gp.lats == gp2.lats)
all(gp.lons == gp2.lons)

# calculate the lons and lats by projecting into the future
rlons2, rlats2 = gp2.reconstruct(535, return_array=True)

# Set up a GeoAxis plot
fig = plt.figure(figsize=(32,24), dpi=600)
ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude = 0))
ax.set_global()
ax.gridlines(color='0.7',linestyle='--', xlocs=np.arange(-180,180,15), ylocs=np.arange(-90,90,15))
ax.set_title("Projected - plate ID not specified")

# Plot shapefile features, subduction zones and MOR boundaries at 0 Ma
gplot.time = 540 # Ma
gplot.plot_coastlines(ax, color='DarkGrey')
sc = ax.scatter(paleodem_A['Longitude'], paleodem_A['Latitude'], color='orange', 
                transform=ccrs.PlateCarree(), label='A', s = 10)
sc = ax.scatter(rlons2, rlats2, color='blue',
                transform=ccrs.PlateCarree(), label='Projected', s= 10)
ax.legend(frameon=False)
plt.show()

# Now it does work....

# I looked at the code and I think the error is in the creation of the features, so recalculate as if plate_ID was absent

# create a 'features' object -> this is also what is created in the gplately.Points function when plate_ID is specified
features = gplately.tools.points_to_features(gp.lons, gp.lats, gp.plate_id)
# run the pygplates.partition_into_plates function to get the 'proper' features -> normally run in gplately.Points when plate_ID == None
partitioned_features = pygplates.partition_into_plates(
  static_polygons, rotation_model, features, reconstruction_time=540, sort_partitioning_plates=None
)
# replace the original features
gp.features = partitioned_features
# calculate new coordinates
rlons3, rlats3 = gp.reconstruct(535, return_array=True)



# Set up a GeoAxis plot
fig = plt.figure(figsize=(32,24), dpi=600)
ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude = 0))
ax.set_global()
ax.gridlines(color='0.7',linestyle='--', xlocs=np.arange(-180,180,15), ylocs=np.arange(-90,90,15))
ax.set_title("Projected - features replaced")

# Plot shapefile features, subduction zones and MOR boundaries at 0 Ma
gplot.time = 540 # Ma
gplot.plot_coastlines(ax, color='DarkGrey')
sc = ax.scatter(paleodem_A['Longitude'], paleodem_A['Latitude'], color='orange', 
                transform=ccrs.PlateCarree(), label='A', s = 10)
sc = ax.scatter(rlons3, rlats3, color='blue',
                transform=ccrs.PlateCarree(), label='Projected', s= 10)
ax.legend(frameon=False)
plt.show()

# and the calculated coordinates are now correct

# however this is problematic if I want to work with pre-assigned plate_IDs, especially if those do not fall on a continental shelve (as those will be changed to 0 with this fix)
2

Yes I think that’s the issue.

When the plate IDs are specified then gplately assumes the initial points are in their present day location (although I might be wrong - see next paragraph). According to the documention (for gplately.Points) the time is what the initial points are reconstructed to (ie, they are not first reverse reconstructed to present day) because it says “The specific geological time (Ma) at which to reconstruct the point data”. However when the plate IDs are not specified then the initial points are assumed to already be in their position at time (not present day). You can see this in the source code where pygplates.partition_into_plates() is called - it will first reconstruct the static polygons to time and then assign plate IDs to the points (which are assumed to already have been reconstructed to time, ie, are already in their position at time) and then reverse reconstruct the points (using the newly assigned plate IDs) to present day (because all features should store present day geometry) - this is documented here.

And when gplately.Points.reconstruct() is called it looks like the initial points are reconstructed from the initial time to a specified reconstruction time, which then calls gplately.PlateReconstruction.reconstruct() that expects the from_time to be zero otherwise it documents throwing NotImplementedError but does not appear to. @mchin and @ben.mather, there might be an issue here? And if non-zero from_time is not yet going to be implemented then perhaps gplately.Points could document that a non-zero time should only be used when plate IDs are not specified, or something like that, not sure what’s the best approach.