Returning geometry of larger plates in PyGplates or GPlately

Hi All

I’m looking for a way to return the geometry of the Plate Polygons from a model not necessarily the static polygons in a model through any of the APIs. As an example below, I’m looking to return the geometry (list of lon/lats) of the larger Farallon plate (highlighted in white) not the Alexander or Wrangellia terranes that are onboard it.

I’ve tried doing this through GPlately, where I make a topological model, get_all_topologies(). As there are a mix of TopologicalClosedPlateBoundary’s and TopologicalNetworks, filtering for just Boundaries leaves large holes in the model where topological networks are. Dissolving the topologies by their plate ID leaves some artifacts where there aren’t plate boundaries.

I’ve tried going through pyGPlates but thats a bit of nightmare iterating through properties.

Lastly, going through the gplates_ws_proxy and using get_topology leaves similar holes where topological networks are.

Any help would be awesome. Its been a few days trying to figure out how to get just large-scale plates and their XY geometry through one of the programming APIs. Hoping to end up with some table of:
Plate ID | Plate Name | Geometry

Thanks in advance
Danny

Hi Danny,

So it sounds like you want the points of the reconstructed plate (at a particular time). In this case you can use pyGPlates to resolve (reconstruct) a topological model to a reconstruction time, and then find the Farallon plate using its feature ID, and then you extract its reconstructed geometry.

So first you would need the feature ID which can be found in GPlates (right-click it and Copy)…

image

Then use something like the following code with pyGPlates:


def find_resolved_plate(topological_snapshot, plate_feature_id):
    """Find a resolved topological boundary in a topological snapshot using the specified feature ID."""
    resolved_plates = topological_snapshot.get_resolved_topologies(
            pygplates.ResolveTopologyType.boundary)
    for resolved_plate in resolved_plates:
        if (resolved_plate.get_feature().get_feature_id().get_string() ==
            plate_feature_id):
            return resolved_plate
    
    return None  # not found

def get_resolved_plate_id_and_name_and_geometry(resolved_plate):
    """Return the plate ID, name and list of boundary (lon, lat) points of a resolved plate."""
    # Get plate ID and name.
    plate_feature = resolved_plate.get_feature()
    plate_name = plate_feature.get_name()
    plate_id = plate_feature.get_reconstruction_plate_id()

    # Get reconstructed plate polygon's points as a list of (lon, lat).
    resolved_plate_polygon = resolved_plate.get_resolved_boundary()
    plate_boundary_list_of_lon_lat = []
    for point in resolved_plate_polygon.get_exterior_ring_points():
        lat, lon = point.to_lat_lon()
        plate_boundary_list_of_lon_lat.append((lon, lat))

   return plate_id, plate_name, plate_boundary_list_of_lon_lat

reconstruction_time = 205
farallon_plate_feature_id = 'GPlates-ec8fe0bf-2fc8-4dd6-8115-04f8d2532e35'

# Create a topological model from topological and rotation features.
topological_model = pygplates.TopologicalModel(topology_features, rotation_model)
# Create a topological snapshot at a specific time.
topological_snapshot = topological_model.topological_snapshot(reconstruction_time)
# Find a particular plate using its feature ID.
farallon_resolved_plate = find_resolved_plate(topological_snapshot, farallon_plate_feature_id)
if farallon_resolved_plate:
    # Get plate ID and name and the reconstructed plate polygon's points as a list of (lon, lat).
    farallon_plate_id, farallon_plate_name, farallon_plate_boundary_lon_lats = get_resolved_plate_id_and_name_and_geometry(farallon_resolved_plate)

…I haven’t tried the above code so there might be a typo in there. But that’s the general idea.

If you wanted to just get all large plates you could filter using the plate area rather than searching by feature IDs. So something like:

def find_large_resolved_plates(topological_snapshot, min_plate_area_steradians):
    """Find resolved topological boundaries in a topological snapshot exceeding a minimum area."""
    large_resolved_plates = []

    resolved_plates = topological_snapshot.get_resolved_topologies(
            pygplates.ResolveTopologyType.boundary)
    for resolved_plate in resolved_plates:
        if (resolved_plate.get_resolved_boundary().get_area() >
            min_plate_area_steradians):
            large_resolved_plates.append(resolved_plate)

    return large_resolved_plates

min_plate_area_square_kms = 20e+6
min_plate_area_steradians = min_plate_area_square_kms / pygplates.Earth.mean_radius_in_kms**2

# Create a topological model from topological and rotation features.
topological_model = pygplates.TopologicalModel(topology_features, rotation_model)
# Reconstruct topologies through a sequence of times.
for time in times:
    # Create a topological snapshot at a specific time.
    topological_snapshot = topological_model.topological_snapshot(time)
    # Filter those rigid plates exceeding a minimum area.
    large_resolved_topological_boundaries = find_large_resolved_plates(topological_snapshot, min_plate_area_steradians)
    for large_resolved_topological_boundary in large_resolved_topological_boundaries:
        # Get plate ID and name and the reconstructed plate polygon's points as a list of (lon, lat).
        plate_id, plate_name, plate_boundary_lon_lats = get_resolved_plate_id_and_name_and_geometry(large_resolved_topological_boundary)

That’s mostly for getting present day geometry. But there are easier functions to use, documented at the top of the Feature class, that mean you don’t have to iterate through properties. For example, to get a regular present-day geometry use get_geometry() or to get a topological geometry use get_topological_geometry(). The latter would retrieve the intersecting lines that form a topological plate polygon (note that you don’t get an actual polygon until the topologies are reconstructed/resolved to a particular reconstruction time which is what the above code does).

A post was split to a new topic: Mapping plate IDs to plate names