I wrote a google sheet and a set of pygplates programs to help streamline some of the process of designing a planet, without having to go through all the work of drawing out a bunch of shapes only to accidentally make an earth clone. You can also use the random points program to generate hot spots, for example.
The sheet above has a bunch of useful features:
- Plate ID Generator
- Generates all the necessary lines for a working plate rotation based on your input. Note, if you set all plates relative to plate ID 1, you can use it to move all your continents into different latitudes at any point in the animation with ease. It also has the helpful side-effect of removing the need to use 3 lines to switch a plate’s relative ID, since every plate has a tie back to one global plate.
- Drift Corrector
- Does what it says on the tin, it corrects the drift caused by moving your plates and not changing the final position at timestamp 0 my. This follows Artifexian’s advice, which I’ll admit I’m not certain is necessary, and could be changed to correct the 0 my line instead.
- CPT generator (with color displays and OKLAB values!)
- Generates colormap files to change how features appear in your simulation. I particularly like using feature type coloring, since it helps visualize what’s happening in a simulation.
- A table of all feature types and their associated possible geometries
- and more…
Using the pygplates programs I’ve written below (click to see the code), you can generate any number of points, a delauney triangulation, and polygons around each point with three-way intersections for each neighbor intersection. I’ve found using this to be helpful for finding new and unique continent shapes.
1. Generating N Random Points
This will generate however many points you want (n_points) and make a new GPlates file with the resulting points.
# generate randomly distributed points
import random
import math
import pygplates
# set a random seed
random.seed(69)
n_points = int(10) # define number of points
# create empty lists and strings for lon and lat
lon_l = []
lat_l = []
coords = []
# generate random lon/lat coordinates in euler space
# using trig here to make sure the distribution is even, also requires valid lon/lat values
for i in range(n_points):
# convert to degrees and into -180 to 180 range
lon = math.degrees(2*math.pi*random.random())-180
lon_l.append(lon)
# convert to degrees and into -90 to 90 range
lat = math.degrees(math.acos(2*random.random()-1))-90
lat_l.append(lat)
# record point in pygplates tuplet format
coord = pygplates.LatLonPoint(lat, lon)
coords.append(coord)
# generate random points in GPlates loadable format
# create feature list
points_f = [] # pyGPlates uses this to construct a list of features with metadata
# create initial ID
id = 100
for i in coords: # use i in coords, meaning for each random point generated above
# define the coordinates for each point
random_geometry = pygplates.PointOnSphere(i)
# create a reconstructable feature
random_feature = pygplates.Feature.create_reconstructable_feature(
pygplates.FeatureType.gpml_unclassified_feature, # define feature type
random_geometry, # define geometry
name='Random Point %i' % id, # assign name
valid_time=(pygplates.GeoTimeInstant.create_distant_past(), pygplates.GeoTimeInstant.create_distant_future()), # assign valid time
reconstruction_plate_id = id # assign plate ID
)
# add feature to feature list
points_f.append(random_feature)
# create readout
name='Random Point %i' % id
print('Created %a at %r, ID: %i' % (name, pygplates.LatLonPoint.to_lat_lon(i), id))
# iterate plate ID
id += 1
# create the feature collection file once every point is in the feature list
points_fc = pygplates.FeatureCollection(points_f)
points_fc.write('random_points.gpml')
2. Connecting all of the dots
This will generate a file with every connection between every point. It’s based on code that finds the nearest points, so you can change max_n to a lower integer if you have a lot of points to compare. You don’t actually need to open neighbors.gpml at any point, since you’ll have python write a more useful visualization in the next step.
import pygplates
# Load some features.
features = pygplates.FeatureCollection('random_points.gpml')
neighbors_f = []
print("Loaded collection has "+str(len(features))+" features")
# Input maximum number of neighbors
max_n = len(features)
# The minimum distance to all features and the nearest feature.
for feature in features:
# Set initial point
point1 = feature.get_geometry()
# # Read out initial point in lat/lon
# print("")
# print("For point "+str(point1.to_lat_lon()))
# print("")
# Initialize neighbors list
neighbors = []
# For each feature, find the distance to the initial point
for feature2 in features:
# Set secondary point
point2 = feature2.get_geometry()
# Find distance between initial and secondary points
d_rad = pygplates.GeometryOnSphere.distance(point1, point2)
# Store distance and points as tuple
tuple = (d_rad, point1.to_lat_lon(), feature.get_reconstruction_plate_id(), point2.to_lat_lon(), feature2.get_reconstruction_plate_id())
# Append data into neighbors list
neighbors.append(tuple)
# # Print distance for each point
# print(str(d_rad)+" from "+str(point2.to_lat_lon()))
# Sort neighbors ascending (by distance, first in tuple)
neighbors.sort()
# Initialize iterator
i = 0
# Initialize nearest neighbor list
nearest_neighbors = []
# Add the first max_n tuples from neighbors to nearest_neighbors
for neighbor in neighbors:
# Check if iterator is less than maximum neighbors
if i <= max_n:
# Append to nearest_neighbors
nearest_neighbors.append(neighbor)
elif i > max_n:
# Do nothing if iterator is greater than maximum neighbors
continue
# Iterate
i = i+1
# Generate polyline features from nearest_neighbors
for neighbor in nearest_neighbors:
distance, i_point, i_id, s_point, s_id = neighbor
# Ignore the entry where distance = 0
if distance == 0:
continue
# Generate a polyline feature from i_point to s_point
elif distance > 0:
neighbor_f = pygplates.Feature.create_reconstructable_feature(
pygplates.FeatureType.gpml_unclassified_feature,
pygplates.PolylineOnSphere([i_point, s_point]),
name = str(round(distance*pygplates.Earth.mean_radius_in_kms))+" km from "+str(s_id),
valid_time = (pygplates.GeoTimeInstant.create_distant_past(), pygplates.GeoTimeInstant.create_distant_future()),
reconstruction_plate_id = i_id,
conjugate_plate_id = s_id
)
neighbors_f.append(neighbor_f)
# Generate and write a feature collection for neighbors
neighbors_fc = pygplates.FeatureCollection(neighbors_f)
neighbors_fc.write('neighbors.gpml')
print("New feature collection generated with "+str(len(neighbors_f)/2)+" connections") # Note that each connection has 2 lines overlapping, which is important later.
3. Delauney triangulation from the chaos
This will go through every intersecting line and remove the shorter of the two lines (so not including precisely overlapping lines)
import pygplates as pg
random_points = pg.FeatureCollection('random_points.gpml')
neighbor_lines = []
neighbors = []
for point in random_points:
for neighbor in random_points:
d = pg.GeometryOnSphere.distance(point.get_geometry(), neighbor.get_geometry())
if d > 0:
neighbors.append((
d,
point.get_geometry().to_lat_lon(),
point.get_reconstruction_plate_id(),
neighbor.get_geometry().to_lat_lon(),
neighbor.get_reconstruction_plate_id()
))
for d, p1, id1, p2, id2 in neighbors:
feature = pg.Feature.create_reconstructable_feature(
pg.FeatureType.gpml_unclassified_feature, # Set feature type
pg.PolylineOnSphere([p1, p2]), # Set geometry
name = str(round(d*pg.Earth.mean_radius_in_kms))+" km from "+str(id2), # Set name
valid_time = (pg.GeoTimeInstant.create_distant_past(), pg.GeoTimeInstant.create_distant_future()), # Set times
reconstruction_plate_id = id1, # Set plate ID
conjugate_plate_id = id2 # Set conjugate plate ID
)
neighbor_lines.append(feature)
# Write the feature collection file
pg.FeatureCollection(neighbor_lines).write('delauney.gpml')
overlappers = []
# Compare each line to every other line to find intersections
for line in neighbor_lines:
len1 = line.get_geometry().get_arc_length()
p1 = line.get_geometry()[0]
p2 = line.get_geometry()[1]
id1 = line.get_reconstruction_plate_id()
id2 = line.get_conjugate_plate_id()
f_id1 = line.get_feature_id()
for neighbor in neighbor_lines:
len2 = neighbor.get_geometry().get_arc_length()
p3 = neighbor.get_geometry()[0]
p4 = neighbor.get_geometry()[1]
id3 = neighbor.get_reconstruction_plate_id()
id4 = neighbor.get_conjugate_plate_id()
f_id2 = neighbor.get_feature_id()
# Ignore shared origins, shared ends, and conjugate copies
if id1 != id3 and id1 != id4 and id2 != id3 and id2 != id4:
intersect = pg.GeometryOnSphere.distance( # Calculate the distance between the lines
pg.PolylineOnSphere([p1, p2]), # represented by each tuple in neighbors[]
pg.PolylineOnSphere([p3, p4]),
geometry1_is_solid=False,
geometry2_is_solid=False
)
# For each intersection, add the longer feature to a list of overlappers
if intersect == 0:
if len1 > len2:
overlappers.append(f_id1)
if len2 > len1:
overlappers.append(f_id2)
network = []
for feature in neighbor_lines:
i = 0
for overlap in overlappers:
if overlap == feature.get_feature_id():
i = 1
if i == 0:
network.append(feature)
pg.FeatureCollection(network).write('delauney.gpml')
4. Making regions (this was very hard to figure out)
Once you have this, you will have a set of polygons centered around each random point generated. Because of how I chose the three-way intersection points, this isn’t a perfect way to do it, but I think that’s part of the charm? At any rate, make sure you open random_points.gpml, midpoints.gpml, and regions.gpml together to see the output.
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('delauney.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
Hope this helps!