Reconstruct regular features

This example shows a couple of different scenarios involving the reconstruction of regular features to geological times.

Regular features exclude topological features, motion path features and flowline features.

Exported reconstructed features to a file 

In this example we reconstruct regular features and export the results to a Shapefile.

Sample code 

import pygplates


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

# Load some features.
features = pygplates.FeatureCollection('features.gpml')

# Reconstruct features to this geological time.
reconstruction_time = 50

# The filename of the exported reconstructed geometries.
# It's a shapefile called 'reconstructed_50Ma.shp'.
export_filename = 'reconstructed_{0}Ma.shp'.format(reconstruction_time)

# Reconstruct the features to the reconstruction time and export them to a shapefile.
pygplates.reconstruct(features, rotation_model, export_filename, reconstruction_time)

Details 

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

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

The reconstructable features are loaded into a pygplates.FeatureCollection.

features = pygplates.FeatureCollection('features.gpml')

The features will be reconstructed to their 50Ma positions.

reconstruction_time = 50

All features are reconstructed to 50Ma using pygplates.reconstruct().

We specify a filename to export the reconstructed geometries to.

pygplates.reconstruct(features, rotation_model, export_filename, reconstruction_time)

Output 

We should now have a file called reconstructed_50Ma.shp containing feature geometries reconstructed to their 50Ma positions.

Calculate distance a feature is reconstructed 

In this example we calculate the distance between a feature geometry’s present day (centroid) location and its reconstructed (centroid) location.

Sample code 

import pygplates


# A function to return the centroid of the geometry (point/multipoint/polyline/polygon).
def get_geometry_centroid(geometry):

    # See if geometry is a polygon.
    try:
        return geometry.get_interior_centroid()
    except AttributeError:
        # Not a polygon so keeping going.
        pass

    # See if geometry is a polyline or multipoint.
    try:
        return geometry.get_centroid()
    except AttributeError:
        # Not a polyline or multipoint so keeping going.
        pass

    # Geometry must be a point - it is already its own centroid.
    return geometry


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

# Load some features.
features = pygplates.FeatureCollection('features.gpml')

# Reconstruct features to this geological time.
reconstruction_time = 50

# Reconstruct the features to the reconstruction time.
reconstructed_feature_geometries = []
pygplates.reconstruct(features, rotation_model, reconstructed_feature_geometries, reconstruction_time)

# Iterate over all reconstructed feature geometries.
for reconstructed_feature_geometry in reconstructed_feature_geometries:

    # Calculate distance between:
    #  - the centroid of the present-day geometry, and
    #  - the centroid of the reconstructed geometry.
    distance_reconstructed = pygplates.GeometryOnSphere.distance(
        get_geometry_centroid(reconstructed_feature_geometry.get_present_day_geometry()),
        get_geometry_centroid(reconstructed_feature_geometry.get_reconstructed_geometry()))

    # Convert distance from radians to Kms.
    distance_reconstructed_in_kms = distance_reconstructed * pygplates.Earth.mean_radius_in_kms

    # Print the associated feature name and plate ID. And print the distance reconstructed.
    print 'Feature: %s' % reconstructed_feature_geometry.get_feature().get_name()
    print '  plate ID: %d' % reconstructed_feature_geometry.get_feature().get_reconstruction_plate_id()
    print '  distance reconstructed: %f kms' % distance_reconstructed_in_kms

Details 

We define a function to return the centroid of a geometry.
We don’t necessarily know whether the geometry is a pygplates.PointOnSphere,
pygplates.MultiPointOnSphere, pygplates.PolylineOnSphere or pygplates.PolygonOnSphere.
Each geometry type requires a different method for obtaining its centroid.
We use the standard Python approach of attempting to use a method and if it fails try something else.
So first we see if it’s a polygon and call pygplates.PolygonOnSphere.get_interior_centroid().
Then we see if it’s a polyline or multipoint - both of which have a get_centroid() method
(pygplates.PolylineOnSphere.get_centroid() and pygplates.MultiPointOnSphere.get_centroid()).
If they all fail then it must be a point geometry so we just return that as the centroid.

def get_geometry_centroid(geometry):
    try:
        return geometry.get_interior_centroid()
    except AttributeError:
        pass
    try:
        return geometry.get_centroid()
    except AttributeError:
        pass
    return geometry

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

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

The reconstructable features are loaded into a pygplates.FeatureCollection.

features = pygplates.FeatureCollection('features.gpml')

The features will be reconstructed to their 50Ma positions.

reconstruction_time = 50

All features are reconstructed to 50Ma using pygplates.reconstruct().

We specify a list for reconstructed_feature_geometries instead of a filename so that we can query the reconstructed geometries easily.

reconstructed_feature_geometries = []
pygplates.reconstruct(features, rotation_model, reconstructed_feature_geometries, reconstruction_time)

We use our get_geometry_centroid() function to find the centroid of the present day and reconstructed geometries.

We use the pygplates.GeometryOnSphere.distance() function to calculate the shortest distance between the two centroids and convert it to kilometres using pygplates.Earth.

distance_reconstructed = pygplates.GeometryOnSphere.distance(
    get_geometry_centroid(reconstructed_feature_geometry.get_present_day_geometry()),
    get_geometry_centroid(reconstructed_feature_geometry.get_reconstructed_geometry()))
distance_reconstructed_in_kms = distance_reconstructed * pygplates.Earth.mean_radius_in_kms

Output 

Feature: Pacific
  plate ID: 982
  distance reconstructed: 3815.013838 kms
Feature: Marie Byrd Land
  plate ID: 804
  distance reconstructed: 514.440695 kms
Feature: Pacific
  plate ID: 901
  distance reconstructed: 3795.781009 kms
Feature: Pacific
  plate ID: 901
  distance reconstructed: 3786.206123 kms
Feature: Pacific
  plate ID: 901
  distance reconstructed: 3786.068477 kms
Feature: Pacific
  plate ID: 901
  distance reconstructed: 3785.868706 kms
Feature: Pacific
  plate ID: 901
  distance reconstructed: 3785.465344 kms
Feature: Pacific
  plate ID: 901
  distance reconstructed: 3788.422368 kms
Feature: Pacific
  plate ID: 901
  distance reconstructed: 3790.540180 kms
Feature: Pacific
  plate ID: 901
  distance reconstructed: 3554.951168 kms
Feature: Pacific
  plate ID: 901
  distance reconstructed: 3553.133934 kms
Feature: Northwest Africa
  plate ID: 714
  distance reconstructed: 643.521413 kms

...