Hi GPlates community!
I am trying to calculate the average velocity of each plate in pygplates. Ideally this would be output in a list so that I can plot it in a bar chart against the name of each plate. I’m not sure how best to go about doing this so would appreciate any help!
Thanks,
Edward
Hi Edward,
The following pyGPlates code prints out each plate’s name and average velocity magnitude. It uses a static distribution of points across the globe and determines which dynamic plates they fall within at each reconstruction time. Then the average velocity of each plate is the average of the velocities at those points that fall within that plate.
Note that I’m simply using a uniform distribution of points in latitude/longitude space, and since that’s not uniform on the 3D sphere I’ve added a cosine(latitude) weighting factor in the average to give an area-weighted average (where the area around each point varies as cosine of its latitude). It would probably be better to use a distribution that is uniform on the 3D sphere (such as HEALPix) but the final averages shouldn’t differ too much if you use a high enough density of points.
Also note that the code only works with non-deforming topologies (dynamic plate polygons). The next release of pyGPlates (early next year) will support queries on deforming topologies, as well as making it easier to calculate velocities in general (the following code is quite low-level) and faster to calculate them (the part in the following code that finds which plate each point is inside is the bottleneck).
Regards,
John
import math
import os.path
import pygplates
import sys
rotation_filenames = [...]
topology_filenames = [...]
# Load rotation files into a rotation model.
rotation_model = pygplates.RotationModel(rotation_filenames)
# Load the topology files (we only want to do this once, not for every reconstruction time).
topology_feature_collections = []
for topology_filename in topology_filenames:
topology_feature_collections.append( pygplates.FeatureCollection(topology_filename) )
# Spacing of point samples in latitude/longitude space.
grid_spacing_degrees = 1.0
num_latitudes = int(math.floor(180.0 / grid_spacing_degrees))
num_longitudes = int(math.floor(360.0 / grid_spacing_degrees))
# Create the sample points.
points = []
for lat_index in range(num_latitudes):
lat = -90 + (lat_index + 0.5) * grid_spacing_degrees
for lon_index in range(num_longitudes):
lon = -180 + (lon_index + 0.5) * grid_spacing_degrees
points.append(pygplates.PointOnSphere(lat, lon))
velocity_delta_time = 1
oldest_reconstruction_time = 10
# Iterate over reconstruction times from 0 to 'oldest_reconstruction_time'.
for reconstruction_time in range(0, oldest_reconstruction_time + 1):
print('Time:', reconstruction_time)
# Resolved topologies at current reconstruction time.
resolved_topologies = []
pygplates.resolve_topologies(topology_feature_collections, rotation_model, resolved_topologies, reconstruction_time)
# Create a dict mapping resolved topologies to the points contained within them.
# Start with an empty list of points for each resolved plate.
points_in_resolved_topologies = {resolved_topology : [] for resolved_topology in resolved_topologies}
# Find the resolved topology plate that each point is contained within.
# This is quite slow. PlateTectonicTools has 'points_in_polygons.find_polygons()' which speeds this up a lot.
# And a future version of pyGPlates will have a similar approach.
for point in points:
for resolved_topology in resolved_topologies:
if resolved_topology.get_resolved_boundary().is_point_in_polygon(point):
points_in_resolved_topologies[resolved_topology].append(point)
break
# Iterate over the resolved topologies (and each one has a list of points contained within it).
for resolved_topology, points_in_resolved_topology in points_in_resolved_topologies.items():
# Skip current plate if there are no points in it for some reason
# (eg, if plate is small and sampling density of points too low).
if not points_in_resolved_topology:
continue
# Get the resolved plate ID.
plate_id = resolved_topology.get_feature().get_reconstruction_plate_id()
# Get the stage rotation of partitioning plate from 'time + velocity_delta_time' to 'time'.
plate_stage_rotation = rotation_model.get_rotation(reconstruction_time, plate_id, reconstruction_time + velocity_delta_time)
# Calculate velocities at the points in the current plate.
plate_velocity_vectors = pygplates.calculate_velocities(
points_in_resolved_topology,
plate_stage_rotation,
velocity_delta_time,
pygplates.VelocityUnits.cms_per_yr)
# Calculate average velocity magnitude.
#
# Since we're using a non-uniform lat/lon distribution of points we can't just sum up the
# velocities and divide by the number of velocities. Instead need to weight by area.
# The area varies with latitude according to cosine(latitude). Max at equator, min at poles.
sum_weighted_velocity_magnitudes = 0.0
sum_weights = 0.0
for index in range(len(points_in_resolved_topology)):
lat, lon = points_in_resolved_topology[index].to_lat_lon()
velocity = plate_velocity_vectors[index]
weight = math.cos(math.radians(lat))
sum_weighted_velocity_magnitudes += weight * velocity.get_magnitude()
sum_weights += weight
average_plate_velocity = sum_weighted_velocity_magnitudes / sum_weights
# Print plate name and average velocity (cms/yr).
#
# This can instead be written to a text file.
print(resolved_topology.get_feature().get_name(), average_plate_velocity)