Import geometries and assign plate IDs

This example shows how to import points or polylines from a text file and assign plate IDs to them. The resulting feature data is then saved to a file format that pyGPlates can load directly.

Note

Importing features is different than loading features. When feature data is not in a file format that pyGPlates can load then it needs to be imported into a GPlates-compatible file format. This usually involves more than just importing the geometry data. For example, plate IDs need to be assigned so that pyGPlates can reconstruct the feature data to geological times. And other feature metadata such as name and description should also be assigned.

Import points from a text file and assign plate IDs 

This example:

reads some point locations from a text file,
assigns a plate ID and valid time period to each point, and
writes the assigned points to a GPML file ready for use in GPlates.

Sample code 

import pygplates


# Load one or more rotation files into a rotation model.
rotation_model = pygplates.RotationModel('rotations.rot')

# The static polygons file will be used to assign plate IDs and valid time periods.
static_polygons_filename = 'static_polygons.gpml'

# Filename of the input points file that we will read.
input_points_filename = 'input_points.txt'

# Filename of the output points file that we will write.
# Note that it is a GPML file (has extension '.gpml').
# This enables it to be read into GPlates.
output_points_filename = 'output_points.gpml'

# Parse the input points text file containing a lon/lat point per line.
# Longitude/latitude is the order that GMT files ('.xy', '.gmt') use.
# Note that pyGPlates specifies points in the opposite (latitude/longitude) order.
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

        # Split the line into strings (separated by whitespace).
        line_string_list = line.split()

        # Need at least two strings per line (for latitude and longitude).
        if len(line_string_list) < 2:
            print 'Line %d: Ignoring point - line does not have at least two white-space separated strings.' % line_number
            continue

        # Attempt to convert each string into a floating-point number.
        try:
            # Use GMT (lon/lat) order.
            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.
        # Note that pyGPlates uses the opposite (lat/lon) order to GMT (lon/lat).
        input_points.append(pygplates.PointOnSphere(lat, lon))

# Create a feature for each point we read from the input file.
point_features = []
for point in input_points:

    # Create an unclassified feature.
    point_feature = pygplates.Feature()

    # Set the feature's geometry to the input point read from the text file.
    point_feature.set_geometry(point)

    point_features.append(point_feature)

# Use the static polygons to assign plate IDs and valid time periods.
# Each point feature is partitioned into one of the static polygons and assigned its
# reconstruction plate ID and valid time period.
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)

Input 

An example input text file (in longitude/latitude order) looks like:

-79.747867      -1.444159
-79.786712      -1.654002
-79.872547      -2.221801
-79.858122      -6.951201
-78.850008      -9.359851
-76.020448      -13.798207
-75.549659      -14.315297
-75.411320      -14.456342
-74.335501      -15.543422
-72.539796      -17.187214
-71.922547      -17.773935
-71.381735      -18.373316
-70.979182      -18.850190
-70.786266      -19.126329
-70.571175      -19.417365
-70.343507      -19.716224
-70.280285      -19.858811
-70.107565      -20.531859
-70.059697      -22.248895

Details 

The rotations are loaded from a rotation file into a pygplates.RotationModel.

rotation_model = pygplates.RotationModel('rotations.rot')

The points will be read from an input text file and written to an output GPML file.

input_points_filename = 'input_points.txt'
output_points_filename = 'output_points.gpml'

The input points file is opened and read line-by-line.

with open(input_points_filename, 'r') as input_points_file:
    for line_number, line in enumerate(input_points_file):

Each line contains a latitude string and a longitude string.
We attempt to convert them to floating-point numbers.
If that fails then we catch the ValueError exception that Python raises and ignore that line in the file.

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

An unclassified feature is created for each point we read from the input file.
Leaving the feature type empty in pygplates.Feature()
defaults to a feature type of pygplates.FeatureType.gpml_unclassified_feature.
Ideally we should pick a specific feature type such as
pygplates.FeatureType.gpml_hot_spot,
perhaps reading it from the input file (as an extra column).
And ideally we should also import extra metadata such as feature name
and description.

point_feature = pygplates.Feature()

The point geometry is set on the point feature using pygplates.Feature.set_geometry().

If we don’t do this then the feature cannot be used in spatial calculations and will not display on the globe in GPlates.

point_feature.set_geometry(point)

Each point feature is partitioned into one of the static polygons and assigned its
reconstruction plate ID and valid time period using pygplates.partition_into_plates().
The static polygons have global coverage at present day (the default reconstruction time
for pygplates.partition_into_plates()) and should therefore partition all the input points.
We also explicitly specify the argument properties_to_copy to assign both the reconstruction plate
ID and valid time period (the default is just to assign the reconstruction plate ID).

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])

Finally we put the list of assigned features into a pygplates.FeatureCollection so that we can write them out to a file using pygplates.FeatureCollection.write().

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

Load points from a GMT file and assign plate IDs 

This example is similiar to Import points from a text file and assign plate IDs except it takes advantage of the ability of pyGPlates to load a GMT file to avoid having to manually parse a text file line-by-line.

Sample code 

# Load a GMT file (instead of manually reading a '.txt' file line-by-line).
point_features = pygplates.FeatureCollection('input_points.gmt')

# Use the static polygons to assign plate IDs and valid time periods.
# Each point feature is partitioned into one of the static polygons and assigned its
# reconstruction plate ID and valid time period.
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)

Input 

An example input text file (in longitude/latitude order) looks like:

-79.747867      -1.444159
-79.786712      -1.654002
-79.872547      -2.221801
-79.858122      -6.951201
-78.850008      -9.359851
-76.020448      -13.798207
-75.549659      -14.315297
-75.411320      -14.456342
-74.335501      -15.543422
-72.539796      -17.187214
-71.922547      -17.773935
-71.381735      -18.373316
-70.979182      -18.850190
-70.786266      -19.126329
-70.571175      -19.417365
-70.343507      -19.716224
-70.280285      -19.858811
-70.107565      -20.531859
-70.059697      -22.248895

Details 

Since GPlates can directly load GMT '.gmt' files, an alternative is to
change the filename extension, of your text file, to '.gmt'. The feature metadata will be missing
from your text file so only the geometry data will get loaded, but this achieves the same effect
as the above example.
As with the previous example, there can be more than two numbers per line, but only the first two
are used (as longitude and latitude) - note that if you were to load your '.gmt' file into GPlates
the extra data would cause it to give a warning about flattening 2.5D to 2D.
Note that, as with the previous example, the data should be in GMT (longitude/latitude) order.

point_features = pygplates.FeatureCollection('input_points.gmt')

The rest of the sample code is similar to Import points from a text file and assign plate IDs.

Import polylines from a text file and assign plate IDs 

This example is similiar to Import points from a text file and assign plate IDs except it imports polylines instead of points.

Sample code 

import pygplates


# A function that create a polyline feature from some points and adds to a list of features.
def add_polyline_feature_from_points(polyline_features, points, line_number):

    # If have no points then nothing to do.
    if not points:
        return

    # Need at least two points for a polyline.
    if len(points) >= 2:
        polyline = pygplates.PolylineOnSphere(points)

        polyline_feature = pygplates.Feature() # 'unclassified' feature
        polyline_feature.set_geometry(polyline)

        polyline_features.append(polyline_feature)

    # If only one point then emit warning.
    else:
        print 'Line %d: Ignoring polyline - polyline has only one point.' % (line_number-1)

    # Clear points list for next feature.
    # Clear the list in-place so that all references to the list see an empty list.
    del points[:]


# Load one or more rotation files into a rotation model.
rotation_model = pygplates.RotationModel('rotations.rot')

# The static polygons file will be used to assign plate IDs and valid time periods.
static_polygons_filename = 'static_polygons.gpml'

# Filename of the input polylines file that we will read.
input_polylines_filename = 'input_polylines.txt'

# Filename of the output polylines file that we will write.
# Note that it is a GPML file (has extension '.gpml').
# This enables it to be read into GPlates.
output_polylines_filename = 'output_polylines.gpml'

# Parse the input polylines text file containing groups of lon/lat points per line.
# Longitude/latitude is the order that GMT files ('.xy', '.gmt') use.
# Note that pyGPlates specifies points in the opposite (latitude/longitude) order.
polyline_features = []
polyline_points = []
with open(input_polylines_filename, 'r') as input_polylines_file:
    for line_number, line in enumerate(input_polylines_file):

        # Make line number 1-based instead of 0-based.
        line_number = line_number + 1

        # See if line begins with '>'.
        # This is was separates groups of points into polylines.
        if line.strip().startswith('>'):

            # Generate the previous polyline feature if we have two or more points.
            add_polyline_feature_from_points(polyline_features, polyline_points, line_number)

            # Skip to next line.
            continue

        # Split the line into strings (separated by whitespace).
        line_string_list = line.split()

        # Need at least two strings per line (for latitude and longitude).
        if len(line_string_list) < 2:
            print 'Line %d: Ignoring point - line does not have at least two white-space separated strings.' % line_number
            continue

        # Attempt to convert each string into a floating-point number.
        try:
            # Use GMT (lon/lat) order.
            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.
        # Note that pyGPlates uses the opposite (lat/lon) order to GMT (lon/lat).
        polyline_points.append(pygplates.PointOnSphere(lat, lon))

    # If we have any points leftover then generate the last polyline feature.
    # This happens if last line does not start with '>'.
    add_polyline_feature_from_points(polyline_features, polyline_points, line_number)

# Use the static polygons to assign plate IDs and valid time periods.
# Each polyline feature is partitioned into one or more of the static polygons and assigned their
# reconstruction plate IDs and valid time periods.
assigned_polyline_features = pygplates.partition_into_plates(
    static_polygons_filename,
    rotation_model,
    polyline_features,
    properties_to_copy = [
        pygplates.PartitionProperty.reconstruction_plate_id,
        pygplates.PartitionProperty.valid_time_period])

# Write the assigned polyline features to the output GPML file (ready for use in GPlates).
assigned_polyline_feature_collection = pygplates.FeatureCollection(assigned_polyline_features)
assigned_polyline_feature_collection.write(output_polylines_filename)

Input 

An example input text file (in longitude/latitude order) containing three polylines looks like:

>
  -79.747867      -1.444159
  -79.786712      -1.654002
  -79.872547      -2.221801
  -79.858122      -6.951201
  -78.850008      -9.359851
  -76.020448      -13.798207
>
  -75.549659      -14.315297
  -75.411320      -14.456342
  -74.335501      -15.543422
  -72.539796      -17.187214
  -71.922547      -17.773935
  -71.381735      -18.373316
  -70.979182      -18.850190
  -70.786266      -19.126329
>
  -70.571175      -19.417365
  -70.343507      -19.716224
  -70.280285      -19.858811
  -70.107565      -20.531859
  -70.059697      -22.248895

Note

The > symbol is used to group points into polylines.

Details 

The rotations are loaded from a rotation file into a pygplates.RotationModel.

rotation_model = pygplates.RotationModel('rotations.rot')

The polylines will be read from an input text file and written to an output GPML file.

input_polylines_filename = 'input_polylines.txt'
output_polylines_filename = 'output_polylines.gpml'

The input polylines file is opened and read line-by-line.

with open(input_polylines_filename, 'r') as input_polylines_file:
    for line_number, line in enumerate(input_polylines_file):

If a line begins with a '>' character then it separates those points in lines before it into one polyline and those points in lines after it into another polyline.

Here the points in prior lines are used to create a new polyline feature.

if line.strip().startswith('>'):
    add_polyline_feature_from_points(polyline_features, polyline_points, line_number)
    continue

Each line contains a latitude string and a longitude string.
We attempt to convert them to floating-point numbers.
If that fails then we catch the ValueError exception that Python raises and ignore that line in the file.

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

Keep track of the points for the current polyline so we can create the polyline once we’ve reached the last point (of the current polyline).

polyline_points.append(pygplates.PointOnSphere(lat, lon))

Then function add_polyline_feature_from_points creates a polyline feature from a list of points and adds it to a list of polyline features.

def add_polyline_feature_from_points(polyline_features, points, line_number):
    ...

If there are at least two points (required for a polyline) then a pygplates.PolylineOnSphere geometry is created from the points.

if len(points) >= 2:
    polyline = pygplates.PolylineOnSphere(points)

An unclassified feature is created for each polyline we read from the input file.
Leaving the feature type empty in pygplates.Feature()
defaults to a feature type of pygplates.FeatureType.gpml_unclassified_feature.
Ideally we should pick a specific feature type such as
pygplates.FeatureType.gpml_subduction_zone.
And ideally we should also import extra metadata such as feature name
and description.

polyline_feature = pygplates.Feature()

The polyline geometry is set on the polyline feature using pygplates.Feature.set_geometry().

If we don’t do this then the feature cannot be used in spatial calculations and will not display on the globe in GPlates.

polyline_feature.set_geometry(polyline)

Each polyline feature is partitioned into one or more of the static polygons and assigned their
reconstruction plate IDs and valid time periods using pygplates.partition_into_plates().
Note that the default argument for the partition_method parameter is pygplates.PartitionMethod.split_into_plates
which results in two polylines being generated (returned) when one polyline intersects two static polygons (plates).
The static polygons have global coverage at present day (the default reconstruction time
for pygplates.partition_into_plates()) and should therefore partition all the input polylines.
We also explicitly specify the argument properties_to_copy to assign both the reconstruction plate
ID and valid time period (the default is just to assign the reconstruction plate ID).

assigned_polyline_features = pygplates.partition_into_plates(
    static_polygons_filename,
    rotation_model,
    polyline_features,
    properties_to_copy = [
        pygplates.PartitionProperty.reconstruction_plate_id,
        pygplates.PartitionProperty.valid_time_period])

Finally we put the list of assigned features into a pygplates.FeatureCollection so that we can write them out to a file using pygplates.FeatureCollection.write().

assigned_polyline_feature_collection = pygplates.FeatureCollection(assigned_polyline_features)
assigned_polyline_feature_collection.write(output_polylines_filename)

Load polylines from a GMT file and assign plate IDs 

This example is similiar to Import polylines from a text file and assign plate IDs except it takes advantage of the ability of pyGPlates to load a GMT file to avoid having to manually parse a text file line-by-line.

Sample code 

# Load a GMT file (instead of manually reading a '.txt' file line-by-line).
polyline_features = pygplates.FeatureCollection('input_polylines.gmt')

# Use the static polygons to assign plate IDs and valid time periods.
# Each polyline feature is partitioned into one or more of the static polygons and assigned their
# reconstruction plate IDs and valid time periods.
assigned_polyline_features = pygplates.partition_into_plates(
    static_polygons_filename,
    rotation_model,
    polyline_features,
    properties_to_copy = [
        pygplates.PartitionProperty.reconstruction_plate_id,
        pygplates.PartitionProperty.valid_time_period])

# Write the assigned polyline features to the output GPML file (ready for use in GPlates).
assigned_polyline_feature_collection = pygplates.FeatureCollection(assigned_polyline_features)
assigned_polyline_feature_collection.write(output_polylines_filename)

Input 

An example input text file (in longitude/latitude order) containing three polylines looks like:

>
  -79.747867      -1.444159
  -79.786712      -1.654002
  -79.872547      -2.221801
  -79.858122      -6.951201
  -78.850008      -9.359851
  -76.020448      -13.798207
>
  -75.549659      -14.315297
  -75.411320      -14.456342
  -74.335501      -15.543422
  -72.539796      -17.187214
  -71.922547      -17.773935
  -71.381735      -18.373316
  -70.979182      -18.850190
  -70.786266      -19.126329
>
  -70.571175      -19.417365
  -70.343507      -19.716224
  -70.280285      -19.858811
  -70.107565      -20.531859
  -70.059697      -22.248895

Note

The > symbol is used by GMT to group points into polylines.

Details 

Since GPlates can directly load GMT '.gmt' files, an alternative is to
change the filename extension, of your text file, to '.gmt'. The feature metadata will be missing
from your text file so only the geometry data will get loaded, but this achieves the same effect
as the above example.
As with the previous example, there can be more than two numbers per line, but only the first two
are used (as longitude and latitude) - note that if you were to load your '.gmt' file into GPlates
the extra data would cause it to give a warning about flattening 2.5D to 2D.
Note that, as with the previous example, the data should be in GMT (longitude/latitude) order.

polyline_features = pygplates.FeatureCollection('input_polylines.gmt')

The rest of the sample code is similar to Import polylines from a text file and assign plate IDs.