Find the total length of ridges and subduction zones

This example resolves topological plate polygons (and deforming networks) and determines:

  • the total length of all mid-ocean ridges, and
  • the total length of all subduction zones

...over a series of geological times.

Sample code

import pygplates


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

# Load the topological plate polygon features (can also include deforming networks).
topology_features = pygplates.FeatureCollection('topologies.gpml')

# Our geological times will be from 0Ma to 'num_time_steps' Ma (inclusive) in 1 My intervals.
num_time_steps = 140

# 'time' = 0, 1, 2, ... , 140
for time in range(num_time_steps + 1):

    # Resolve our topological plate polygons (and deforming networks) to the current 'time'.
    # We generate both the resolved topology boundaries and the boundary sections between them.
    resolved_topologies = []
    shared_boundary_sections = []
    pygplates.resolve_topologies(topology_features, rotation_model, resolved_topologies, time, shared_boundary_sections)

    # We will accumulate the total ridge and subduction zone lengths for the current 'time'.
    total_ridge_length = 0
    total_subduction_zone_length = 0

    # Iterate over the shared boundary sections.
    for shared_boundary_section in shared_boundary_sections:

        # Skip sections that are not ridges or subduction zones.
        if (shared_boundary_section.get_feature().get_feature_type() != pygplates.FeatureType.gpml_subduction_zone and
            shared_boundary_section.get_feature().get_feature_type() != pygplates.FeatureType.gpml_mid_ocean_ridge):
            continue

        # Iterate over the shared sub-segments to accumulate their lengths.
        shared_sub_segments_length = 0
        for shared_sub_segment in shared_boundary_section.get_shared_sub_segments():

            # Each sub-segment has a polyline with a length.
            shared_sub_segments_length += shared_sub_segment.get_resolved_geometry().get_arc_length()

        # The shared sub-segments contribute either to the ridges or to the subduction zones.
        if shared_boundary_section.get_feature().get_feature_type() == pygplates.FeatureType.gpml_mid_ocean_ridge:
            total_ridge_length += shared_sub_segments_length
        else:
            total_subduction_zone_length += shared_sub_segments_length

    # The lengths are for a unit-length sphere so we must multiple by the Earth's radius.
    total_ridge_length_in_kms = total_ridge_length * pygplates.Earth.mean_radius_in_kms
    total_subduction_zone_length_in_kms = total_subduction_zone_length * pygplates.Earth.mean_radius_in_kms

    print "At time %dMa, total ridge length is %f kms and total subduction zone length is %f kms." % (
            time, total_ridge_length_in_kms, total_subduction_zone_length_in_kms)

Details

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

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

The topological features are loaded into a pygplates.FeatureCollection.

topology_features = pygplates.FeatureCollection('topologies.gpml')
The topological features are resolved to the current time using pygplates.resolve_topologies().
By default both pygplates.ResolvedTopologicalBoundary (used for dynamic plate polygons) and pygplates.ResolvedTopologicalNetwork (used for deforming regions) are appended to the list resolved_topologies.
Additionally the resolved topological sections are appended to the list shared_boundary_sections.
resolved_topologies = []
shared_boundary_sections = []
pygplates.resolve_topologies(topology_features, rotation_model, resolved_topologies, time, shared_boundary_sections)
The resolved topological sections are actually what we’re interested in because they contain no duplicate sub-segments.
If we were to iterate over the resolved topologies and their sub-segments, as we do in the Find average area and subducting boundary proportion of topologies sample code, then those sub-segments would have been counted twice (since two adjacent plate polygons will both have sub-segments at the same shared boundary).

The feature type of each topological section is checked to see if it’s a ridge or subduction zone feature type and ignored if it’s neither.

if (shared_boundary_section.get_feature().get_feature_type() != pygplates.FeatureType.gpml_subduction_zone and
    shared_boundary_section.get_feature().get_feature_type() != pygplates.FeatureType.gpml_mid_ocean_ridge):
    continue
Not all parts of a topological section feature’s geometry contribute to the boundaries of topologies.
Little bits at the ends get clipped off.
The parts that do contribute can be found using pygplates.ResolvedTopologicalSection.get_shared_sub_segments().
So we iterate over these and accumulate the lengths of each sub-segment obtained with pygplates.PolylineOnSphere.get_arc_length().
shared_sub_segments_length = 0
for shared_sub_segment in shared_boundary_section.get_shared_sub_segments():
    shared_sub_segments_length += shared_sub_segment.get_resolved_geometry().get_arc_length()

The lengths are for a unit-length sphere so we must multiple by the Earth’s radius (see pygplates.Earth).

total_ridge_length_in_kms = total_ridge_length * pygplates.Earth.mean_radius_in_kms
total_subduction_zone_length_in_kms = total_subduction_zone_length * pygplates.Earth.mean_radius_in_kms

Finally the results for the current ‘time’ are printed.

print "At time %dMa, total ridge length is %f kms and total subduction zone length is %f kms." % (
        time, total_ridge_length_in_kms, total_subduction_zone_length_in_kms)

Output

At time 0Ma, total ridge length is 87002.773452 kms and total subduction zone length is 63502.688936 kms.
At time 1Ma, total ridge length is 87018.115101 kms and total subduction zone length is 63229.149473 kms.
At time 2Ma, total ridge length is 87041.183740 kms and total subduction zone length is 62003.392960 kms.
At time 3Ma, total ridge length is 87156.095568 kms and total subduction zone length is 61475.263778 kms.
At time 4Ma, total ridge length is 89792.644317 kms and total subduction zone length is 61149.051087 kms.
At time 5Ma, total ridge length is 89856.487644 kms and total subduction zone length is 60915.010934 kms.
At time 6Ma, total ridge length is 102897.926344 kms and total subduction zone length is 62442.122395 kms.
At time 7Ma, total ridge length is 102805.357344 kms and total subduction zone length is 62170.240868 kms.
At time 8Ma, total ridge length is 104766.806279 kms and total subduction zone length is 61901.033731 kms.