It would be helpful to be able to set points to move relative to three plate IDs. In addition to this, pygplates should get a function to generate a three-way midpoint between three points. There are a lot of uses for this feature, like places where three oceanic plates diverge, generating complex point-based networks for orogenies, and figuring out how different subduction zones interact.
This math paper describes ways to calculate the multiple different center points (vertex-median, areal, etc.) Only adding one of them would add a lot of value, adding more than one would give us even more freedom of movement. In theory, it would be possible to add n-plate relative motion, but that might be gratuitous.
I wrote some code to generate n random points, then all the lines that connect to every other point, then to delete the longer line of every intersection (delauney triangulation). I’ll include that
and this in another post, since it’s a lot more useful for the worldbuilders. From there, I wrote this monster of a program that takes the output of the delauney triangulation and adds vertex median midpoints to generate regions with each random point at its center, resulting in a randomly generated global map of plates that’s very useful for starting up a worldbuilding project.
Although the code wasn’t written for my suggestion directly, I hope it’s helpful if you decide to implement it!
Also… I’m not a programmer. any code improvement suggestions would be much appreciated.
Reveal pygplates code
import pygplates as pg
import numpy as np
import math
# load the feature collections file
random_points = pg.FeatureCollection('random_points.gpml')
# Load the neighbor lines feature collection
lines = pg.FeatureCollection('neighbors.gpml')
# Create an empty list for the midpoint features
midpoints_f = []
# Define a sorter for the cluster list
def sort(cluster):
return(sorted(cluster, key = lambda x: float(x[0]), reverse = False))
# Find each point's cluster of lines and generate a region topology
for point in random_points:
# Get the point ID for comparison, and set a datum to the north pole
point_id = point.get_reconstruction_plate_id() # Get plate ID
datum = pg.GreatCircleArc(point.get_geometry(),pg.LatLonPoint.north_pole) # Create datum arc
datum_axis = datum.get_rotation_axis_lat_lon() # Get datum axis
print(point_id) # Print for log
# Clear/generate a cluster list, which we'll polar sort for the topology generation
cluster = []
# For each line in the neighbors feature collection, find the lines with matching plate IDs.
for line in lines:
line_id = line.get_reconstruction_plate_id()
# If the line originates from the point (shares its ID)
if line_id == point_id:
conj_line_id = line.get_conjugate_plate_id()
# Find the conjugate point to create an arc
for p in random_points:
p_id = p.get_reconstruction_plate_id()
# If match is found, generate an arc from center to conjugate and compare its axis to datum.
if p_id == conj_line_id:
arc = pg.GreatCircleArc(point.get_geometry(),p.get_geometry()) # Create arc
arc_axis = arc.get_rotation_axis_lat_lon() # Get axis
arc_d = (datum_axis[1]+180)-(arc_axis[1]+180) # Calculate difference from datum axis
cluster.append((arc_d, arc, p))
# Do nothing for non-matching cases
# Do nothing for non-matching cases
cluster = sort(cluster)
i = 1
for set in cluster:
# Find the next set in the cluster for generating the vertex median
if i < len(cluster):
next = cluster[i]
# If pulling from the last set, use the first set
elif i == len(cluster):
next = cluster[0]
print(set[0], "from datum to", set[2].get_reconstruction_plate_id(), "next is", next[2].get_reconstruction_plate_id())
# Set the midpoint for the current arc
set_mid = set[1].get_arc_point(0.5)
# Generate features at the midpoint
midpoint_f = pg.Feature.create_reconstructable_feature(
pg.FeatureType.gpml_unclassified_feature,
set_mid,
name = "",
valid_time = (pg.GeoTimeInstant.create_distant_past(), pg.GeoTimeInstant.create_distant_future()),
reconstruction_plate_id = point.get_reconstruction_plate_id(),
conjugate_plate_id = set[2].get_reconstruction_plate_id()
)
hs_midpoint_f = pg.Feature.create_tectonic_section(
pg.FeatureType.gpml_unclassified_feature,
set_mid,
name = "",
valid_time = (pg.GeoTimeInstant.create_distant_past(), pg.GeoTimeInstant.create_distant_future()),
left_plate = point.get_reconstruction_plate_id(),
right_plate = set[2].get_reconstruction_plate_id(),
reconstruction_method = 'HalfStageRotationVersion2'
)
midpoints_f.append(midpoint_f)
midpoints_f.append(hs_midpoint_f)
# Set the geometry to find the vertex median intersection
set_p = set[2].get_geometry().to_lat_lon_point_list()
next_mid = next[1].get_arc_point(0.5).to_lat_lon_point_list()
next_p = next[2].get_geometry().to_lat_lon_point_list()
set_mid = set[1].get_arc_point(0.5).to_lat_lon_point_list()
# Define points in great circle 1
p1_lat1 = set_p[0].get_latitude()
p1_long1 = set_p[0].get_longitude()
p1_lat2 = next_mid[0].get_latitude()
p1_long2 = next_mid[0].get_longitude()
# Define points in great circle 2
p2_lat1 = next_p[0].get_latitude()
p2_long1 = next_p[0].get_longitude()
p2_lat2 = set_mid[0].get_latitude()
p2_long2 = set_mid[0].get_longitude()
# Convert points in great circle 1, degrees to radians
p1_lat1_rad = ((math.pi * p1_lat1) / 180.0)
p1_long1_rad = ((math.pi * p1_long1) / 180.0)
p1_lat2_rad = ((math.pi * p1_lat2) / 180.0)
p1_long2_rad = ((math.pi * p1_long2) / 180.0)
# Convert points in great circle 2, degrees to radians
p2_lat1_rad = ((math.pi * p2_lat1) / 180.0)
p2_long1_rad = ((math.pi * p2_long1) / 180.0)
p2_lat2_rad = ((math.pi * p2_lat2) / 180.0)
p2_long2_rad = ((math.pi * p2_long2) / 180.0)
# Put in polar coordinates
x1 = math.cos(p1_lat1_rad) * math.cos(p1_long1_rad)
y1 = math.cos(p1_lat1_rad) * math.sin(p1_long1_rad)
z1 = math.sin(p1_lat1_rad)
x2 = math.cos(p1_lat2_rad) * math.cos(p1_long2_rad)
y2 = math.cos(p1_lat2_rad) * math.sin(p1_long2_rad)
z2 = math.sin(p1_lat2_rad)
cx1 = math.cos(p2_lat1_rad) * math.cos(p2_long1_rad)
cy1 = math.cos(p2_lat1_rad) * math.sin(p2_long1_rad)
cz1 = math.sin(p2_lat1_rad)
cx2 = math.cos(p2_lat2_rad) * math.cos(p2_long2_rad)
cy2 = math.cos(p2_lat2_rad) * math.sin(p2_long2_rad)
cz2 = math.sin(p2_lat2_rad)
# Get normal to planes containing great circles
# np.cross product of vector to each point from the origin
N1 = np.cross([x1, y1, z1], [x2, y2, z2])
N2 = np.cross([cx1, cy1, cz1], [cx2, cy2, cz2])
# Find line of intersection between two planes
L = np.cross(N1, N2)
# Find two intersection points
X1 = L / np.sqrt(L[0]**2 + L[1]**2 + L[2]**2)
X2 = -X1
i_lat1 = math.asin(X1[2]) * 180./np.pi
i_long1 = math.atan2(X1[1], X1[0]) * 180./np.pi
i_lat2 = math.asin(X2[2]) * 180./np.pi
i_long2 = math.atan2(X2[1], X2[0]) * 180./np.pi
# If the point being tested is above the equator, use the second intersection
if point.get_geometry().to_lat_lon_point_list()[0].get_latitude() > 0:
vertex_median = pg.PointOnSphere(i_lat2, i_long2)
else:
vertex_median = pg.PointOnSphere(i_lat1, i_long1)
# Generate a vertex point
vm_f = pg.Feature.create_reconstructable_feature(
pg.FeatureType.gpml_unclassified_feature,
vertex_median,
name = str("vertex median of "+str(point.get_reconstruction_plate_id())+" between "+str(set[2].get_reconstruction_plate_id())+" and "+str(next[2].get_reconstruction_plate_id())),
valid_time = (pg.GeoTimeInstant.create_distant_past(), pg.GeoTimeInstant.create_distant_future()),
reconstruction_plate_id = point.get_reconstruction_plate_id(),
conjugate_plate_id = set[2].get_reconstruction_plate_id()
)
midpoints_f.append(vm_f)
# Iterate
i = i+1
midpoints_fc = pg.FeatureCollection(midpoints_f)
midpoints_fc.write('midpoints.gpml')
def find_referenced_features(
features,
referenced_feature_id_strings):
referenced_feature_features = []
for referenced_feature_id_string in referenced_feature_id_strings:
for feature in features:
if feature.get_feature_id().get_string() == referenced_feature_id_string:
referenced_feature_features.append(feature)
break
return referenced_feature_features
def create_topological_sections(
referenced_features,
topological_geometry_type):
topological_sections = []
for referenced_feature in referenced_features:
topological_section = pg.GpmlTopologicalSection.create(
referenced_feature,
topological_geometry_type=topological_geometry_type)
if topological_section:
topological_sections.append(topological_section)
return topological_sections
def create_topological_polygon_feature(
features,
topological_polygon_referenced_feature_id_strings,
topological_polygon_feature_type):
topological_polygon_referenced_features = find_referenced_features(
features,
topological_polygon_referenced_feature_id_strings)
topological_polygon_sections = create_topological_sections(
topological_polygon_referenced_features,
pg.GpmlTopologicalPolygon)
topological_polygon_feature = pg.Feature.create_topological_feature(
topological_polygon_feature_type,
pg.GpmlTopologicalPolygon(topological_polygon_sections))
return topological_polygon_feature
# The topologies that will be created from each point's cluster
topologies = []
for point in random_points:
# Clear/create the cluster list
cluster = []
origin_id = point.get_reconstruction_plate_id()
datum = pg.GreatCircleArc(point.get_geometry(),pg.LatLonPoint.north_pole) # Create datum arc
datum_axis = datum.get_rotation_axis_lat_lon() # Get datum axis
# Compare the axis of every vertex sharing the origin's plate ID to the datum
for vertex in midpoints_fc:
vertex_id = vertex.get_reconstruction_plate_id()
if vertex_id == origin_id:
arc = pg.GreatCircleArc(point.get_geometry(),vertex.get_geometry()) # Create arc
arc_axis = arc.get_rotation_axis_lat_lon() # Get axis
arc_d = (datum_axis[1]+180)-(arc_axis[1]+180) # Calculate difference from datum axis
cluster.append((arc_d, arc, vertex))
# Sort the list of queried features
cluster = sort(cluster)
# The list of geometries that will be turned into a topology
region_polygon = []
for set in cluster:
region_polygon.append(set[2].get_feature_id().get_string()) # This is the format the function needs
region_f = create_topological_polygon_feature(
midpoints_fc, # Feature collection to pull sections from
region_polygon, # List of feature IDs, in order
pg.FeatureType.gpml_topological_closed_plate_boundary)
topologies.append(region_f) # Add the region boundary to topologies
pg.FeatureCollection(topologies).write('regions.gpml') # Write the topologies file