Is there a way to combine the geometry of a line topology (especially cutting off the fringes) with the properties (plate ids) of the original features? We can make a line topology out of multiple features. An example of such a situation is where there are two mid-ocean ridges more or less in line, but between different sets of plates, such as between the North and South Atlantic (or even more complex triple junction configurations). If the topology itself is made one mid-ocean ridge feature it can have only two plate ids, but we need at least four.
Hi Thomas,
So you have two mid-ocean ridges joined end-to-end into a single topological line. But each ridge feature has a different pair of conjugate plate IDs. And so the problem becomes what plate IDs do you specify in the topological line feature.
I would suggest not specifying any plate IDs since it’s usually not meaningful to specify plate IDs for topological lines or networks (since they are typically deforming). The exception is topological boundaries since these are non-deforming and so a plate ID specifies the movement of the rigid interior crust (eg, calculating surface velocities). In fact, when you create a topology in GPlates, it will only prompt you for a plate ID if the topology is a boundary.
You can create your topological line with feature type MidOceanRidge. It won’t ask for the conjugate plate IDs, but you can still add them if you want. However they are not necessary and, if specified, are not used to determine the motion of the topological line.
So regardless of whether you specify plate IDs or not, your topological line should still move as expected. In other words, the two original mid-ocean ridges will move as they normally do (using their plate IDs) and, at each reconstruction time, get joined/intersected into a single resolved mid-ocean ridge line.
You can then add your topological MidOceanRidge feature to the two adjoining plates. If you have a triple junction then you’ll probably need two separate line features (eg, MidOceanRidge added to two plates and SubductionZone added to three plates).
Regards,
John
Hi John,
thank you very much for your explanation. Including some sketches and explanation of what I would like to do, makes my problem clearer. We want to calculate isochrons with pyGPlates. This is a sketch of an initial situation with two ridges between four plates moving at different velocities.
After some time, the area where the mid ocean ridges overlap becomes a mess because isochrons are positioned over each other as illustrated in this sketch.
When the overlapping parts of the ridges would be cut off before creating the isochrons it would become less of a mess. Do you know a way to do that? I was thinking to convert lines to a line topology, doing the cut off and convert back to the lines (without the parts that were cut off when making the line topology). If this isn’t possible I guess I could use another tool like the pyGPlates distance to determine the intersection of the two lines or something from shapely.
Hi Thomas,
That’s interesting, I see your original intent more clearly now. So you want to cut off the overlapping parts of the mid-ocean ridges before they are used to create isochrons. And each created isochron gets the conjugate plate IDs from the ridge used to create it.
In this case you can still use a topological line to intersect the two mid-ocean ridges and cut off the overlapping parts. Then iterate over the two clipped ridges and create isochrons from them. This assumes the ridges actually intersect. If they don’t then they’ll just be joined together (end of one to the start of the other).
I’ve written some pyGPlates code to show how to do this (but haven’t run or tested it).
rotation_model = pygplates.RotationModel(...)
reconstruction_time = 100
mid_ocean_ridge_feature_1 = ...
mid_ocean_ridge_feature_2 = ...
# Create a mid-ocean ridge topological line feature referencing two regular
# mid-ocean ridge features.
#
# NOTE: Requires pyGPlates revision 24.
mid_ocean_ridge_topological_line_feature = pygplates.Feature.create_topological_feature(
pygplates.FeatureType.gpml_mid_ocean_ridge,
pygplates.GpmlTopologicalLine([
pygplates.GpmlTopologicalSection.create(
mid_ocean_ridge_feature_1,
topological_geometry_type=pygplates.GpmlTopologicalLine),
pygplates.GpmlTopologicalSection.create(
mid_ocean_ridge_feature_2,
topological_geometry_type=pygplates.GpmlTopologicalLine)]),
valid_time=(float('inf'), 0.0))
# Resolve the mid-ocean ridge topological line at the reconstruction time.
resolved_topological_lines = []
pygplates.resolve_topologies(
[mid_ocean_ridge_topological_line_feature,
mid_ocean_ridge_feature_1,
mid_ocean_ridge_feature_2],
rotation_model,
resolved_topological_lines,
reconstruction_time,
# We're only interested in topological *lines*...
resolve_topology_types=pygplates.ResolveTopologyType.line)
# There's only one topological line...
resolved_topological_line = resolved_topological_lines[0]
# Get the two clipped reconstructed mid-ocean ridge lines.
line_sub_segments = resolved_topological_line.get_line_sub_segments()
clipped_mid_ocean_ridge_features = [line_sub_segment.get_resolved_feature()
for line_sub_segment in line_sub_segments]
# Then create isochrons from these two clipped ridges.
#
# The following is similar to the code at:
# http://www.gplates.org/docs/pygplates/sample-code/pygplates_create_conjugate_isochrons_from_ridge.html
#
isochron_features = []
for clipped_mid_ocean_ridge_feature in clipped_mid_ocean_ridge_features:
# The geometry has already been reconstructed to 'reconstruction_time' due
# to pygplates.resolve_topologies() above. We don't need to reconstruct it here.
# Normally a feature contains present-day geometry, so this might be a little surprising.
clipped_reconstructed_mid_ocean_ridge_geometry = clipped_mid_ocean_ridge_feature.get_geometry()
# Create the left and right isochrons.
# Since they are conjugates they have swapped left and right plate IDs.
# And reverse reconstruct the mid-ocean ridge geometries to present day.
left_isochron_feature = pygplates.Feature.create_reconstructable_feature(
pygplates.FeatureType.gpml_isochron,
clipped_reconstructed_mid_ocean_ridge_geometry,
name = clipped_mid_ocean_ridge_feature.get_name(None),
description = clipped_mid_ocean_ridge_feature.get_description(None),
valid_time = (reconstruction_time, 0),
reconstruction_plate_id = clipped_mid_ocean_ridge_feature.get_left_plate(),
conjugate_plate_id = clipped_mid_ocean_ridge_feature.get_right_plate(),
# Feature should store present day geometry, so reverse reconstruct from 'reconstruction_time'...
reverse_reconstruct = (rotation_model, reconstruction_time))
right_isochron_feature = pygplates.Feature.create_reconstructable_feature(
pygplates.FeatureType.gpml_isochron,
clipped_reconstructed_mid_ocean_ridge_geometry,
name = clipped_mid_ocean_ridge_feature.get_name(None),
description = clipped_mid_ocean_ridge_feature.get_description(None),
valid_time = (reconstruction_time, 0),
reconstruction_plate_id = clipped_mid_ocean_ridge_feature.get_right_plate(),
conjugate_plate_id = clipped_mid_ocean_ridge_feature.get_left_plate(),
# Feature should store present day geometry, so reverse reconstruct from 'reconstruction_time'...
reverse_reconstruct = (rotation_model, reconstruction_time))
isochron_features.append(left_isochron_feature)
isochron_features.append(right_isochron_feature)
# Write isochrons to a file, so can load into GPlates.
pygplates.FeatureCollection(isochron_features).write('isochrons.gpml')
Regards,
John
Hi John,
that is awesome! More than I could ask for 
Seems only the create_topological_feature is not available in pyGPlates revision 18, but that part can be done in GPlates quite easily. I will test it and let you know!
Have a good weekend! Regards,
Thomas
Hi John,
Trying to get this to work, there are a few things I run into.
One is: how to define which part of the line feature to cut off and which part to preserve. In GPlates it seems that the part that is clicked is preserved. Is there a way in pyGPlates to define this? It seems to depend on the order in which the features are provided to create_topological_feature. With two lines there are only two options for the order in which you input them [1,2] or [2,1], while there are four possible combinations of segments 1a, 1b, 2a, 2b (with a and b being segments on either side of the intersection of lines 1 and 2).
Another problem is that resolve_topologies returns an empty list, although the mid_ocean_ridge_topological_line_feature appears okay when I save it to a gpml. Do you happen to have a list of common things that go wrong when using this method?
Best regards,
Thomas
Hi Thomas,
Try adding reverse_order=True or reverse_order=False to pygplates.GpmlTopologicalSection.create(). Since there’s two calls to it in the above code, there are the 4 combinations you mention.
By the way, here’s the documentation for that function (since it’s not publicly released yet; we’re releasing this very soon though, as one final beta release to support Python 3 before pyGPlates 1.0 officially comes out). Note the part circled in red…
Oh yeah, sorry I forgot to specify the valid_time parameter when creating the mid-ocean ridge topological line. I’ve corrected that in the code above (my previous post) where valid_time=(float('inf'), 0.0) is added to specify that the ridge lives from distant past float('inf') to present day 0.0. You can change to a more suitable time range if you like. When valid_time is not specified, the feature does not appear to exist at any time (and so resolving topologies returns an empty list).
Regards,
John
Hi John,
I am still trying to get this to work. This is my procedure:
- make the topological line feature in GPlates (as it appears easier to control what we are making);
- save it as a DeformingRegionEdge (to be able to subsequently extract it easily from the file, as there are no other features of this type) in the same gpml file as the original lines;
- in pyGPlates create a FeatureCollection by reading the gpml file consisting of the original lines and the topology;
- loop through the features until there is a DeformingRegionEdge;
- input that feature into the code above (here, only copied the explanatory headers):
# Resolve the mid-ocean ridge topological line at the reconstruction time.,
# There's only one topological line...,
# Get the two clipped reconstructed mid-ocean ridge lines.
This crashes as resolved_topological_lines is an empty list.
It appears to me that the geometries of the line features might be not connected to the line topology as that is the only feature without a geometry.
I have sent you the gpml, rot and Python files. Could you have a look and put me on the right track?
Thanks in advance. Best regards,
Thomas
Hi Thomas,
Your GPML file looks fine - it has the topological line feature connected to two regular line features. Unfortunately your Python file was scrubbed for security reasons, so I couldn’t look at it, however I think the issue is in my original example code where I failed to include the two regular line features, along with the topological line feature, when resolving the topologies. So I’ve adjusted the post above with the following change…
pygplates.resolve_topologies(
[mid_ocean_ridge_topological_line_feature,
mid_ocean_ridge_feature_1,
mid_ocean_ridge_feature_2],
...
…that includes all three features, since mid_ocean_ridge_topological_line_feature references mid_ocean_ridge_feature_1 and mid_ocean_ridge_feature_2, and so they all must be included.
By the way, so that’s for pyGPlates but I should also note that in GPlates, when you build the topological line (from two regular intersecting lines), it chooses the longest clipped part of each line. If that’s not what you want then pyGPlates will give you more choice, as noted above. Or you can just edit the gpml:reverseOrder flags in the GPML file produced by GPlates.
Regards,
John
Thanks John, that puts me back on track.
But now for the next step: as there are multiple topologies in my model, I would like to create lists of the TopologicalSection with it’s associated targetFeatures and use these lists as input for resolving in a next step. So far I have been able to determine whether a feature is a topological section with this code:
for feature in feature_collection:
for p in feature:
if "GpmlTopologicalSection" in str(p.get_value()):
topological_line_feature = feature
But I cannot find how to get the lines out.
Is there a method available in pyGPlates to access the elements of the gpml:centerLineOf property to extract the associated line features? My work around solution would be to manually add the plate ids of the line features to the topology, but that is a job computers are usually better at.
Regards,
Thomas
Hi Thomas,
So have you already created the topologies in your model at this stage (and you just want to resolve them) ? If so, have you created the topologies in GPlates, or are you now using pyGPlates for that?
If you’ve already got your topologies then you can just throw all your features (topological and regular features) at pygplates.resolve_topologies. I only passed the mid-ocean ridge topology and its two referenced features as an example. Any regular features not required for resolving topologies are simply ignored. So you can do something like the following (assumes you’ve created topologies in GPlates and are loading them into pyGPlates via a GPML file)…
all_features = pygplates.FeatureCollection('my_features.gpml')
resolved_topologies = []
pygplates.resolve_topologies(
all_features,
rotation_model,
resolved_topologies,
reconstruction_time,
...)
for resolved_topology in resolved_topologies:
...
…where you’d then extract the reconstructed clipped line features from each resolved_topology as shown in a previous post (I’m assuming each resolved_topology is a topological line mid-ocean ridge).
Regards,
John