Calculate velocities by plate ID¶

This example calculates velocities at all points in geometries of a collection of features using the plate IDs of those features and a rotation model.

Sample code¶

import pygplates

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

# Load the features that contain the geometries we will calculate velocities at.
# Calling them 'domain' features since using them as input to velocities (but can be any type of feature).
domain_features = pygplates.FeatureCollection('features.gpml')

# Calculate velocities at 10Ma using a delta time interval of 1Ma.
reconstruction_time = 10
delta_time = 1

# All reconstructed geometry points and associated (magnitude, azimuth, inclination) velocities.
all_reconstructed_points = []
all_velocities = []

# Iterate over all geometries in all domain features and calculate velocities at each of their points.
for domain_feature in domain_features:

# We need the feature's plate ID to get the equivalent stage rotation of that tectonic plate.
domain_plate_id = domain_feature.get_reconstruction_plate_id()

# Get the rotation of plate 'domain_plate_id' from present day (0Ma) to 'reconstruction_time'.
equivalent_total_rotation = rotation_model.get_rotation(reconstruction_time, domain_plate_id)

# Get the rotation of plate 'domain_plate_id' from 'reconstruction_time + delta_time' to 'reconstruction_time'.
equivalent_stage_rotation = rotation_model.get_rotation(
reconstruction_time, domain_plate_id, reconstruction_time + delta_time)

# A feature usually has a single geometry but it could have more - iterate over them all.
for geometry in domain_feature.get_geometries():

# Reconstruct the geometry to 'reconstruction_time'.
reconstructed_geometry = equivalent_total_rotation * geometry
reconstructed_points = reconstructed_geometry.get_points()

# Calculate velocities at the reconstructed geometry points.
# This is from 'reconstruction_time + delta_time' to 'reconstruction_time' on plate 'domain_plate_id'.
velocity_vectors = pygplates.calculate_velocities(reconstructed_points, equivalent_stage_rotation, delta_time)

# Convert global 3D velocity vectors to local (magnitude, azimuth, inclination) tuples (one tuple per point).
velocities = pygplates.LocalCartesian.convert_from_geocentric_to_magnitude_azimuth_inclination(
reconstructed_points, velocity_vectors)

# Append results for the current geometry to the final results.
all_reconstructed_points.extend(reconstructed_points)
all_velocities.extend(velocities)


Details¶

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

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


The features to calculate velocities at are loaded into a pygplates.FeatureCollection. They can be any type of feature as long as they have a reconstruction plate ID (and of course some geometry).

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

The velocities are calculated at geometries reconstructed to their 10Ma positions.
And the velocities are calculated over a 1My interval from 11Ma to 10Ma.
reconstruction_time = 10
delta_time = 1


pygplates.RotationModel enables to calculate both the rotation from present day to 10Ma of a particular tectonic plate relative to the anchor plate (defaults to zero):

equivalent_total_rotation = rotation_model.get_rotation(reconstruction_time, domain_plate_id)


…and the stage rotation from 11Ma to 10Ma:

equivalent_stage_rotation = rotation_model.get_rotation(
reconstruction_time, domain_plate_id, reconstruction_time + delta_time)

A pygplates.Feature usually contains a single geometry property but sometimes it contains more.
Actually domain_feature.get_geometries() is just a convenient alternative to domain_feature.get_geometry(property_return=PropertyReturn.all).
for geometry in domain_feature.get_geometries():


The geometries extracted from features are in present day coordinates and need to be reconstructed to their 10Ma positions.

reconstructed_geometry = equivalent_total_rotation * geometry

reconstructed_points = reconstructed_geometry.get_points()

The velocities are calculated at the reconstructed geometry positions (10Ma) using the stage rotation.
This returns a list of pygplates.Vector3D (one global cartesian velocity vector per geometry point).
velocity_vectors = pygplates.calculate_velocities(reconstructed_points, equivalent_stage_rotation, delta_time)

If the velocities need to be in local (magnitude, azimuth, inclination) coordinates then the global cartesian vectors can be converted using pygplates.LocalCartesian.convert_from_geocentric_to_magnitude_azimuth_inclination().
Note that each point in reconstructed_points determines a separate local coordinate system. For example, the velocity azimuth is relative to North as viewed from a particular point position.
velocities = pygplates.LocalCartesian.convert_from_geocentric_to_magnitude_azimuth_inclination(
reconstructed_points, velocity_vectors)

Finally we add the reconstructed points and velocities to two large lists for all features.
all_reconstructed_points.extend(reconstructed_points)
all_velocities.extend(velocities)
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 0, "mathjax_asciimath": 1, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0, "math_score": 0.5906901955604553, "perplexity": 7359.518539617965}, "config": {"markdown_headings": false, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "warc_path": "s3://commoncrawl/crawl-data/CC-MAIN-2021-10/segments/1614178362513.50/warc/CC-MAIN-20210301121225-20210301151225-00361.warc.gz"}