Reconstruct crustal thickness and tectonic subsidence

This example reconstructs (and deforms) points through topological plate polygons (and deforming networks) to find the following quantities over a time period spanning a past geological time to present day:

• the reconstructed positions of the initial points,

• their crustal stretching factors, and

• their tectonic subsidences.

Sample code

import pygplates


# Create a topological model from our topological features (can be plate polygons and/or deforming networks)
# and rotation file(s).
topological_model = pygplates.TopologicalModel('topologies.gpml', 'rotations.rot')

# Our reconstruction will span from 250Ma to present day in 1 Myr intervals.
initial_time = 250
time_increment = 1

# Reconstruct the initial points through the topological model from initial time to present day.
reconstructed_time_span = topological_model.reconstruct_geometry(
        initial_points,
        initial_time,
        time_increment=time_increment)

# Get the history of reconstructed positions, crustal stretching and tectonic subsidence of the initial points from initial time to present day.
for time in range(initial_time, -1, -time_increment):  # initial_time, initial_time - time_increment, ..., 0
    # Reconstructed positions at the current 'time'.
    reconstructed_points = reconstructed_time_span.get_geometry_points(time)

    # In which resolved topologies are the reconstructed positions located at the current 'time'.
    reconstructed_topology_point_locations = reconstructed_time_span.get_topology_point_locations(time)

    # The crustal stretching factor and tectonic subsidence (in kms) for each point at the current 'time'.
    reconstructed_scalar_values = reconstructed_time_span.get_scalar_values(time)
    reconstructed_crustal_stretching_factors = reconstructed_scalar_values[pygplates.ScalarType.gpml_crustal_stretching_factor]
    reconstructed_tectonic_subsidences = reconstructed_scalar_values[pygplates.ScalarType.gpml_tectonic_subsidence]

Details

First create a topological model from topological features and rotation files.

The topological features can be plate polygons and/or deforming networks.

More than one file containing topological features can be specified here, however we’re only specifying one file.

Also note that more than one rotation file (or even a single pygplates.RotationModel) can be specified here, however we’re only specifying a single rotation file.

topological_model = pygplates.TopologicalModel('topologies.gpml', 'rotations.rot')

Next we get our topological model to reconstruct (and deform) our initial_points from their position at initial_time to present day in 1 Myr intervals using pygplates.TopologicalModel.reconstruct_geometry(). This returns a reconstructed geometry time span containing a history of the reconstructed point positions and their associated scalar values, such as crustal stretching factor and tectonic subsidence, that change over time when passing through deforming networks. Since we did not specify the initial_scalars parameter, the initial crustal stretching factors and initial tectonic subsidences will have default values (of 1.0 and 0.0 respectively) indicating no initial stretching and sea-level subsidence. And since we did not specify the deactivate_points parameter, the default deactivation is used which matches GPlates (when using ‘reconstruct by topologies’ in a green visual layer). This means any initial points on oceanic crust could get subducted and disappear (get deactivated).

reconstructed_time_span = topological_model.reconstruct_geometry(
        initial_points,
        initial_time,
        time_increment=time_increment)

Get the reconstructed positions at the current time. Note that some of the initial points can be deactivated, in which case the number of points in reconstructed_points could be less than initial_points.

reconstructed_points = reconstructed_time_span.get_geometry_points(time)

Query in which resolved topologies the reconstructed positions are located at the current time. Note that the number of values in reconstructed_topology_point_locations matches the number of points in reconstructed_points.

reconstructed_topology_point_locations = reconstructed_time_span.get_topology_point_locations(time)

Extract the reconstructed/deformed crustal stretching factor and tectonic subsidence (in kms) for each point in reconstructed_topology_point_locations. Note that pygplates.ReconstructedGeometryTimeSpan.get_scalar_values() returns a Python dictionary mapping scalar types to their scalar values (unless a specific scalar type is specified). From that dictionary we then extract the crustal stretching factor and tectonic subsidence values. Also note that these builtin scalar types are always available (even when no initial scalar values are provided).

reconstructed_scalar_values = reconstructed_time_span.get_scalar_values(time)
reconstructed_crustal_stretching_factors = reconstructed_scalar_values[pygplates.ScalarType.gpml_crustal_stretching_factor]
reconstructed_tectonic_subsidences = reconstructed_scalar_values[pygplates.ScalarType.gpml_tectonic_subsidence]