Create conjugate isochrons from a ridge

This example creates a conjugate pair of isochrons from a mid-ocean ridge at each specified geological time in a series.

Sample code

import pygplates


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

# Load the mid-ocean ridge features.
ridge_features = pygplates.FeatureCollection('ridges.gpml')

# The times at which to create isochrons.
isochron_creation_times = [40, 30, 20, 10, 0]

# We'll store the created isochrons here - later we'll write it to a file.
isochron_feature_collection = pygplates.FeatureCollection()

# Iterate over the ridge features.
for ridge_feature in ridge_features:

    # Get the ridge left and right plate ids, and time of appearance.
    left_plate_id = ridge_feature.get_left_plate()
    right_plate_id = ridge_feature.get_right_plate()
    time_of_appearance, time_of_disappearance = ridge_feature.get_valid_time()

    # Iterate over our list of creation times for the left/right isochrons.
    for isochron_creation_time in isochron_creation_times:

        # If creation time is later than ridge birth time then we can create an isochron.
        if isochron_creation_time < time_of_appearance:

            # Reconstruct the mid-ocean ridge to isochron creation time.
            # The ridge geometry will be in the same position as the left/right isochrons at that time.
            reconstructed_ridges = []
            pygplates.reconstruct(ridge_feature, rotation_model, reconstructed_ridges, isochron_creation_time)

            # Get the isochron geometry from the ridge reconstruction.
            # This is the geometry at 'isochron_creation_time' (not present day).
            isochron_geometry_at_creation_time = [reconstructed_ridge.get_reconstructed_geometry()
                    for reconstructed_ridge in reconstructed_ridges]

            # Create the left and right isochrons.
            # Since they are conjugates they have swapped left and right plate IDs.
            # And reverse reconstruct the mid-ocean ridge geometries to present day.
            left_isochron_feature = pygplates.Feature.create_reconstructable_feature(
                    pygplates.FeatureType.gpml_isochron,
                    isochron_geometry_at_creation_time,
                    name = ridge_feature.get_name(None),
                    description = ridge_feature.get_description(None),
                    valid_time = (isochron_creation_time, 0),
                    reconstruction_plate_id = left_plate_id,
                    conjugate_plate_id = right_plate_id,
                    reverse_reconstruct = (rotation_model, isochron_creation_time))
            right_isochron_feature = pygplates.Feature.create_reconstructable_feature(
                    pygplates.FeatureType.gpml_isochron,
                    isochron_geometry_at_creation_time,
                    name = ridge_feature.get_name(None),
                    description = ridge_feature.get_description(None),
                    valid_time = (isochron_creation_time, 0),
                    reconstruction_plate_id = right_plate_id,
                    conjugate_plate_id = left_plate_id,
                    reverse_reconstruct = (rotation_model, isochron_creation_time))

            # Add isochrons to feature collection.
            isochron_feature_collection.add(left_isochron_feature)
            isochron_feature_collection.add(right_isochron_feature)

# Write the isochrons to a new file.
isochron_feature_collection.write('isochrons.gpml')

Details

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

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

The ridge features are loaded into a pygplates.FeatureCollection.

ridge_features = pygplates.FeatureCollection('ridges.gpml')

The plate IDs and time period are obtained using pygplates.Feature.get_left_plate(), pygplates.Feature.get_right_plate() and pygplates.Feature.get_valid_time().

left_plate_id = ridge_feature.get_left_plate()
right_plate_id = ridge_feature.get_right_plate()
time_of_appearance, time_of_disappearance = ridge_feature.get_valid_time()

Smaller time values are closer to present day (younger).

if isochron_creation_time < time_of_appearance:

The ridges are reconstructed to their locations at time ‘isochron_creation_time’ using pygplates.reconstruct().

reconstructed_ridges = []
pygplates.reconstruct(ridge_feature, rotation_model, reconstructed_ridges, isochron_creation_time)

A Python list comprehension is used to build a list of pygplates.GeometryOnSphere from a list of pygplates.ReconstructedFeatureGeometry.

isochron_geometry_at_creation_time = [reconstructed_ridge.get_reconstructed_geometry()
        for reconstructed_ridge in reconstructed_ridges]

Isochron features are created using pygplates.Feature.create_reconstructable_feature().

left_isochron_feature = pygplates.Feature.create_reconstructable_feature(
        pygplates.FeatureType.gpml_isochron,
        isochron_geometry_at_creation_time,
        name = ridge_feature.get_name(None),
        description = ridge_feature.get_description(None),
        valid_time = (isochron_creation_time, 0),
        reconstruction_plate_id = left_plate_id,
        conjugate_plate_id = right_plate_id,
        reverse_reconstruct = (rotation_model, isochron_creation_time))

The reverse_reconstruct parameter is needed because all features must store their geometry in present day coordinates which means reverse reconstructing from isochron_creation_time to present day using the rotation model.

Note

The use of None in, for example, ridge_feature.get_name(None) results in a name property only getting created if the ridge feature has a name.

And finally the isochrons are saved to a new file using pygplates.FeatureCollection.write().

isochron_feature_collection.write('isochrons.gpml')

Advanced

If we want to be a bit more robust then we can check that our ridge features are actually ridges and we can make sure they contain left/right plate IDs and a time of appearance/disappearance:

...

# Iterate over the ridge features.
for ridge_feature in ridge_features:

    # Ignore anything that's not a mid-ocean ridge.
    if ridge_feature.get_feature_type() != pygplates.FeatureType.gpml_mid_ocean_ridge:
        continue

    # Get the ridge left and right plate ids, and time of appearance.
    # We don't need to specify 'None', but if we do then it allows us to test if the ridge feature
    # is missing plate IDs or begin/end time period.
    left_plate_id = ridge_feature.get_left_plate(None)
    right_plate_id = ridge_feature.get_right_plate(None)
    valid_time = ridge_feature.get_valid_time(None)

    # Ignore mid-ocean ridges that don't have a left and right plate id and time of appearance.
    if (left_plate_id is None or
        right_plate_id is None or
        valid_time is None):
        continue

    # Extract time of appearance/disappearance from the tuple.
    time_of_appearance, time_of_disappearance = valid_time

    ...

By specifying None in:

left_plate_id = ridge_feature.get_left_plate(None)
right_plate_id = ridge_feature.get_right_plate(None)
valid_time = ridge_feature.get_valid_time(None)
...we will get None returned to us if the feature property (eg, left plate ID) is missing in the ridge feature.
If we didn’t specify None then a default value would be returned if a property was missing. For get_left_plate() and get_right_plate() this is plate ID 0 and for get_valid_time() this is a time period from distant past to distant future.

Table Of Contents

Previous topic

Find overriding plate of closest subducting line

Next topic

Split an isochron into ridge and transform segments

This Page