Reconstruction differs between GUI and pygplates

Q&A
1

Hi all,

I am having this weird issue when trying to reconstruct a coastline back in time, at 300 Ma, using the Merdith 2021 dataset (i.e., using 1000_0_rotfile_Merdith_et_al).

When looking at a series of coordinates comprising a coastline at present using the GUI:
27.5034 ; 66.6876
25.8708 ; 66.4792
25.466 ; 66.7494
25.2607 ; 66.6729

the reconstructed coordinates (again, using the GUI) I get are:
-35.6073 ; 59.9325
-36.4288 ; 58.1729
-36.8636 ; 57.9425
-36.9276 ; 57.6838

when trying the very same reconstruction using pygplates + gplately, I’m getting:
-35.60729953 ; 59.93246794
-48.43720088 ; 51.61909315
-48.86747549 ; 51.32369867
-48.92629734 ; 51.00663505

Note how the first coordinate matches the reconstruction with GUI, but the remaining 3 coordinates do not.

Here’s the Python 3 code to reproduce the issue:

import pygplates
import gplately

gdownload = gplately.download.DataServer("Merdith2021")
rotation_model, topology_features, static_polygons = gdownload.get_plate_reconstruction_files()
model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)
coastlines, continents, _ = gdownload.get_topology_geometries()
gplot = gplately.PlotTopologies(
            model, coastlines=coastlines, continents=continents)
gplot.time = 0

gpts1 = gplately.Points(model, [66.68764182400014, 66.47917004564209, 66.74944466750146, 66.67292015518123], [27.50343290900005, 25.870788285673104, 25.46597337481638, 25.260725144855893])
rec_lngs1, rec_lats1 = gpts1.reconstruct(time=300.0, return_array=True)

print("look, below, the difference between reconstruction in Python vs. using the GUI")
print(rec_lats1)
# first matches, last three don't
# -35.60729953, -48.43720088, -48.86747549, -48.92629734

print(rec_lngs1)
# first matches, last three don't
# 59.93246794, 51.61909315, 51.32369867, 51.00663505

I am pretty new to GPlates and this could totally be some dumb oversight on my end. (Let’s hope it’s actually that!)

Thanks so much for any help you can provide!
cheers,
Fabio

2

Answering my own question, I think the problem is with gplately perhaps.

# ... code for initializing coastline_geometry_obj, rotation_model, etc ...
reconstructed_feature_geometries = list()
pygplates.reconstruct(coastline_geometry_obj, rotation_model, reconstructed_feature_geometries, 300)
rec_feat = reconstructed_feature_geometries[0].get_reconstructed_geometry()
rec_geo_cl_points = [p.to_lat_lon() for p in rec_feat]

print("See below the right reconstructions, which suggests the issue is with gplately?")
print(rec_geo_cl_points)

# [(-35.607299527307546, 59.93246793546048), (-36.42878142628977, 58.17287923711001), (-36.863566268945554, 57.942455584993006), (-36.927622155807114, 57.68375084121882)]

Anyway, if anyone has seen this issue before or has any thoughts, please do share with me.

Thanks so much!
cheers,
Fabio

3

Hi Fabio,

It might be that gplately is assigning plate IDs to your points (since your gplately.Points call does not specify the plate_id parameter). And in your second example (using pygplates) your coastline_geometry_obj could already have plate IDs (but they might be different than what gplately is assigning using the Merdith2021 static polygons).

Can you show how you created coastline_geometry_obj?

4

Hi John,

You turned out to be exactly right.

I did not think that specifying the ID was strictly necessary because gplately automatically figures out which plate ID to use in most – but not all! – cases. There are even tutorials where the authors don’t specify plate IDs. It works fine if a point sits smack in the middle of a plate.

If a point sits right at the boundary of two plates, gplately anchors the point on one or the other more or less randomly. I could not find any regularity in how gplately does that.

In any case, problem solved!

Thanks a lot!
Fabio

5

Hi Fabio,

Yes it’s nice that gplately automatically assigns the plate ID like that.

Underneath it’s using pygplates via its partition_into_plates function. And by default that function will sort by plate ID (because SortPartitioningPlates.by_partition_type_then_plate_id is the default sorting flag). So I think if a point sits right at the boundary (and hits both plates) then the highest plate ID is used.

I say “hits both plates” because, depending on how the static polygons were constructed, there could be a tiny sliver of a gap between two adjacent polygons - so it’s possible that a point could miss both plates, hit one or the other plate, or hit both plates (based on numerical tolerance - there’s a tiny expanded tolerance used in the internal point-in-polygon test).