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

Hi There GPlates Folks

I’ve been continuing on with some other parts of this but am still struggling with getting the types of features I want from models. I don’t know if I’m pulling the wrong type of topological class. In the first post of this chain, it seemed like I was after TopologicalClosedPlateBoundaries in GPlates image at the top. In the image below I’ve tried to show the features plate features I’m after in GPlates and how I’ve tried reconstructing them through pyGlates and GPlately

The code suggested at the top returns topological boundaries from the Muller2019 model but this removes topological networks creating large holes in North America and South America (middle image). Returning both networks and boundaries and dissolving plateID does a pretty good job but still leaves remnant networks that have different plateID’s from the plates that the sit on (exampld of Yucatan and Scotia Plates).

I’m hoping to do this at any time step in any model without much GIS QC on overlapping polygons or which ones should be combined. Let me know if you think this is possible. Thanks again for your time.

If your code selects via feature type, then try including TopologicalNetwork (in addition to TopologicalClosedPlateBoundary). That should close the gaps (where the deforming networks are) in the Muller2019 model - for example in your middle screenshot (purple).

Although I think you’re probably already doing that since you mentioned the following…

…in which case, yes, having the same plate ID (for a deforming network and the plate it sits on) would work in cases like South America, but not in other cases (where you observe artifacts in your yellow image).

I’m not sure the best way to merge/fill them (when the plate IDs differ). One possible option is to link the two plate IDs (as candidates for merging). For example, the Yucatan network has plate ID 2004 and it sits on rigid plate 101 - if you look in the rotation file they are linked (also see image below) - and this could be detected/automated with pyGPlates.

Also it might help to exclude the inactive networks so that there’s less networks to merge. In the Muller2019 model you exclude the “Inactive_Meshes_and_Topologies” file. These inactive networks won’t leave holes underneath when not present.

Just so you know, later models (like what’s in GPlates 2.5), have deforming networks that overlay the rigid plates. So if you remove the network it won’t leave a hole underneath.

Perfect thank you sir!