pyGPlates: restore points to paleomap

Q&A
1

Hi, I am new to pyGPlates and have read the online document for a while… I have a list of lat and long points and want to restore the positions back to a certain geological time (e.g., 100 Ma). I know I need to use : pygplates.reconstruct(features, rotation_model, export_filename, reconstruction_time)

But how to create a reconstructable feature from lat and long (csv or xlsx data, point data) at first place? I see the examples online are mostly for coastline? Any existing workflow I can follow?

Thanks much!

2

Hi Wenrong,

There’s a script in PlateTectonicTools called convert_xy_to_gplates.py that can create pygplates features from text files containing rows longitude and latitude (or vice versa).

By default it creates point features but you can create polyline or polygon features (using lines beginning with >), and you can optionally import metadata (eg, a line containing > FeatureType = Coastline defines a Coastline feature).

Once you’ve installed PlateTectonicTools you can run python -m ptt.convert_xy_to_gplates --help to see the overview and options of that specific script.

The above tool is basically an extension of this pyGPlates sample code to import from a text file (note that the sample code also shows how to assign plate IDs to the imported geometry so it can be reconstructed).

Regards,
John

4

Thanks John! Following the example import geometries, I am able to restore the present-day points to their paleo-positions and generate associated gpmlz and shape files.

Another native question is how to append properties associated with points and restore them back together with the lat and long to a geological time? For example, I have SiO2 wt% for some plutons having lat and long and I want to see how the SiO2 was spatially distributed at 100 Ma? The examples on pyGPlates only teach how to restore lat and long …

The ptt.convert_xy_to_gplates is probably an overkill for my stage… anyway thanks for the info !

5

Great! That was quick.

You can add SiO2 as an extra column in your text file (after lon and lat) and then set that as a property on each point feature that you’re creating.

Have a look at this thread for a discussion on how to do this. Since it’s not a standard GPlates property there are some issues that that thread talks about (like being able to load these imported features into GPlates). But essentially you could do something like:

feature.add(
    pygplates.PropertyName.create_gpml('SiO2'),
    pygplates.XsDouble(SiO2_value),
    pygplates.VerifyInformationModel.no)

…and that would be fine if you stay in pyGPlates land where you could later get back the value from the reconstructed feature (after having reconstructed with reconstruct using the group_with_feature option) with something like:

SiO2_value = reconstructed_feature.get_value(
    pygplates.PropertyName.create_gpml('SiO2')).get_double()

However, as mentioned in the thread, loading into GPlates would be a problem, in which case the thread shows how to use a Shapefile attribute (eg, feature.set_shapefile_attribute('SiO2', SiO2_value)).

6

Hi John, since I will deal with reconstructed points in QGIS, I used the set_shapefile_attribute function. After trail and errors, I think it works.

I notice if I do a regular point reconstruction (without adding an additional shapefile attribute), the resulted reconstructed shapefile has ~15 attributes including Anchor, Time, Reconfile1, PlateID…Speard_Asy, Import_Age, etc.

After I used set_shapefile_attribute, I only have 4 attributes: Anchor, Time, Reconfile1 and SiO2. Missing other attributes do not bother me for now since SiO2 are tied with points… I may put the set_shapefile_attribute in a wrong place…

The input text file (input_points_new.txt) has three columns: lat, long and SiO2

Here is my code:

import pygplates
rotation_model = pygplates.RotationModel(‘rotations.rot’)
static_polygons_filename = ‘static_polygons.gpmlz’
input_points_filename = ‘input_points_new.txt’
output_points_filename = ‘output_points.gpmlz’

Parse the input points text file containing a lon/lat point per line.

input_points = []
with open(input_points_filename, ‘r’) as input_points_file:
for line_number, line in enumerate(input_points_file):
# Make line number 1-based instead of 0-based.
line_number = line_number + 1
line_string_list = line.split()
# Attempt to convert each string into a floating-point number.
try:
lon = float(line_string_list[0])
lat = float(line_string_list[1])
except ValueError:
print (‘Line %d: Ignoring point - cannot read lon/lat values.’ % line_number)
continue
# Create a pyGPlates point from the latitude and longitude, and add it to our list of points.
input_points.append(pygplates.PointOnSphere(lat, lon))

get SiO2 values from txt as a list

SiO2_value_list = []
file = open(‘input_points_new.txt’, ‘r’)
lines = file.readlines()
for line in lines:
a = line.split()
x = a[2]
SiO2_value_list.append(x)
file.close()

convert SiO2 values from list to float

SiO2_value_floats = []
for item in SiO2_value_list:
SiO2_value_floats.append(float(item))

Create a feature for each point we read from the input file.

i=0 # counter for SiO2 value float
point_features = []
for point in input_points:
point_feature = pygplates.Feature()

# add additional features to shapefile
point_feature.set_shapefile_attribute('SiO2',SiO2_value_floats[i])

point_feature.set_geometry(point)

point_features.append(point_feature)

i=i+1

Use the static polygons to assign plate IDs and valid time periods.

assigned_point_features = pygplates.partition_into_plates(
static_polygons_filename,
rotation_model,
point_features,
properties_to_copy = [
pygplates.PartitionProperty.reconstruction_plate_id,
pygplates.PartitionProperty.valid_time_period])

Write the assigned point features to the output GPML file (ready for use in GPlates).

assigned_point_feature_collection = pygplates.FeatureCollection(assigned_point_features)
assigned_point_feature_collection.write(output_points_filename)

#%% Restore points to a geological time based on rotation model

rotation_model = pygplates.RotationModel(‘rotations.rot’)
features = pygplates.FeatureCollection(‘output_points.gpmlz’)
reconstruction_time = 100
export_filename = ‘reconstructed_points_{0}Ma.shp’.format(reconstruction_time)
pygplates.reconstruct(features, rotation_model, export_filename, reconstruction_time)

7

Great!

What you are doing is correct. The difference in shapefile attributes is confusing and is mostly due to the fact that pyGPlates is re-using some code in the GPlates export pipeline.

To explain the difference…

Firstly, the common behaviour is to set ANCHOR, TIME and RECONFILE1 (to the anchor plate, reconstruction time and rotation file). This is why they are seen in both cases (ie, when you reconstruct with and without shapefile attributes).

Secondly, the different behaviour depends on whether the features have shapefile attributes or not. When they do have attributes (eg, after calling set_shapefile_attribute()) then only those attributes are written. When they don’t have attributes then the standard shapefile attributes (that GPlates uses to map shapefile attributes to GPlates properties, eg, PLATEID1, FROMAGE, TOAGE, etc) are written. The assumption here being that the standard shapefile attributes are always present - this is usually true when you load a Shapefile in GPlates but as you are finding out this is not necessary true for pyGPlates (especially when creating your own features rather than loading them from a Shapefile).

This whole area concerning Shapefile mapping is being overhauled for reasons like those addressed in this thread. But what you’re doing is fine and you only need the SiO2 attribute so this doesn’t affect that.