From c23f3c889c6efbb66d9c5aa8df335fe863660e88 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Wed, 29 Oct 2025 14:01:41 +1100 Subject: [PATCH 01/17] An improved version of collision detection for seafloor gridding. This naturally deactivates seed points at convergent plate boundaries without having to compare velocity deltas and distances to boundaries with empirical thresholds. --- gplately/lib/reconstruct_by_topologies.py | 358 +++++++++++++++++++++- gplately/oceans.py | 4 +- 2 files changed, 359 insertions(+), 3 deletions(-) diff --git a/gplately/lib/reconstruct_by_topologies.py b/gplately/lib/reconstruct_by_topologies.py index f7c9c7ab..4d207703 100644 --- a/gplately/lib/reconstruct_by_topologies.py +++ b/gplately/lib/reconstruct_by_topologies.py @@ -842,6 +842,361 @@ def _find_resolved_topologies_containing_points(self): ) +class _ReconstructByTopologiesImplV2(_ReconstructByTopologies): + """An improved version of _ReconstructByTopologiesImpl. + + This uses a different approach to collision detection of seed points + (to determine if they get subducted going forward in time and hence deactivated). + Previously we detected if the velocity delta of a seed point exceeded a threshold (when that point + transitioned from one plate to another) - and the point had to be close enough to the boundary. + Instead, we now take the plate containing a seed point at the current time and resolve it at the next time + (ie, current time plus time step) and see if its resolved shape still contains that seed point + (reconstructed to the next time). If it does then the seed point remains active (otherwise it's deactivated). + Note that the *same* plate is resolved at the current and next times. Usually when you resolve topologies + at a different time you will get different plates (eg, there could be plate merges or splits). + So, to resolve the *same* plate, some fenagling is required to ensure that a resolved plate boundary at the + current time can also be resolved at the next time. The end result is that any boundary around a plate that + converges will naturally deactivate seed points near that boundary (without having to use thresholds, etc). + And this doesn't require plate boundaries to be labelled as converging (ie, subduction zones). + """ + + def __init__( + self, + rotation_features_or_model, + topology_features, + reconstruction_begin_time, + reconstruction_end_time, + reconstruction_time_interval, + points, + *, + point_begin_times: Union[np.ndarray, list, None] = None, + point_end_times: Union[np.ndarray, list, None] = None, + ): + + super().__init__( + reconstruction_begin_time, + reconstruction_end_time, + reconstruction_time_interval, + points, + point_begin_times, + point_end_times, + ) + + self.rotation_model = pygplates.RotationModel(rotation_features_or_model) # type: ignore + + # Turn topology data into a list of features (if not already). + self.topology_features = pygplates.FeaturesFunctionArgument( # type: ignore + topology_features + ).get_features() + + def _begin_reconstruction_impl(self): + # Activate any original points whose begin time is equal to (or older than) the newly updated current time. + # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. + self._activate_deactivate_points_using_their_begin_end_times() + + def _reconstruct_to_next_time_impl(self): + + current_time = self.get_current_time() + # Positive/negative time step means reconstructing backward/forward in time. + next_time = current_time + self.reconstruction_time_step + + curr_topological_snapshot = pygplates.TopologicalSnapshot( # type: ignore + self.topology_features, self.rotation_model, current_time + ) + curr_resolved_topologies = curr_topological_snapshot.get_resolved_topologies() + + # Some of 'curr_points' will be None, so 'curr_active_points' contains only the active (not None) + # points, and 'curr_active_points_indices' is the same length as 'curr_active_points' but indexes + # into 'curr_points' so we can quickly find which point in 'curr_points' is associated with a + # particular point in 'curr_active_points'. + curr_active_points = [] + curr_active_points_indices = [] + for point_index, curr_point in enumerate(self.curr_points): + if curr_point is None: + # Current point is not currently active, so we cannot reconstruct it to the next point. + self.next_points[point_index] = None + continue + curr_active_points.append(curr_point) + curr_active_points_indices.append(point_index) + + # For each active current point find the resolved topology containing it. + curr_active_point_resolved_topologies = ( + _ptt.utils.points_in_polygons.find_polygons( + curr_active_points, + [ + resolved_topology.get_resolved_boundary() + for resolved_topology in curr_resolved_topologies + ], + curr_resolved_topologies, + ) + ) + + # List of all unique resolved topologies that the currently active points fall within (in the current time step). + curr_unique_resolved_topologies = [] + # List of indices into 'curr_unique_resolved_topologies' for each currently active point. + active_point_resolved_topology_indices = [] + # Map of current resolved topologies to their index into 'curr_unique_resolved_topologies'. + map_curr_resolved_topology_to_index = {} + + # Iterate over the resolved topologies containing all currently active points and get a list of the unique resolved topologies. + # + # Note: A resolved topology for a currently active point can be None if it fell outside all resolved topologies. + # In this case we'll just deactivate the point. Previously we would keep it active but just not reconstruct it + # (ie, the current and next positions would be the same). However the point might end up in weird locations. + # So it's best to just remove it. And this shouldn't happen very often at all for topologies with global coverage. + curr_active_point_index = 0 + while curr_active_point_index < len(curr_active_points): + + curr_active_point_resolved_topology = curr_active_point_resolved_topologies[ + curr_active_point_index + ] + # See if current active point fell outside all current resolved topologies. + if curr_active_point_resolved_topology is None: + # Index into the active and inactive points 'self.curr_points'. + curr_point_index = curr_active_points_indices[curr_active_point_index] + + # Current point is currently active but it fell outside all resolved topologies. + # So we deactivate it. + self.next_points[curr_point_index] = None + + # Also remove evidence that the current point is active. + del curr_active_points[curr_active_point_index] + del curr_active_points_indices[curr_active_point_index] + del curr_active_point_resolved_topologies[curr_active_point_index] + + # Continue to next active point. + # + # Note: We don't increment the active point index. + # We just removed an active point which essentially does the same thing. + continue + + # If resolved topology has not been encountered yet then add it to our list of unique resolved topologies. + if ( + curr_active_point_resolved_topology + not in map_curr_resolved_topology_to_index + ): + curr_active_point_resolved_topology_index = len( + curr_unique_resolved_topologies + ) + curr_unique_resolved_topologies.append( + curr_active_point_resolved_topology + ) + map_curr_resolved_topology_to_index[ + curr_active_point_resolved_topology + ] = curr_active_point_resolved_topology_index + + # The index of resolved topology (containing current point) into 'curr_unique_resolved_topologies'. + active_point_resolved_topology_indices.append( + map_curr_resolved_topology_to_index[curr_active_point_resolved_topology] + ) + + # Increment to the next active point. + curr_active_point_index += 1 + + # All topology features and their referenced topological section features resolved at the next time step. + next_topology_features = [] + next_topological_section_features = [] + + # Map of each current topological section feature to its corresponding *next* topological section feature. + map_curr_to_next_topological_section_feature = {} + + # For each current resolved topology, create a corresponding *next* topological boundary feature and + # all the topological section features referenced by them. + for curr_resolved_topology in curr_unique_resolved_topologies: + + # Iterate over the boundary sub-segments of the current resolved topology and create a new topological + # section feature for each one (if its topological section feature hasn't been encountered yet). + next_topological_section_property_values = [] + for ( + curr_boundary_sub_segment + ) in curr_resolved_topology.get_boundary_sub_segments(): # type: ignore + + curr_topological_section_recon_geom = ( + curr_boundary_sub_segment.get_topological_section() + ) + # Topological section is either a ReconstructedFeatureGeometry or a ResolvedTopologicalLine. + if isinstance( + curr_topological_section_recon_geom, pygplates.ReconstructedFeatureGeometry # type: ignore + ): + curr_topological_section_feature = ( + curr_topological_section_recon_geom.get_feature() + ) + # If the current topological section feature hasn't been encountered yet then create an associated + # *next* topological section feature that is either the same as the current topological section feature + # (if it's still valid at the *next* time) or a clone of the current topological section feature with + # the valid time range modified to be valid at the *next* time. + if ( + curr_topological_section_feature + not in map_curr_to_next_topological_section_feature + ): + # If feature is not valid at 'next_time' then clone the feature and set a valid time that includes 'next_time'. + if not curr_topological_section_feature.is_valid_at_time( + next_time + ): + next_topological_section_feature = ( + curr_topological_section_feature.clone() + ) + next_topological_section_feature.set_valid_time( + pygplates.GeoTimeInstant.create_distant_past(), 0.0 # type: ignore + ) + else: + next_topological_section_feature = ( + curr_topological_section_feature + ) + + next_topological_section_features.append( + next_topological_section_feature + ) + + # Map the current to the next (topological section feature). + map_curr_to_next_topological_section_feature[ + curr_topological_section_feature + ] = next_topological_section_feature + + # Get the next topological section feature associated with the current one. + next_topological_section_feature = ( + map_curr_to_next_topological_section_feature[ + curr_topological_section_feature + ] + ) + + # Create the *next* topological section associated with the current one. + next_topological_section_property_value = pygplates.GpmlTopologicalSection.create( # type: ignore + next_topological_section_feature, + geometry_property_name=curr_topological_section_recon_geom.get_property().get_name(), + reverse_order=curr_boundary_sub_segment.was_geometry_reversed_in_topology(), + topological_geometry_type=pygplates.GpmlTopologicalPolygon, # type: ignore + ) + if next_topological_section_property_value: + next_topological_section_property_values.append( + next_topological_section_property_value + ) + + else: + # TODO: Handle topological sections that are topological lines. + pass + + # Create a topological boundary feature that is essentially the current topological boundary + # resolved to the *next* time step. + # + # Note: If the current topology is a deforming network we still create a topological closed plate boundary + # (which is normally used for rigid plates) because we only need to detect if a seed point + # reconstructed to the next time step is contained within the network's *boundary*. + # + # Note: 'valid_time' argument is not specified which means valid for all time. + # This works for us since we only need it to be valid at the *next* time step. + next_topology_feature = pygplates.Feature.create_topological_feature( # type: ignore + pygplates.FeatureType.gpml_topological_closed_plate_boundary, # type: ignore + pygplates.GpmlTopologicalPolygon( # type: ignore + next_topological_section_property_values + ), + ) + + next_topology_features.append(next_topology_feature) + + # Resolved the topologies to the *next* time step. + next_topological_snapshot = pygplates.TopologicalSnapshot( # type: ignore + next_topology_features + next_topological_section_features, + self.rotation_model, + next_time, + ) + next_resolved_topologies = next_topological_snapshot.get_resolved_topologies() + + # Map the next topology features to their indices into 'next_topology_features'. + # + # This will also be the indices into 'curr_unique_resolved_topologies' and 'next_unique_resolved_topologies'. + map_next_topology_feature_to_index = { + next_topology_feature: index + for index, next_topology_feature in enumerate(next_topology_features) + } + # + # List of all topologies (that currently active points fall within) but resolved to the *next* time step. + # + # Note: This is a one-to-one mapping with 'curr_unique_resolved_topologies'. + next_unique_resolved_topologies = [None] * len(curr_unique_resolved_topologies) + # + # Store each next resolved topology in 'next_unique_resolved_topologies'. + # + # Each resolved topology in 'curr_unique_resolved_topologies' will be associated with + # the same index into 'next_unique_resolved_topologies'. + for next_resolved_topology in next_resolved_topologies: + next_topology_feature_index = map_next_topology_feature_to_index.get( + next_resolved_topology.get_feature() + ) + # It's possible that the next topology could not get resolved at the next time step. + # This generally shouldn't happen though. + if next_topology_feature_index is not None: + next_unique_resolved_topologies[next_topology_feature_index] = ( + next_resolved_topology + ) + + # Iterate over all currently active points to reconstruct them to the next time step and + # see if they've been consumed by a convergent plate boundary. + for curr_active_point_index, curr_active_point in enumerate(curr_active_points): + + # Index into the active and inactive points 'self.curr_points'. + curr_point_index = curr_active_points_indices[curr_active_point_index] + + # Index into 'curr_unique_resolved_topologies' and 'next_unique_resolved_topologies' for the currently active point. + unique_resolved_topologies_index = active_point_resolved_topology_indices[ + curr_active_point_index + ] + + # Resolved topology containing the currently active point at the current time step. + curr_active_point_resolved_topology = curr_unique_resolved_topologies[ + unique_resolved_topologies_index + ] + + # Topology containing the currently active point at the current but resolved to the next time step. + next_active_point_resolved_topology = next_unique_resolved_topologies[ + unique_resolved_topologies_index + ] + if next_active_point_resolved_topology is None: + # The topology (containing the currently active point) could not be resolved to the next time step. + # So we deactivate it. This generally shouldn't happen though. + self.next_points[curr_point_index] = None + continue + + # Reconstruct the currently active point from its position at current time to its position at the next time step. + next_active_point = curr_active_point_resolved_topology.reconstruct_point( + curr_active_point, next_time + ) + + # See if the location (of the current active point) reconstructed to the *next* time step + # is inside the current topology resolved to the *next* time step. + # + # If it is then is was not subducting going forward in time (or consumed by a mid-ocean ridge going backward in time). + if next_active_point_resolved_topology.get_point_location( + next_active_point + ).not_located_in_resolved_topology(): + self.next_points[curr_point_index] = None + continue + + # The current active point was not consumed by a plate boundary. + self.next_points[curr_point_index] = next_active_point + + # + # Set up for next loop iteration. + # + # Rotate previous, current and next point arrays. + # The new previous will be the old current. + # The new current will be the old next. + # The new next will be the old previous (but values are ignored and overridden in next time step; just re-using its memory). + self.prev_points, self.curr_points, self.next_points = ( + self.curr_points, + self.next_points, + self.prev_points, + ) + # + # Move the current time to the next time. + self.current_time_index += 1 + current_time = self.get_current_time() + + # Activate any original points whose begin time is equal to (or older than) the newly updated current time. + # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. + self._activate_deactivate_points_using_their_begin_end_times() + + class _ReconstructByTopologicalModelImpl(_ReconstructByTopologies): """Similar to ``_ReconstructByTopologiesImpl`` except uses the ``pygplates.TopologicalModel`` class to reconstruct seed points. @@ -990,7 +1345,8 @@ def _reconstruct_to_next_time_impl(self): reconstructed_time_span = self.topological_model.reconstruct_geometry( current_active_points, initial_time=current_time, - youngest_time=next_time, + oldest_time=max(current_time, next_time), + youngest_time=min(current_time, next_time), time_increment=self.reconstruction_time_interval, # must be positive deactivate_points=self.detect_collisions, ) diff --git a/gplately/oceans.py b/gplately/oceans.py index 91cbb317..724a32d4 100644 --- a/gplately/oceans.py +++ b/gplately/oceans.py @@ -33,6 +33,7 @@ _ContinentCollision, _DefaultCollision, _ReconstructByTopologiesImpl, + _ReconstructByTopologiesImplV2, _ReconstructByTopologicalModelImpl, ) from .lib.reconstruct_continents import ReconstructContinents @@ -1408,7 +1409,7 @@ def _build_and_reconstruct_ocean_seed_points( detect_collisions=collision_spec, ) else: - topology_reconstruction = _ReconstructByTopologiesImpl( + topology_reconstruction = _ReconstructByTopologiesImplV2( self.plate_reconstruction.rotation_model, self.plate_reconstruction.topology_features, from_time, @@ -1416,7 +1417,6 @@ def _build_and_reconstruct_ocean_seed_points( self._ridge_time_step, points, point_begin_times=appearance_times, - detect_collisions=collision_spec, ) # Initialise the reconstruction. From 892553d8bf9f46bed3967545dee1bce7f5582179 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Wed, 29 Oct 2025 22:43:13 +1100 Subject: [PATCH 02/17] Support topological lines in new seafloor gridding collision detection. Topological lines can be sections of topological boundaries/networks. New collision detection was recently added in 521b6fe. --- gplately/lib/reconstruct_by_topologies.py | 231 ++++++++++++++++------ 1 file changed, 170 insertions(+), 61 deletions(-) diff --git a/gplately/lib/reconstruct_by_topologies.py b/gplately/lib/reconstruct_by_topologies.py index 4d207703..c645f79b 100644 --- a/gplately/lib/reconstruct_by_topologies.py +++ b/gplately/lib/reconstruct_by_topologies.py @@ -993,8 +993,10 @@ def _reconstruct_to_next_time_impl(self): # Increment to the next active point. curr_active_point_index += 1 - # All topology features and their referenced topological section features resolved at the next time step. - next_topology_features = [] + # All topology *boundary* features and their referenced topological *section* features resolved at the next time step. + # + # Note: Topological *section* features can include topological lines. + next_topological_boundary_features = [] next_topological_section_features = [] # Map of each current topological section feature to its corresponding *next* topological section feature. @@ -1006,78 +1008,184 @@ def _reconstruct_to_next_time_impl(self): # Iterate over the boundary sub-segments of the current resolved topology and create a new topological # section feature for each one (if its topological section feature hasn't been encountered yet). - next_topological_section_property_values = [] + next_boundary_section_property_values = [] for ( curr_boundary_sub_segment ) in curr_resolved_topology.get_boundary_sub_segments(): # type: ignore - curr_topological_section_recon_geom = ( - curr_boundary_sub_segment.get_topological_section() - ) - # Topological section is either a ReconstructedFeatureGeometry or a ResolvedTopologicalLine. - if isinstance( - curr_topological_section_recon_geom, pygplates.ReconstructedFeatureGeometry # type: ignore + curr_boundary_section_feature = curr_boundary_sub_segment.get_feature() + # If the current boundary section feature hasn't been encountered yet then create an associated + # *next* boundary section feature that is either the same as the current boundary section feature + # (if it's still valid at the *next* time) or a clone of the current boundary section feature with + # the valid time range modified to be valid at the *next* time. + if ( + curr_boundary_section_feature + not in map_curr_to_next_topological_section_feature ): - curr_topological_section_feature = ( - curr_topological_section_recon_geom.get_feature() + curr_boundary_section_recon_geom = ( + curr_boundary_sub_segment.get_topological_section() ) - # If the current topological section feature hasn't been encountered yet then create an associated - # *next* topological section feature that is either the same as the current topological section feature - # (if it's still valid at the *next* time) or a clone of the current topological section feature with - # the valid time range modified to be valid at the *next* time. - if ( - curr_topological_section_feature - not in map_curr_to_next_topological_section_feature + # Boundary section is either a ReconstructedFeatureGeometry or a ResolvedTopologicalLine. + # + # If it's a ResolvedTopologicalLine then we need to iterate through its sub-segments and create a + # new line section for each one, and finally create a new topological line feature from them. + if isinstance( + curr_boundary_section_recon_geom, pygplates.ResolvedTopologicalLine # type: ignore ): + + # Iterate over the sub-segments of the current resolved topological line and create a new line + # section feature for each one (if its topological section feature hasn't been encountered yet). + next_line_section_property_values = [] + for ( + curr_line_sub_segment + ) in curr_boundary_section_recon_geom.get_line_sub_segments(): # type: ignore + curr_line_section_feature = ( + curr_line_sub_segment.get_feature() + ) + # If the current line section feature hasn't been encountered yet then create an associated + # *next* line section feature that is either the same as the current line section feature + # (if it's still valid at the *next* time) or a clone of the current line section feature with + # the valid time range modified to be valid at the *next* time. + if ( + curr_line_section_feature + not in map_curr_to_next_topological_section_feature + ): + # Line section must be a ReconstructedFeatureGeometry (ie, can't be ResolvedTopologicalLine) + # since a topological line cannot have sections that are also topological lines. + + # If feature is not valid at 'next_time' then clone the feature and set a valid time that includes 'next_time'. + if not curr_line_section_feature.is_valid_at_time( + next_time + ): + next_line_section_feature = ( + curr_line_section_feature.clone() + ) + next_line_section_feature.set_valid_time( + pygplates.GeoTimeInstant.create_distant_past(), 0.0 # type: ignore + ) + else: + next_line_section_feature = ( + curr_line_section_feature + ) + + next_topological_section_features.append( + next_line_section_feature + ) + + # Map the current to the next (line section feature). + map_curr_to_next_topological_section_feature[ + curr_line_section_feature + ] = next_line_section_feature + + # Get the next line section feature associated with the current one. + next_line_section_feature = ( + map_curr_to_next_topological_section_feature[ + curr_line_section_feature + ] + ) + + # The property name of the geometry referenced by the next line section. + next_line_section_geometry_property_name = ( + curr_line_sub_segment.get_topological_section() + .get_property() + .get_name() + ) + # Whether to reverse the geometry referenced by the next line section. + next_line_section_reverse_order = ( + curr_line_sub_segment.was_geometry_reversed_in_topology() + ) + + # Create the *next* line section associated with the current one. + next_line_section_property_value = pygplates.GpmlTopologicalSection.create( # type: ignore + next_line_section_feature, + geometry_property_name=next_line_section_geometry_property_name, + reverse_order=next_line_section_reverse_order, + topological_geometry_type=pygplates.GpmlTopologicalLine, # type: ignore + ) + if next_line_section_property_value: + next_line_section_property_values.append( + next_line_section_property_value + ) + + # Create a topological line feature that enables the current topological line + # to be resolved to the *next* time step. + # + # Note: 'valid_time' argument is not specified which means valid for all time. + # This works for us since we only need it to be valid at the *next* time step. + next_line_feature = pygplates.Feature( # type: ignore + pygplates.FeatureType.gpml_unclassified_feature # type: ignore + ) + next_line_feature.set_topological_geometry( + pygplates.GpmlTopologicalLine( # type: ignore + next_line_section_property_values + ), + # Use the same geometry property name as the current boundary sub-segment so + # that it can be found by the topological polygon(s) that will reference it... + curr_boundary_sub_segment.get_topological_section() + .get_property() + .get_name(), + ) + + next_boundary_section_feature = next_line_feature + + else: # boundary section is a ReconstructedFeatureGeometry ... + # If feature is not valid at 'next_time' then clone the feature and set a valid time that includes 'next_time'. - if not curr_topological_section_feature.is_valid_at_time( + if not curr_boundary_section_feature.is_valid_at_time( next_time ): - next_topological_section_feature = ( - curr_topological_section_feature.clone() + next_boundary_section_feature = ( + curr_boundary_section_feature.clone() ) - next_topological_section_feature.set_valid_time( + next_boundary_section_feature.set_valid_time( pygplates.GeoTimeInstant.create_distant_past(), 0.0 # type: ignore ) else: - next_topological_section_feature = ( - curr_topological_section_feature + next_boundary_section_feature = ( + curr_boundary_section_feature ) - next_topological_section_features.append( - next_topological_section_feature - ) + next_topological_section_features.append( + next_boundary_section_feature + ) - # Map the current to the next (topological section feature). - map_curr_to_next_topological_section_feature[ - curr_topological_section_feature - ] = next_topological_section_feature + # Map the current to the next (boundary section feature). + map_curr_to_next_topological_section_feature[ + curr_boundary_section_feature + ] = next_boundary_section_feature - # Get the next topological section feature associated with the current one. - next_topological_section_feature = ( - map_curr_to_next_topological_section_feature[ - curr_topological_section_feature - ] - ) + # Get the next boundary section feature associated with the current one. + next_boundary_section_feature = ( + map_curr_to_next_topological_section_feature[ + curr_boundary_section_feature + ] + ) - # Create the *next* topological section associated with the current one. - next_topological_section_property_value = pygplates.GpmlTopologicalSection.create( # type: ignore - next_topological_section_feature, - geometry_property_name=curr_topological_section_recon_geom.get_property().get_name(), - reverse_order=curr_boundary_sub_segment.was_geometry_reversed_in_topology(), - topological_geometry_type=pygplates.GpmlTopologicalPolygon, # type: ignore - ) - if next_topological_section_property_value: - next_topological_section_property_values.append( - next_topological_section_property_value - ) + # The property name of the geometry referenced by the next boundary section. + next_boundary_section_geometry_property_name = ( + curr_boundary_sub_segment.get_topological_section() + .get_property() + .get_name() + ) + # Whether to reverse the geometry referenced by the next boundary section. + next_boundary_section_reverse_order = ( + curr_boundary_sub_segment.was_geometry_reversed_in_topology() + ) - else: - # TODO: Handle topological sections that are topological lines. - pass + # Create the *next* boundary section associated with the current one. + next_boundary_section_property_value = pygplates.GpmlTopologicalSection.create( # type: ignore + next_boundary_section_feature, + geometry_property_name=next_boundary_section_geometry_property_name, + reverse_order=next_boundary_section_reverse_order, + topological_geometry_type=pygplates.GpmlTopologicalPolygon, # type: ignore + ) + if next_boundary_section_property_value: + next_boundary_section_property_values.append( + next_boundary_section_property_value + ) - # Create a topological boundary feature that is essentially the current topological boundary - # resolved to the *next* time step. + # Create a topological boundary feature that enables the current topological boundary + # to be resolved to the *next* time step. # # Note: If the current topology is a deforming network we still create a topological closed plate boundary # (which is normally used for rigid plates) because we only need to detect if a seed point @@ -1085,29 +1193,30 @@ def _reconstruct_to_next_time_impl(self): # # Note: 'valid_time' argument is not specified which means valid for all time. # This works for us since we only need it to be valid at the *next* time step. - next_topology_feature = pygplates.Feature.create_topological_feature( # type: ignore + next_boundary_feature = pygplates.Feature.create_topological_feature( # type: ignore pygplates.FeatureType.gpml_topological_closed_plate_boundary, # type: ignore pygplates.GpmlTopologicalPolygon( # type: ignore - next_topological_section_property_values + next_boundary_section_property_values ), ) - - next_topology_features.append(next_topology_feature) + next_topological_boundary_features.append(next_boundary_feature) # Resolved the topologies to the *next* time step. next_topological_snapshot = pygplates.TopologicalSnapshot( # type: ignore - next_topology_features + next_topological_section_features, + next_topological_section_features + next_topological_boundary_features, self.rotation_model, next_time, ) next_resolved_topologies = next_topological_snapshot.get_resolved_topologies() - # Map the next topology features to their indices into 'next_topology_features'. + # Map the next topology *boundary* features to their indices into 'next_topological_boundary_features'. # # This will also be the indices into 'curr_unique_resolved_topologies' and 'next_unique_resolved_topologies'. map_next_topology_feature_to_index = { next_topology_feature: index - for index, next_topology_feature in enumerate(next_topology_features) + for index, next_topology_feature in enumerate( + next_topological_boundary_features + ) } # # List of all topologies (that currently active points fall within) but resolved to the *next* time step. @@ -1165,7 +1274,7 @@ def _reconstruct_to_next_time_impl(self): # See if the location (of the current active point) reconstructed to the *next* time step # is inside the current topology resolved to the *next* time step. # - # If it is then is was not subducting going forward in time (or consumed by a mid-ocean ridge going backward in time). + # If it is outside then it is subducting going forward in time (or consumed by a mid-ocean ridge going backward in time). if next_active_point_resolved_topology.get_point_location( next_active_point ).not_located_in_resolved_topology(): From e14a7c234733cea600da7613d5063b0e00882cad Mon Sep 17 00:00:00 2001 From: John Cannon Date: Tue, 11 Nov 2025 12:54:24 +1100 Subject: [PATCH 03/17] Move creation of *next* topology features into its own class. --- gplately/lib/reconstruct_by_topologies.py | 321 ++++++++++++---------- 1 file changed, 170 insertions(+), 151 deletions(-) diff --git a/gplately/lib/reconstruct_by_topologies.py b/gplately/lib/reconstruct_by_topologies.py index c645f79b..3d069dae 100644 --- a/gplately/lib/reconstruct_by_topologies.py +++ b/gplately/lib/reconstruct_by_topologies.py @@ -931,19 +931,15 @@ def _reconstruct_to_next_time_impl(self): ) ) - # List of all unique resolved topologies that the currently active points fall within (in the current time step). - curr_unique_resolved_topologies = [] - # List of indices into 'curr_unique_resolved_topologies' for each currently active point. - active_point_resolved_topology_indices = [] - # Map of current resolved topologies to their index into 'curr_unique_resolved_topologies'. - map_curr_resolved_topology_to_index = {} + # Map of current resolved topologies to the active points contained within them (their active point indices). + map_curr_resolved_topology_to_active_point_indices = {} # Iterate over the resolved topologies containing all currently active points and get a list of the unique resolved topologies. # # Note: A resolved topology for a currently active point can be None if it fell outside all resolved topologies. # In this case we'll just deactivate the point. Previously we would keep it active but just not reconstruct it # (ie, the current and next positions would be the same). However the point might end up in weird locations. - # So it's best to just remove it. And this shouldn't happen very often at all for topologies with global coverage. + # So it's best to just remove it. And this shouldn't happen very often at all (for topologies with global coverage). curr_active_point_index = 0 while curr_active_point_index < len(curr_active_points): @@ -970,41 +966,162 @@ def _reconstruct_to_next_time_impl(self): # We just removed an active point which essentially does the same thing. continue - # If resolved topology has not been encountered yet then add it to our list of unique resolved topologies. + # If resolved topology has not been encountered yet then give it an empty list of active point indices. if ( curr_active_point_resolved_topology - not in map_curr_resolved_topology_to_index + not in map_curr_resolved_topology_to_active_point_indices ): - curr_active_point_resolved_topology_index = len( - curr_unique_resolved_topologies - ) - curr_unique_resolved_topologies.append( - curr_active_point_resolved_topology - ) - map_curr_resolved_topology_to_index[ + map_curr_resolved_topology_to_active_point_indices[ curr_active_point_resolved_topology - ] = curr_active_point_resolved_topology_index + ] = [] - # The index of resolved topology (containing current point) into 'curr_unique_resolved_topologies'. - active_point_resolved_topology_indices.append( - map_curr_resolved_topology_to_index[curr_active_point_resolved_topology] - ) + # Add the index of the currently active point to the resolved topology containing it. + map_curr_resolved_topology_to_active_point_indices[ + curr_active_point_resolved_topology + ].append(curr_active_point_index) # Increment to the next active point. curr_active_point_index += 1 - # All topology *boundary* features and their referenced topological *section* features resolved at the next time step. - # - # Note: Topological *section* features can include topological lines. - next_topological_boundary_features = [] - next_topological_section_features = [] - - # Map of each current topological section feature to its corresponding *next* topological section feature. - map_curr_to_next_topological_section_feature = {} + # The current resolved topologies that contain active points. + curr_active_resolved_topologies = list( + map_curr_resolved_topology_to_active_point_indices.keys() + ) + # The *next* topology features and their topological section features. + next_topological_features = self._NextTopologicalFeatures(next_time) # For each current resolved topology, create a corresponding *next* topological boundary feature and # all the topological section features referenced by them. - for curr_resolved_topology in curr_unique_resolved_topologies: + for curr_resolved_topology in curr_active_resolved_topologies: + next_topological_features.add_current_topology(curr_resolved_topology) + + # Resolve the topologies to the *next* time step. + next_topological_snapshot = pygplates.TopologicalSnapshot( # type: ignore + next_topological_features.next_topological_section_features + + next_topological_features.next_topological_boundary_features, + self.rotation_model, + next_time, + ) + + # First, map the *next* *boundary* features to their indices into 'next_topological_boundary_features'. + # + # This will also be the indices into 'curr_active_resolved_topologies' (due to order of insertion). + map_next_topology_boundary_feature_to_index = { + next_topology_feature: index + for index, next_topology_feature in enumerate( + next_topological_features.next_topological_boundary_features + ) + } + # + # Then, map each current resolved topology to its next resolved topology. + map_curr_to_next_resolved_topologies = {} + for ( + next_resolved_topology + ) in next_topological_snapshot.get_resolved_topologies(): + index = map_next_topology_boundary_feature_to_index.get( + next_resolved_topology.get_feature() + ) + # It's possible that the next topology could not get resolved at the next time step. + # This generally shouldn't happen though. + if index is not None: + curr_resolved_topology = curr_active_resolved_topologies[index] + map_curr_to_next_resolved_topologies[curr_resolved_topology] = ( + next_resolved_topology + ) + + # Iterate over all currently active points to reconstruct them to the next time step and + # see if they've been consumed by a convergent plate boundary. + # + # Note: Active points are grouped by the current resolved topology containing them. + for ( + curr_resolved_topology, + curr_resolved_topology_active_point_indices, + ) in map_curr_resolved_topology_to_active_point_indices.items(): + + # The next resolved topology associated with the current resolved topology. + next_resolved_topology = map_curr_to_next_resolved_topologies.get( + curr_resolved_topology + ) + if next_resolved_topology is None: + # The current topology could not be resolved to the next time step. + # So we deactivate all points that it contains. This generally shouldn't happen though. + for ( + curr_active_point_index + ) in curr_resolved_topology_active_point_indices: + curr_point_index = curr_active_points_indices[ + curr_active_point_index + ] + self.next_points[curr_point_index] = None + continue + + # Iterate over the currently active points contained by the current resolved topology. + for curr_active_point_index in curr_resolved_topology_active_point_indices: + + # Reconstruct the currently active point from its position at current time to its position at the next time step. + curr_active_point = curr_active_points[curr_active_point_index] + next_active_point = curr_resolved_topology.reconstruct_point( + curr_active_point, next_time + ) + + # Index into the active and inactive points 'self.curr_points'. + curr_point_index = curr_active_points_indices[curr_active_point_index] + + # See if the location (of the current active point) reconstructed to the *next* time step + # is inside the current topology resolved to the *next* time step. + # + # If it is outside then it is subducting going forward in time (or consumed by a mid-ocean ridge going backward in time). + # It doesn't necessarily have to happen at a subduction zone (or mid-ocean ridge). It can be any part of a plate boundary + # that is convergent forward in time (or divergent forward in time; which is convergent backward in time). + if next_resolved_topology.get_point_location( + next_active_point + ).not_located_in_resolved_topology(): + self.next_points[curr_point_index] = None + continue + + # The current active point was not consumed by a plate boundary. + self.next_points[curr_point_index] = next_active_point + + # + # Set up for next loop iteration. + # + # Rotate previous, current and next point arrays. + # The new previous will be the old current. + # The new current will be the old next. + # The new next will be the old previous (but values are ignored and overridden in next time step; just re-using its memory). + self.prev_points, self.curr_points, self.next_points = ( + self.curr_points, + self.next_points, + self.prev_points, + ) + # + # Move the current time to the next time. + self.current_time_index += 1 + current_time = self.get_current_time() + + # Activate any original points whose begin time is equal to (or older than) the newly updated current time. + # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. + self._activate_deactivate_points_using_their_begin_end_times() + + class _NextTopologicalFeatures(object): + """The *next* topology boundary features and their topological section features. + + These are essentially the same as the *current* topology features (and their topological section features) but + with their valid time periods modified (if needed) to include the *next* time. + """ + + def __init__(self, next_time): + + self.next_time = next_time + # All topology *boundary* features and their referenced topological *section* features ready to be resolved at the next time. + # + # Note: Topological *section* features can include topological lines. + self.next_topological_boundary_features = [] + self.next_topological_section_features = [] + + # Map of each current topological section feature to its corresponding *next* topological section feature. + self._map_curr_to_next_topological_section_feature = {} + + def add_current_topology(self, curr_resolved_topology): # Iterate over the boundary sub-segments of the current resolved topology and create a new topological # section feature for each one (if its topological section feature hasn't been encountered yet). @@ -1020,7 +1137,7 @@ def _reconstruct_to_next_time_impl(self): # the valid time range modified to be valid at the *next* time. if ( curr_boundary_section_feature - not in map_curr_to_next_topological_section_feature + not in self._map_curr_to_next_topological_section_feature ): curr_boundary_section_recon_geom = ( curr_boundary_sub_segment.get_topological_section() @@ -1048,14 +1165,14 @@ def _reconstruct_to_next_time_impl(self): # the valid time range modified to be valid at the *next* time. if ( curr_line_section_feature - not in map_curr_to_next_topological_section_feature + not in self._map_curr_to_next_topological_section_feature ): # Line section must be a ReconstructedFeatureGeometry (ie, can't be ResolvedTopologicalLine) # since a topological line cannot have sections that are also topological lines. # If feature is not valid at 'next_time' then clone the feature and set a valid time that includes 'next_time'. if not curr_line_section_feature.is_valid_at_time( - next_time + self.next_time ): next_line_section_feature = ( curr_line_section_feature.clone() @@ -1068,18 +1185,18 @@ def _reconstruct_to_next_time_impl(self): curr_line_section_feature ) - next_topological_section_features.append( + self.next_topological_section_features.append( next_line_section_feature ) # Map the current to the next (line section feature). - map_curr_to_next_topological_section_feature[ + self._map_curr_to_next_topological_section_feature[ curr_line_section_feature ] = next_line_section_feature # Get the next line section feature associated with the current one. next_line_section_feature = ( - map_curr_to_next_topological_section_feature[ + self._map_curr_to_next_topological_section_feature[ curr_line_section_feature ] ) @@ -1102,10 +1219,13 @@ def _reconstruct_to_next_time_impl(self): reverse_order=next_line_section_reverse_order, topological_geometry_type=pygplates.GpmlTopologicalLine, # type: ignore ) - if next_line_section_property_value: - next_line_section_property_values.append( - next_line_section_property_value + if not next_line_section_property_value: + raise RuntimeError( + "Unable to create topological section for a topological line." ) + next_line_section_property_values.append( + next_line_section_property_value + ) # Create a topological line feature that enables the current topological line # to be resolved to the *next* time step. @@ -1132,7 +1252,7 @@ def _reconstruct_to_next_time_impl(self): # If feature is not valid at 'next_time' then clone the feature and set a valid time that includes 'next_time'. if not curr_boundary_section_feature.is_valid_at_time( - next_time + self.next_time ): next_boundary_section_feature = ( curr_boundary_section_feature.clone() @@ -1145,18 +1265,18 @@ def _reconstruct_to_next_time_impl(self): curr_boundary_section_feature ) - next_topological_section_features.append( + self.next_topological_section_features.append( next_boundary_section_feature ) # Map the current to the next (boundary section feature). - map_curr_to_next_topological_section_feature[ + self._map_curr_to_next_topological_section_feature[ curr_boundary_section_feature ] = next_boundary_section_feature # Get the next boundary section feature associated with the current one. next_boundary_section_feature = ( - map_curr_to_next_topological_section_feature[ + self._map_curr_to_next_topological_section_feature[ curr_boundary_section_feature ] ) @@ -1179,10 +1299,13 @@ def _reconstruct_to_next_time_impl(self): reverse_order=next_boundary_section_reverse_order, topological_geometry_type=pygplates.GpmlTopologicalPolygon, # type: ignore ) - if next_boundary_section_property_value: - next_boundary_section_property_values.append( - next_boundary_section_property_value + if not next_boundary_section_property_value: + raise RuntimeError( + "Unable to create topological section for a topological polygon." ) + next_boundary_section_property_values.append( + next_boundary_section_property_value + ) # Create a topological boundary feature that enables the current topological boundary # to be resolved to the *next* time step. @@ -1199,111 +1322,7 @@ def _reconstruct_to_next_time_impl(self): next_boundary_section_property_values ), ) - next_topological_boundary_features.append(next_boundary_feature) - - # Resolved the topologies to the *next* time step. - next_topological_snapshot = pygplates.TopologicalSnapshot( # type: ignore - next_topological_section_features + next_topological_boundary_features, - self.rotation_model, - next_time, - ) - next_resolved_topologies = next_topological_snapshot.get_resolved_topologies() - - # Map the next topology *boundary* features to their indices into 'next_topological_boundary_features'. - # - # This will also be the indices into 'curr_unique_resolved_topologies' and 'next_unique_resolved_topologies'. - map_next_topology_feature_to_index = { - next_topology_feature: index - for index, next_topology_feature in enumerate( - next_topological_boundary_features - ) - } - # - # List of all topologies (that currently active points fall within) but resolved to the *next* time step. - # - # Note: This is a one-to-one mapping with 'curr_unique_resolved_topologies'. - next_unique_resolved_topologies = [None] * len(curr_unique_resolved_topologies) - # - # Store each next resolved topology in 'next_unique_resolved_topologies'. - # - # Each resolved topology in 'curr_unique_resolved_topologies' will be associated with - # the same index into 'next_unique_resolved_topologies'. - for next_resolved_topology in next_resolved_topologies: - next_topology_feature_index = map_next_topology_feature_to_index.get( - next_resolved_topology.get_feature() - ) - # It's possible that the next topology could not get resolved at the next time step. - # This generally shouldn't happen though. - if next_topology_feature_index is not None: - next_unique_resolved_topologies[next_topology_feature_index] = ( - next_resolved_topology - ) - - # Iterate over all currently active points to reconstruct them to the next time step and - # see if they've been consumed by a convergent plate boundary. - for curr_active_point_index, curr_active_point in enumerate(curr_active_points): - - # Index into the active and inactive points 'self.curr_points'. - curr_point_index = curr_active_points_indices[curr_active_point_index] - - # Index into 'curr_unique_resolved_topologies' and 'next_unique_resolved_topologies' for the currently active point. - unique_resolved_topologies_index = active_point_resolved_topology_indices[ - curr_active_point_index - ] - - # Resolved topology containing the currently active point at the current time step. - curr_active_point_resolved_topology = curr_unique_resolved_topologies[ - unique_resolved_topologies_index - ] - - # Topology containing the currently active point at the current but resolved to the next time step. - next_active_point_resolved_topology = next_unique_resolved_topologies[ - unique_resolved_topologies_index - ] - if next_active_point_resolved_topology is None: - # The topology (containing the currently active point) could not be resolved to the next time step. - # So we deactivate it. This generally shouldn't happen though. - self.next_points[curr_point_index] = None - continue - - # Reconstruct the currently active point from its position at current time to its position at the next time step. - next_active_point = curr_active_point_resolved_topology.reconstruct_point( - curr_active_point, next_time - ) - - # See if the location (of the current active point) reconstructed to the *next* time step - # is inside the current topology resolved to the *next* time step. - # - # If it is outside then it is subducting going forward in time (or consumed by a mid-ocean ridge going backward in time). - if next_active_point_resolved_topology.get_point_location( - next_active_point - ).not_located_in_resolved_topology(): - self.next_points[curr_point_index] = None - continue - - # The current active point was not consumed by a plate boundary. - self.next_points[curr_point_index] = next_active_point - - # - # Set up for next loop iteration. - # - # Rotate previous, current and next point arrays. - # The new previous will be the old current. - # The new current will be the old next. - # The new next will be the old previous (but values are ignored and overridden in next time step; just re-using its memory). - self.prev_points, self.curr_points, self.next_points = ( - self.curr_points, - self.next_points, - self.prev_points, - ) - # - # Move the current time to the next time. - self.current_time_index += 1 - current_time = self.get_current_time() - - # Activate any original points whose begin time is equal to (or older than) the newly updated current time. - # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. - self._activate_deactivate_points_using_their_begin_end_times() + self.next_topological_boundary_features.append(next_boundary_feature) class _ReconstructByTopologicalModelImpl(_ReconstructByTopologies): From 714f1892192eb45e5ee137b3a1cc08ba73b534b6 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Fri, 14 Nov 2025 16:47:12 +1100 Subject: [PATCH 04/17] Detect consistent topoloogy boundary intersections in seafloor gridding. Resolving the same current topological boundary/network to the next time might generate different number of intersections between boundary sections. We now detect this so we can generate a topology boundary polygon using a different approach when this happens. --- gplately/lib/reconstruct_by_topologies.py | 304 +++++++++++++++++++++- 1 file changed, 302 insertions(+), 2 deletions(-) diff --git a/gplately/lib/reconstruct_by_topologies.py b/gplately/lib/reconstruct_by_topologies.py index 3d069dae..7a371b53 100644 --- a/gplately/lib/reconstruct_by_topologies.py +++ b/gplately/lib/reconstruct_by_topologies.py @@ -1054,6 +1054,14 @@ def _reconstruct_to_next_time_impl(self): self.next_points[curr_point_index] = None continue + # Get the boundary polygon of the next resolved topology. + # + # If the next resolved topology is consistent with the current resolved topology then we're fine, + # otherwise we need to replace the boundary of the next resolved topology with something more consistent. + next_resolved_topology_boundary = self._NextTopologicalBoundary( + curr_resolved_topology, next_resolved_topology + ).get_next_resolved_topology_boundary() + # Iterate over the currently active points contained by the current resolved topology. for curr_active_point_index in curr_resolved_topology_active_point_indices: @@ -1072,9 +1080,9 @@ def _reconstruct_to_next_time_impl(self): # If it is outside then it is subducting going forward in time (or consumed by a mid-ocean ridge going backward in time). # It doesn't necessarily have to happen at a subduction zone (or mid-ocean ridge). It can be any part of a plate boundary # that is convergent forward in time (or divergent forward in time; which is convergent backward in time). - if next_resolved_topology.get_point_location( + if not next_resolved_topology_boundary.is_point_in_polygon( next_active_point - ).not_located_in_resolved_topology(): + ): self.next_points[curr_point_index] = None continue @@ -1324,6 +1332,298 @@ def add_current_topology(self, curr_resolved_topology): ) self.next_topological_boundary_features.append(next_boundary_feature) + class _NextTopologicalBoundary(object): + """The next resolved topology boundary. + + If the next resolved topology is NOT consistent with the current resolved topology then we need to replace + the boundary of the next topology with the boundary of the current topology moved to the next time. + This involves taking the *clipped* boundary sub-segments of the current topology, moving them to the next time, and + joining them together. This is not as good as actually resolving the *unclipped* boundary sections of the + current topology at the next time, which would produce a more accurate plate boundary for the next resolved topology. + But this case shouldn't happen very often. + + The next resolved topology is not consistent with the current resolved topology if its boundary sections don't have + the same number of intersections with neighbouring sections. For example, a boundary section of the current resolved topology + intersects its neighbour once (as expected) but the same boundary section of the next resolved topology (which is really + just the same topology resolved to the next time) intersects the same neighbour twice. In this case, the next resolved topology + would likely have an unexpectedly different plate boundary shape than the current resolved topology, which might deactivate + a bunch of points that it shouldn't. This happens because pyGPlates only considers the first intersection (if there are two or more). + Usually the topological model is built such that each neighbouring section only intersects once. However, when you resolve + that same topology at a different time (even if it's just 1 Myr different) then the intersections may not be what the + model builder intended (especially if the current topology's time period does not include the *next* time). + """ + + # Distance from a polyline-polyline intersection such that if we clipped that much distance off the polylines + # around that intersection then the clipped polylines are unlikely to intersect (unless they're almost parallel). + # + # We are a little conversative on this (ie, a little big) to ensure we don't get too close. + INTERSECTION_THRESHOLD_RADIANS = 1e-4 # approx 600 metres + + # Maximum number of intersections to test between two intersecting polylines. + # + # If two polylines reach the max number of intersections then we've likely detected a partial overlap + # between the two polylines (where a segment from each polyline are on top of each other and + # hence essentially have an infinite number of intersections). + MAX_NUM_INTERSECTIONS = 5 + + def __init__(self, curr_resolved_topology, next_resolved_topology): + + self.curr_resolved_topology = curr_resolved_topology + self.next_resolved_topology = next_resolved_topology + + def get_next_resolved_topology_boundary(self): + + if self._are_current_and_next_resolved_topologies_consistent(): + return self.next_resolved_topology.get_resolved_boundary() + + # TODO: Replace this with taking the *clipped* boundary sub-segments of the current topology, + # moving them to the next time, and joining them together. + return self.next_resolved_topology.get_resolved_boundary() + + def _are_current_and_next_resolved_topologies_consistent(self): + + curr_boundary_sub_segments = ( + self.curr_resolved_topology.get_boundary_sub_segments() + ) + next_boundary_sub_segments = ( + self.next_resolved_topology.get_boundary_sub_segments() + ) + if len(curr_boundary_sub_segments) != len(next_boundary_sub_segments): + raise RuntimeError( + "Current and next topologies have a different number of boundary sub-segments." + ) + num_boundary_sub_segments = len(curr_boundary_sub_segments) + + # If there's only 1 sub-segment then it cannot intersect with itself (but it can still be converted to a polygon). + # + # Note: If there's 2 sub-segments then they can intersect each other twice (to form a closed polygon). + # We will still process them to make sure they stay consistent (same number of intersections at current and next times). + if num_boundary_sub_segments <= 1: + return True + + # Compare number of neighbour intersections for the current and next boundary sub-segments. + for boundary_sub_segment_index in range(num_boundary_sub_segments): + curr_boundary_sub_segment = curr_boundary_sub_segments[ + boundary_sub_segment_index + ] + next_boundary_sub_segment = next_boundary_sub_segments[ + boundary_sub_segment_index + ] + + prev_boundary_sub_segment_index = boundary_sub_segment_index - 1 + if prev_boundary_sub_segment_index < 0: + prev_boundary_sub_segment_index += num_boundary_sub_segments + + # Previous neighbour sub-segment (for the current and next resolved topologies). + prev_curr_boundary_sub_segment = curr_boundary_sub_segments[ + prev_boundary_sub_segment_index + ] + prev_next_boundary_sub_segment = next_boundary_sub_segments[ + prev_boundary_sub_segment_index + ] + + # Number of neighbour intersections of sub-segment (for the current and next resolved topologies). + curr_num_intersections = self._get_num_intersections( + prev_curr_boundary_sub_segment.get_topological_section_geometry(), + curr_boundary_sub_segment.get_topological_section_geometry(), + ) + next_num_intersections = self._get_num_intersections( + prev_next_boundary_sub_segment.get_topological_section_geometry(), + next_boundary_sub_segment.get_topological_section_geometry(), + ) + + # If the number of intersections differ then the current and next resolved topologies are inconsistent. + if curr_num_intersections != next_num_intersections: + return False + + return True + + def _get_num_intersections(self, polyline1, polyline2, num_intersections=0): + + distance_tuple = pygplates.GeometryOnSphere.distance( # type: ignore + polyline1, + polyline2, + # We're only detecting zero distance (an intersection), so as an optimisation we + # don't need to calculate an accurate non-zero distance (instead just return None)... + distance_threshold_radians=self.INTERSECTION_THRESHOLD_RADIANS, + return_closest_positions=True, + return_closest_indices=True, + ) + + if distance_tuple is not None: + ( + distance, + polyline1_intersection, + polyline2_intersection, + polyline1_segment_index, + polyline2_segment_index, + ) = distance_tuple + + if distance == 0.0: + # We found an intersection (between 'polyline1' and 'polyline2'). + num_intersections += 1 + + # If we reached the max number of intersections then we've likely detected a partial overlap + # between the two polylines (where a segment from each polyline are on top of each other and + # hence essentially have an infinite number of intersections). + # + # In this case we'll just return the maximum number of intersections. This means if both the + # polylines overlap at both the current and next times then they'll return the same number + # of intersections (max) and compare equal and hence be considered consistent. + if num_intersections == self.MAX_NUM_INTERSECTIONS: + return num_intersections + + # Split 'polyline1' into two polylines at its intersection point. + polyline1a = self._pre_split_polyline( + polyline1, polyline1_segment_index, polyline1_intersection + ) + polyline1b = self._post_split_polyline( + polyline1, polyline1_segment_index, polyline1_intersection + ) + + # Split 'polyline2' into two polylines at its intersection point. + polyline2a = self._pre_split_polyline( + polyline2, polyline2_segment_index, polyline2_intersection + ) + polyline2b = self._post_split_polyline( + polyline2, polyline2_segment_index, polyline2_intersection + ) + + # Test for intersections among the four combinations of split polylines. + # + # Note: It's possible there may be less than four split polylines if + # either (or both) polylines were intersected at one of their end points. + if polyline1a: + if polyline2a: + num_intersections = self._get_num_intersections( + polyline1a, polyline2a, num_intersections + ) + if num_intersections == self.MAX_NUM_INTERSECTIONS: + return num_intersections + if polyline2b: + num_intersections = self._get_num_intersections( + polyline1a, polyline2b, num_intersections + ) + if num_intersections == self.MAX_NUM_INTERSECTIONS: + return num_intersections + if polyline1b: + if polyline2a: + num_intersections = self._get_num_intersections( + polyline1b, polyline2a, num_intersections + ) + if num_intersections == self.MAX_NUM_INTERSECTIONS: + return num_intersections + if polyline2b: + num_intersections = self._get_num_intersections( + polyline1b, polyline2b, num_intersections + ) + if num_intersections == self.MAX_NUM_INTERSECTIONS: + return num_intersections + + return num_intersections + + def _pre_split_polyline( + self, polyline, polyline_segment_index, polyline_intersection + ): + # Get the points *before* the intersection (including the intersection). + # This will have at least two points (needed for creating a polyline). + pre_split_points = list(polyline[: polyline_segment_index + 1]) + pre_split_points.append(polyline_intersection) + + # Remove any zero length segments next to the intersected end of the polyline. + # Because we've already detected the intersection and don't want it detected in subsequent intersection tests. + while True: + pre_split_last_segment = pygplates.GreatCircleArc( # type: ignore + pre_split_points[-2], + pre_split_points[-1], + ) + if not pre_split_last_segment.is_zero_length(): + break + + # Remove the last segment (by removing last point). + del pre_split_points[-1] + if len(pre_split_points) < 2: + # All segments before the intersection are zero length. + # So there's no pre-split polyline. + return None + + # We've encountered a non-zero length segment. Let's shorten it such that it will no longer intersect the intersection point. + # Because we've already detected the intersection and don't want it detected in subsequent intersection tests. + if ( + self.INTERSECTION_THRESHOLD_RADIANS + < pre_split_last_segment.get_arc_length() + ): + # Last segment is *longer* than the intersection threshold, so shorten it by that amount. + # Replace last point with rotated point. + pre_split_points[-1] = ( + pygplates.FiniteRotation( # type: ignore + pre_split_last_segment.get_rotation_axis(), + -self.INTERSECTION_THRESHOLD_RADIANS, # negative rotates away from segment *end* point + ) + * pre_split_points[-1] + ) + else: + # Last segment is *shorter* than the intersection threshold. But we know it's not a zero length segment, + # so it should be far enough away from the intersection point. So just remove the segment (by removing last point). + del pre_split_points[-1] + if len(pre_split_points) < 2: + # There are no segments before the intersection. + # So there's no pre-split polyline. + return None + + return pygplates.PolylineOnSphere(pre_split_points) # type: ignore + + def _post_split_polyline( + self, polyline, polyline_segment_index, polyline_intersection + ): + # Get the points *after* the intersection (including the intersection). + # This will have at least two points (needed for creating a polyline). + post_split_points = [polyline_intersection] + post_split_points.extend(polyline[polyline_segment_index + 1 :]) + + # Remove any zero length segments next to the intersected start of the polyline. + # Because we've already detected the intersection and don't want it detected in subsequent intersection tests. + while True: + post_split_first_segment = pygplates.GreatCircleArc( # type: ignore + post_split_points[0], + post_split_points[1], + ) + if not post_split_first_segment.is_zero_length(): + break + + # Remove the first segment (by removing first point). + del post_split_points[0] + if len(post_split_points) < 2: + # All segments after the intersection are zero length. + # So there's no post-split polyline. + return None + + # We've encountered a non-zero length segment. Let's shorten it such that it will no longer intersect the intersection point. + # Because we've already detected the intersection and don't want it detected in subsequent intersection tests. + if ( + self.INTERSECTION_THRESHOLD_RADIANS + < post_split_first_segment.get_arc_length() + ): + # First segment is *longer* than the intersection threshold, so shorten it by that amount. + # Replace first point with rotated point. + post_split_points[0] = ( + pygplates.FiniteRotation( # type: ignore + post_split_first_segment.get_rotation_axis(), + self.INTERSECTION_THRESHOLD_RADIANS, + ) + * post_split_points[0] + ) + else: + # First segment is *shorter* than the intersection threshold. But we know it's not a zero length segment, + # so it should be far enough away from the intersection point. So just remove the segment (by removing first point). + del post_split_points[0] + if len(post_split_points) < 2: + # There are no segments after the intersection. + # So there's no post-split polyline. + return None + + return pygplates.PolylineOnSphere(post_split_points) # type: ignore + class _ReconstructByTopologicalModelImpl(_ReconstructByTopologies): """Similar to ``_ReconstructByTopologiesImpl`` except uses the ``pygplates.TopologicalModel`` class to reconstruct seed points. From 81314e0e043c88f008b2c0e2cfb9b42819aac471 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Fri, 21 Nov 2025 15:01:29 +1100 Subject: [PATCH 05/17] Return different resolved boundary when inconsistency detected. Instead of returning the boundary of the next resolved topology, return the boundary sub-segments of the current resolved topology moved to the next time and joined together. --- gplately/lib/reconstruct_by_topologies.py | 136 +++++++++++++++++++++- 1 file changed, 131 insertions(+), 5 deletions(-) diff --git a/gplately/lib/reconstruct_by_topologies.py b/gplately/lib/reconstruct_by_topologies.py index 7a371b53..3b07e49d 100644 --- a/gplately/lib/reconstruct_by_topologies.py +++ b/gplately/lib/reconstruct_by_topologies.py @@ -1059,7 +1059,9 @@ def _reconstruct_to_next_time_impl(self): # If the next resolved topology is consistent with the current resolved topology then we're fine, # otherwise we need to replace the boundary of the next resolved topology with something more consistent. next_resolved_topology_boundary = self._NextTopologicalBoundary( - curr_resolved_topology, next_resolved_topology + curr_resolved_topology, + next_resolved_topology, + self.rotation_model, ).get_next_resolved_topology_boundary() # Iterate over the currently active points contained by the current resolved topology. @@ -1366,19 +1368,143 @@ class _NextTopologicalBoundary(object): # hence essentially have an infinite number of intersections). MAX_NUM_INTERSECTIONS = 5 - def __init__(self, curr_resolved_topology, next_resolved_topology): + def __init__( + self, curr_resolved_topology, next_resolved_topology, rotation_model + ): self.curr_resolved_topology = curr_resolved_topology self.next_resolved_topology = next_resolved_topology + self.curr_time = curr_resolved_topology.get_reconstruction_time() + self.next_time = next_resolved_topology.get_reconstruction_time() + self.rotation_model = rotation_model def get_next_resolved_topology_boundary(self): + # If the current and next resolved topologies are consistent then + # just return the resolved boundary of the next resolved topology. if self._are_current_and_next_resolved_topologies_consistent(): return self.next_resolved_topology.get_resolved_boundary() - # TODO: Replace this with taking the *clipped* boundary sub-segments of the current topology, - # moving them to the next time, and joining them together. - return self.next_resolved_topology.get_resolved_boundary() + # Instead of returning the boundary of the next resolved topology, return the boundary + # sub-segments of the current resolved topology moved to the next time and joined together. + return self._get_next_boundary_from_current_boundary_sub_segments() + + def _get_next_boundary_from_current_boundary_sub_segments(self): + + # + # Instead of the resolved boundary of the *next* resolved topology we'll: + # - take the *clipped* boundary sub-segments of the *current* resolved topology, + # - move them to the next time, and + # - join them together into a polygon boundary. + # + + # Collect the current boundary sub-segment features to use for the next boundary. + next_boundary_sub_segment_features = [] + next_boundary_sub_segment_reversals = [] + for ( + curr_boundary_sub_segment + ) in self.curr_resolved_topology.get_boundary_sub_segments(): + # See if the current boundary sub-segment is from a topological *line*. + # + # If it is then we need to iterate over its sub-segments. This is because we need features that are *reconstructable*. + # Topological lines are not reconstructable (because they're topologies linking reconstructable features). + # Only the sub-segments of topological lines are reconstructable. + sub_sub_segments = curr_boundary_sub_segment.get_sub_segments() + if sub_sub_segments: + curr_boundary_sub_segment_reversal = ( + curr_boundary_sub_segment.was_geometry_reversed_in_topology() + ) + # Each sub-segment of the resolved topological *line* will contribute to the current boundary sub-segment. + # + # Note: Only a sub-section of the full resolved topological *line* will contribute + # to the boundary of the current resolved topological boundary/network. + # Hence not all the sub-segments of the resolved topological *line* will contribute + # to the boundary of the current resolved topological boundary/network. + # + # Note: We need to reverse the order of sub-sub-segments if the resolved topological *line* is reversed + # in the current resolved boundary. + next_boundary_sub_segment_features.extend( + sub_sub_segment.get_resolved_feature() + for sub_sub_segment in ( + reversed(sub_sub_segments) + if curr_boundary_sub_segment_reversal + else sub_sub_segments + ) + ) + # If a sub-sub-segment is reversed in the resolved topological *line* which is, in turn, reversed in the + # current resolved boundary then the sub-sub-segment is NOT reversed in the current resolved boundary. + next_boundary_sub_segment_reversals.extend( + curr_boundary_sub_segment_reversal + ^ sub_sub_segment.was_geometry_reversed_in_topology() + for sub_sub_segment in ( + reversed(sub_sub_segments) + if curr_boundary_sub_segment_reversal + else sub_sub_segments + ) + ) + else: + next_boundary_sub_segment_features.append( + curr_boundary_sub_segment.get_resolved_feature() + ) + next_boundary_sub_segment_reversals.append( + curr_boundary_sub_segment.was_geometry_reversed_in_topology() + ) + + # The sub-segment resolved features are already reconstructed to 'curr_time'. + # So we need to reverse reconstruct them to present day (before we can reconstruct them to 'next_time'). + # + # Note: We don't need to clone these features before modifying them because they were already essentially cloned + # when 'ResolvedTopologicalSubSegment.get_resolved_feature()' was called. + pygplates.reverse_reconstruct( # type:ignore + next_boundary_sub_segment_features, + self.rotation_model, + self.curr_time, + ) + + # The resolved sub-segment features are valid at the 'curr_time' but not necessarily at 'next_time'. + # If not valid at 'next_time' then set a valid time that includes 'next_time'. + for feature in next_boundary_sub_segment_features: + if not feature.is_valid_at_time(self.next_time): + feature.set_valid_time( + pygplates.GeoTimeInstant.create_distant_past(), 0.0 # type: ignore + ) + + # Reconstruct to 'next_time' (from present day). + next_boundary_sub_segment_reconstructed_geometries = ( + pygplates.ReconstructSnapshot( # type:ignore + next_boundary_sub_segment_features, + self.rotation_model, + self.next_time, + ).get_reconstructed_geometries( + same_order_as_reconstructable_features=True + ) + ) + + # The order and length of reconstructed geometries should match that of the features. + if len(next_boundary_sub_segment_reconstructed_geometries) != len( + next_boundary_sub_segment_features + ): + raise RuntimeError( + "Expected number of reconstructed geometries to match number of features." + ) + + # Iterate over the next boundary sub-segments and join them to form a polygon boundary. + next_boundary_polygon_points = [] + for ( + feature_index, + next_boundary_sub_segment_reconstructed_geometry, + ) in enumerate(next_boundary_sub_segment_reconstructed_geometries): + next_boundary_polygon_points.extend( + reversed( + next_boundary_sub_segment_reconstructed_geometry.get_reconstructed_geometry_points() + ) + if next_boundary_sub_segment_reversals[feature_index] + else next_boundary_sub_segment_reconstructed_geometry.get_reconstructed_geometry_points() + ) + + return pygplates.PolygonOnSphere( # type:ignore + next_boundary_polygon_points + ) def _are_current_and_next_resolved_topologies_consistent(self): From 3861c1b334a123a8f630feee1ad4325ca21a5b82 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Fri, 21 Nov 2025 20:00:58 +1100 Subject: [PATCH 06/17] Replace only inconsistent boundary sub-segments (instead of all). The consistent sub-segments (resolved at 'next' time) can be kept. --- gplately/lib/reconstruct_by_topologies.py | 346 ++++++++++++++-------- 1 file changed, 226 insertions(+), 120 deletions(-) diff --git a/gplately/lib/reconstruct_by_topologies.py b/gplately/lib/reconstruct_by_topologies.py index 3b07e49d..8d38280b 100644 --- a/gplately/lib/reconstruct_by_topologies.py +++ b/gplately/lib/reconstruct_by_topologies.py @@ -1338,11 +1338,13 @@ class _NextTopologicalBoundary(object): """The next resolved topology boundary. If the next resolved topology is NOT consistent with the current resolved topology then we need to replace - the boundary of the next topology with the boundary of the current topology moved to the next time. - This involves taking the *clipped* boundary sub-segments of the current topology, moving them to the next time, and - joining them together. This is not as good as actually resolving the *unclipped* boundary sections of the - current topology at the next time, which would produce a more accurate plate boundary for the next resolved topology. - But this case shouldn't happen very often. + replace any inconsistent boundary sub-segments of the next resolved topology with the equivalent ones in the + current resolved topology (but moved to the next time). For each inconsistent boundary sub-segment, this involves + taking the boundary sub-segment of the *current* resolved topology and moving it to the *next* time, and using that + (instead of the boundary sub-segment of the *next* resolved topology). This is not as ideal as actually resolving the + all the *unclipped* boundary sections of the current topology at the *next* time, which would produce a more accurate + plate boundary for the next resolved topology. But the inconsistency forces us to take an alternative approach. + However this case shouldn't happen very often. The next resolved topology is not consistent with the current resolved topology if its boundary sections don't have the same number of intersections with neighbouring sections. For example, a boundary section of the current resolved topology @@ -1374,6 +1376,22 @@ def __init__( self.curr_resolved_topology = curr_resolved_topology self.next_resolved_topology = next_resolved_topology + + self.curr_boundary_sub_segments = ( + curr_resolved_topology.get_boundary_sub_segments() + ) + self.next_boundary_sub_segments = ( + next_resolved_topology.get_boundary_sub_segments() + ) + + if len(self.curr_boundary_sub_segments) != len( + self.next_boundary_sub_segments + ): + raise RuntimeError( + "Current and next topologies have a different number of boundary sub-segments." + ) + self.num_boundary_sub_segments = len(self.curr_boundary_sub_segments) + self.curr_time = curr_resolved_topology.get_reconstruction_time() self.next_time = next_resolved_topology.get_reconstruction_time() self.rotation_model = rotation_model @@ -1382,169 +1400,247 @@ def get_next_resolved_topology_boundary(self): # If the current and next resolved topologies are consistent then # just return the resolved boundary of the next resolved topology. - if self._are_current_and_next_resolved_topologies_consistent(): + inconsistent_boundary_sub_segment_indices = ( + self._get_inconsistent_boundary_sub_segments() + ) + if not inconsistent_boundary_sub_segment_indices: return self.next_resolved_topology.get_resolved_boundary() - # Instead of returning the boundary of the next resolved topology, return the boundary - # sub-segments of the current resolved topology moved to the next time and joined together. - return self._get_next_boundary_from_current_boundary_sub_segments() + # Instead of returning the boundary of the next resolved topology, replace any inconsistent boundary sub-segments + # of the next resolved topology with the equivalent ones in the current resolved topology (but moved to the next time). + return self._get_next_boundary_from_consistent_boundary_sub_segments( + inconsistent_boundary_sub_segment_indices + ) - def _get_next_boundary_from_current_boundary_sub_segments(self): + def _get_next_boundary_from_consistent_boundary_sub_segments( + self, inconsistent_boundary_sub_segment_indices + ): # - # Instead of the resolved boundary of the *next* resolved topology we'll: - # - take the *clipped* boundary sub-segments of the *current* resolved topology, - # - move them to the next time, and - # - join them together into a polygon boundary. + # Instead of the resolved boundary of the *next* resolved topology we'll replace the inconsistent boundary sub-segments + # of the *next* resolved topology with the equivalent ones in the *current* resolved topology (but moved to the next time). # - # Collect the current boundary sub-segment features to use for the next boundary. - next_boundary_sub_segment_features = [] - next_boundary_sub_segment_reversals = [] - for ( - curr_boundary_sub_segment - ) in self.curr_resolved_topology.get_boundary_sub_segments(): - # See if the current boundary sub-segment is from a topological *line*. - # - # If it is then we need to iterate over its sub-segments. This is because we need features that are *reconstructable*. - # Topological lines are not reconstructable (because they're topologies linking reconstructable features). - # Only the sub-segments of topological lines are reconstructable. - sub_sub_segments = curr_boundary_sub_segment.get_sub_segments() - if sub_sub_segments: - curr_boundary_sub_segment_reversal = ( - curr_boundary_sub_segment.was_geometry_reversed_in_topology() + # Iterate over the next boundary sub-segments and join them to form a polygon boundary. + next_boundary_polygon_points = [] + for boundary_sub_segment_index in range(self.num_boundary_sub_segments): + + # If the boundary sub-segment is consistent (ie, not inconsistent) then just add the + # (potentially reversed) points of the boundary sub-segment of the *next* resolved topology. + if ( + boundary_sub_segment_index + not in inconsistent_boundary_sub_segment_indices + ): + next_boundary_sub_segment = self.next_boundary_sub_segments[ + boundary_sub_segment_index + ] + next_boundary_polygon_points.extend( + reversed( + next_boundary_sub_segment.get_resolved_geometry_points() + ) + if next_boundary_sub_segment.was_geometry_reversed_in_topology() + else next_boundary_sub_segment.get_resolved_geometry_points() ) - # Each sub-segment of the resolved topological *line* will contribute to the current boundary sub-segment. + + continue + + # + # The boundary sub-segment is inconsistent. + # + # So take the boundary sub-segment of the *current* resolved topology and move it to the *next* time, + # and use that (instead of the boundary sub-segment of the *next* resolved topology). + # + + curr_boundary_sub_segment = self.curr_boundary_sub_segments[ + boundary_sub_segment_index + ] + curr_boundary_sub_segment_reversal = ( + curr_boundary_sub_segment.was_geometry_reversed_in_topology() + ) + + # See if the boundary sub-segment is from a topological *line*. + curr_boundary_sub_sub_segments = ( + curr_boundary_sub_segment.get_sub_segments() + ) + if not curr_boundary_sub_sub_segments: # - # Note: Only a sub-section of the full resolved topological *line* will contribute - # to the boundary of the current resolved topological boundary/network. - # Hence not all the sub-segments of the resolved topological *line* will contribute - # to the boundary of the current resolved topological boundary/network. + # Boundary sub-segment is NOT from a topological *line*. # - # Note: We need to reverse the order of sub-sub-segments if the resolved topological *line* is reversed - # in the current resolved boundary. - next_boundary_sub_segment_features.extend( - sub_sub_segment.get_resolved_feature() - for sub_sub_segment in ( - reversed(sub_sub_segments) - if curr_boundary_sub_segment_reversal - else sub_sub_segments - ) + curr_boundary_sub_segment_feature = ( + curr_boundary_sub_segment.get_resolved_feature() + ) + + # The sub-segment resolved feature is already reconstructed to 'curr_time'. + # So we need to reverse reconstruct it to present day (before we can reconstruct it to 'next_time'). + # + # Note: We don't need to clone this feature before modifying it because it was already essentially cloned + # when 'ResolvedTopologicalSubSegment.get_resolved_feature()' was called. + pygplates.reverse_reconstruct( # type:ignore + curr_boundary_sub_segment_feature, + self.rotation_model, + self.curr_time, ) - # If a sub-sub-segment is reversed in the resolved topological *line* which is, in turn, reversed in the - # current resolved boundary then the sub-sub-segment is NOT reversed in the current resolved boundary. - next_boundary_sub_segment_reversals.extend( - curr_boundary_sub_segment_reversal - ^ sub_sub_segment.was_geometry_reversed_in_topology() - for sub_sub_segment in ( - reversed(sub_sub_segments) - if curr_boundary_sub_segment_reversal - else sub_sub_segments + + # The resolved sub-segment feature is valid at the 'curr_time' but not necessarily at 'next_time'. + # If not valid at 'next_time' then set a valid time that includes 'next_time'. + if not curr_boundary_sub_segment_feature.is_valid_at_time( + self.next_time + ): + curr_boundary_sub_segment_feature.set_valid_time( + pygplates.GeoTimeInstant.create_distant_past(), 0.0 # type: ignore ) + + # Reconstruct to 'next_time' (from present day). + next_boundary_sub_segment_reconstructed_geometries = ( + pygplates.ReconstructSnapshot( # type:ignore + curr_boundary_sub_segment_feature, + self.rotation_model, + self.next_time, + ).get_reconstructed_geometries() ) - else: - next_boundary_sub_segment_features.append( - curr_boundary_sub_segment.get_resolved_feature() + + # Should only be one reconstructed geometry. + if len(next_boundary_sub_segment_reconstructed_geometries) != 1: + raise RuntimeError("Expected a single reconstructed geometry.") + next_boundary_sub_segment_reconstructed_geometry = ( + next_boundary_sub_segment_reconstructed_geometries[0] ) - next_boundary_sub_segment_reversals.append( - curr_boundary_sub_segment.was_geometry_reversed_in_topology() + + # Add the (potentially reversed) points of the sub-segment to the next polygon boundary. + next_boundary_polygon_points.extend( + reversed( + next_boundary_sub_segment_reconstructed_geometry.get_reconstructed_geometry_points() + ) + if curr_boundary_sub_segment_reversal + else next_boundary_sub_segment_reconstructed_geometry.get_reconstructed_geometry_points() ) - # The sub-segment resolved features are already reconstructed to 'curr_time'. - # So we need to reverse reconstruct them to present day (before we can reconstruct them to 'next_time'). - # - # Note: We don't need to clone these features before modifying them because they were already essentially cloned - # when 'ResolvedTopologicalSubSegment.get_resolved_feature()' was called. - pygplates.reverse_reconstruct( # type:ignore - next_boundary_sub_segment_features, - self.rotation_model, - self.curr_time, - ) + continue + + # + # The boundary sub-segment IS from a topological *line*, so we need to iterate over its sub-segments. + # This is because we need features that are *reconstructable*. + # Topological lines are not reconstructable (because they're topologies linking reconstructable features). + # Only the sub-segments of topological lines are reconstructable. + # + # Each sub-segment of the resolved topological *line* will contribute to the boundary sub-segment. + # + # Note: Only a sub-section of the full resolved topological *line* will contribute + # to the boundary of the current resolved topological boundary/network. + # Hence not all the sub-segments of the resolved topological *line* will contribute + # to the boundary of the current resolved topological boundary/network. + # - # The resolved sub-segment features are valid at the 'curr_time' but not necessarily at 'next_time'. - # If not valid at 'next_time' then set a valid time that includes 'next_time'. - for feature in next_boundary_sub_segment_features: - if not feature.is_valid_at_time(self.next_time): - feature.set_valid_time( - pygplates.GeoTimeInstant.create_distant_past(), 0.0 # type: ignore + # Collect the sub-sub-segment features to use for the next boundary. + # + # Note: We need to reverse the order of sub-sub-segments if the resolved topological *line* is reversed + # in the current resolved boundary. + curr_boundary_sub_sub_segment_features = [ + sub_sub_segment.get_resolved_feature() + for sub_sub_segment in ( + reversed(curr_boundary_sub_sub_segments) + if curr_boundary_sub_segment_reversal + else curr_boundary_sub_sub_segments + ) + ] + # If a sub-sub-segment is reversed in the resolved topological *line* which is, in turn, reversed in the + # current resolved boundary then the sub-sub-segment is NOT reversed in the current resolved boundary. + curr_boundary_sub_sub_segment_reversals = [ + curr_boundary_sub_segment_reversal + ^ sub_sub_segment.was_geometry_reversed_in_topology() + for sub_sub_segment in ( + reversed(curr_boundary_sub_sub_segments) + if curr_boundary_sub_segment_reversal + else curr_boundary_sub_sub_segments ) + ] + + # The sub-segment resolved features are already reconstructed to 'curr_time'. + # So we need to reverse reconstruct them to present day (before we can reconstruct them to 'next_time'). + # + # Note: We don't need to clone these features before modifying them because they were already essentially cloned + # when 'ResolvedTopologicalSubSegment.get_resolved_feature()' was called. + pygplates.reverse_reconstruct( # type:ignore + curr_boundary_sub_sub_segment_features, + self.rotation_model, + self.curr_time, + ) - # Reconstruct to 'next_time' (from present day). - next_boundary_sub_segment_reconstructed_geometries = ( - pygplates.ReconstructSnapshot( # type:ignore - next_boundary_sub_segment_features, + # The resolved sub-segment features are valid at the 'curr_time' but not necessarily at 'next_time'. + # If not valid at 'next_time' then set a valid time that includes 'next_time'. + for feature in curr_boundary_sub_sub_segment_features: + if not feature.is_valid_at_time(self.next_time): + feature.set_valid_time( + pygplates.GeoTimeInstant.create_distant_past(), 0.0 # type: ignore + ) + + # Reconstruct to 'next_time' (from present day). + next_boundary_sub_sub_segment_reconstructed_geometries = pygplates.ReconstructSnapshot( # type:ignore + curr_boundary_sub_sub_segment_features, self.rotation_model, self.next_time, ).get_reconstructed_geometries( + # Make sure the order is retained... same_order_as_reconstructable_features=True ) - ) - # The order and length of reconstructed geometries should match that of the features. - if len(next_boundary_sub_segment_reconstructed_geometries) != len( - next_boundary_sub_segment_features - ): - raise RuntimeError( - "Expected number of reconstructed geometries to match number of features." - ) + # The order and length of reconstructed geometries should match that of the features. + if len(next_boundary_sub_sub_segment_reconstructed_geometries) != len( + curr_boundary_sub_sub_segment_features + ): + raise RuntimeError( + "Expected number of reconstructed geometries to match number of features." + ) - # Iterate over the next boundary sub-segments and join them to form a polygon boundary. - next_boundary_polygon_points = [] - for ( - feature_index, - next_boundary_sub_segment_reconstructed_geometry, - ) in enumerate(next_boundary_sub_segment_reconstructed_geometries): - next_boundary_polygon_points.extend( - reversed( - next_boundary_sub_segment_reconstructed_geometry.get_reconstructed_geometry_points() + # Add the (potentially reversed) points of the sub-segments to the next polygon boundary. + for ( + feature_index, + sub_sub_segment_reconstructed_geometry, + ) in enumerate(next_boundary_sub_sub_segment_reconstructed_geometries): + next_boundary_polygon_points.extend( + reversed( + sub_sub_segment_reconstructed_geometry.get_reconstructed_geometry_points() + ) + if curr_boundary_sub_sub_segment_reversals[feature_index] + else sub_sub_segment_reconstructed_geometry.get_reconstructed_geometry_points() ) - if next_boundary_sub_segment_reversals[feature_index] - else next_boundary_sub_segment_reconstructed_geometry.get_reconstructed_geometry_points() - ) + # Join the (potentially reversed) points of the boundary sub-segments and join them together into a polygon boundary. + # Each sub-segment is either from the *current* or *next* resolved topology. return pygplates.PolygonOnSphere( # type:ignore next_boundary_polygon_points ) - def _are_current_and_next_resolved_topologies_consistent(self): - - curr_boundary_sub_segments = ( - self.curr_resolved_topology.get_boundary_sub_segments() - ) - next_boundary_sub_segments = ( - self.next_resolved_topology.get_boundary_sub_segments() - ) - if len(curr_boundary_sub_segments) != len(next_boundary_sub_segments): - raise RuntimeError( - "Current and next topologies have a different number of boundary sub-segments." - ) - num_boundary_sub_segments = len(curr_boundary_sub_segments) + def _get_inconsistent_boundary_sub_segments(self): + """Returns a set of indices of any boundary sub-segments that are inconsistent, otherwise None.""" # If there's only 1 sub-segment then it cannot intersect with itself (but it can still be converted to a polygon). # # Note: If there's 2 sub-segments then they can intersect each other twice (to form a closed polygon). # We will still process them to make sure they stay consistent (same number of intersections at current and next times). - if num_boundary_sub_segments <= 1: - return True + if self.num_boundary_sub_segments <= 1: + return None + + inconsistent_boundary_sub_segment_indices = None # Compare number of neighbour intersections for the current and next boundary sub-segments. - for boundary_sub_segment_index in range(num_boundary_sub_segments): - curr_boundary_sub_segment = curr_boundary_sub_segments[ + for boundary_sub_segment_index in range(self.num_boundary_sub_segments): + curr_boundary_sub_segment = self.curr_boundary_sub_segments[ boundary_sub_segment_index ] - next_boundary_sub_segment = next_boundary_sub_segments[ + next_boundary_sub_segment = self.next_boundary_sub_segments[ boundary_sub_segment_index ] prev_boundary_sub_segment_index = boundary_sub_segment_index - 1 if prev_boundary_sub_segment_index < 0: - prev_boundary_sub_segment_index += num_boundary_sub_segments + prev_boundary_sub_segment_index += self.num_boundary_sub_segments # Previous neighbour sub-segment (for the current and next resolved topologies). - prev_curr_boundary_sub_segment = curr_boundary_sub_segments[ + prev_curr_boundary_sub_segment = self.curr_boundary_sub_segments[ prev_boundary_sub_segment_index ] - prev_next_boundary_sub_segment = next_boundary_sub_segments[ + prev_next_boundary_sub_segment = self.next_boundary_sub_segments[ prev_boundary_sub_segment_index ] @@ -1560,9 +1656,19 @@ def _are_current_and_next_resolved_topologies_consistent(self): # If the number of intersections differ then the current and next resolved topologies are inconsistent. if curr_num_intersections != next_num_intersections: - return False + # We only create a set when we actually need one + # (most often we won't since most resolved topologies are consistent). + if inconsistent_boundary_sub_segment_indices is None: + inconsistent_boundary_sub_segment_indices = set() + # Add the previous and current sub-segments affected by the intersection(s). + inconsistent_boundary_sub_segment_indices.add( + prev_boundary_sub_segment_index + ) + inconsistent_boundary_sub_segment_indices.add( + boundary_sub_segment_index + ) - return True + return inconsistent_boundary_sub_segment_indices def _get_num_intersections(self, polyline1, polyline2, num_intersections=0): From 7e85186acc9738eb63ce239e4d61bd6176746108 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Mon, 24 Nov 2025 14:17:38 +1100 Subject: [PATCH 07/17] Restore continent collision detection (essentially removed in 0582168). And we can now sample all points together, rather than one by one. --- gplately/lib/reconstruct_by_topologies.py | 55 +++++++++++++++++++++++ gplately/oceans.py | 2 + 2 files changed, 57 insertions(+) diff --git a/gplately/lib/reconstruct_by_topologies.py b/gplately/lib/reconstruct_by_topologies.py index 8d38280b..e6e161b8 100644 --- a/gplately/lib/reconstruct_by_topologies.py +++ b/gplately/lib/reconstruct_by_topologies.py @@ -867,11 +867,17 @@ def __init__( reconstruction_begin_time, reconstruction_end_time, reconstruction_time_interval, + continent_mask_filepath_format, points, *, point_begin_times: Union[np.ndarray, list, None] = None, point_end_times: Union[np.ndarray, list, None] = None, ): + """ + continent_mask_filepath_format: str + The format of the file path of the continental mask grids that is converted to a + file path using ``continent_mask_filepath_format.format(time)``. + """ super().__init__( reconstruction_begin_time, @@ -889,6 +895,8 @@ def __init__( topology_features ).get_features() + self.continent_mask_filepath_format = continent_mask_filepath_format + def _begin_reconstruction_impl(self): # Activate any original points whose begin time is equal to (or older than) the newly updated current time. # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. @@ -1091,6 +1099,9 @@ def _reconstruct_to_next_time_impl(self): # The current active point was not consumed by a plate boundary. self.next_points[curr_point_index] = next_active_point + # See if any 'next' points collide with the continents and deactivate those that do. + self._detect_continent_collisions(next_time) + # # Set up for next loop iteration. # @@ -1112,6 +1123,50 @@ def _reconstruct_to_next_time_impl(self): # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. self._activate_deactivate_points_using_their_begin_end_times() + def _detect_continent_collisions(self, next_time): + # Get those 'next' points that are active. + next_active_points = [] + next_active_points_indices = [] + for point_index, next_point in enumerate(self.next_points): + if next_point: + next_active_points.append(next_point) + next_active_points_indices.append(point_index) + + # Convert points to lat/lon. + points_lat = np.empty(len(next_active_points)) + points_lon = np.empty(len(next_active_points)) + for active_point_index, point in enumerate(next_active_points): + point_lat, point_lon = point.to_lat_lon() + points_lat[active_point_index] = point_lat + points_lon[active_point_index] = point_lon + + # Read the continent mask at 'next_time'. + continent_mask_filepath = self.continent_mask_filepath_format.format(next_time) + gridZ, gridX, gridY = read_netcdf_grid( + continent_mask_filepath, return_grids=True + ) + ni, nj = gridZ.shape + xmin = np.nanmin(gridX) + xmax = np.nanmax(gridX) + ymin = np.nanmin(gridY) + ymax = np.nanmax(gridY) + + # Sample continent mask grid, which is one over continents and zero over oceans. + points_i = (ni - 1) * ((points_lat - ymin) / (ymax - ymin)) + points_j = (nj - 1) * ((points_lon - xmin) / (xmax - xmin)) + points_i_uint = np.rint(points_i).astype(np.uint) + points_j_uint = np.rint(points_j).astype(np.uint) + try: + mask_values = gridZ[points_i_uint, points_j_uint] + except IndexError: + points_i = np.clip(np.rint(points_i), 0, ni - 1).astype(np.int_) + points_j = np.clip(np.rint(points_j), 0, nj - 1).astype(np.int_) + mask_values = gridZ[points_i, points_j] + + # Deactivate any points that sampled inside the continent mask. + for active_point_index in np.where(mask_values >= 0.5)[0]: + self.next_points[next_active_points_indices[active_point_index]] = None + class _NextTopologicalFeatures(object): """The *next* topology boundary features and their topological section features. diff --git a/gplately/oceans.py b/gplately/oceans.py index 724a32d4..348fb959 100644 --- a/gplately/oceans.py +++ b/gplately/oceans.py @@ -1415,6 +1415,8 @@ def _build_and_reconstruct_ocean_seed_points( from_time, self._min_time, self._ridge_time_step, + # This filename string should not have a time formatted into it - this is taken care of later... + self.continent_mask_filepath, points, point_begin_times=appearance_times, ) From 321f91fd3f9adeb8aeb76d6161eb66e0c3772b74 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Mon, 24 Nov 2025 20:25:52 +1100 Subject: [PATCH 08/17] Remove original ReconstructByTopologies implementation. And remove implementation based on TopologicalModel since it's no longer needed (since we now use a different collision detection algorithm and also the new implementation can reconstruct points that are in deforming networks, which the original implementation could not). --- gplately/lib/reconstruct_by_topologies.py | 1013 ++------------------- gplately/oceans.py | 95 +- 2 files changed, 106 insertions(+), 1002 deletions(-) diff --git a/gplately/lib/reconstruct_by_topologies.py b/gplately/lib/reconstruct_by_topologies.py index e6e161b8..c4c05876 100644 --- a/gplately/lib/reconstruct_by_topologies.py +++ b/gplately/lib/reconstruct_by_topologies.py @@ -14,7 +14,6 @@ # with this program; if not, write to Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. # -from abc import ABC, abstractmethod import logging import math import os @@ -29,401 +28,49 @@ logger = logging.getLogger("gplately") -class _DefaultCollision(object): - """ - Default collision detection function class (the function is the '__call__' method). - """ - - DEFAULT_GLOBAL_COLLISION_PARAMETERS = (7.0, 10.0) - """ - Default collision parameters for all feature types. - - This is a 2-tuple of (threshold velocity delta in kms/my, threshold distance to boundary per My in kms/my): - Here we default to the same constants used internally in GPlates 2.0 (ie, 7.0 and 10.0). - """ - - def __init__( - self, - global_collision_parameters=DEFAULT_GLOBAL_COLLISION_PARAMETERS, - feature_specific_collision_parameters=None, - ): - """ - global_collision_parameters: The collision parameters to use for any feature type not specified in 'feature_specific_collision_parameters'. - Should be a 2-tuple of (threshold velocity delta in kms/my, threshold distance to boundary per My in kms/my). - The first threshold parameter means: - A point that transitions from one plate to another can disappear if the change in velocity exceeds this threshold. - The second threshold parameter means: - Only those transitioning points exceeding the threshold velocity delta and that are close enough to a plate boundary can disappear. - The distance is proportional to the relative velocity (change in velocity), plus a constant offset based on the threshold distance to boundary - to account for plate boundaries that change shape significantly from one time step to the next - (note that some boundaries are meant to do this and others are a result of digitisation). - The actual distance threshold used is (threshold_distance_to_boundary + relative_velocity) * time_interval - Defaults to parameters used in GPlates 2.0, if not specified. - - feature_specific_collision_parameters: Optional sequence of collision parameters specific to feature types. - If specified then should be a sequence of 2-tuples, with each 2-tuple specifying (feature_type, collision_parameters). - And where each 'collision_parameters' is a 2-tuple of (threshold velocity delta in kms/my, threshold distance to boundary per My in kms/my). - See 'global_collision_parameters' for details on these thresholds. - Any feature type not specified here defaults to using 'global_collision_parameters'. - """ - - # Convert list of (feature_type, collision_parameters) tuples to a dictionary. - if feature_specific_collision_parameters: - self.feature_specific_collision_parameters = dict( - feature_specific_collision_parameters - ) - else: - self.feature_specific_collision_parameters = dict() - # Fallback for any feature type not specified in the optional feature-specific list. - self.global_collision_parameters = global_collision_parameters - - # Used to improve performance by caching velocity stage rotations in a dict (for a specific reconstruction time). - self.velocity_stage_rotation_dict = {} - self.velocity_stage_rotation_time = None - - def __call__( - self, - rotation_model, - time, - reconstruction_time_interval, - prev_point, - curr_point, - prev_topology_plate_id, - prev_resolved_plate_boundary, - curr_topology_plate_id, - curr_resolved_plate_boundary, - ): - """ - Returns True if a collision was detected. - - If transitioning from a rigid plate to another rigid plate with a different plate ID then - calculate the difference in velocities and continue testing as follows - (otherwise, if there's no transition, then the point is still active)... - - If the velocity difference is below a threshold then we assume the previous plate was split, - or two plates joined. In this case the point has not subducted (forward in time) or - been consumed by a mid-ocean (backward in time) and hence is still active. - - If the velocity difference is large enough then we see if the distance of the *previous* position - to the polygon boundary (of rigid plate containing it) exceeds a threshold. - If the distance exceeds the threshold then the point is far enough away from the boundary that it - cannot be subducted or consumed by it and hence the point is still active. - However if the point is close enough then we assume the point was subducted/consumed - (remember that the point switched plate IDs). - Also note that the threshold distance increases according to the velocity difference to account for fast - moving points (that would otherwise tunnel through the boundary and accrete onto the other plate). - The reason for testing the distance from the *previous* point, and not from the *current* point, is: - - (i) A topological boundary may *appear* near the current point (such as a plate split at the current time) - and we don't want that split to consume the current point regardless of the velocity difference. - It won't get consumed because the *previous* point was not near a boundary (because before split happened). - If the velocity difference is large enough then it might cause the current point to transition to the - adjacent split plate in the *next* time step (and that's when it should get consumed, not in the current time step). - An example of this is a mid-ocean ridge suddenly appearing (going forward in time). - - (ii) A topological boundary may *disappear* near the current point (such as a plate merge at the current time) - and we want that merge to consume the current point if the velocity difference is large enough. - In this case the *previous* point is near a boundary (because before plate merged) and hence can be - consumed (provided velocity difference is large enough). And since the boundary existed in the previous - time step, it will affect position of the current point (and whether it gets consumed or not). - An example of this is a mid-ocean ridge suddenly disappearing (going backward in time). - - ...note that items (i) and (ii) above apply both going forward and backward in time. - """ - - # See if a collision occurred. - if ( - curr_topology_plate_id != prev_topology_plate_id - and prev_topology_plate_id is not None - and curr_topology_plate_id is not None - ): - # - # Speed up by caching velocity stage rotations in a dict. - # - if time != self.velocity_stage_rotation_time: - # We've just switched to a new time so clear the cache. - # - # We only cache stage rotations for a specific time. - # We only really need to cache different plate IDs at the same 'time', so this avoids caching for all times - # (which would also require including 'time' in the key) and using memory unnecessarily. - self.velocity_stage_rotation_dict.clear() - self.velocity_stage_rotation_time = time - prev_location_velocity_stage_rotation = ( - self.velocity_stage_rotation_dict.get(prev_topology_plate_id) - ) - if not prev_location_velocity_stage_rotation: - prev_location_velocity_stage_rotation = rotation_model.get_rotation( - time + 1, prev_topology_plate_id, time - ) - self.velocity_stage_rotation_dict[prev_topology_plate_id] = ( - prev_location_velocity_stage_rotation - ) - curr_location_velocity_stage_rotation = ( - self.velocity_stage_rotation_dict.get(curr_topology_plate_id) - ) - if not curr_location_velocity_stage_rotation: - curr_location_velocity_stage_rotation = rotation_model.get_rotation( - time + 1, curr_topology_plate_id, time - ) - self.velocity_stage_rotation_dict[curr_topology_plate_id] = ( - curr_location_velocity_stage_rotation - ) - - # Note that even though the current point is not inside the previous boundary (because different plate ID), we can still - # calculate a velocity using its plate ID (because we really should use the same point in our velocity comparison). - prev_location_velocity = pygplates.calculate_velocities( # type: ignore - (curr_point,), - prev_location_velocity_stage_rotation, - 1, - pygplates.VelocityUnits.kms_per_my, - )[0] - curr_location_velocity = pygplates.calculate_velocities( # type: ignore - (curr_point,), - curr_location_velocity_stage_rotation, - 1, - pygplates.VelocityUnits.kms_per_my, - )[0] - - delta_velocity = curr_location_velocity - prev_location_velocity - delta_velocity_magnitude = delta_velocity.get_magnitude() - - # If we have feature-specific collision parameters then iterate over the boundary sub-segments of the *previous* topological boundary - # and test proximity to each sub-segment individually (with sub-segment feature type specific collision parameters). - # Otherwise just test proximity to the entire boundary polygon using the global collision parameters. - if self.feature_specific_collision_parameters: - for ( - prev_boundary_sub_segment - ) in prev_resolved_plate_boundary.get_boundary_sub_segments(): - # Use feature-specific collision parameters if found (falling back to global collision parameters). - ( - threshold_velocity_delta, - threshold_distance_to_boundary_per_my, - ) = self.feature_specific_collision_parameters.get( - prev_boundary_sub_segment.get_feature().get_feature_type(), - # Default to global collision parameters if no collision parameters specified for sub-segment's feature type... - self.global_collision_parameters, - ) - - # Since each feature type could use different collision parameters we must use the current boundary sub-segment instead of the boundary polygon. - if self._detect_collision_using_collision_parameters( - reconstruction_time_interval, - delta_velocity_magnitude, - prev_point, - prev_boundary_sub_segment.get_resolved_geometry(), - threshold_velocity_delta, - threshold_distance_to_boundary_per_my, - ): - # Detected a collision. - return True - else: - # No feature-specific collision parameters so use global fallback. - ( - threshold_velocity_delta, - threshold_distance_to_boundary_per_my, - ) = self.global_collision_parameters - - # Since all feature types use the same collision parameters we can use the boundary polygon instead of iterating over its sub-segments. - if self._detect_collision_using_collision_parameters( - reconstruction_time_interval, - delta_velocity_magnitude, - prev_point, - prev_resolved_plate_boundary.get_resolved_boundary(), - threshold_velocity_delta, - threshold_distance_to_boundary_per_my, - ): - # Detected a collision. - return True - - return False - - def _detect_collision_using_collision_parameters( - self, - reconstruction_time_interval, - delta_velocity_magnitude, - prev_point, - prev_boundary_geometry, - threshold_velocity_delta, - threshold_distance_to_boundary_per_my, - ): - if delta_velocity_magnitude > threshold_velocity_delta: - # Add the minimum distance threshold to the delta velocity threshold. - # The delta velocity threshold only allows those points that are close enough to the boundary to reach - # it given their current relative velocity. - # The minimum distance threshold accounts for sudden changes in the shape of a plate boundary - # which are no supposed to represent a new or shifted boundary but are just a result of the topology - # builder/user digitising a new boundary line that differs noticeably from that of the previous time period. - distance_threshold_radians = ( - (threshold_distance_to_boundary_per_my + delta_velocity_magnitude) - * reconstruction_time_interval - / pygplates.Earth.equatorial_radius_in_kms - ) - distance_threshold_radians = min(distance_threshold_radians, math.pi) - distance_threshold_radians = max(distance_threshold_radians, 0.0) - - distance = pygplates.GeometryOnSphere.distance( - prev_point, - prev_boundary_geometry, - distance_threshold_radians=float(distance_threshold_radians), - ) - if distance is not None: - # Detected a collision. - return True - - return False - - -_DEFAULT_COLLISION = _DefaultCollision() - +class ReconstructByTopologies(object): + """Reconstruct points using topologies and deactivates those that subduct at converging plate boundaries. -class _ContinentCollision(object): - """ - Continental collision detection function class (the function is the '__call__' method). + This uses a new approach to collision detection of seed points + (to determine if they get subducted going forward in time and hence deactivated). + Previously we detected if the velocity delta of a seed point exceeded a threshold (when that point + transitioned from one plate to another) - and the point had to be close enough to the boundary. + Instead, we now take the plate containing a seed point at the current time and resolve it at the next time + (ie, current time plus time step) and see if its resolved shape still contains that seed point + (reconstructed to the next time). If it does then the seed point remains active (otherwise it's deactivated). + Note that the *same* plate is resolved at the current and next times. Usually when you resolve topologies + at a different time you will get different plates (eg, there could be plate merges or splits). + So, to resolve the *same* plate, some fenagling is required to ensure that a resolved plate boundary at the + current time can also be resolved at the next time. The end result is that any boundary around a plate that + converges will naturally deactivate seed points near that boundary (without having to use thresholds, etc). + And this doesn't require plate boundaries to be labelled as converging (ie, subduction zones). """ def __init__( self, - grd_output_dir, - chain_collision_detection=_DEFAULT_COLLISION, - verbose=False, - ): - """ - grd_output_dir: The directory containing the continental grids. - - chain_collision_detection: Another collision detection class/function to reference if we find no collision. - If None then no collision detection is chained. Defaults to the default collision detection. - - verbose: Print progress messages - """ - - self.grd_output_dir = grd_output_dir - self.chain_collision_detection = chain_collision_detection - self.verbose = verbose - - # Load a new grid each time the reconstruction time changes. - self.grid_time = None - - @property - def grid_time(self): - return self._grid_time - - @grid_time.setter - def grid_time(self, time): - if time is None: - self._grid_time = time - else: - filename = "{:s}".format(self.grd_output_dir.format(time)) - if self.verbose: - print( - "Points masked against grid: {0}".format(os.path.basename(filename)) - ) - gridZ, gridX, gridY = read_netcdf_grid(filename, return_grids=True) - self.gridZ = gridZ - self.ni, self.nj = gridZ.shape - self.xmin = np.nanmin(gridX) - self.xmax = np.nanmax(gridX) - self.ymin = np.nanmin(gridY) - self.ymax = np.nanmax(gridY) - self._grid_time = float(time) - - def __call__( - self, - rotation_model, - time, - reconstruction_time_interval, - prev_point, - curr_point, - prev_topology_plate_id, - prev_resolved_plate_boundary, - curr_topology_plate_id, - curr_resolved_plate_boundary, - ): - """ - Returns True if a collision with a continent was detected, or returns result of - chained collision detection if 'self.chain_collision_detection' is not None. - """ - # Load the grid for the current time if encountering a new time. - if time != self.grid_time: - self.grid_time = time - self.continent_deletion_count = 0 - - # Sample mask grid, which is one over continents and zero over oceans. - point_lat, point_lon = curr_point.to_lat_lon() - point_i = (self.ni - 1) * ((point_lat - self.ymin) / (self.ymax - self.ymin)) - point_j = (self.nj - 1) * ((point_lon - self.xmin) / (self.xmax - self.xmin)) - point_i_uint = np.rint(point_i).astype(np.uint) - point_j_uint = np.rint(point_j).astype(np.uint) - try: - mask_value = self.gridZ[point_i_uint, point_j_uint] - except IndexError: - point_i = np.clip(np.rint(point_i), 0, self.ni - 1).astype(np.int_) - point_j = np.clip(np.rint(point_j), 0, self.nj - 1).astype(np.int_) - mask_value = self.gridZ[point_i, point_j] - if mask_value >= 0.5: - # Detected a collision. - self.continent_deletion_count += 1 - return True - - # We didn't find a collision, so ask the chained collision detection if it did (if we have anything chained). - if self.chain_collision_detection: - return self.chain_collision_detection( - rotation_model, - time, - reconstruction_time_interval, - prev_point, - curr_point, - prev_topology_plate_id, - prev_resolved_plate_boundary, - curr_topology_plate_id, - curr_resolved_plate_boundary, - ) - - return False - - -class _ReconstructByTopologies(ABC): - """Reconstruct geometries using topologies. Currently only points are supported.""" - - @abstractmethod - def _begin_reconstruction_impl(self): - """Called from 'begin_reconstruction()' and implemented in derived class. - - Derived class should call 'self._activate_deactivate_points_using_their_begin_end_times()' to introduce points - whose begin time is equal to (or older than) the initial time. - - Derived class can optionally do anything else needed before starting reconstruction over the time intervals. - """ - pass # This is an abstract method, no implementation here. - - @abstractmethod - def _reconstruct_to_next_time_impl(self): - """Called from 'reconstruct_to_next_time()' and implemented in derived class. - - Derived class should reconstruct the current active points in 'self.curr_points' to the next points 'self.next_points' from - the current time 'self.get_current_time()' to the next time 'self.get_current_time() + self.reconstruction_time_step'. - And collision detected for points should deactivate those points (ie, set None in 'self.next_points'). - - Derived class then should cycle the three point arrays in preparation for the next time interval. - Essentially do this: - self.prev_points, self.curr_points, self.next_points = self.curr_points, self.next_points, self.prev_points - - Derived class then should increment the current time index (to update the current time to the next time). - Essentially do this: - self.current_time_index += 1 - - Finally, derived class needs to call 'self._activate_deactivate_points_using_their_begin_end_times()' to introduce points - whose begin time is the new current time, and remove points whose end time is the new current time. - """ - pass # This is an abstract method, no implementation here. - - def __init__( - self, + rotation_features_or_model, + topology_features, reconstruction_begin_time, reconstruction_end_time, reconstruction_time_interval, points, + *, point_begin_times: Union[np.ndarray, list, None] = None, point_end_times: Union[np.ndarray, list, None] = None, + continent_mask_filepath_format=None, ): + """ + continent_mask_filepath_format: str + The format of the file path of the continental mask grids that is converted to a + file path using ``continent_mask_filepath_format.format(time)``. + """ + + self.rotation_model = pygplates.RotationModel(rotation_features_or_model) # type: ignore + + # Turn topology data into a list of features (if not already). + self.topology_features = pygplates.FeaturesFunctionArgument( # type: ignore + topology_features + ).get_features() # Set up an array of reconstruction times covering the reconstruction time span. self.reconstruction_begin_time = reconstruction_begin_time @@ -485,6 +132,8 @@ def __init__( "Length of 'point_end_times' must match length of 'points'." ) + self.continent_mask_filepath_format = continent_mask_filepath_format + def reconstruct(self): # Initialise the reconstruction. self.begin_reconstruction() @@ -510,8 +159,9 @@ def begin_reconstruction(self): self.point_has_been_activated = np.zeros(self.num_points, dtype=bool) self.num_activated_points = 0 - # Call derived class implementation. - self._begin_reconstruction_impl() + # Activate any original points whose begin time is equal to (or older than) the newly updated current time. + # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. + self._activate_deactivate_points_using_their_begin_end_times() def get_current_time_index(self): return self.current_time_index @@ -540,368 +190,12 @@ def reconstruct_to_next_time(self): if self.num_activated_points == self.num_points and not any(self.curr_points): return False - # Call derived class implementation. + # Call the main implementation. self._reconstruct_to_next_time_impl() # We successfully reconstructed to the next time. return True - def _activate_deactivate_points_using_their_begin_end_times(self): - current_time = self.get_current_time() - - # Iterate over all points and activate/deactivate as necessary depending on each point's valid time range. - for point_index in range(self.num_points): - if self.curr_points[point_index] is None: - if not self.point_has_been_activated[point_index]: - # Point is not active and has never been activated, so see if can activate it. - if ( - current_time <= self.point_begin_times[point_index] - and current_time >= self.point_end_times[point_index] - ): - # The initial point is assumed to be the position at the current time - # which is typically the point's begin time (approximately). - # But it could be the beginning of the reconstruction time span (specified in constructor) - # if that falls in the middle of the point's valid time range - in this case the - # initial point position is assumed to be in a position that is some time *after* - # it appeared (at its begin time) - and this can happen, for example, if you have a - # uniform grids of points at some intermediate time and want to see how they - # reconstruct to either a younger or older time (remembering that points can - # be subducted forward in time and consumed back into a mid-ocean ridge going - # backward in time). - self.curr_points[point_index] = self.points[point_index] - self.point_has_been_activated[point_index] = True - self.num_activated_points += 1 - else: - # Point is active, so see if can deactivate it. - if not ( - current_time <= self.point_begin_times[point_index] - and current_time >= self.point_end_times[point_index] - ): - self.curr_points[point_index] = None - - -class _ReconstructByTopologiesImpl(_ReconstructByTopologies): - """Reconstruct geometries using topologies. Currently only points are supported.""" - - use_plate_partitioner = False - """If the use_plate_partitioner is True then use pygplates.PlatePartitioner to partition points, - otherwise use faster points_in_polygons.find_polygons().""" - - def __init__( - self, - rotation_features_or_model, - topology_features, - reconstruction_begin_time, - reconstruction_end_time, - reconstruction_time_interval, - points, - *, - point_begin_times: Union[np.ndarray, list, None] = None, - point_end_times: Union[np.ndarray, list, None] = None, - point_plate_ids: Union[np.ndarray, list, None] = None, - detect_collisions=_DEFAULT_COLLISION, - ): - """ - rotation_features_or_model: Rotation model or feature collection(s), or list of features, or filename(s). - - topology_features: Topology feature collection(s), or list of features, or filename(s) or any combination of those. - - detect_collisions: Collision detection function, or None. Defaults to _DEFAULT_COLLISION. - """ - - super().__init__( - reconstruction_begin_time, - reconstruction_end_time, - reconstruction_time_interval, - points, - point_begin_times, - point_end_times, - ) - - self.rotation_model = pygplates.RotationModel(rotation_features_or_model) - - # Turn topology data into a list of features (if not already). - self.topology_features = pygplates.FeaturesFunctionArgument( - topology_features - ).get_features() - - # Use the specified point plate IDs if provided (otherwise use '0'). - # These plate IDs are only used when a point falls outside all resolved topologies during a time step. - if point_plate_ids is None: - self.point_plate_ids = np.zeros(self.num_points, dtype=int) - else: - # Make sure numpy array (if not already). - self.point_plate_ids = np.asarray(point_plate_ids) - if len(self.point_plate_ids) != self.num_points: - raise ValueError( - "Length of 'point_plate_ids' must match length of 'points'." - ) - - self.detect_collisions = detect_collisions - - def _begin_reconstruction_impl(self): - # Set up topology arrays (corresponding to active/inactive points at same indices as the point arrays). - self.prev_topology_plate_ids = [None] * self.num_points - self.curr_topology_plate_ids = [None] * self.num_points - self.prev_resolved_plate_boundaries = [None] * self.num_points - self.curr_resolved_plate_boundaries = [None] * self.num_points - - # Activate any original points whose begin time is equal to (or older than) the newly updated current time. - # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. - # - # Note: This should be done before calling 'self._find_resolved_topologies_containing_points()' so that new - # points just appearing at their begin time will have associated topologies (that contain them). - self._activate_deactivate_points_using_their_begin_end_times() - - self._find_resolved_topologies_containing_points() - - def _reconstruct_to_next_time_impl(self): - # Cache stage rotations by plate ID. - # Use different dicts since using different rotation models and time steps, etc. - reconstruct_stage_rotation_dict = {} - - current_time = self.get_current_time() - - # Iterate over all points to reconstruct them to the next time step. - for point_index in range(self.num_points): - curr_point = self.curr_points[point_index] - if curr_point is None: - # Current point is not currently active. - # So we cannot reconstruct to next time. - self.next_points[point_index] = None - continue - - # Get plate ID of resolved topology containing current point - # (this was determined in last call to '_find_resolved_topologies_containing_points()'). - curr_plate_id = self.curr_topology_plate_ids[point_index] - if curr_plate_id is None: - # Current point is currently active but it fell outside all resolved polygons. - # So instead we just reconstruct using its plate ID (that was manually assigned by the user/caller). - curr_plate_id = self.point_plate_ids[point_index] - - # Get the stage rotation that will move the point from where it is at the current time to its - # location at the next time step, based on the plate id that contains the point at the current time. - - # Speed up by caching stage rotations in a dict. - stage_rotation = reconstruct_stage_rotation_dict.get(curr_plate_id) - if not stage_rotation: - stage_rotation = self.rotation_model.get_rotation( - # Positive/negative time step means reconstructing backward/forward in time. - current_time + self.reconstruction_time_step, - curr_plate_id, - current_time, - ) - reconstruct_stage_rotation_dict[curr_plate_id] = stage_rotation - - # Use the stage rotation to reconstruct the tracked point from position at current time - # to position at the next time step. - self.next_points[point_index] = stage_rotation * curr_point - - # - # Set up for next loop iteration. - # - # Rotate previous, current and next point arrays. - # The new previous will be the old current. - # The new current will be the old next. - # The new next will be the old previous (but values are ignored and overridden in next time step; just re-using its memory). - self.prev_points, self.curr_points, self.next_points = ( - self.curr_points, - self.next_points, - self.prev_points, - ) - # - # Swap previous and current topology arrays. - # The new previous will be the old current. - # The new current will be the old previous (but values are ignored and overridden in next time step; just re-using its memory). - self.prev_topology_plate_ids, self.curr_topology_plate_ids = ( - self.curr_topology_plate_ids, - self.prev_topology_plate_ids, - ) - self.prev_resolved_plate_boundaries, self.curr_resolved_plate_boundaries = ( - self.curr_resolved_plate_boundaries, - self.prev_resolved_plate_boundaries, - ) - # - # Move the current time to the next time. - self.current_time_index += 1 - current_time = self.get_current_time() - - # Activate any original points whose begin time is equal to (or older than) the newly updated current time. - # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. - # - # Note: This should be done before calling 'self._find_resolved_topologies_containing_points()' so that new - # points just appearing at their begin time will have associated topologies (that contain them). - self._activate_deactivate_points_using_their_begin_end_times() - - self._find_resolved_topologies_containing_points() - - # Iterate over all points to detect collisions. - if self.detect_collisions: - for point_index in range(self.num_points): - prev_point = self.prev_points[point_index] - curr_point = self.curr_points[point_index] - if prev_point is None or curr_point is None: - # If current point is not currently active then no need to detect a collision for it (to deactivate it). - # Also previous point might just have been activated now, at end of current time step, and hence - # not active at beginning of time step. - continue - - # Get plate IDs of resolved topology containing previous and current point - # (this was determined in last call to '_find_resolved_topologies_containing_points()'). - # - # Note that could be None, so the collision detection needs to handle that. - prev_plate_id = self.prev_topology_plate_ids[point_index] - curr_plate_id = self.curr_topology_plate_ids[point_index] - - # Detect collisions at the end of the current time step since we need previous, and current, points and topologies. - # De-activate point (in 'curr_points') if subducted (forward in time) or consumed back into MOR (backward in time). - if self.detect_collisions( - self.rotation_model, - current_time, - self.reconstruction_time_interval, - prev_point, - curr_point, - prev_plate_id, - self.prev_resolved_plate_boundaries[point_index], - curr_plate_id, - self.curr_resolved_plate_boundaries[point_index], - ): - # An inactive point in 'curr_points' becomes None. - # It may have been reconstructed from the previous time step to a valid position - # but now we override that result as inactive. - self.curr_points[point_index] = None - - def _find_resolved_topologies_containing_points(self): - current_time = self.get_current_time() - - # Resolve the plate polygons for the current time. - resolved_topologies = [] - pygplates.resolve_topologies( # type: ignore - self.topology_features, - self.rotation_model, - resolved_topologies, - current_time, - ) - - if _ReconstructByTopologiesImpl.use_plate_partitioner: - # Create a plate partitioner from the resolved polygons. - plate_partitioner = pygplates.PlatePartitioner( - resolved_topologies, self.rotation_model - ) - else: - # Some of 'curr_points' will be None so 'curr_valid_points' contains only the valid (not None) - # points, and 'curr_valid_points_indices' is the same length as 'curr_points' but indexes into - # 'curr_valid_points' so we can quickly find which point (and hence which resolved topology) - # in 'curr_valid_points' is associated with the a particular point in 'curr_points'. - curr_valid_points = [] - curr_valid_points_indices = [None] * self.num_points - for point_index, curr_point in enumerate(self.curr_points): - if curr_point is not None: - curr_valid_points_indices[point_index] = len(curr_valid_points) # type: ignore - curr_valid_points.append(curr_point) - # For each valid current point find the resolved topology containing it. - resolved_topologies_containing_curr_valid_points = ( - _ptt.utils.points_in_polygons.find_polygons( - curr_valid_points, - [ - resolved_topology.get_resolved_boundary() - for resolved_topology in resolved_topologies - ], - resolved_topologies, - ) - ) - - # Iterate over all points. - for point_index, curr_point in enumerate(self.curr_points): - if curr_point is None: - # Current point is not currently active - so skip it. - self.curr_topology_plate_ids[point_index] = None - self.curr_resolved_plate_boundaries[point_index] = None - continue - - # Find the plate id of the polygon that contains 'curr_point'. - if _ReconstructByTopologiesImpl.use_plate_partitioner: - curr_polygon = plate_partitioner.partition_point(curr_point) # type: ignore - else: - curr_polygon = resolved_topologies_containing_curr_valid_points[ # type: ignore - # Index back into 'curr_valid_points' and hence also into - # 'resolved_topologies_containing_curr_valid_points'. - curr_valid_points_indices[point_index] # type: ignore - ] # type: ignore - self.curr_resolved_plate_boundaries[point_index] = curr_polygon - - # If the polygon is None, that means (presumably) that it fell into a crack between - # topologies. So it will be skipped and thrown away from future iterations. - if curr_polygon is None: - self.curr_topology_plate_ids[point_index] = None - continue - - # Set the plate ID of resolved topology containing current point. - self.curr_topology_plate_ids[point_index] = ( - curr_polygon.get_feature().get_reconstruction_plate_id() - ) - - -class _ReconstructByTopologiesImplV2(_ReconstructByTopologies): - """An improved version of _ReconstructByTopologiesImpl. - - This uses a different approach to collision detection of seed points - (to determine if they get subducted going forward in time and hence deactivated). - Previously we detected if the velocity delta of a seed point exceeded a threshold (when that point - transitioned from one plate to another) - and the point had to be close enough to the boundary. - Instead, we now take the plate containing a seed point at the current time and resolve it at the next time - (ie, current time plus time step) and see if its resolved shape still contains that seed point - (reconstructed to the next time). If it does then the seed point remains active (otherwise it's deactivated). - Note that the *same* plate is resolved at the current and next times. Usually when you resolve topologies - at a different time you will get different plates (eg, there could be plate merges or splits). - So, to resolve the *same* plate, some fenagling is required to ensure that a resolved plate boundary at the - current time can also be resolved at the next time. The end result is that any boundary around a plate that - converges will naturally deactivate seed points near that boundary (without having to use thresholds, etc). - And this doesn't require plate boundaries to be labelled as converging (ie, subduction zones). - """ - - def __init__( - self, - rotation_features_or_model, - topology_features, - reconstruction_begin_time, - reconstruction_end_time, - reconstruction_time_interval, - continent_mask_filepath_format, - points, - *, - point_begin_times: Union[np.ndarray, list, None] = None, - point_end_times: Union[np.ndarray, list, None] = None, - ): - """ - continent_mask_filepath_format: str - The format of the file path of the continental mask grids that is converted to a - file path using ``continent_mask_filepath_format.format(time)``. - """ - - super().__init__( - reconstruction_begin_time, - reconstruction_end_time, - reconstruction_time_interval, - points, - point_begin_times, - point_end_times, - ) - - self.rotation_model = pygplates.RotationModel(rotation_features_or_model) # type: ignore - - # Turn topology data into a list of features (if not already). - self.topology_features = pygplates.FeaturesFunctionArgument( # type: ignore - topology_features - ).get_features() - - self.continent_mask_filepath_format = continent_mask_filepath_format - - def _begin_reconstruction_impl(self): - # Activate any original points whose begin time is equal to (or older than) the newly updated current time. - # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. - self._activate_deactivate_points_using_their_begin_end_times() - def _reconstruct_to_next_time_impl(self): current_time = self.get_current_time() @@ -1123,7 +417,44 @@ def _reconstruct_to_next_time_impl(self): # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. self._activate_deactivate_points_using_their_begin_end_times() + def _activate_deactivate_points_using_their_begin_end_times(self): + current_time = self.get_current_time() + + # Iterate over all points and activate/deactivate as necessary depending on each point's valid time range. + for point_index in range(self.num_points): + if self.curr_points[point_index] is None: + if not self.point_has_been_activated[point_index]: + # Point is not active and has never been activated, so see if can activate it. + if ( + current_time <= self.point_begin_times[point_index] + and current_time >= self.point_end_times[point_index] + ): + # The initial point is assumed to be the position at the current time + # which is typically the point's begin time (approximately). + # But it could be the beginning of the reconstruction time span (specified in constructor) + # if that falls in the middle of the point's valid time range - in this case the + # initial point position is assumed to be in a position that is some time *after* + # it appeared (at its begin time) - and this can happen, for example, if you have a + # uniform grids of points at some intermediate time and want to see how they + # reconstruct to either a younger or older time (remembering that points can + # be subducted forward in time and consumed back into a mid-ocean ridge going + # backward in time). + self.curr_points[point_index] = self.points[point_index] + self.point_has_been_activated[point_index] = True + self.num_activated_points += 1 + else: + # Point is active, so see if can deactivate it. + if not ( + current_time <= self.point_begin_times[point_index] + and current_time >= self.point_end_times[point_index] + ): + self.curr_points[point_index] = None + def _detect_continent_collisions(self, next_time): + # If no continental collision detection requested then just return. + if self.continent_mask_filepath_format is None: + return + # Get those 'next' points that are active. next_active_points = [] next_active_points_indices = [] @@ -1912,191 +1243,6 @@ def _post_split_polyline( return pygplates.PolylineOnSphere(post_split_points) # type: ignore -class _ReconstructByTopologicalModelImpl(_ReconstructByTopologies): - """Similar to ``_ReconstructByTopologiesImpl`` except uses the ``pygplates.TopologicalModel`` class to reconstruct seed points. - - This is currently just a transition towards using ``pygplates.TopologicalModel``. - But note that this still uses Python code, like in class ``_DefaultCollision``, to do the collision detection - (as opposed to collision detection built into pyGPlates, which needs to be updated to better support GPlately). - Also it's not as fast as ``_ReconstructByTopologiesImpl`` which has faster point-in-polygon testing (since it uses a spatial tree). - - So, for the time being it's probably still better to use class ``_ReconstructByTopologiesImpl`` instead. - """ - - class CollisionDelegator(pygplates.ReconstructedGeometryTimeSpan.DeactivatePoints): - """Delegates collision detection from pyGPlates (pygplates.TopologicalModel) to the collision classes accessed by ``_ReconstructByTopologies`` - (like class ``_DefaultCollision``). - - Ultimately there should be less (or no) need to delegate once the default pyGPlates collision handler (``DefaultDeactivatePoints``) supports the - funcionality of classes like ``_DefaultCollision``. Although ``_ContinentCollision`` might still need to be handled outside pyGPlates (ie, delegated to). - """ - - def __init__( - self, detect_collisions, rotation_model, reconstruction_time_interval - ): - super().__init__() - self.detect_collisions = detect_collisions - self.rotation_model = rotation_model - self.reconstruction_time_interval = reconstruction_time_interval - - def deactivate( - self, - prev_point, - prev_location, - prev_time, - curr_point, - curr_location, - curr_time, - ): - # Get plate IDs of resolved topology containing previous and current point. - # - # Note that could be None, so the collision detection needs to handle that. - prev_plate_id = None - curr_plate_id = None - - # See if previous point is located in a resolved rigid plate or deforming network (or neither). - prev_resolved_topology = prev_location.located_in_resolved_boundary() - if not prev_resolved_topology: - prev_resolved_topology = prev_location.located_in_resolved_network() - if prev_resolved_topology: - prev_plate_id = ( - prev_resolved_topology.get_feature().get_reconstruction_plate_id() - ) - - # See if current point is located in a resolved rigid plate or deforming network (or neither). - curr_resolved_topology = curr_location.located_in_resolved_boundary() - if not curr_resolved_topology: - curr_resolved_topology = curr_location.located_in_resolved_network() - if curr_resolved_topology: - curr_plate_id = ( - curr_resolved_topology.get_feature().get_reconstruction_plate_id() - ) - - # Delegate to the Python implementation of collision detection (eg, classes '_DefaultCollision' and '_ContinentCollision'). - return self.detect_collisions( - self.rotation_model, - curr_time, - self.reconstruction_time_interval, - prev_point, - curr_point, - prev_plate_id, - prev_resolved_topology, - curr_plate_id, - curr_resolved_topology, - ) - - def __init__( - self, - rotation_features_or_model, - topology_features, - reconstruction_begin_time, - reconstruction_end_time, - reconstruction_time_interval, - points, - *, - point_begin_times: Union[np.ndarray, list, None] = None, - point_end_times: Union[np.ndarray, list, None] = None, - detect_collisions=_DEFAULT_COLLISION, - ): - """ - rotation_features_or_model: Rotation model or feature collection(s), or list of features, or filename(s). - - topology_features: Topology feature collection(s), or list of features, or filename(s) or any combination of those. - - detect_collisions: Collision detection function, or None. Defaults to _DEFAULT_COLLISION. - """ - - super().__init__( - reconstruction_begin_time, - reconstruction_end_time, - reconstruction_time_interval, - points, - point_begin_times, - point_end_times, - ) - - self.topological_model = pygplates.TopologicalModel( - topology_features, - rotation_features_or_model, - # Only really need to cache 2 topological snapshots since we progressively move through time (and hence never return to previous times). - # Might only need to cache 1 topological snapshot (not sure). Best to be safe with 2... - topological_snapshot_cache_size=2, - ) - - self.detect_collisions = self.CollisionDelegator( - detect_collisions, - self.topological_model.get_rotation_model(), - self.reconstruction_time_interval, - ) - - def _begin_reconstruction_impl(self): - # Activate any original points whose begin time is equal to (or older than) the newly updated current time. - # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. - self._activate_deactivate_points_using_their_begin_end_times() - - def _reconstruct_to_next_time_impl(self): - - current_time = self.get_current_time() - next_time = current_time + self.reconstruction_time_step - - # Find the active points. - # These will be reconstructed to the next time step. - current_active_points = [] - current_active_point_indices = [] - for point_index in range(self.num_points): - current_point = self.curr_points[point_index] - if current_point is None: - # Current point is not currently active. - # So we cannot reconstruct to next time. - self.next_points[point_index] = None - continue - - # Keep track of the currently active points. - # And their original indices into ALL points (active and inactive). - current_active_points.append(current_point) - current_active_point_indices.append(point_index) - - # Reconstruct active points to the next time step. - reconstructed_time_span = self.topological_model.reconstruct_geometry( - current_active_points, - initial_time=current_time, - oldest_time=max(current_time, next_time), - youngest_time=min(current_time, next_time), - time_increment=self.reconstruction_time_interval, # must be positive - deactivate_points=self.detect_collisions, - ) - - # Store the next points back to their original locations in ALL points (active and inactive). - next_points = reconstructed_time_span.get_geometry_points( - next_time, return_inactive_points=True - ) - for next_point_index, next_point in enumerate(next_points): - # Get index into ALL points (active and inactive). - point_index = current_active_point_indices[next_point_index] - # Update next point (note that this can be None if a collision was detected). - self.next_points[point_index] = next_point - - # - # Set up for next loop iteration. - # - # Rotate previous, current and next point arrays. - # The new previous will be the old current. - # The new current will be the old next. - # The new next will be the old previous (but values are ignored and overridden in next time step; just re-using its memory). - self.prev_points, self.curr_points, self.next_points = ( - self.curr_points, - self.next_points, - self.prev_points, - ) - # - # Move the current time to the next time. - self.current_time_index += 1 - - # Activate any original points whose begin time is equal to (or older than) the newly updated current time. - # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. - self._activate_deactivate_points_using_their_begin_end_times() - - def reconstruct_points_by_topologies( rotation_features_or_model, topology_features, @@ -2104,14 +1250,14 @@ def reconstruct_points_by_topologies( reconstruction_end_time, reconstruction_time_interval, points, + *, point_begin_times=None, point_end_times=None, - point_plate_ids=None, - detect_collisions=_DEFAULT_COLLISION, + continent_mask_filepath_format=None, ): - """Reconstruct points using the topological polygons.""" + """Reconstruct points using topologies.""" - topology_reconstruction = _ReconstructByTopologiesImpl( + topology_reconstruction = ReconstructByTopologies( rotation_features_or_model, topology_features, reconstruction_begin_time, @@ -2120,8 +1266,7 @@ def reconstruct_points_by_topologies( points, point_begin_times=point_begin_times, point_end_times=point_end_times, - point_plate_ids=point_plate_ids, - detect_collisions=detect_collisions, + continent_mask_filepath_format=continent_mask_filepath_format, ) return topology_reconstruction.reconstruct() diff --git a/gplately/oceans.py b/gplately/oceans.py index 348fb959..ac1e0834 100644 --- a/gplately/oceans.py +++ b/gplately/oceans.py @@ -30,11 +30,7 @@ from . import grids, tools from .lib.reconstruct_by_topologies import ( - _ContinentCollision, - _DefaultCollision, - _ReconstructByTopologiesImpl, - _ReconstructByTopologiesImplV2, - _ReconstructByTopologicalModelImpl, + ReconstructByTopologies, ) from .lib.reconstruct_continents import ReconstructContinents from .parallel import get_num_cpus @@ -1146,22 +1142,23 @@ def _calc_shifted_mor_point_and_spreading_rate( return shifted_mor_points, np.array(point_spreading_rates) + def reconstruct_by_topological_model(self): + """Alias for :meth:`reconstruct_by_topologies`. + + Introduced in version ``2.0``, this method used to use `pygplates.TopologicalModel`_ class to reconstruct seed points. + Now it's just an alias for :meth:`reconstruct_by_topologies` and is deprecated. + + .. deprecated:: 2.1 + + Use :meth:`reconstruct_by_topologies` instead. + """ + self.reconstruct_by_topologies() + def reconstruct_by_topologies(self): """Obtain all active ocean seed points which are points that have not been consumed at subduction zones or have not collided with continental polygons. Active points' latitudes, longitues, seafloor ages, spreading rates and all other general z-values are saved to a gridding input file (.npz). """ - self._reconstruct_by_topologies_impl(use_topological_model=False) - - def reconstruct_by_topological_model(self): - """Use `pygplates.TopologicalModel`_ class to reconstruct seed points. - This method is an alternative to :meth:`reconstruct_by_topologies()` which uses Python code to do the reconstruction. - - .. _pygplates.TopologicalModel: https://www.gplates.org/docs/pygplates/generated/pygplates.topologicalmodel - """ - self._reconstruct_by_topologies_impl(use_topological_model=True) - - def _reconstruct_by_topologies_impl(self, use_topological_model): logger.info("Preparing to reconstruct ocean seed points using topologies...") @@ -1212,15 +1209,12 @@ def _reconstruct_by_topologies_impl(self, use_topological_model): partial( _build_and_reconstruct_ocean_seed_points_parallel, seafloor_grid=self, - use_topological_model=use_topological_model, ), self._times, ) else: for time in self._times: - self._build_and_reconstruct_ocean_seed_points( - time, use_topological_model=use_topological_model - ) + self._build_and_reconstruct_ocean_seed_points(time) logger.info( "Aggregating reconstructed ocean seed point data for gridding input..." @@ -1306,9 +1300,7 @@ def _reconstruct_by_topologies_impl(self, use_topological_model): all_reconstructed_seed_point_data ) - def _build_and_reconstruct_ocean_seed_points( - self, from_time, use_topological_model - ): + def _build_and_reconstruct_ocean_seed_points(self, from_time): """Creates ocean seed points at `from_time` and reconstructs them to `min_time`. Ocean seed points can be either the *initial* ocean seed points at `from_time=max_time` or @@ -1378,49 +1370,19 @@ def _build_and_reconstruct_ocean_seed_points( # Begin the reconstruction by topology process. # - # Specify the default collision detection region as subduction zones - default_collision = _DefaultCollision( - feature_specific_collision_parameters=[ - ( - pygplates.FeatureType.gpml_subduction_zone, - self.subduction_collision_parameters, - ) - ] - ) - # In addition to the default subduction detection, also detect continental collisions - collision_spec = _ContinentCollision( - # This filename string should not have a time formatted into it - this is - # taken care of later. - self.continent_mask_filepath, - default_collision, - verbose=False, + # Call the reconstruct-by-topologies object. + topology_reconstruction = ReconstructByTopologies( + self.plate_reconstruction.rotation_model, + self.plate_reconstruction.topology_features, + from_time, + self._min_time, + self._ridge_time_step, + points, + point_begin_times=appearance_times, + # This filename string should not have a time formatted into it - this is taken care of later... + continent_mask_filepath_format=self.continent_mask_filepath, ) - # Call the reconstruct by topologies object - if use_topological_model: - topology_reconstruction = _ReconstructByTopologicalModelImpl( - self.plate_reconstruction.rotation_model, - self.plate_reconstruction.topology_features, - from_time, - self._min_time, - self._ridge_time_step, - points, - point_begin_times=appearance_times, - detect_collisions=collision_spec, - ) - else: - topology_reconstruction = _ReconstructByTopologiesImplV2( - self.plate_reconstruction.rotation_model, - self.plate_reconstruction.topology_features, - from_time, - self._min_time, - self._ridge_time_step, - # This filename string should not have a time formatted into it - this is taken care of later... - self.continent_mask_filepath, - points, - point_begin_times=appearance_times, - ) - # Initialise the reconstruction. topology_reconstruction.begin_reconstruction() @@ -2014,9 +1976,6 @@ def _generate_debug_files_containing_reconstructed_ocean_seed_point_data_paralle def _build_and_reconstruct_ocean_seed_points_parallel( time, seafloor_grid, - use_topological_model, ): # Dispatch from parallel helper function to SeafloorGrid class method. - seafloor_grid._build_and_reconstruct_ocean_seed_points( - time, use_topological_model=use_topological_model - ) + seafloor_grid._build_and_reconstruct_ocean_seed_points(time) From 6ef25d0e62092dea8a380d42ca1ef10d497a940a Mon Sep 17 00:00:00 2001 From: John Cannon Date: Tue, 25 Nov 2025 17:12:24 +1100 Subject: [PATCH 09/17] Use a boolean array to track active seafloor seed points. And no longer return all the points (active and inactive). Now returns only the active points and their indices into the original points. --- gplately/lib/reconstruct_by_topologies.py | 293 +++++++++++----------- gplately/oceans.py | 57 ++--- 2 files changed, 167 insertions(+), 183 deletions(-) diff --git a/gplately/lib/reconstruct_by_topologies.py b/gplately/lib/reconstruct_by_topologies.py index c4c05876..bda55f54 100644 --- a/gplately/lib/reconstruct_by_topologies.py +++ b/gplately/lib/reconstruct_by_topologies.py @@ -79,33 +79,33 @@ def __init__( raise ValueError("'reconstruction_time_interval' must be positive.") # Reconstruction can go forward or backward in time. if self.reconstruction_begin_time > self.reconstruction_end_time: - self.reconstruction_time_step = -reconstruction_time_interval + self._reconstruction_time_step = -reconstruction_time_interval else: - self.reconstruction_time_step = reconstruction_time_interval + self._reconstruction_time_step = reconstruction_time_interval # Get number of times including end time if time span is a multiple of time step. # The '1' is because, for example, 2 time intervals is 3 times. # The '1e-6' deals with limited floating-point precision, eg, we want (3.0 - 0.0) / 1.0 to be 3.0 and not 2.999999 (which gets truncated to 2). - self.num_times = 1 + int( + self._num_times = 1 + int( math.floor( 1e-6 + float(self.reconstruction_end_time - self.reconstruction_begin_time) - / self.reconstruction_time_step + / self._reconstruction_time_step ) ) # It's possible the time step is larger than the time span, in which case we change it to equal the time span # unless the reconstruction begin and end times are equal (in which case there'll be only one reconstruction snapshot). # This guarantees there'll be at least one time step (which has two times; one at either end of interval). if ( - self.num_times == 1 + self._num_times == 1 and self.reconstruction_end_time != self.reconstruction_begin_time ): - self.num_times = 2 - self.reconstruction_time_step = ( + self._num_times = 2 + self._reconstruction_time_step = ( self.reconstruction_end_time - self.reconstruction_begin_time ) - self.reconstruction_time_interval = math.fabs(self.reconstruction_time_step) + self._reconstruction_time_interval = math.fabs(self._reconstruction_time_step) - self.last_time_index = self.num_times - 1 + self._last_time_index = self._num_times - 1 self.points = points self.num_points = len(points) @@ -132,7 +132,7 @@ def __init__( "Length of 'point_end_times' must match length of 'points'." ) - self.continent_mask_filepath_format = continent_mask_filepath_format + self._continent_mask_filepath_format = continent_mask_filepath_format def reconstruct(self): # Initialise the reconstruction. @@ -144,55 +144,79 @@ def reconstruct(self): while self.reconstruct_to_next_time(): pass - return self.get_active_current_points() + return self.get_current_active_points() def begin_reconstruction(self): - self.current_time_index = 0 + self._current_time_index = 0 - # Set up point arrays. - # Store active and inactive points here (inactive points have None in corresponding entries). - self.prev_points = [None] * self.num_points - self.curr_points = [None] * self.num_points - self.next_points = [None] * self.num_points + # Store active and inactive points here. + # + # NOTE: Deactivated points will NOT get reset to None - we'll just ignore them. + self._all_current_points = [None] * self.num_points + + # A boolean array indicating which points (in 'self._all_current_points') are active. + self._point_is_currently_active = np.zeros(self.num_points, dtype=bool) # Each point can only get activated once (after deactivation it cannot be reactivated). - self.point_has_been_activated = np.zeros(self.num_points, dtype=bool) - self.num_activated_points = 0 + self._point_has_been_activated = np.zeros(self.num_points, dtype=bool) # Activate any original points whose begin time is equal to (or older than) the newly updated current time. # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. self._activate_deactivate_points_using_their_begin_end_times() + # Deactivate any newly actived points that are masked by the continents. + self._deactivate_continent_masked_points() + def get_current_time_index(self): - return self.current_time_index + return self._current_time_index def get_current_time(self): return ( self.reconstruction_begin_time - + self.current_time_index * self.reconstruction_time_step + + self._current_time_index * self._reconstruction_time_step ) - def get_all_current_points(self): - return self.curr_points + def get_current_active_points(self): + """Return those points that are currently active.""" + return [ + self._all_current_points[point_index] + for point_index in self.get_current_active_point_indices() + ] - def get_active_current_points(self): - # Return only the active points (the ones that are not None). - return [point for point in self.get_all_current_points() if point is not None] + def get_current_active_point_indices(self): + """Return the indices of the currently active points (into the original points).""" + return np.where(self._point_is_currently_active)[0] def reconstruct_to_next_time(self): # If we're at the last time then there is no next time to reconstruct to. - if self.current_time_index == self.last_time_index: + if self._current_time_index == self._last_time_index: return False # If all points have been previously activated, but none are currently active then we're finished. # This means all points have entered their valid time range *and* either exited their time range or # have been deactivated (subducted forward in time or consumed by MOR backward in time). - if self.num_activated_points == self.num_points and not any(self.curr_points): + if np.all(self._point_has_been_activated) and not np.any( + self._point_is_currently_active + ): return False # Call the main implementation. self._reconstruct_to_next_time_impl() + # Move the current time to the next time. + self._current_time_index += 1 + + # Activate any original points whose begin time is equal to (or older than) the newly updated current time. + # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. + self._activate_deactivate_points_using_their_begin_end_times() + + # Deactivate any active points that are masked by the continents. + # + # This includes deactivating any newly activated points. + # + # Note: This is done at the newly updated current time. + self._deactivate_continent_masked_points() + # We successfully reconstructed to the next time. return True @@ -200,26 +224,19 @@ def _reconstruct_to_next_time_impl(self): current_time = self.get_current_time() # Positive/negative time step means reconstructing backward/forward in time. - next_time = current_time + self.reconstruction_time_step + next_time = current_time + self._reconstruction_time_step curr_topological_snapshot = pygplates.TopologicalSnapshot( # type: ignore self.topology_features, self.rotation_model, current_time ) curr_resolved_topologies = curr_topological_snapshot.get_resolved_topologies() - # Some of 'curr_points' will be None, so 'curr_active_points' contains only the active (not None) - # points, and 'curr_active_points_indices' is the same length as 'curr_active_points' but indexes - # into 'curr_points' so we can quickly find which point in 'curr_points' is associated with a - # particular point in 'curr_active_points'. - curr_active_points = [] - curr_active_points_indices = [] - for point_index, curr_point in enumerate(self.curr_points): - if curr_point is None: - # Current point is not currently active, so we cannot reconstruct it to the next point. - self.next_points[point_index] = None - continue - curr_active_points.append(curr_point) - curr_active_points_indices.append(point_index) + # Get the currently active points and their indices (into the original points). + curr_active_point_indices = self.get_current_active_point_indices() + curr_active_points = [ + self._all_current_points[point_index] + for point_index in curr_active_point_indices + ] # For each active current point find the resolved topology containing it. curr_active_point_resolved_topologies = ( @@ -233,7 +250,7 @@ def _reconstruct_to_next_time_impl(self): ) ) - # Map of current resolved topologies to the active points contained within them (their active point indices). + # Map of current resolved topologies to the active points contained within them (their point indices). map_curr_resolved_topology_to_active_point_indices = {} # Iterate over the resolved topologies containing all currently active points and get a list of the unique resolved topologies. @@ -242,33 +259,23 @@ def _reconstruct_to_next_time_impl(self): # In this case we'll just deactivate the point. Previously we would keep it active but just not reconstruct it # (ie, the current and next positions would be the same). However the point might end up in weird locations. # So it's best to just remove it. And this shouldn't happen very often at all (for topologies with global coverage). - curr_active_point_index = 0 - while curr_active_point_index < len(curr_active_points): + for index, curr_active_point_resolved_topology in enumerate( + curr_active_point_resolved_topologies + ): + + # Index into the active and inactive points 'self._all_current_points'. + curr_point_index = curr_active_point_indices[index] - curr_active_point_resolved_topology = curr_active_point_resolved_topologies[ - curr_active_point_index - ] # See if current active point fell outside all current resolved topologies. if curr_active_point_resolved_topology is None: - # Index into the active and inactive points 'self.curr_points'. - curr_point_index = curr_active_points_indices[curr_active_point_index] - # Current point is currently active but it fell outside all resolved topologies. # So we deactivate it. - self.next_points[curr_point_index] = None - - # Also remove evidence that the current point is active. - del curr_active_points[curr_active_point_index] - del curr_active_points_indices[curr_active_point_index] - del curr_active_point_resolved_topologies[curr_active_point_index] + self._point_is_currently_active[curr_point_index] = False # Continue to next active point. - # - # Note: We don't increment the active point index. - # We just removed an active point which essentially does the same thing. continue - # If resolved topology has not been encountered yet then give it an empty list of active point indices. + # If resolved topology has not been encountered yet then give it an empty list of point indices. if ( curr_active_point_resolved_topology not in map_curr_resolved_topology_to_active_point_indices @@ -280,10 +287,12 @@ def _reconstruct_to_next_time_impl(self): # Add the index of the currently active point to the resolved topology containing it. map_curr_resolved_topology_to_active_point_indices[ curr_active_point_resolved_topology - ].append(curr_active_point_index) + ].append(curr_point_index) - # Increment to the next active point. - curr_active_point_index += 1 + # Prevent usage since might no longer be representative of the active points because + # we might've just deactivated some points that fell outside all resolved topologies. + del curr_active_points + del curr_active_point_indices # The current resolved topologies that contain active points. curr_active_resolved_topologies = list( @@ -347,13 +356,9 @@ def _reconstruct_to_next_time_impl(self): if next_resolved_topology is None: # The current topology could not be resolved to the next time step. # So we deactivate all points that it contains. This generally shouldn't happen though. - for ( - curr_active_point_index - ) in curr_resolved_topology_active_point_indices: - curr_point_index = curr_active_points_indices[ - curr_active_point_index - ] - self.next_points[curr_point_index] = None + self._point_is_currently_active[ + curr_resolved_topology_active_point_indices + ] = False continue # Get the boundary polygon of the next resolved topology. @@ -367,17 +372,14 @@ def _reconstruct_to_next_time_impl(self): ).get_next_resolved_topology_boundary() # Iterate over the currently active points contained by the current resolved topology. - for curr_active_point_index in curr_resolved_topology_active_point_indices: + for point_index in curr_resolved_topology_active_point_indices: # Reconstruct the currently active point from its position at current time to its position at the next time step. - curr_active_point = curr_active_points[curr_active_point_index] + curr_active_point = self._all_current_points[point_index] next_active_point = curr_resolved_topology.reconstruct_point( curr_active_point, next_time ) - # Index into the active and inactive points 'self.curr_points'. - curr_point_index = curr_active_points_indices[curr_active_point_index] - # See if the location (of the current active point) reconstructed to the *next* time step # is inside the current topology resolved to the *next* time step. # @@ -387,92 +389,80 @@ def _reconstruct_to_next_time_impl(self): if not next_resolved_topology_boundary.is_point_in_polygon( next_active_point ): - self.next_points[curr_point_index] = None + self._point_is_currently_active[point_index] = False continue # The current active point was not consumed by a plate boundary. - self.next_points[curr_point_index] = next_active_point - - # See if any 'next' points collide with the continents and deactivate those that do. - self._detect_continent_collisions(next_time) - - # - # Set up for next loop iteration. - # - # Rotate previous, current and next point arrays. - # The new previous will be the old current. - # The new current will be the old next. - # The new next will be the old previous (but values are ignored and overridden in next time step; just re-using its memory). - self.prev_points, self.curr_points, self.next_points = ( - self.curr_points, - self.next_points, - self.prev_points, - ) - # - # Move the current time to the next time. - self.current_time_index += 1 - current_time = self.get_current_time() - - # Activate any original points whose begin time is equal to (or older than) the newly updated current time. - # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. - self._activate_deactivate_points_using_their_begin_end_times() + # So update its position. + self._all_current_points[point_index] = next_active_point def _activate_deactivate_points_using_their_begin_end_times(self): current_time = self.get_current_time() - # Iterate over all points and activate/deactivate as necessary depending on each point's valid time range. - for point_index in range(self.num_points): - if self.curr_points[point_index] is None: - if not self.point_has_been_activated[point_index]: - # Point is not active and has never been activated, so see if can activate it. - if ( - current_time <= self.point_begin_times[point_index] - and current_time >= self.point_end_times[point_index] - ): - # The initial point is assumed to be the position at the current time - # which is typically the point's begin time (approximately). - # But it could be the beginning of the reconstruction time span (specified in constructor) - # if that falls in the middle of the point's valid time range - in this case the - # initial point position is assumed to be in a position that is some time *after* - # it appeared (at its begin time) - and this can happen, for example, if you have a - # uniform grids of points at some intermediate time and want to see how they - # reconstruct to either a younger or older time (remembering that points can - # be subducted forward in time and consumed back into a mid-ocean ridge going - # backward in time). - self.curr_points[point_index] = self.points[point_index] - self.point_has_been_activated[point_index] = True - self.num_activated_points += 1 - else: - # Point is active, so see if can deactivate it. - if not ( - current_time <= self.point_begin_times[point_index] - and current_time >= self.point_end_times[point_index] - ): - self.curr_points[point_index] = None + # Get indices of points that are not active and have never been activated. + curr_inactive_and_not_yet_activated_point_indices = np.where( + ~self._point_is_currently_active & ~self._point_has_been_activated + )[0] + # See which of those points we can activate. + activate_point_indices = curr_inactive_and_not_yet_activated_point_indices[ + ( + current_time + <= self.point_begin_times[ + curr_inactive_and_not_yet_activated_point_indices + ] + ) + & ( + current_time + >= self.point_end_times[ + curr_inactive_and_not_yet_activated_point_indices + ] + ) + ] + # Activate those points. + self._point_is_currently_active[activate_point_indices] = True + # Copy original point into currently active points array. + for point_index in activate_point_indices: + self._all_current_points[point_index] = self.points[point_index] + # Mark those points as having been activated. + self._point_has_been_activated[activate_point_indices] = True + + # Get indices of points that are active (and hence must have been previously activated). + curr_active_point_indices = np.where(self._point_is_currently_active)[0] + # See which of those points we can deactivate. + deactivate_point_indices = curr_active_point_indices[ + ~( + (current_time <= self.point_begin_times[curr_active_point_indices]) + & (current_time >= self.point_end_times[curr_active_point_indices]) + ) + ] + # Deactivate those points. + self._point_is_currently_active[deactivate_point_indices] = False + # Leave deactivated points in the currently active points array (we'll just ignore them). - def _detect_continent_collisions(self, next_time): + def _deactivate_continent_masked_points(self): # If no continental collision detection requested then just return. - if self.continent_mask_filepath_format is None: + if self._continent_mask_filepath_format is None: return - # Get those 'next' points that are active. - next_active_points = [] - next_active_points_indices = [] - for point_index, next_point in enumerate(self.next_points): - if next_point: - next_active_points.append(next_point) - next_active_points_indices.append(point_index) + # Get the currently active points and their indices (into the original points). + curr_active_point_indices = self.get_current_active_point_indices() + curr_active_points = [ + self._all_current_points[point_index] + for point_index in curr_active_point_indices + ] # Convert points to lat/lon. - points_lat = np.empty(len(next_active_points)) - points_lon = np.empty(len(next_active_points)) - for active_point_index, point in enumerate(next_active_points): + curr_points_lat = np.empty(len(curr_active_points)) + curr_points_lon = np.empty(len(curr_active_points)) + for index, point in enumerate(curr_active_points): point_lat, point_lon = point.to_lat_lon() - points_lat[active_point_index] = point_lat - points_lon[active_point_index] = point_lon + curr_points_lat[index] = point_lat + curr_points_lon[index] = point_lon - # Read the continent mask at 'next_time'. - continent_mask_filepath = self.continent_mask_filepath_format.format(next_time) + # Read the continent mask at the current time. + continent_mask_filepath = self._continent_mask_filepath_format.format( + self.get_current_time() + ) gridZ, gridX, gridY = read_netcdf_grid( continent_mask_filepath, return_grids=True ) @@ -483,8 +473,8 @@ def _detect_continent_collisions(self, next_time): ymax = np.nanmax(gridY) # Sample continent mask grid, which is one over continents and zero over oceans. - points_i = (ni - 1) * ((points_lat - ymin) / (ymax - ymin)) - points_j = (nj - 1) * ((points_lon - xmin) / (xmax - xmin)) + points_i = (ni - 1) * ((curr_points_lat - ymin) / (ymax - ymin)) + points_j = (nj - 1) * ((curr_points_lon - xmin) / (xmax - xmin)) points_i_uint = np.rint(points_i).astype(np.uint) points_j_uint = np.rint(points_j).astype(np.uint) try: @@ -495,8 +485,9 @@ def _detect_continent_collisions(self, next_time): mask_values = gridZ[points_i, points_j] # Deactivate any points that sampled inside the continent mask. - for active_point_index in np.where(mask_values >= 0.5)[0]: - self.next_points[next_active_points_indices[active_point_index]] = None + self._point_is_currently_active[ + curr_active_point_indices[np.where(mask_values >= 0.5)[0]] + ] = False class _NextTopologicalFeatures(object): """The *next* topology boundary features and their topological section features. diff --git a/gplately/oceans.py b/gplately/oceans.py index ac1e0834..b2a562f9 100644 --- a/gplately/oceans.py +++ b/gplately/oceans.py @@ -1396,48 +1396,41 @@ def _build_and_reconstruct_ocean_seed_points(self, from_time): # This is the same size as 'points', which is either the initial ocean points at `time=max_time` or the MOR seed points at `time`. # Here 'current_points', despite being the same size, differs in that any *inactive* points at the current time # are None and the *active* points are at their reconstructed position at the current time. - current_points = topology_reconstruction.get_all_current_points() - - # Get the indices of the currently *active* points into all the points. - # - # Note: Points that are None represent currently *inactive* points. - # These are points that have either: - # 1) not been activated yet because the current time is older than their appearance time, or - # 2) have been deactivated because the current time is younger than their dissappearance time, or - # 3) have been deactivated through collision detection (ie, collided with a continent or trench). - # However, note that (2) does not apply here because the points currently don't have a hardwired disappearance time - # (it's implicitly '-inf'). Instead we rely on collision detection to deactivate points. - # And note that (1) does not apply here either, because the initial ocean seed points all exist at 'time=max_time' and - # all MOR seed points created at 'time' also exist at 'time'. - # - active_point_indices = np.fromiter( - ( - point_index - for point_index, point in enumerate(current_points) - if point is not None - ), - dtype=np.int32, - ) + active_points = topology_reconstruction.get_current_active_points() # logger.debug( - # f"At {time:.2f} Ma, {len(active_point_indices)} of the {len(points)} points created at {from_time:.2f} Ma are still active." + # f"At {time:.2f} Ma, {len(active_points)} of the {len(points)} points created at {from_time:.2f} Ma are still active." # ) # Store the reconstructed seed point data if any seed points are currently active. - if len(active_point_indices) > 0: + if active_points: # Latitudes and longitudes of the active points. # Store as 32-bit floating-point to save memory/disk-space. - active_latitudes = np.empty(len(active_point_indices), dtype=np.float32) - active_longitudes = np.empty( - len(active_point_indices), dtype=np.float32 + active_latitudes = np.empty(len(active_points), dtype=np.float32) + active_longitudes = np.empty(len(active_points), dtype=np.float32) + for index, active_point in enumerate(active_points): + active_latitudes[index], active_longitudes[index] = ( + active_point.to_lat_lon() + ) + + # Get the indices of the currently *active* points into all the points (active and inactive). + # + # Note: Points that are None represent currently *inactive* points. + # These are points that have either: + # 1) not been activated yet because the current time is older than their appearance time, or + # 2) have been deactivated because the current time is younger than their dissappearance time, or + # 3) have been deactivated through collision detection (ie, collided with a continent or trench). + # However, note that (2) does not apply here because the points currently don't have a hardwired disappearance time + # (it's implicitly '-inf'). Instead we rely on collision detection to deactivate points. + # And note that (1) does not really apply here either, because the initial ocean seed points all exist at the + # start time 'from_time=max_time' and all MOR seed points created at 'from_time' also exist at 'from_time'. + # + active_point_indices = ( + topology_reconstruction.get_current_active_point_indices() ) - for index, point_index in enumerate(active_point_indices): - active_latitudes[index], active_longitudes[index] = current_points[ - point_index - ].to_lat_lon() - # Store the active reconstructed seed point locations their active point indices for the current time. + # Store the active reconstructed seed point locations and their active point indices for the current time. # # Use a time "index" (associated with the current time - in the range `min_time <= time <= from_time`). time_index = topology_reconstruction.get_current_time_index() From d90ca40cecf8251f0a92bc79f5a0fa6f54c6e823 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Wed, 26 Nov 2025 10:39:38 +1100 Subject: [PATCH 10/17] Get topological snapshot from PlateReconstruction. More chance of reusing cached resolved topologies. --- gplately/lib/reconstruct_by_topologies.py | 28 ++++++++--------------- gplately/oceans.py | 3 +-- 2 files changed, 11 insertions(+), 20 deletions(-) diff --git a/gplately/lib/reconstruct_by_topologies.py b/gplately/lib/reconstruct_by_topologies.py index bda55f54..1c3e8c4f 100644 --- a/gplately/lib/reconstruct_by_topologies.py +++ b/gplately/lib/reconstruct_by_topologies.py @@ -48,8 +48,7 @@ class ReconstructByTopologies(object): def __init__( self, - rotation_features_or_model, - topology_features, + plate_reconstruction, reconstruction_begin_time, reconstruction_end_time, reconstruction_time_interval, @@ -65,12 +64,7 @@ def __init__( file path using ``continent_mask_filepath_format.format(time)``. """ - self.rotation_model = pygplates.RotationModel(rotation_features_or_model) # type: ignore - - # Turn topology data into a list of features (if not already). - self.topology_features = pygplates.FeaturesFunctionArgument( # type: ignore - topology_features - ).get_features() + self.plate_reconstruction = plate_reconstruction # Set up an array of reconstruction times covering the reconstruction time span. self.reconstruction_begin_time = reconstruction_begin_time @@ -226,10 +220,10 @@ def _reconstruct_to_next_time_impl(self): # Positive/negative time step means reconstructing backward/forward in time. next_time = current_time + self._reconstruction_time_step - curr_topological_snapshot = pygplates.TopologicalSnapshot( # type: ignore - self.topology_features, self.rotation_model, current_time - ) - curr_resolved_topologies = curr_topological_snapshot.get_resolved_topologies() + # Resolved the topologies at the current time. + curr_resolved_topologies = self.plate_reconstruction.topological_snapshot( + current_time, include_topological_slab_boundaries=False + ).get_resolved_topologies() # Get the currently active points and their indices (into the original points). curr_active_point_indices = self.get_current_active_point_indices() @@ -310,7 +304,7 @@ def _reconstruct_to_next_time_impl(self): next_topological_snapshot = pygplates.TopologicalSnapshot( # type: ignore next_topological_features.next_topological_section_features + next_topological_features.next_topological_boundary_features, - self.rotation_model, + self.plate_reconstruction.rotation_model, next_time, ) @@ -368,7 +362,7 @@ def _reconstruct_to_next_time_impl(self): next_resolved_topology_boundary = self._NextTopologicalBoundary( curr_resolved_topology, next_resolved_topology, - self.rotation_model, + self.plate_reconstruction.rotation_model, ).get_next_resolved_topology_boundary() # Iterate over the currently active points contained by the current resolved topology. @@ -1235,8 +1229,7 @@ def _post_split_polyline( def reconstruct_points_by_topologies( - rotation_features_or_model, - topology_features, + plate_reconstruction, reconstruction_begin_time, reconstruction_end_time, reconstruction_time_interval, @@ -1249,8 +1242,7 @@ def reconstruct_points_by_topologies( """Reconstruct points using topologies.""" topology_reconstruction = ReconstructByTopologies( - rotation_features_or_model, - topology_features, + plate_reconstruction, reconstruction_begin_time, reconstruction_end_time, reconstruction_time_interval, diff --git a/gplately/oceans.py b/gplately/oceans.py index b2a562f9..5f192938 100644 --- a/gplately/oceans.py +++ b/gplately/oceans.py @@ -1372,8 +1372,7 @@ def _build_and_reconstruct_ocean_seed_points(self, from_time): # Call the reconstruct-by-topologies object. topology_reconstruction = ReconstructByTopologies( - self.plate_reconstruction.rotation_model, - self.plate_reconstruction.topology_features, + self.plate_reconstruction, from_time, self._min_time, self._ridge_time_step, From a608a465c0191c8c08354dde88dca7b2e063dc99 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Wed, 26 Nov 2025 12:11:40 +1100 Subject: [PATCH 11/17] Prefer resolved networks over resolved plates (former overlays latter). --- gplately/lib/reconstruct_by_topologies.py | 58 +++++++++++++++++++---- 1 file changed, 48 insertions(+), 10 deletions(-) diff --git a/gplately/lib/reconstruct_by_topologies.py b/gplately/lib/reconstruct_by_topologies.py index 1c3e8c4f..b3f88a80 100644 --- a/gplately/lib/reconstruct_by_topologies.py +++ b/gplately/lib/reconstruct_by_topologies.py @@ -221,9 +221,17 @@ def _reconstruct_to_next_time_impl(self): next_time = current_time + self._reconstruction_time_step # Resolved the topologies at the current time. - curr_resolved_topologies = self.plate_reconstruction.topological_snapshot( + curr_topological_snapshot = self.plate_reconstruction.topological_snapshot( current_time, include_topological_slab_boundaries=False - ).get_resolved_topologies() + ) + + # Get the resolved networks and rigid plates. + curr_resolved_networks = curr_topological_snapshot.get_resolved_topologies( + pygplates.ResolveTopologyType.network # type:ignore + ) + curr_resolved_boundaries = curr_topological_snapshot.get_resolved_topologies( + pygplates.ResolveTopologyType.boundary # type:ignore + ) # Get the currently active points and their indices (into the original points). curr_active_point_indices = self.get_current_active_point_indices() @@ -232,18 +240,40 @@ def _reconstruct_to_next_time_impl(self): for point_index in curr_active_point_indices ] - # For each active current point find the resolved topology containing it. - curr_active_point_resolved_topologies = ( - _ptt.utils.points_in_polygons.find_polygons( + # For each active current point find the resolved topology network or rigid plate containing it. + curr_active_points_spatial_tree = ( + _ptt.utils.points_spatial_tree.PointsSpatialTree(curr_active_points) + ) + curr_active_point_resolved_networks = ( + _ptt.utils.points_in_polygons.find_polygons_using_points_spatial_tree( + curr_active_points, + curr_active_points_spatial_tree, + [ + resolved_topology.get_resolved_boundary() + for resolved_topology in curr_resolved_networks + ], + curr_resolved_networks, + ) + ) + curr_active_point_resolved_boundaries = ( + _ptt.utils.points_in_polygons.find_polygons_using_points_spatial_tree( curr_active_points, + curr_active_points_spatial_tree, [ resolved_topology.get_resolved_boundary() - for resolved_topology in curr_resolved_topologies + for resolved_topology in curr_resolved_boundaries ], - curr_resolved_topologies, + curr_resolved_boundaries, ) ) + if len(curr_active_point_resolved_networks) != len( + curr_active_point_resolved_boundaries + ): + raise RuntimeError( + "Unexpected length mismatch of arrays of resolved topologies containing points." + ) + # Map of current resolved topologies to the active points contained within them (their point indices). map_curr_resolved_topology_to_active_point_indices = {} @@ -253,13 +283,21 @@ def _reconstruct_to_next_time_impl(self): # In this case we'll just deactivate the point. Previously we would keep it active but just not reconstruct it # (ie, the current and next positions would be the same). However the point might end up in weird locations. # So it's best to just remove it. And this shouldn't happen very often at all (for topologies with global coverage). - for index, curr_active_point_resolved_topology in enumerate( - curr_active_point_resolved_topologies - ): + for index in range(len(curr_active_points)): # Index into the active and inactive points 'self._all_current_points'. curr_point_index = curr_active_point_indices[index] + # If point is inside a resolved *network* then prefer that, otherwise see if it's inside a rigid plate. + # This is because deforming networks can overlay rigid plates. + curr_active_point_resolved_topology = curr_active_point_resolved_networks[ + index + ] + if curr_active_point_resolved_topology is None: + curr_active_point_resolved_topology = ( + curr_active_point_resolved_boundaries[index] + ) + # See if current active point fell outside all current resolved topologies. if curr_active_point_resolved_topology is None: # Current point is currently active but it fell outside all resolved topologies. From 4fb6a683eeb4cbb10355a22589f0199baf3fa7b7 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Thu, 27 Nov 2025 23:12:43 +1100 Subject: [PATCH 12/17] No longer need subduction collision parameters in seafloor gridding. Due to improved collision detection added in 0582168. --- Notebooks/10-SeafloorGrids.ipynb | 21 +++------------------ gplately/oceans.py | 20 +++++++++----------- 2 files changed, 12 insertions(+), 29 deletions(-) diff --git a/Notebooks/10-SeafloorGrids.ipynb b/Notebooks/10-SeafloorGrids.ipynb index b74e6aac..a0662cb9 100644 --- a/Notebooks/10-SeafloorGrids.ipynb +++ b/Notebooks/10-SeafloorGrids.ipynb @@ -301,19 +301,8 @@ "\n", "One way an oceanic point can be deleted from the ocean basin at a certain timestep is if it is collides into continental crust. At each timestep continental crust is represented as a **continental mask** which is a binary grid that partitions continental regions from oceanic regions. By default, continental masks are generated from the reconstructed continental polygons (where the continent polygons features are specified with ``continent_polygon_features``). But they can also be generated using [continent contouring](https://github.com/EarthByte/continent-contouring) (by specifying ``use_continent_contouring=True``) or you can explicitly provide them as your own continent masks (by specifying them with ``continent_mask_filename``).\n", "\n", - "### Subduction parameters\n", - "Another way an oceanic point can be deleted from the ocean basin at a certain timestep is if it is approaching a plate boundary such that a **velocity and displacement test** is passed:\n", - "\n", - "#### Velocity test\n", - "First, given the trajectory of a point at the next time step, a point will cross a subducting plate boundary (and thus will be deleted) if the difference between its velocity on its current plate AND the velocity it will have on the other plate is greater than this velocity difference is higher than a specified **threshold delta velocity in kms/Myr**. \n", - "\n", - "#### Displacement test\n", - "If the proximity of the point's previous time position to the plate boundary it is approaching is higher than a set distance threshold (**threshold distance to boundary per Myr in kms/Myr**), then the point is far enough away from the boundary that it cannot be subducted or consumed by it, and \n", - "hence the point is still active.\n", - "\n", - "The parameter needed to encapsulate these thresholds is:\n", - "* **`subduction_collision_parameters`**: This a tuple with two elements, by default (5.0, 10.0), for the (threshold velocity delta in kms/my, threshold distance to boundary per My in kms/my)\n", - "***" + "### Subduction collision\n", + "Another way an oceanic point can be deleted from the ocean basin at a certain timestep is if it is approaching a convergent plate boundary. This case is now handled automatically without specifying any parameters." ] }, { @@ -323,8 +312,7 @@ "metadata": {}, "outputs": [], "source": [ - "continent_polygon_features=continents\n", - "subduction_collision_parameters=(5.0, 10.0)" + "continent_polygon_features=continents" ] }, { @@ -392,9 +380,6 @@ " refinement_levels = refinement_levels,\n", " initial_ocean_mean_spreading_rate = initial_ocean_mean_spreading_rate,\n", " \n", - " # Subduction parameters\n", - " subduction_collision_parameters = subduction_collision_parameters,\n", - " \n", " # Methodology parameters\n", " resume_from_checkpoints = resume_from_checkpoints,\n", "\n", diff --git a/gplately/oceans.py b/gplately/oceans.py index 5f192938..754f6fcc 100644 --- a/gplately/oceans.py +++ b/gplately/oceans.py @@ -104,13 +104,6 @@ class SeafloorGrid(object): Reconstruction by topologies involves determining which points are active and inactive (collided with a continent or subducted at a trench) for each reconstruction time step. This is done using a hidden object in :class:`PlateReconstruction` called ``ReconstructByTopologies``. - If an ocean point with a certain velocity on one plate ID transitions into another rigid plate ID at another timestep (with another velocity), - the velocity difference between both plates is calculated. The point may have subducted/collided with a continent if this velocity difference - is higher than a specified velocity threshold (which can be controlled with ``subduction_collision_parameters``). To ascertain whether the point - should be deactivated, a displacement test is conducted. If the proximity of the point's previous time position to the polygon boundary it is - approaching is higher than a set distance threshold, then the point is far enough away from the boundary that it cannot be subducted or consumed by it, - and hence the point is still active. Otherwise, it is deemed inactive and deleted from the ocean basin mesh. - With each reconstruction time step, points from mid-ocean ridges (which have more accurate spreading rates and attributed valid times) will spread across the ocean floor. Eventually, points will be pushed into continental boundaries or subduction zones, where they are deleted. Ideally, all initial ocean points (from the Stripy icosahedral mesh) should be deleted over time. However, not all will be deleted - such points typically reside near continental boundaries. @@ -193,6 +186,15 @@ def wrapper(self, *args, **kwargs): ) kwargs["continent_polygon_features"] = continent_polygon_features + # + # Handle the deprecated 'subduction_collision_parameters' argument. + # + if kwargs.pop("subduction_collision_parameters", None): + warnings.warn( + "`subduction_collision_parameters` keyword argument has been deprecated, it is no longer used", + DeprecationWarning, + ) + # # Handle the deprecated 'zval_names' argument. # @@ -220,7 +222,6 @@ def __init__( ridge_sampling: float = 0.5, extent: Tuple = (-180, 180, -90, 90), grid_spacing: float = 0.1, - subduction_collision_parameters=(5.0, 10.0), initial_ocean_mean_spreading_rate: float = 75.0, resume_from_checkpoints=False, continent_polygon_features=None, @@ -258,8 +259,6 @@ def __init__( grid_spacing : float, default=0.1 The degree spacing/interval with which to space grid points across all masking and final grids. If ``grid_spacing`` is provided, all grids will use it. If not, ``grid_spacing`` defaults to 0.1. - subduction_collision_parameters : len-2 tuple of float, default=(5.0, 10.0) - A 2-tuple of (threshold velocity delta in kms/my, threshold distance to boundary per My in kms/my) initial_ocean_mean_spreading_rate : float, default=75. A spreading rate to uniformly allocate to points that define the initial ocean basin. These points will have inaccurate ages, but most of them will be phased @@ -331,7 +330,6 @@ def __init__( # Topological parameters self.refinement_levels = refinement_levels self.ridge_sampling = ridge_sampling - self.subduction_collision_parameters = subduction_collision_parameters self.initial_ocean_mean_spreading_rate = initial_ocean_mean_spreading_rate # Gridding parameters From cee759a86a663b464550d806e5ea64d65521f0b3 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Wed, 3 Jun 2026 13:00:29 +1000 Subject: [PATCH 13/17] Renormalize line endings of Notebooks/10-SeafloorGrids.ipynb. So that subsequent diffs don't show the whole file as changed. --- Notebooks/10-SeafloorGrids.ipynb | 1268 +++++++++++++++--------------- 1 file changed, 634 insertions(+), 634 deletions(-) diff --git a/Notebooks/10-SeafloorGrids.ipynb b/Notebooks/10-SeafloorGrids.ipynb index a0662cb9..25828d85 100644 --- a/Notebooks/10-SeafloorGrids.ipynb +++ b/Notebooks/10-SeafloorGrids.ipynb @@ -1,634 +1,634 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "6a6c21c7", - "metadata": {}, - "source": [ - "## 10 - Seafloor Grids\n", - "\n", - "An adaptation of [agegrid-01](https://github.com/siwill22/agegrid-0.1) written by Simon Williams, Nicky Wright and John Cannon for gridding general z-values onto seafloor basin points using GPlately.\n", - "\n", - "This notebook demonstrates how to generate seafloor age and spreading-rate grids through time. The figure below, generated by this notebook, shows seafloor age at 65 Ma. You may adjust the time parameter to generate grids at different geological times.\n", - "\n", - "![Figure: seafloor_age_65Ma.png](https://gplates.github.io/gplately/latest/sphinx/html/_images/seafloor_age_65Ma.png)\n", - "\n", - "### Citation:\n", - "Simon Williams, Nicky M. Wright, John Cannon, Nicolas Flament, R. Dietmar Müller, Reconstructing seafloor age distributions in lost ocean basins, Geoscience Frontiers, Volume 12, Issue 2, 2021, Pages 769-780, ISSN 1674-9871,\n", - "https://doi.org/10.1016/j.gsf.2020.06.004." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "db118eb4", - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "#\n", - "# Setting GPLATELY_DEBUG to a value above 100 will create the following debug files\n", - "# during seafloor gridding (that can be loaded into GPlates):\n", - "#\n", - "# - A separate debug file for each time (from 'max_time' to 'min_time) containing the seed points\n", - "# reconstructed from that time to 'min_time'.\n", - "# * The debug file with \"initial_ocean_basin\" in the filename contains reconstructions of the initial ocean basin seed points, and\n", - "# * each debug file with \"mid_ocean_ridge\" in the filename contains reconstructions of the seed points created along mid-ocean ridges\n", - "# that formed at that time specified in the filename.\n", - "# - One big debug file containing the reconstructions of ALL seed points created at ALL times.\n", - "#\n", - "# These files are located in the 'debug' sub-directory of the save directory.\n", - "#\n", - "#os.environ[\"GPLATELY_DEBUG\"] = \"200\"\n", - "\n", - "import gplately\n", - "\n", - "import pygplates\n", - "import glob\n", - "import matplotlib.pyplot as plt\n", - "import cartopy.crs as ccrs\n", - "from plate_model_manager import PlateModelManager" - ] - }, - { - "cell_type": "markdown", - "id": "9a4bddee", - "metadata": {}, - "source": [ - "### Define a rotation model, topology features and continents for the `PlateReconstruction` model\n", - "There are two ways to do this. To use local files, set `use_local_files = True`. To use `PlateModelManager`, set `use_local_files = False`." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b571bf94", - "metadata": {}, - "outputs": [], - "source": [ - "use_local_files = False" - ] - }, - { - "cell_type": "markdown", - "id": "4936baad", - "metadata": {}, - "source": [ - "#### 1) Manually pointing to files" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e5319664", - "metadata": {}, - "outputs": [], - "source": [ - "if use_local_files:\n", - " # Method 1: manually point to files\n", - " # you can download the model files at \n", - " # https://www.earthbyte.org/webdav/ftp/Data_Collections/Merdith_etal_2021_ESR/SM2-Merdith_et_al_1_Ga_reconstruction_v1.2.4.zip\n", - " input_directory = \"./SM2-Merdith_et_al_1_Ga_reconstruction_v1.2.4\"\n", - " \n", - " rotation_model = glob.glob(os.path.join(input_directory, '*.rot'))\n", - " static_polygons = input_directory+\"/shapes_static_polygons_Merdith_et_al.gpml\"\n", - " topology_features = [\n", - " input_directory+\"/250-0_plate_boundaries_Merdith_et_al.gpml\",\n", - " input_directory+\"/410-250_plate_boundaries_Merdith_et_al.gpml\",\n", - " input_directory+\"/1000-410-Convergence_Merdith_et_al.gpml\",\n", - " input_directory+\"/1000-410-Divergence_Merdith_et_al.gpml\",\n", - " input_directory+\"/1000-410-Topologies_Merdith_et_al.gpml\",\n", - " input_directory+\"/1000-410-Transforms_Merdith_et_al.gpml\",\n", - " input_directory+\"/TopologyBuildingBlocks_Merdith_et_al.gpml\"\n", - " ]\n", - " \n", - " continents = input_directory+\"/shapes_continents.gpml\"\n", - " coastlines = input_directory+\"/shapes_coastlines_Merdith_et_al_v2.gpmlz\"\n", - " COBs = None\n", - "\n", - " from pathlib import Path\n", - " if not Path(input_directory).is_dir():\n", - " raise FileNotFoundError(f\"The input directory does not exist: {input_directory}. Ensure your files are in the correct locations or download the example input files at https://www.earthbyte.org/webdav/ftp/Data_Collections/Merdith_etal_2021_ESR/SM2-Merdith_et_al_1_Ga_reconstruction_v1.2.4.zip\")\n", - " " - ] - }, - { - "cell_type": "markdown", - "id": "1519da7c", - "metadata": {}, - "source": [ - "#### 2) Using `PlateModelManager`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3ebf9727", - "metadata": {}, - "outputs": [], - "source": [ - "if not use_local_files:\n", - " # Method 2: use PlateModelManager\n", - " pm_manager = PlateModelManager()\n", - " plate_model = pm_manager.get_model(\"Merdith2021\", data_dir=\"plate-model-repo\")\n", - " \n", - " rotation_model = plate_model.get_rotation_model()\n", - " topology_features = plate_model.get_topologies()\n", - " static_polygons = plate_model.get_static_polygons()\n", - " \n", - " continents = plate_model.get_layer('ContinentalPolygons') " - ] - }, - { - "cell_type": "markdown", - "id": "3560443e", - "metadata": {}, - "source": [ - "## The SeafloorGrid object\n", - "...is a collection of methods to generate seafloor grids.\n", - "\n", - "The critical input parameters are:\n", - "\n", - "### Plate model parameters\n", - "* **`plate_reconstruction`**: The gplately `PlateReconstruction` object defined in the cell below. This object is a collection of methods for calculating plate tectonic stats through geological time.\n", - "* **`plot_topologies`**: The gplately `PlotTopologies` object defined in the cell below. This object is a collection of methods for resolving geological features needed for gridding to a certain geological time.\n", - "\n", - "#### Define the `PlateReconstruction` and `PlotTopologies` objects" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "12bec53e", - "metadata": {}, - "outputs": [], - "source": [ - "model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)\n", - "\n", - "gplot = gplately.PlotTopologies(model, continents=continents)" - ] - }, - { - "cell_type": "markdown", - "id": "57d6e8b2", - "metadata": {}, - "source": [ - "### Time parameters\n", - "* **`max_time`**: The first time step to start generating age grids. This is the time step where we define certain **gridding initial conditions**, which will be explained below.\n", - "* **`min_time`**: The final step to generate age grids. This is the time when recursive reconstructions starting from `max_time` stop.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "17add2d7", - "metadata": {}, - "outputs": [], - "source": [ - "max_time = 410.\n", - "min_time = 65." - ] - }, - { - "cell_type": "markdown", - "id": "703d2315", - "metadata": {}, - "source": [ - "### Ridge resolution parameters\n", - "With each reconstruction time step, mid-ocean ridge segments emerge and spread across the ocean floor. In `gplately`, these ridges are lines, or a tessellation of infinitesimal points. The spatio-temporal resolution of these ridge points can be controlled by two parameters:\n", - "* **`ridge_time_step`**: The \"delta time\" or time increment between successive resolutions of ridges, and hence successive grids. By default this is 1Myr, so grids are produced every millionth year.\n", - "* **`ridge_sampling`**: This controls the geographical resolution (in degrees) with which ridge lines are partitioned into points. The larger this is, the larger the spacing between successive ridge points, and hence the smaller the number of ridge points at each timestep. By default this is 0.5 degrees, so points are spaced 0.5 degrees apart. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b959a7b4", - "metadata": {}, - "outputs": [], - "source": [ - "ridge_time_step = 5. # We set the time step to 5 Myr to reduce the notebook’s runtime.\n", - "ridge_sampling = 0.5" - ] - }, - { - "cell_type": "markdown", - "id": "c1c20f7d", - "metadata": {}, - "source": [ - "***\n", - "### Grid resolution parameters\n", - "The ridge resolution parameters mentioned above take care of the spatio-temporal positioning of ridge features/points. However, these ridge points need to be interpolated onto regular grids - separate parameters are needed for these grids. \n", - "\n", - "#### For the `max_time` initial condition\n", - "At `max_time`, i.e. 410Ma, the `PlateReconstruction` model has not been initialised at 411Ma and any time(s) before that to sculpt the geological history before `max_time`. In terms of the workflow, all geological history before `max_time` is unknown. Thus, the global seafloor point distriubtion, and each point's spreading rate and age must be manually defined. \n", - "\n", - "The initial point distribution is an icosahedral point mesh, made using [`stripy`](https://github.com/underworldcode/stripy/tree/294354c00dd72e085a018e69c345d9353c6fafef).\n", - "\n", - "* **`refinement_levels`**: A unitless integer that controls the number of points in the `max_time` icosahedral point mesh.\n", - " \n", - " 6 is the default. Any higher will result in a finer mesh.\n", - " \n", - " \n", - "* **`initial_ocean_mean_spreading_rate`** (in units of mm/yr or km/myr): Since the geological history before and at `max_time` is unknown, we will need to manually define the spreading rates and ages of all ocean points. This is manually set to a uniform spreading rate of 75 (mm/yr or km/myr). Each point's age is equal to its proximity to the nearest mid ocean ridge (assuming that ridge is the source of the point) divided by half this uniform spreading rate (half to account for spreading to the left and right of a ridge). \n", - "\n", - "#### For the interpolated regular grids\n", - "The regular grid on which all data is interpolated has a resolution that can be controlled by:\n", - "\n", - "* **`grid_spacing`**: The degree spacing between successive nodes in the grid. By default, this is 0.1 degrees. **Acceptable degree spacings are 0.1, 0.25, 0.5, 0.75, or 1 degree(s)** because these allow cleanly divisible grid nodes across global latitude and longitude extents. Anything greater than 1 will be rounded to 1, and anything between these acceptable spacings will be rounded to the nearest acceptable spacing. \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5ec00f36", - "metadata": {}, - "outputs": [], - "source": [ - "refinement_levels = 6\n", - "initial_ocean_mean_spreading_rate = 50.\n", - "\n", - "# Gridding parameter\n", - "grid_spacing = 0.25\n", - "\n", - "extent=[-180,180,-90,90]" - ] - }, - { - "cell_type": "markdown", - "id": "8464bfb4", - "metadata": {}, - "source": [ - "***\n", - "### File saving and naming parameters\n", - "When grids are produced, they are saved to:\n", - "* **`save_directory`**: A string to a directory that must exist already.\n", - "\n", - "These files are named according to:\n", - "* **`file_collection`**: A string used to help with naming all output files. \n", - "***" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "106b61f5", - "metadata": {}, - "outputs": [], - "source": [ - "# continent masks, initial ocean seed points, and gridding input files are kept here\n", - "output_parent_directory = os.path.join(\n", - " \"NotebookFiles\",\n", - " \"Notebook10\",\n", - ")\n", - "save_directory = os.path.join(\n", - " output_parent_directory,\n", - " \"seafloor_grid_output\",\n", - ")\n", - "print(f\"The ouput files created by this notebook will be saved in folder {save_directory}.\" )\n", - "\n", - "# A string to help name files according to a plate model \"Merdith2021\"\n", - "file_collection = \"Merdith2021\"\n" - ] - }, - { - "cell_type": "markdown", - "id": "a328fa6d", - "metadata": {}, - "source": [ - "***\n", - "### Continent collision\n", - "\n", - "One way an oceanic point can be deleted from the ocean basin at a certain timestep is if it is collides into continental crust. At each timestep continental crust is represented as a **continental mask** which is a binary grid that partitions continental regions from oceanic regions. By default, continental masks are generated from the reconstructed continental polygons (where the continent polygons features are specified with ``continent_polygon_features``). But they can also be generated using [continent contouring](https://github.com/EarthByte/continent-contouring) (by specifying ``use_continent_contouring=True``) or you can explicitly provide them as your own continent masks (by specifying them with ``continent_mask_filename``).\n", - "\n", - "### Subduction collision\n", - "Another way an oceanic point can be deleted from the ocean basin at a certain timestep is if it is approaching a convergent plate boundary. This case is now handled automatically without specifying any parameters." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "68450bf6", - "metadata": {}, - "outputs": [], - "source": [ - "continent_polygon_features=continents" - ] - }, - { - "cell_type": "markdown", - "id": "bd978b09", - "metadata": {}, - "source": [ - "***\n", - "### Methodology parameters\n", - "* **`resume_from_checkpoints`**: Gridding can take a while (depending on the number of CPUs used). If this is changed to `True`, after gridding was interrupted, then rerunning the gridding cell will resume from where it left off.\n", - "* **`nprocs`**: The number of CPUs to use for parts of the code that are parallelized. The default is `-2` which uses all CPUs except one (to keep the system responsive).\n", - "***" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b1a6fbee", - "metadata": {}, - "outputs": [], - "source": [ - "# Methodology parameters\n", - "resume_from_checkpoints = False\n", - "nprocs = -2" - ] - }, - { - "cell_type": "markdown", - "id": "5038c7e6", - "metadata": {}, - "source": [ - "Use all parameters to define `SeafloorGrid`:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ed5fc9e9", - "metadata": {}, - "outputs": [], - "source": [ - "# The SeafloorGrid object with all aforementioned parameters\n", - "seafloorgrid = gplately.SeafloorGrid(\n", - " \n", - " plate_reconstruction = model, \n", - " \n", - " # Time parameters\n", - " max_time = max_time,\n", - " min_time = min_time,\n", - " \n", - " # Ridge tessellation parameters\n", - " ridge_time_step = ridge_time_step,\n", - " ridge_sampling = ridge_sampling,\n", - " \n", - " # Gridding parameters\n", - " grid_spacing = grid_spacing,\n", - "\n", - " extent = extent,\n", - " \n", - " # Naming parameters\n", - " save_directory = save_directory,\n", - " file_collection = file_collection,\n", - " \n", - " # Initial condition parameters\n", - " refinement_levels = refinement_levels,\n", - " initial_ocean_mean_spreading_rate = initial_ocean_mean_spreading_rate,\n", - " \n", - " # Methodology parameters\n", - " resume_from_checkpoints = resume_from_checkpoints,\n", - "\n", - " # Continental polygons.\n", - " continent_polygon_features = continent_polygon_features,\n", - "\n", - " nprocs = nprocs\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "ab3cbf87", - "metadata": {}, - "source": [ - "## How to use `SeafloorGrid`\n", - "\n", - "### 1) Run `seafloorgrid.reconstruct_by_topologies()`\n", - "Running `seafloorgrid.reconstruct_by_topologies()` prepares all ocean basin points and their seafloor ages and spreading rates for gridding, and reconstructs all active points per timestep to form the grids per timestep.\n", - "\n", - "#### Continental masks\n", - "A netCDF4 continent mask is produced for each _time_ from `max_time` to `min_time`.\n", - "\n", - "#### Initial ocean basin\n", - "At `max_time`, an initial ocean seed point icosahedral mesh fills the ocean basin. Each initial seed point is allocated an age based on its distance from the nearest mid-ocean ridge. And these are reconstructed to `min_time` to produce an `.npz` file (containing reconstructed seed point data for seed points created at `max_time`).\n", - "\n", - "#### Mid-ocean ridges\n", - "At each successive timestep (`ridge_time_step`), new seed points emerge from spreading ridges, and they are allocated their own spreading rates. The seed points emerging at a particular _time_ are reconstructed to `min_time` to produce an `.npz` file (containing reconstructed seed point data for seed points created at that _time_). This results in one `.npz` file for each _time_ from `max_time - ridge_time_step` to `min_time`.\n", - "\n", - "#### Seafloor age and spreading rate\n", - "After all seed points (initial and mid-ocean ridge) have been reconstructed to `min_time`, then one `.npz` file is generated for each _time_ from `max_time` to `min_time`. These files contain the gridding input that will later be used to generate the seafloor age and spreading rate grids.\n", - "\n", - "#### Resuming an interrupted gridding run\n", - "*If the `resume_from_checkpoints` parameter is passed as `True`, after `seafloorgrid.reconstruct_by_topologies()` was interrupted, then you can rerun the cell. The workflow will pick up from where it left off (provided the continent mask files and reconstructed seed point data files have not been erased).*\n", - "\n", - "To overwrite all files in `save_directory`, pass `resume_from_checkpoints` as `False` (this is the default behaviour).\n", - "\n", - "#### Reconstruct by topologies\n", - "The `ReconstructByTopologies` object (written by Simon Williams, Nicky Wright, and John Cannon) handles the reconstruction of seed points. `ReconstructByTopologies` (RBT) identifies active points on the ocean basin per timestep. It works as follows:\n", - "\n", - "If an ocean point on one plate ID transitions into another rigid plate ID at the next timestep, RBT calculates the point's velocity difference between both plates. The point **may** have subducted/collided with a continent at this boundary if this velocity difference is higher than a set velocity threshold. To ascertain whether the point should indeed be deactivated, a second test is conducted: RBT checks the previous time position of the point and calculates this point’s proximity to the boundary of the plate ID polygon it is approaching. If this distance is higher than a set distance threshold, then the point is far enough away from the boundary that it cannot be subducted or consumed by it and hence the point is still active. Else, it is deactivated/deleted.\n", - "\n", - "Once all active points and their z-values are identified, they are written to the gridding input file (.npz) for that timestep.\n", - "\n", - "Below is an example of the initial condition: the ocean basin is populated with an `intial_mean_ocean_spreading_rate` of 50 mm/yr at a `max_time` of 410Ma. Reconstruction over 10 Myr to 400 Ma sees new seed points emerging from ridges with their own plate-model-ascribed spreading rates.\n", - "\n", - "![init_condition](https://raw.githubusercontent.com/GPlates/gplately/master/Notebooks/NotebookFiles/Notebook10/Merdith2021_initial_sr_conditions.png)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "143ea077", - "metadata": {}, - "outputs": [], - "source": [ - "import datetime\n", - "import time\n", - "\n", - "start = time.time()\n", - "seafloorgrid.reconstruct_by_topologies()\n", - "end = time.time()\n", - "duration = datetime.timedelta(seconds=end - start)\n", - "print(\"Duration: {}\".format(duration))" - ] - }, - { - "cell_type": "markdown", - "id": "a0a0ca60", - "metadata": {}, - "source": [ - "### 2) Run `seafloorgrid.lat_lon_z_to_netCDF` to write grids to netCDF\n", - "Calling `seafloorgrid.lat_lon_z_to_netCDF` grids one set of z-data per latitude-longitude pair from each timestep's gridding input file (produced in `seafloorgrid.reconstruct_by_topologies()`). Grids are in netCDF format.\n", - "\n", - "The desired z-data to grid is identified using a `zval_name`. \n", - "For example, seafloor age grids can be produced using `SEAFLOOR_AGE`, and spreading rate grids are `SPREADING_RATE`. \n", - "\n", - "Use `unmasked = True` to output both the masked and unmasked versions of the grids." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "891183f3", - "metadata": {}, - "outputs": [], - "source": [ - "seafloorgrid.lat_lon_z_to_netCDF(\"SEAFLOOR_AGE\", unmasked=True)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "d51a41a5", - "metadata": {}, - "outputs": [], - "source": [ - "seafloorgrid.lat_lon_z_to_netCDF(\"SPREADING_RATE\", unmasked=True)" - ] - }, - { - "cell_type": "markdown", - "id": "8213b46a", - "metadata": {}, - "source": [ - "### Plotting a sample age grid and spreading rate grid\n", - "Read one netCDF grid using GPlately's `Raster` object from `grids`, and plot it using the `PlotTopologies` object.\n", - "\n", - "Notice the evolution of seafloor spreading rate from the initial value set with `initial_ocean_mean_spreading_rate`. Eventually, this initial uniform spreading rate will be phased out with sufficient recursive reconstruction." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "06a92ce9-ecea-49b9-9c99-c33438e75379", - "metadata": {}, - "outputs": [], - "source": [ - "time = min_time\n", - "# set the plot_unmasked_grid variable to True if you would like to see the unmasked age grid.\n", - "plot_unmasked_grid = False\n", - "\n", - "# age grids\n", - "if not plot_unmasked_grid:\n", - " agegrid_filename = f\"{seafloorgrid.save_directory}/SEAFLOOR_AGE/{seafloorgrid.file_collection}_SEAFLOOR_AGE_grid_{time:0.2f}Ma.nc\"\n", - "else:\n", - " agegrid_filename = f\"{seafloorgrid.save_directory}/SEAFLOOR_AGE/{seafloorgrid.file_collection}_SEAFLOOR_AGE_grid_unmasked_{time:0.2f}Ma.nc\"\n", - "\n", - "age_grid = gplately.grids.Raster(agegrid_filename)\n", - "\n", - "# Prepare plots\n", - "fig = plt.figure(figsize=(8,4), dpi=240, linewidth=1)\n", - "ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))\n", - "\n", - "gplot.time = time\n", - "plt.title(f\"Seafloor Age at {int(time)} Ma (Model:{seafloorgrid.file_collection})\")\n", - "\n", - "agegrid_cmap = \"YlGnBu\"\n", - "agegrid_cpt_file = \"NotebookFiles/agegrid.cpt\"\n", - "if not os.path.isfile(agegrid_cpt_file):\n", - " print(f\"The file {agegrid_cpt_file} does not exist locally. It can be found at https://raw.githubusercontent.com/GPlates/gplately/refs/heads/master/Notebooks/NotebookFiles/agegrid.cpt\")\n", - " print(f\"The '{agegrid_cmap}' will be used to plot the age grid.\")\n", - " agegrid_cpt_file = None\n", - "if agegrid_cpt_file:\n", - " try:\n", - " from gplately.mapping.gmt_cpt import get_cmap_from_gmt_cpt\n", - " agegrid_cmap=get_cmap_from_gmt_cpt(agegrid_cpt_file)\n", - " except ImportError as e:\n", - " print(\"To use GMT cpt file, ensure you are using GPlately version 2.1 or higher.\")\n", - " print(e)\n", - "\n", - "im = gplot.plot_grid(\n", - " ax, \n", - " age_grid.data, \n", - " cmap = agegrid_cmap,\n", - " vmin = 0, \n", - " vmax =410,\n", - ")\n", - "gplot.plot_ridges(ax,linewidth=1)\n", - "if not plot_unmasked_grid:\n", - " gplot.plot_continents(ax, edgecolor=\"grey\", facecolor=\"lightgrey\", linewidth=0.5)\n", - "plt.colorbar(im, label=\"Seafloor Age (Myr)\", shrink=0.5, pad=0.05)\n", - "\n", - "# Save figure\n", - "plt.savefig(\n", - " f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_age_{int(time)}Ma.png\",\n", - " bbox_inches = 'tight',\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "56e03e39", - "metadata": {}, - "outputs": [], - "source": [ - "time = min_time\n", - "# set the plot_unmasked_grid variable to True if you would like to see the unmasked spreading rate grid.\n", - "plot_unmasked_grid = False\n", - "\n", - "# spreading rate grids\n", - "if not plot_unmasked_grid:\n", - " srgrid_filename = f\"{seafloorgrid.save_directory}/SPREADING_RATE/{seafloorgrid.file_collection}_SPREADING_RATE_grid_{time:0.2f}Ma.nc\"\n", - "else:\n", - " srgrid_filename = f\"{seafloorgrid.save_directory}/SPREADING_RATE/{seafloorgrid.file_collection}_SPREADING_RATE_grid_unmasked_{time:0.2f}Ma.nc\"\n", - "\n", - "sr_grid = gplately.grids.Raster(srgrid_filename)\n", - "\n", - "# Prepare plots\n", - "fig = plt.figure(figsize=(8,4), dpi=240, linewidth=1)\n", - "ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))\n", - "\n", - "gplot.time = time\n", - "plt.title(f\"Seafloor Spreading Rate at {int(time)} Ma (Model:{seafloorgrid.file_collection})\")\n", - "\n", - "spreading_rate_cmap = \"magma\"\n", - "spreading_rate_cpt_file = \"NotebookFiles/spreading_full_rate.cpt\"\n", - "if not os.path.isfile(spreading_rate_cpt_file):\n", - " print(f\"The file {spreading_rate_cpt_file} does not exist locally. It can be found at https://raw.githubusercontent.com/GPlates/gplately/refs/heads/master/Notebooks/NotebookFiles/spreading_full_rate.cpt\")\n", - " print(f\"The '{spreading_rate_cmap}' will be used to plot the spreading rate grid.\")\n", - " spreading_rate_cpt_file = None\n", - "if spreading_rate_cpt_file:\n", - " try:\n", - " from gplately.mapping.gmt_cpt import get_cmap_from_gmt_cpt\n", - " spreading_rate_cmap=get_cmap_from_gmt_cpt(spreading_rate_cpt_file)\n", - " except ImportError as e:\n", - " print(\"To use GMT cpt file, ensure you are using GPlately version 2.1 or higher.\")\n", - " print(e)\n", - "\n", - "im = gplot.plot_grid(\n", - " ax, \n", - " sr_grid.data, \n", - " cmap = spreading_rate_cmap,\n", - " vmin = 0, \n", - " vmax =200,\n", - ")\n", - "gplot.plot_ridges(ax, linewidth=1)\n", - "if not plot_unmasked_grid:\n", - " gplot.plot_continents(ax, edgecolor=\"grey\", facecolor=\"lightgrey\", linewidth=0.5)\n", - "plt.colorbar(im, label=\"Spreading rate (mm/yr)\", shrink=0.5, pad=0.05, ticks=[0,50,100,150,200])\n", - "\n", - "# Save figure\n", - "plt.savefig(\n", - " f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_spreading_rate_{int(time)}Ma.png\",\n", - " bbox_inches = 'tight',\n", - ")" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.11.11" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} +{ + "cells": [ + { + "cell_type": "markdown", + "id": "6a6c21c7", + "metadata": {}, + "source": [ + "## 10 - Seafloor Grids\n", + "\n", + "An adaptation of [agegrid-01](https://github.com/siwill22/agegrid-0.1) written by Simon Williams, Nicky Wright and John Cannon for gridding general z-values onto seafloor basin points using GPlately.\n", + "\n", + "This notebook demonstrates how to generate seafloor age and spreading-rate grids through time. The figure below, generated by this notebook, shows seafloor age at 65 Ma. You may adjust the time parameter to generate grids at different geological times.\n", + "\n", + "![Figure: seafloor_age_65Ma.png](https://gplates.github.io/gplately/latest/sphinx/html/_images/seafloor_age_65Ma.png)\n", + "\n", + "### Citation:\n", + "Simon Williams, Nicky M. Wright, John Cannon, Nicolas Flament, R. Dietmar Müller, Reconstructing seafloor age distributions in lost ocean basins, Geoscience Frontiers, Volume 12, Issue 2, 2021, Pages 769-780, ISSN 1674-9871,\n", + "https://doi.org/10.1016/j.gsf.2020.06.004." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db118eb4", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "#\n", + "# Setting GPLATELY_DEBUG to a value above 100 will create the following debug files\n", + "# during seafloor gridding (that can be loaded into GPlates):\n", + "#\n", + "# - A separate debug file for each time (from 'max_time' to 'min_time) containing the seed points\n", + "# reconstructed from that time to 'min_time'.\n", + "# * The debug file with \"initial_ocean_basin\" in the filename contains reconstructions of the initial ocean basin seed points, and\n", + "# * each debug file with \"mid_ocean_ridge\" in the filename contains reconstructions of the seed points created along mid-ocean ridges\n", + "# that formed at that time specified in the filename.\n", + "# - One big debug file containing the reconstructions of ALL seed points created at ALL times.\n", + "#\n", + "# These files are located in the 'debug' sub-directory of the save directory.\n", + "#\n", + "#os.environ[\"GPLATELY_DEBUG\"] = \"200\"\n", + "\n", + "import gplately\n", + "\n", + "import pygplates\n", + "import glob\n", + "import matplotlib.pyplot as plt\n", + "import cartopy.crs as ccrs\n", + "from plate_model_manager import PlateModelManager" + ] + }, + { + "cell_type": "markdown", + "id": "9a4bddee", + "metadata": {}, + "source": [ + "### Define a rotation model, topology features and continents for the `PlateReconstruction` model\n", + "There are two ways to do this. To use local files, set `use_local_files = True`. To use `PlateModelManager`, set `use_local_files = False`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b571bf94", + "metadata": {}, + "outputs": [], + "source": [ + "use_local_files = False" + ] + }, + { + "cell_type": "markdown", + "id": "4936baad", + "metadata": {}, + "source": [ + "#### 1) Manually pointing to files" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e5319664", + "metadata": {}, + "outputs": [], + "source": [ + "if use_local_files:\n", + " # Method 1: manually point to files\n", + " # you can download the model files at \n", + " # https://www.earthbyte.org/webdav/ftp/Data_Collections/Merdith_etal_2021_ESR/SM2-Merdith_et_al_1_Ga_reconstruction_v1.2.4.zip\n", + " input_directory = \"./SM2-Merdith_et_al_1_Ga_reconstruction_v1.2.4\"\n", + " \n", + " rotation_model = glob.glob(os.path.join(input_directory, '*.rot'))\n", + " static_polygons = input_directory+\"/shapes_static_polygons_Merdith_et_al.gpml\"\n", + " topology_features = [\n", + " input_directory+\"/250-0_plate_boundaries_Merdith_et_al.gpml\",\n", + " input_directory+\"/410-250_plate_boundaries_Merdith_et_al.gpml\",\n", + " input_directory+\"/1000-410-Convergence_Merdith_et_al.gpml\",\n", + " input_directory+\"/1000-410-Divergence_Merdith_et_al.gpml\",\n", + " input_directory+\"/1000-410-Topologies_Merdith_et_al.gpml\",\n", + " input_directory+\"/1000-410-Transforms_Merdith_et_al.gpml\",\n", + " input_directory+\"/TopologyBuildingBlocks_Merdith_et_al.gpml\"\n", + " ]\n", + " \n", + " continents = input_directory+\"/shapes_continents.gpml\"\n", + " coastlines = input_directory+\"/shapes_coastlines_Merdith_et_al_v2.gpmlz\"\n", + " COBs = None\n", + "\n", + " from pathlib import Path\n", + " if not Path(input_directory).is_dir():\n", + " raise FileNotFoundError(f\"The input directory does not exist: {input_directory}. Ensure your files are in the correct locations or download the example input files at https://www.earthbyte.org/webdav/ftp/Data_Collections/Merdith_etal_2021_ESR/SM2-Merdith_et_al_1_Ga_reconstruction_v1.2.4.zip\")\n", + " " + ] + }, + { + "cell_type": "markdown", + "id": "1519da7c", + "metadata": {}, + "source": [ + "#### 2) Using `PlateModelManager`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3ebf9727", + "metadata": {}, + "outputs": [], + "source": [ + "if not use_local_files:\n", + " # Method 2: use PlateModelManager\n", + " pm_manager = PlateModelManager()\n", + " plate_model = pm_manager.get_model(\"Merdith2021\", data_dir=\"plate-model-repo\")\n", + " \n", + " rotation_model = plate_model.get_rotation_model()\n", + " topology_features = plate_model.get_topologies()\n", + " static_polygons = plate_model.get_static_polygons()\n", + " \n", + " continents = plate_model.get_layer('ContinentalPolygons') " + ] + }, + { + "cell_type": "markdown", + "id": "3560443e", + "metadata": {}, + "source": [ + "## The SeafloorGrid object\n", + "...is a collection of methods to generate seafloor grids.\n", + "\n", + "The critical input parameters are:\n", + "\n", + "### Plate model parameters\n", + "* **`plate_reconstruction`**: The gplately `PlateReconstruction` object defined in the cell below. This object is a collection of methods for calculating plate tectonic stats through geological time.\n", + "* **`plot_topologies`**: The gplately `PlotTopologies` object defined in the cell below. This object is a collection of methods for resolving geological features needed for gridding to a certain geological time.\n", + "\n", + "#### Define the `PlateReconstruction` and `PlotTopologies` objects" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12bec53e", + "metadata": {}, + "outputs": [], + "source": [ + "model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)\n", + "\n", + "gplot = gplately.PlotTopologies(model, continents=continents)" + ] + }, + { + "cell_type": "markdown", + "id": "57d6e8b2", + "metadata": {}, + "source": [ + "### Time parameters\n", + "* **`max_time`**: The first time step to start generating age grids. This is the time step where we define certain **gridding initial conditions**, which will be explained below.\n", + "* **`min_time`**: The final step to generate age grids. This is the time when recursive reconstructions starting from `max_time` stop.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17add2d7", + "metadata": {}, + "outputs": [], + "source": [ + "max_time = 410.\n", + "min_time = 65." + ] + }, + { + "cell_type": "markdown", + "id": "703d2315", + "metadata": {}, + "source": [ + "### Ridge resolution parameters\n", + "With each reconstruction time step, mid-ocean ridge segments emerge and spread across the ocean floor. In `gplately`, these ridges are lines, or a tessellation of infinitesimal points. The spatio-temporal resolution of these ridge points can be controlled by two parameters:\n", + "* **`ridge_time_step`**: The \"delta time\" or time increment between successive resolutions of ridges, and hence successive grids. By default this is 1Myr, so grids are produced every millionth year.\n", + "* **`ridge_sampling`**: This controls the geographical resolution (in degrees) with which ridge lines are partitioned into points. The larger this is, the larger the spacing between successive ridge points, and hence the smaller the number of ridge points at each timestep. By default this is 0.5 degrees, so points are spaced 0.5 degrees apart. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b959a7b4", + "metadata": {}, + "outputs": [], + "source": [ + "ridge_time_step = 5. # We set the time step to 5 Myr to reduce the notebook’s runtime.\n", + "ridge_sampling = 0.5" + ] + }, + { + "cell_type": "markdown", + "id": "c1c20f7d", + "metadata": {}, + "source": [ + "***\n", + "### Grid resolution parameters\n", + "The ridge resolution parameters mentioned above take care of the spatio-temporal positioning of ridge features/points. However, these ridge points need to be interpolated onto regular grids - separate parameters are needed for these grids. \n", + "\n", + "#### For the `max_time` initial condition\n", + "At `max_time`, i.e. 410Ma, the `PlateReconstruction` model has not been initialised at 411Ma and any time(s) before that to sculpt the geological history before `max_time`. In terms of the workflow, all geological history before `max_time` is unknown. Thus, the global seafloor point distriubtion, and each point's spreading rate and age must be manually defined. \n", + "\n", + "The initial point distribution is an icosahedral point mesh, made using [`stripy`](https://github.com/underworldcode/stripy/tree/294354c00dd72e085a018e69c345d9353c6fafef).\n", + "\n", + "* **`refinement_levels`**: A unitless integer that controls the number of points in the `max_time` icosahedral point mesh.\n", + " \n", + " 6 is the default. Any higher will result in a finer mesh.\n", + " \n", + " \n", + "* **`initial_ocean_mean_spreading_rate`** (in units of mm/yr or km/myr): Since the geological history before and at `max_time` is unknown, we will need to manually define the spreading rates and ages of all ocean points. This is manually set to a uniform spreading rate of 75 (mm/yr or km/myr). Each point's age is equal to its proximity to the nearest mid ocean ridge (assuming that ridge is the source of the point) divided by half this uniform spreading rate (half to account for spreading to the left and right of a ridge). \n", + "\n", + "#### For the interpolated regular grids\n", + "The regular grid on which all data is interpolated has a resolution that can be controlled by:\n", + "\n", + "* **`grid_spacing`**: The degree spacing between successive nodes in the grid. By default, this is 0.1 degrees. **Acceptable degree spacings are 0.1, 0.25, 0.5, 0.75, or 1 degree(s)** because these allow cleanly divisible grid nodes across global latitude and longitude extents. Anything greater than 1 will be rounded to 1, and anything between these acceptable spacings will be rounded to the nearest acceptable spacing. \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5ec00f36", + "metadata": {}, + "outputs": [], + "source": [ + "refinement_levels = 6\n", + "initial_ocean_mean_spreading_rate = 50.\n", + "\n", + "# Gridding parameter\n", + "grid_spacing = 0.25\n", + "\n", + "extent=[-180,180,-90,90]" + ] + }, + { + "cell_type": "markdown", + "id": "8464bfb4", + "metadata": {}, + "source": [ + "***\n", + "### File saving and naming parameters\n", + "When grids are produced, they are saved to:\n", + "* **`save_directory`**: A string to a directory that must exist already.\n", + "\n", + "These files are named according to:\n", + "* **`file_collection`**: A string used to help with naming all output files. \n", + "***" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "106b61f5", + "metadata": {}, + "outputs": [], + "source": [ + "# continent masks, initial ocean seed points, and gridding input files are kept here\n", + "output_parent_directory = os.path.join(\n", + " \"NotebookFiles\",\n", + " \"Notebook10\",\n", + ")\n", + "save_directory = os.path.join(\n", + " output_parent_directory,\n", + " \"seafloor_grid_output\",\n", + ")\n", + "print(f\"The ouput files created by this notebook will be saved in folder {save_directory}.\" )\n", + "\n", + "# A string to help name files according to a plate model \"Merdith2021\"\n", + "file_collection = \"Merdith2021\"\n" + ] + }, + { + "cell_type": "markdown", + "id": "a328fa6d", + "metadata": {}, + "source": [ + "***\n", + "### Continent collision\n", + "\n", + "One way an oceanic point can be deleted from the ocean basin at a certain timestep is if it is collides into continental crust. At each timestep continental crust is represented as a **continental mask** which is a binary grid that partitions continental regions from oceanic regions. By default, continental masks are generated from the reconstructed continental polygons (where the continent polygons features are specified with ``continent_polygon_features``). But they can also be generated using [continent contouring](https://github.com/EarthByte/continent-contouring) (by specifying ``use_continent_contouring=True``) or you can explicitly provide them as your own continent masks (by specifying them with ``continent_mask_filename``).\n", + "\n", + "### Subduction collision\n", + "Another way an oceanic point can be deleted from the ocean basin at a certain timestep is if it is approaching a convergent plate boundary. This case is now handled automatically without specifying any parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "68450bf6", + "metadata": {}, + "outputs": [], + "source": [ + "continent_polygon_features=continents" + ] + }, + { + "cell_type": "markdown", + "id": "bd978b09", + "metadata": {}, + "source": [ + "***\n", + "### Methodology parameters\n", + "* **`resume_from_checkpoints`**: Gridding can take a while (depending on the number of CPUs used). If this is changed to `True`, after gridding was interrupted, then rerunning the gridding cell will resume from where it left off.\n", + "* **`nprocs`**: The number of CPUs to use for parts of the code that are parallelized. The default is `-2` which uses all CPUs except one (to keep the system responsive).\n", + "***" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b1a6fbee", + "metadata": {}, + "outputs": [], + "source": [ + "# Methodology parameters\n", + "resume_from_checkpoints = False\n", + "nprocs = -2" + ] + }, + { + "cell_type": "markdown", + "id": "5038c7e6", + "metadata": {}, + "source": [ + "Use all parameters to define `SeafloorGrid`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ed5fc9e9", + "metadata": {}, + "outputs": [], + "source": [ + "# The SeafloorGrid object with all aforementioned parameters\n", + "seafloorgrid = gplately.SeafloorGrid(\n", + " \n", + " plate_reconstruction = model, \n", + " \n", + " # Time parameters\n", + " max_time = max_time,\n", + " min_time = min_time,\n", + " \n", + " # Ridge tessellation parameters\n", + " ridge_time_step = ridge_time_step,\n", + " ridge_sampling = ridge_sampling,\n", + " \n", + " # Gridding parameters\n", + " grid_spacing = grid_spacing,\n", + "\n", + " extent = extent,\n", + " \n", + " # Naming parameters\n", + " save_directory = save_directory,\n", + " file_collection = file_collection,\n", + " \n", + " # Initial condition parameters\n", + " refinement_levels = refinement_levels,\n", + " initial_ocean_mean_spreading_rate = initial_ocean_mean_spreading_rate,\n", + " \n", + " # Methodology parameters\n", + " resume_from_checkpoints = resume_from_checkpoints,\n", + "\n", + " # Continental polygons.\n", + " continent_polygon_features = continent_polygon_features,\n", + "\n", + " nprocs = nprocs\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "ab3cbf87", + "metadata": {}, + "source": [ + "## How to use `SeafloorGrid`\n", + "\n", + "### 1) Run `seafloorgrid.reconstruct_by_topologies()`\n", + "Running `seafloorgrid.reconstruct_by_topologies()` prepares all ocean basin points and their seafloor ages and spreading rates for gridding, and reconstructs all active points per timestep to form the grids per timestep.\n", + "\n", + "#### Continental masks\n", + "A netCDF4 continent mask is produced for each _time_ from `max_time` to `min_time`.\n", + "\n", + "#### Initial ocean basin\n", + "At `max_time`, an initial ocean seed point icosahedral mesh fills the ocean basin. Each initial seed point is allocated an age based on its distance from the nearest mid-ocean ridge. And these are reconstructed to `min_time` to produce an `.npz` file (containing reconstructed seed point data for seed points created at `max_time`).\n", + "\n", + "#### Mid-ocean ridges\n", + "At each successive timestep (`ridge_time_step`), new seed points emerge from spreading ridges, and they are allocated their own spreading rates. The seed points emerging at a particular _time_ are reconstructed to `min_time` to produce an `.npz` file (containing reconstructed seed point data for seed points created at that _time_). This results in one `.npz` file for each _time_ from `max_time - ridge_time_step` to `min_time`.\n", + "\n", + "#### Seafloor age and spreading rate\n", + "After all seed points (initial and mid-ocean ridge) have been reconstructed to `min_time`, then one `.npz` file is generated for each _time_ from `max_time` to `min_time`. These files contain the gridding input that will later be used to generate the seafloor age and spreading rate grids.\n", + "\n", + "#### Resuming an interrupted gridding run\n", + "*If the `resume_from_checkpoints` parameter is passed as `True`, after `seafloorgrid.reconstruct_by_topologies()` was interrupted, then you can rerun the cell. The workflow will pick up from where it left off (provided the continent mask files and reconstructed seed point data files have not been erased).*\n", + "\n", + "To overwrite all files in `save_directory`, pass `resume_from_checkpoints` as `False` (this is the default behaviour).\n", + "\n", + "#### Reconstruct by topologies\n", + "The `ReconstructByTopologies` object (written by Simon Williams, Nicky Wright, and John Cannon) handles the reconstruction of seed points. `ReconstructByTopologies` (RBT) identifies active points on the ocean basin per timestep. It works as follows:\n", + "\n", + "If an ocean point on one plate ID transitions into another rigid plate ID at the next timestep, RBT calculates the point's velocity difference between both plates. The point **may** have subducted/collided with a continent at this boundary if this velocity difference is higher than a set velocity threshold. To ascertain whether the point should indeed be deactivated, a second test is conducted: RBT checks the previous time position of the point and calculates this point’s proximity to the boundary of the plate ID polygon it is approaching. If this distance is higher than a set distance threshold, then the point is far enough away from the boundary that it cannot be subducted or consumed by it and hence the point is still active. Else, it is deactivated/deleted.\n", + "\n", + "Once all active points and their z-values are identified, they are written to the gridding input file (.npz) for that timestep.\n", + "\n", + "Below is an example of the initial condition: the ocean basin is populated with an `intial_mean_ocean_spreading_rate` of 50 mm/yr at a `max_time` of 410Ma. Reconstruction over 10 Myr to 400 Ma sees new seed points emerging from ridges with their own plate-model-ascribed spreading rates.\n", + "\n", + "![init_condition](https://raw.githubusercontent.com/GPlates/gplately/master/Notebooks/NotebookFiles/Notebook10/Merdith2021_initial_sr_conditions.png)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "143ea077", + "metadata": {}, + "outputs": [], + "source": [ + "import datetime\n", + "import time\n", + "\n", + "start = time.time()\n", + "seafloorgrid.reconstruct_by_topologies()\n", + "end = time.time()\n", + "duration = datetime.timedelta(seconds=end - start)\n", + "print(\"Duration: {}\".format(duration))" + ] + }, + { + "cell_type": "markdown", + "id": "a0a0ca60", + "metadata": {}, + "source": [ + "### 2) Run `seafloorgrid.lat_lon_z_to_netCDF` to write grids to netCDF\n", + "Calling `seafloorgrid.lat_lon_z_to_netCDF` grids one set of z-data per latitude-longitude pair from each timestep's gridding input file (produced in `seafloorgrid.reconstruct_by_topologies()`). Grids are in netCDF format.\n", + "\n", + "The desired z-data to grid is identified using a `zval_name`. \n", + "For example, seafloor age grids can be produced using `SEAFLOOR_AGE`, and spreading rate grids are `SPREADING_RATE`. \n", + "\n", + "Use `unmasked = True` to output both the masked and unmasked versions of the grids." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "891183f3", + "metadata": {}, + "outputs": [], + "source": [ + "seafloorgrid.lat_lon_z_to_netCDF(\"SEAFLOOR_AGE\", unmasked=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d51a41a5", + "metadata": {}, + "outputs": [], + "source": [ + "seafloorgrid.lat_lon_z_to_netCDF(\"SPREADING_RATE\", unmasked=True)" + ] + }, + { + "cell_type": "markdown", + "id": "8213b46a", + "metadata": {}, + "source": [ + "### Plotting a sample age grid and spreading rate grid\n", + "Read one netCDF grid using GPlately's `Raster` object from `grids`, and plot it using the `PlotTopologies` object.\n", + "\n", + "Notice the evolution of seafloor spreading rate from the initial value set with `initial_ocean_mean_spreading_rate`. Eventually, this initial uniform spreading rate will be phased out with sufficient recursive reconstruction." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "06a92ce9-ecea-49b9-9c99-c33438e75379", + "metadata": {}, + "outputs": [], + "source": [ + "time = min_time\n", + "# set the plot_unmasked_grid variable to True if you would like to see the unmasked age grid.\n", + "plot_unmasked_grid = False\n", + "\n", + "# age grids\n", + "if not plot_unmasked_grid:\n", + " agegrid_filename = f\"{seafloorgrid.save_directory}/SEAFLOOR_AGE/{seafloorgrid.file_collection}_SEAFLOOR_AGE_grid_{time:0.2f}Ma.nc\"\n", + "else:\n", + " agegrid_filename = f\"{seafloorgrid.save_directory}/SEAFLOOR_AGE/{seafloorgrid.file_collection}_SEAFLOOR_AGE_grid_unmasked_{time:0.2f}Ma.nc\"\n", + "\n", + "age_grid = gplately.grids.Raster(agegrid_filename)\n", + "\n", + "# Prepare plots\n", + "fig = plt.figure(figsize=(8,4), dpi=240, linewidth=1)\n", + "ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))\n", + "\n", + "gplot.time = time\n", + "plt.title(f\"Seafloor Age at {int(time)} Ma (Model:{seafloorgrid.file_collection})\")\n", + "\n", + "agegrid_cmap = \"YlGnBu\"\n", + "agegrid_cpt_file = \"NotebookFiles/agegrid.cpt\"\n", + "if not os.path.isfile(agegrid_cpt_file):\n", + " print(f\"The file {agegrid_cpt_file} does not exist locally. It can be found at https://raw.githubusercontent.com/GPlates/gplately/refs/heads/master/Notebooks/NotebookFiles/agegrid.cpt\")\n", + " print(f\"The '{agegrid_cmap}' will be used to plot the age grid.\")\n", + " agegrid_cpt_file = None\n", + "if agegrid_cpt_file:\n", + " try:\n", + " from gplately.mapping.gmt_cpt import get_cmap_from_gmt_cpt\n", + " agegrid_cmap=get_cmap_from_gmt_cpt(agegrid_cpt_file)\n", + " except ImportError as e:\n", + " print(\"To use GMT cpt file, ensure you are using GPlately version 2.1 or higher.\")\n", + " print(e)\n", + "\n", + "im = gplot.plot_grid(\n", + " ax, \n", + " age_grid.data, \n", + " cmap = agegrid_cmap,\n", + " vmin = 0, \n", + " vmax =410,\n", + ")\n", + "gplot.plot_ridges(ax,linewidth=1)\n", + "if not plot_unmasked_grid:\n", + " gplot.plot_continents(ax, edgecolor=\"grey\", facecolor=\"lightgrey\", linewidth=0.5)\n", + "plt.colorbar(im, label=\"Seafloor Age (Myr)\", shrink=0.5, pad=0.05)\n", + "\n", + "# Save figure\n", + "plt.savefig(\n", + " f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_age_{int(time)}Ma.png\",\n", + " bbox_inches = 'tight',\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "56e03e39", + "metadata": {}, + "outputs": [], + "source": [ + "time = min_time\n", + "# set the plot_unmasked_grid variable to True if you would like to see the unmasked spreading rate grid.\n", + "plot_unmasked_grid = False\n", + "\n", + "# spreading rate grids\n", + "if not plot_unmasked_grid:\n", + " srgrid_filename = f\"{seafloorgrid.save_directory}/SPREADING_RATE/{seafloorgrid.file_collection}_SPREADING_RATE_grid_{time:0.2f}Ma.nc\"\n", + "else:\n", + " srgrid_filename = f\"{seafloorgrid.save_directory}/SPREADING_RATE/{seafloorgrid.file_collection}_SPREADING_RATE_grid_unmasked_{time:0.2f}Ma.nc\"\n", + "\n", + "sr_grid = gplately.grids.Raster(srgrid_filename)\n", + "\n", + "# Prepare plots\n", + "fig = plt.figure(figsize=(8,4), dpi=240, linewidth=1)\n", + "ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))\n", + "\n", + "gplot.time = time\n", + "plt.title(f\"Seafloor Spreading Rate at {int(time)} Ma (Model:{seafloorgrid.file_collection})\")\n", + "\n", + "spreading_rate_cmap = \"magma\"\n", + "spreading_rate_cpt_file = \"NotebookFiles/spreading_full_rate.cpt\"\n", + "if not os.path.isfile(spreading_rate_cpt_file):\n", + " print(f\"The file {spreading_rate_cpt_file} does not exist locally. It can be found at https://raw.githubusercontent.com/GPlates/gplately/refs/heads/master/Notebooks/NotebookFiles/spreading_full_rate.cpt\")\n", + " print(f\"The '{spreading_rate_cmap}' will be used to plot the spreading rate grid.\")\n", + " spreading_rate_cpt_file = None\n", + "if spreading_rate_cpt_file:\n", + " try:\n", + " from gplately.mapping.gmt_cpt import get_cmap_from_gmt_cpt\n", + " spreading_rate_cmap=get_cmap_from_gmt_cpt(spreading_rate_cpt_file)\n", + " except ImportError as e:\n", + " print(\"To use GMT cpt file, ensure you are using GPlately version 2.1 or higher.\")\n", + " print(e)\n", + "\n", + "im = gplot.plot_grid(\n", + " ax, \n", + " sr_grid.data, \n", + " cmap = spreading_rate_cmap,\n", + " vmin = 0, \n", + " vmax =200,\n", + ")\n", + "gplot.plot_ridges(ax, linewidth=1)\n", + "if not plot_unmasked_grid:\n", + " gplot.plot_continents(ax, edgecolor=\"grey\", facecolor=\"lightgrey\", linewidth=0.5)\n", + "plt.colorbar(im, label=\"Spreading rate (mm/yr)\", shrink=0.5, pad=0.05, ticks=[0,50,100,150,200])\n", + "\n", + "# Save figure\n", + "plt.savefig(\n", + " f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_spreading_rate_{int(time)}Ma.png\",\n", + " bbox_inches = 'tight',\n", + ")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 6dcdd5ac3b218940af85f3501ec763ec44de304d Mon Sep 17 00:00:00 2001 From: John Cannon Date: Wed, 3 Jun 2026 13:06:54 +1000 Subject: [PATCH 14/17] Set initial seafloor spreading rate to 75 mm/yr to match default. --- Notebooks/10-SeafloorGrids.ipynb | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/Notebooks/10-SeafloorGrids.ipynb b/Notebooks/10-SeafloorGrids.ipynb index 25828d85..b87866ed 100644 --- a/Notebooks/10-SeafloorGrids.ipynb +++ b/Notebooks/10-SeafloorGrids.ipynb @@ -43,7 +43,6 @@ "\n", "import gplately\n", "\n", - "import pygplates\n", "import glob\n", "import matplotlib.pyplot as plt\n", "import cartopy.crs as ccrs\n", @@ -246,7 +245,7 @@ "outputs": [], "source": [ "refinement_levels = 6\n", - "initial_ocean_mean_spreading_rate = 50.\n", + "initial_ocean_mean_spreading_rate = 75.\n", "\n", "# Gridding parameter\n", "grid_spacing = 0.25\n", @@ -357,7 +356,7 @@ "# The SeafloorGrid object with all aforementioned parameters\n", "seafloorgrid = gplately.SeafloorGrid(\n", " \n", - " plate_reconstruction = model, \n", + " model,\n", " \n", " # Time parameters\n", " max_time = max_time,\n", From 3a7e4fac10b666fab86d0f6daec923cc64b140fd Mon Sep 17 00:00:00 2001 From: John Cannon Date: Wed, 3 Jun 2026 23:29:23 +1000 Subject: [PATCH 15/17] Optionally generate seafloor age and spreading rate videos. Also merge separate seafloor age and spreading rate plots into a single function. --- Notebooks/10-SeafloorGrids.ipynb | 296 +++++++++++++++++++++---------- 1 file changed, 199 insertions(+), 97 deletions(-) diff --git a/Notebooks/10-SeafloorGrids.ipynb b/Notebooks/10-SeafloorGrids.ipynb index b87866ed..ef6605bb 100644 --- a/Notebooks/10-SeafloorGrids.ipynb +++ b/Notebooks/10-SeafloorGrids.ipynb @@ -497,115 +497,217 @@ "metadata": {}, "outputs": [], "source": [ - "time = min_time\n", - "# set the plot_unmasked_grid variable to True if you would like to see the unmasked age grid.\n", - "plot_unmasked_grid = False\n", - "\n", - "# age grids\n", - "if not plot_unmasked_grid:\n", - " agegrid_filename = f\"{seafloorgrid.save_directory}/SEAFLOOR_AGE/{seafloorgrid.file_collection}_SEAFLOOR_AGE_grid_{time:0.2f}Ma.nc\"\n", - "else:\n", - " agegrid_filename = f\"{seafloorgrid.save_directory}/SEAFLOOR_AGE/{seafloorgrid.file_collection}_SEAFLOOR_AGE_grid_unmasked_{time:0.2f}Ma.nc\"\n", - "\n", - "age_grid = gplately.grids.Raster(agegrid_filename)\n", - "\n", - "# Prepare plots\n", - "fig = plt.figure(figsize=(8,4), dpi=240, linewidth=1)\n", - "ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))\n", - "\n", - "gplot.time = time\n", - "plt.title(f\"Seafloor Age at {int(time)} Ma (Model:{seafloorgrid.file_collection})\")\n", - "\n", - "agegrid_cmap = \"YlGnBu\"\n", - "agegrid_cpt_file = \"NotebookFiles/agegrid.cpt\"\n", - "if not os.path.isfile(agegrid_cpt_file):\n", - " print(f\"The file {agegrid_cpt_file} does not exist locally. It can be found at https://raw.githubusercontent.com/GPlates/gplately/refs/heads/master/Notebooks/NotebookFiles/agegrid.cpt\")\n", - " print(f\"The '{agegrid_cmap}' will be used to plot the age grid.\")\n", - " agegrid_cpt_file = None\n", - "if agegrid_cpt_file:\n", - " try:\n", - " from gplately.mapping.gmt_cpt import get_cmap_from_gmt_cpt\n", - " agegrid_cmap=get_cmap_from_gmt_cpt(agegrid_cpt_file)\n", - " except ImportError as e:\n", - " print(\"To use GMT cpt file, ensure you are using GPlately version 2.1 or higher.\")\n", - " print(e)\n", - "\n", - "im = gplot.plot_grid(\n", - " ax, \n", - " age_grid.data, \n", - " cmap = agegrid_cmap,\n", - " vmin = 0, \n", - " vmax =410,\n", - ")\n", - "gplot.plot_ridges(ax,linewidth=1)\n", - "if not plot_unmasked_grid:\n", - " gplot.plot_continents(ax, edgecolor=\"grey\", facecolor=\"lightgrey\", linewidth=0.5)\n", - "plt.colorbar(im, label=\"Seafloor Age (Myr)\", shrink=0.5, pad=0.05)\n", - "\n", - "# Save figure\n", - "plt.savefig(\n", - " f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_age_{int(time)}Ma.png\",\n", - " bbox_inches = 'tight',\n", - ")" + "def plot(time, zval_name, save_fig=False):\n", + " \n", + " if zval_name != \"SEAFLOOR_AGE\" and zval_name != \"SPREADING_RATE\":\n", + " raise ValueError(f\"Unknown zval_name '{zval_name}'\")\n", + " \n", + " # set the plot_unmasked_grid variable to True if you would like to see the unmasked age/spreading-rate grid.\n", + " plot_unmasked_grid = False\n", + " \n", + " # age/spreading-rate grids\n", + " if not plot_unmasked_grid:\n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " grid_filename = f\"{seafloorgrid.save_directory}/SEAFLOOR_AGE/{seafloorgrid.file_collection}_SEAFLOOR_AGE_grid_{time:0.2f}Ma.nc\"\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " grid_filename = f\"{seafloorgrid.save_directory}/SPREADING_RATE/{seafloorgrid.file_collection}_SPREADING_RATE_grid_{time:0.2f}Ma.nc\"\n", + " else:\n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " grid_filename = f\"{seafloorgrid.save_directory}/SEAFLOOR_AGE/{seafloorgrid.file_collection}_SEAFLOOR_AGE_grid_unmasked_{time:0.2f}Ma.nc\"\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " grid_filename = f\"{seafloorgrid.save_directory}/SPREADING_RATE/{seafloorgrid.file_collection}_SPREADING_RATE_grid_unmasked_{time:0.2f}Ma.nc\"\n", + " \n", + " grid = gplately.grids.Raster(grid_filename)\n", + " \n", + " # Prepare plots\n", + " fig = plt.figure(figsize=(8,4), dpi=240, linewidth=1)\n", + " ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))\n", + " \n", + " gplot.time = time\n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " plt.title(f\"Seafloor Age at {int(time)} Ma (Model:{seafloorgrid.file_collection})\")\n", + " cmap = \"YlGnBu\"\n", + " cpt_file = \"NotebookFiles/agegrid.cpt\"\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " plt.title(f\"Seafloor Spreading Rate at {int(time)} Ma (Model:{seafloorgrid.file_collection})\")\n", + " cmap = \"magma\"\n", + " cpt_file = \"NotebookFiles/spreading_full_rate.cpt\"\n", + " \n", + " if not os.path.isfile(cpt_file):\n", + " print(f\"The file {cpt_file} does not exist locally. It can be found at https://raw.githubusercontent.com/GPlates/gplately/refs/heads/master/Notebooks/{cpt_file}\")\n", + " print(f\"The '{cmap}' will be used to plot the grid.\")\n", + " cpt_file = None\n", + " if cpt_file:\n", + " try:\n", + " from gplately.mapping.gmt_cpt import get_cmap_from_gmt_cpt\n", + " cmap=get_cmap_from_gmt_cpt(cpt_file)\n", + " except ImportError as e:\n", + " print(\"To use GMT cpt file, ensure you are using GPlately version 2.1 or higher.\")\n", + " print(e)\n", + "\n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " vmin, vmax = 0, 410\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " vmin, vmax = 0, 200\n", + " \n", + " im = gplot.plot_grid(\n", + " ax, \n", + " grid.data, \n", + " cmap = cmap,\n", + " vmin = vmin, \n", + " vmax = vmax,\n", + " )\n", + " gplot.plot_ridges(ax,linewidth=1)\n", + " if not plot_unmasked_grid:\n", + " gplot.plot_continents(ax, edgecolor=\"grey\", facecolor=\"lightgrey\", linewidth=0.5)\n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " plt.colorbar(im, label=\"Seafloor Age (Myr)\", shrink=0.5, pad=0.05, ticks=[0,100,200,300,400])\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " plt.colorbar(im, label=\"Spreading rate (mm/yr)\", shrink=0.5, pad=0.05, ticks=[0,50,100,150,200])\n", + " \n", + " if save_fig:\n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " fig_filename = f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_age_{int(time)}Ma.png\"\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " fig_filename = f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_spreading_rate_{int(time)}Ma.png\"\n", + " \n", + " # Save figure.\n", + " plt.savefig(fig_filename, bbox_inches = 'tight')\n", + " else:\n", + " plt.show()" ] }, { "cell_type": "code", "execution_count": null, - "id": "56e03e39", + "id": "87cd68e6-fb41-42e8-8d3a-4ee317ed002d", "metadata": {}, "outputs": [], "source": [ - "time = min_time\n", - "# set the plot_unmasked_grid variable to True if you would like to see the unmasked spreading rate grid.\n", - "plot_unmasked_grid = False\n", - "\n", - "# spreading rate grids\n", - "if not plot_unmasked_grid:\n", - " srgrid_filename = f\"{seafloorgrid.save_directory}/SPREADING_RATE/{seafloorgrid.file_collection}_SPREADING_RATE_grid_{time:0.2f}Ma.nc\"\n", - "else:\n", - " srgrid_filename = f\"{seafloorgrid.save_directory}/SPREADING_RATE/{seafloorgrid.file_collection}_SPREADING_RATE_grid_unmasked_{time:0.2f}Ma.nc\"\n", + "plot(min_time, zval_name=\"SEAFLOOR_AGE\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "08979858-bd12-4cf5-accd-b372405ecc21", + "metadata": {}, + "outputs": [], + "source": [ + "plot(min_time, zval_name=\"SPREADING_RATE\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "38d323c9-772e-4151-9f0d-66065b6f4a73", + "metadata": {}, + "outputs": [], + "source": [ + "def save_figures(\n", + " zval_name,\n", + " min_time=min_time,\n", + " max_time=max_time,\n", + " time_step=ridge_time_step):\n", + " \n", + " from joblib import Parallel, delayed\n", + " import numpy as np\n", "\n", - "sr_grid = gplately.grids.Raster(srgrid_filename)\n", + " time_range = np.arange(max_time, min_time-0.5*time_step, -time_step)\n", + " \n", + " use_parallel = True\n", + " if use_parallel:\n", + " # Produce plots in a parallel routine - a progress bar and time taken will be shown.\n", + " parallel_plots = Parallel(n_jobs=nprocs, verbose=1)(\n", + " delayed(plot) (time, zval_name, save_fig=True) for time in time_range)\n", + " \n", + " else:\n", + " for time in time_range:\n", + " plot(time, save_fig=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eb5be45a-2725-47cd-be7c-f512a25574ca", + "metadata": {}, + "outputs": [], + "source": [ + "def save_video(\n", + " zval_name,\n", + " min_time=min_time,\n", + " max_time=max_time,\n", + " time_step=ridge_time_step):\n", + " \n", + " try:\n", + " import moviepy.editor as mpy # moviepy 1.x\n", + " except ImportError:\n", + " import moviepy as mpy # moviepy 2.x\n", + " import numpy as np\n", "\n", - "# Prepare plots\n", - "fig = plt.figure(figsize=(8,4), dpi=240, linewidth=1)\n", - "ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))\n", + " if zval_name != \"SEAFLOOR_AGE\" and zval_name != \"SPREADING_RATE\":\n", + " raise ValueError(f\"Unknown zval_name '{zval_name}'\")\n", + " \n", + " time_range = np.arange(max_time, min_time-0.5*time_step, -time_step)\n", + " \n", + " frame_list = []\n", + " for time in time_range:\n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " fig_filename = f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_age_{int(time)}Ma.png\"\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " fig_filename = f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_spreading_rate_{int(time)}Ma.png\"\n", + " \n", + " frame_list.append(fig_filename)\n", + " \n", + " clip = mpy.ImageSequenceClip(frame_list, fps=20)\n", + " \n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " video_filename = f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_age.mp4\"\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " video_filename = f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_spreading_rate.mp4\"\n", + " \n", + " clip.write_videofile(\n", + " video_filename,\n", + " fps=20,\n", + " codec=\"libx264\",\n", + " bitrate=\"8000k\",\n", + " audio=False,\n", + " logger=None,\n", + " ffmpeg_params=[\n", + " \"-vf\",\n", + " \"pad=ceil(iw/2)*2:ceil(ih/2)*2\",\n", + " \"-pix_fmt\",\n", + " \"yuv420p\",\n", + " ],\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e585b0cf-cd2f-4ead-b19b-4f926ea3139c", + "metadata": {}, + "outputs": [], + "source": [ + "generate_seafloor_age_video = False\n", "\n", - "gplot.time = time\n", - "plt.title(f\"Seafloor Spreading Rate at {int(time)} Ma (Model:{seafloorgrid.file_collection})\")\n", + "if generate_seafloor_age_video:\n", + " save_figures(zval_name=\"SEAFLOOR_AGE\")\n", + " save_video(zval_name=\"SEAFLOOR_AGE\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3202b8cc-d134-46ea-b1db-4fc747adaec8", + "metadata": {}, + "outputs": [], + "source": [ + "generate_spreading_rate_video = False\n", "\n", - "spreading_rate_cmap = \"magma\"\n", - "spreading_rate_cpt_file = \"NotebookFiles/spreading_full_rate.cpt\"\n", - "if not os.path.isfile(spreading_rate_cpt_file):\n", - " print(f\"The file {spreading_rate_cpt_file} does not exist locally. It can be found at https://raw.githubusercontent.com/GPlates/gplately/refs/heads/master/Notebooks/NotebookFiles/spreading_full_rate.cpt\")\n", - " print(f\"The '{spreading_rate_cmap}' will be used to plot the spreading rate grid.\")\n", - " spreading_rate_cpt_file = None\n", - "if spreading_rate_cpt_file:\n", - " try:\n", - " from gplately.mapping.gmt_cpt import get_cmap_from_gmt_cpt\n", - " spreading_rate_cmap=get_cmap_from_gmt_cpt(spreading_rate_cpt_file)\n", - " except ImportError as e:\n", - " print(\"To use GMT cpt file, ensure you are using GPlately version 2.1 or higher.\")\n", - " print(e)\n", - "\n", - "im = gplot.plot_grid(\n", - " ax, \n", - " sr_grid.data, \n", - " cmap = spreading_rate_cmap,\n", - " vmin = 0, \n", - " vmax =200,\n", - ")\n", - "gplot.plot_ridges(ax, linewidth=1)\n", - "if not plot_unmasked_grid:\n", - " gplot.plot_continents(ax, edgecolor=\"grey\", facecolor=\"lightgrey\", linewidth=0.5)\n", - "plt.colorbar(im, label=\"Spreading rate (mm/yr)\", shrink=0.5, pad=0.05, ticks=[0,50,100,150,200])\n", - "\n", - "# Save figure\n", - "plt.savefig(\n", - " f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_spreading_rate_{int(time)}Ma.png\",\n", - " bbox_inches = 'tight',\n", - ")" + "if generate_spreading_rate_video:\n", + " save_figures(zval_name=\"SPREADING_RATE\")\n", + " save_video(zval_name=\"SPREADING_RATE\")" ] } ], From cd678ef21cf2b16571550c3be42e1ac9507f0b08 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Thu, 25 Jun 2026 22:43:30 +1000 Subject: [PATCH 16/17] Some minor suggestions from claude's code review. --- gplately/lib/reconstruct_by_topologies.py | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/gplately/lib/reconstruct_by_topologies.py b/gplately/lib/reconstruct_by_topologies.py index b3f88a80..b032785d 100644 --- a/gplately/lib/reconstruct_by_topologies.py +++ b/gplately/lib/reconstruct_by_topologies.py @@ -227,10 +227,10 @@ def _reconstruct_to_next_time_impl(self): # Get the resolved networks and rigid plates. curr_resolved_networks = curr_topological_snapshot.get_resolved_topologies( - pygplates.ResolveTopologyType.network # type:ignore + pygplates.ResolveTopologyType.network # type: ignore ) curr_resolved_boundaries = curr_topological_snapshot.get_resolved_topologies( - pygplates.ResolveTopologyType.boundary # type:ignore + pygplates.ResolveTopologyType.boundary # type: ignore ) # Get the currently active points and their indices (into the original points). @@ -408,6 +408,8 @@ def _reconstruct_to_next_time_impl(self): # Reconstruct the currently active point from its position at current time to its position at the next time step. curr_active_point = self._all_current_points[point_index] + # This should not return None because the currently active point has previously been determined to be inside + # the current resolved topology (in the previous time step). next_active_point = curr_resolved_topology.reconstruct_point( curr_active_point, next_time ) @@ -478,6 +480,9 @@ def _deactivate_continent_masked_points(self): # Get the currently active points and their indices (into the original points). curr_active_point_indices = self.get_current_active_point_indices() + if len(curr_active_point_indices) == 0: + # There are no active points so nothing to do. + return curr_active_points = [ self._all_current_points[point_index] for point_index in curr_active_point_indices @@ -797,7 +802,7 @@ def __init__( self.next_boundary_sub_segments ): raise RuntimeError( - "Current and next topologies have a different number of boundary sub-segments." + "Expected current and next topologies have the same number of boundary sub-segments." ) self.num_boundary_sub_segments = len(self.curr_boundary_sub_segments) @@ -884,7 +889,7 @@ def _get_next_boundary_from_consistent_boundary_sub_segments( # # Note: We don't need to clone this feature before modifying it because it was already essentially cloned # when 'ResolvedTopologicalSubSegment.get_resolved_feature()' was called. - pygplates.reverse_reconstruct( # type:ignore + pygplates.reverse_reconstruct( # type: ignore curr_boundary_sub_segment_feature, self.rotation_model, self.curr_time, @@ -901,7 +906,7 @@ def _get_next_boundary_from_consistent_boundary_sub_segments( # Reconstruct to 'next_time' (from present day). next_boundary_sub_segment_reconstructed_geometries = ( - pygplates.ReconstructSnapshot( # type:ignore + pygplates.ReconstructSnapshot( # type: ignore curr_boundary_sub_segment_feature, self.rotation_model, self.next_time, @@ -969,7 +974,7 @@ def _get_next_boundary_from_consistent_boundary_sub_segments( # # Note: We don't need to clone these features before modifying them because they were already essentially cloned # when 'ResolvedTopologicalSubSegment.get_resolved_feature()' was called. - pygplates.reverse_reconstruct( # type:ignore + pygplates.reverse_reconstruct( # type: ignore curr_boundary_sub_sub_segment_features, self.rotation_model, self.curr_time, @@ -984,7 +989,7 @@ def _get_next_boundary_from_consistent_boundary_sub_segments( ) # Reconstruct to 'next_time' (from present day). - next_boundary_sub_sub_segment_reconstructed_geometries = pygplates.ReconstructSnapshot( # type:ignore + next_boundary_sub_sub_segment_reconstructed_geometries = pygplates.ReconstructSnapshot( # type: ignore curr_boundary_sub_sub_segment_features, self.rotation_model, self.next_time, @@ -1016,7 +1021,7 @@ def _get_next_boundary_from_consistent_boundary_sub_segments( # Join the (potentially reversed) points of the boundary sub-segments and join them together into a polygon boundary. # Each sub-segment is either from the *current* or *next* resolved topology. - return pygplates.PolygonOnSphere( # type:ignore + return pygplates.PolygonOnSphere( # type: ignore next_boundary_polygon_points ) From 4cd01092192c9b101fc2c6a6504752d085b76a8c Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Thu, 25 Jun 2026 12:56:19 +0000 Subject: [PATCH 17/17] ci: regenerate notebooks --- .../14-RuleBasedGPMLProcessingPipeline.ipynb | 50 +++++++++---------- Notebooks/Examples/plot_map_with_pygmt.ipynb | 6 +-- 2 files changed, 28 insertions(+), 28 deletions(-) diff --git a/Notebooks/14-RuleBasedGPMLProcessingPipeline.ipynb b/Notebooks/14-RuleBasedGPMLProcessingPipeline.ipynb index ff30d24a..071da249 100644 --- a/Notebooks/14-RuleBasedGPMLProcessingPipeline.ipynb +++ b/Notebooks/14-RuleBasedGPMLProcessingPipeline.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "3d561d78", + "id": "b2b60c22", "metadata": {}, "source": [ "\n", @@ -19,7 +19,7 @@ }, { "cell_type": "markdown", - "id": "4ca19dbe", + "id": "c24098e8", "metadata": {}, "source": [ "\n", @@ -30,7 +30,7 @@ }, { "cell_type": "markdown", - "id": "743d5977", + "id": "2938cc71", "metadata": { "lines_to_next_cell": 0 }, @@ -41,7 +41,7 @@ { "cell_type": "code", "execution_count": null, - "id": "680ba013", + "id": "40352ed2", "metadata": {}, "outputs": [], "source": [ @@ -118,7 +118,7 @@ }, { "cell_type": "markdown", - "id": "b6bff13f", + "id": "68999076", "metadata": {}, "source": [ "### Search and filter feature collection with pre-defined filters\n", @@ -155,7 +155,7 @@ { "cell_type": "code", "execution_count": null, - "id": "7b9d374f", + "id": "4cb62faa", "metadata": {}, "outputs": [], "source": [ @@ -363,7 +363,7 @@ }, { "cell_type": "markdown", - "id": "45bd7b67", + "id": "41d6296c", "metadata": { "lines_to_next_cell": 0 }, @@ -382,7 +382,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d87cad2f", + "id": "f7e8ca10", "metadata": {}, "outputs": [], "source": [ @@ -418,7 +418,7 @@ }, { "cell_type": "markdown", - "id": "d8fabe5d", + "id": "bbe1f6e0", "metadata": { "lines_to_next_cell": 2 }, @@ -435,7 +435,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16182778", + "id": "b7097090", "metadata": { "lines_to_next_cell": 0 }, @@ -468,7 +468,7 @@ }, { "cell_type": "markdown", - "id": "cbb8ee58", + "id": "c50cec99", "metadata": {}, "source": [ "#### Find Laurussia topological plate and extract the reference features for investigation\n", @@ -484,7 +484,7 @@ { "cell_type": "code", "execution_count": null, - "id": "54860ec7", + "id": "614a6375", "metadata": {}, "outputs": [], "source": [ @@ -515,7 +515,7 @@ }, { "cell_type": "markdown", - "id": "b4fb459c", + "id": "dd52b0c4", "metadata": { "lines_to_next_cell": 0 }, @@ -530,7 +530,7 @@ { "cell_type": "code", "execution_count": null, - "id": "b49c55b4", + "id": "b752f2a9", "metadata": { "lines_to_next_cell": 0 }, @@ -561,7 +561,7 @@ }, { "cell_type": "markdown", - "id": "d5125a0d", + "id": "de192aa5", "metadata": {}, "source": [ "#### Topological reference map\n", @@ -577,7 +577,7 @@ { "cell_type": "code", "execution_count": null, - "id": "fd649426", + "id": "c2b8685f", "metadata": {}, "outputs": [], "source": [ @@ -595,7 +595,7 @@ }, { "cell_type": "markdown", - "id": "ac0440dc", + "id": "10a7e434", "metadata": { "lines_to_next_cell": 2 }, @@ -610,7 +610,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6835cda5", + "id": "b9ca5681", "metadata": {}, "outputs": [], "source": [ @@ -641,7 +641,7 @@ }, { "cell_type": "markdown", - "id": "eca479ff", + "id": "0be467ad", "metadata": {}, "source": [ "if you open the icosahedron_mesh_5.gpmlz file in GPlates, you will see the vertices of the icosahedron mesh,\n", @@ -655,7 +655,7 @@ }, { "cell_type": "markdown", - "id": "a512cac1", + "id": "4e0d6899", "metadata": { "lines_to_next_cell": 2 }, @@ -666,7 +666,7 @@ { "cell_type": "code", "execution_count": null, - "id": "22b6779a", + "id": "320a638d", "metadata": {}, "outputs": [], "source": [ @@ -694,7 +694,7 @@ }, { "cell_type": "markdown", - "id": "4937648e", + "id": "f75e060f", "metadata": {}, "source": [ "If you open icosahedron_vertices_in_region.gpmlz in GPlates, you will see the vertices of the icosahedron mesh that\n", @@ -709,7 +709,7 @@ }, { "cell_type": "markdown", - "id": "469147d9", + "id": "32031886", "metadata": { "lines_to_next_cell": 2 }, @@ -720,7 +720,7 @@ { "cell_type": "code", "execution_count": null, - "id": "194b9e5e", + "id": "e22c0d6f", "metadata": { "lines_to_next_cell": 0 }, @@ -766,7 +766,7 @@ }, { "cell_type": "markdown", - "id": "5f3e7dec", + "id": "f6dcfcc8", "metadata": {}, "source": [ "If you open icosahedron_vertices_within_australia.gpmlz in GPlates, you will see the vertices of the icosahedron mesh that\n", diff --git a/Notebooks/Examples/plot_map_with_pygmt.ipynb b/Notebooks/Examples/plot_map_with_pygmt.ipynb index bf1c073a..b4c17e18 100644 --- a/Notebooks/Examples/plot_map_with_pygmt.ipynb +++ b/Notebooks/Examples/plot_map_with_pygmt.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "a1c4908b", + "id": "32cfe43d", "metadata": {}, "source": [ "This notebook demonstrates how to use the PyGMT plotting maps in GPlately.\n", @@ -14,7 +14,7 @@ }, { "cell_type": "markdown", - "id": "fe001ac0", + "id": "4eb5307d", "metadata": {}, "source": [ "⚠️ This notebook is generated from plot_map_with_pygmt.py using the command `jupytext --to notebook Notebooks/Examples/plot_map_with_pygmt.py -o Notebooks/Examples/plot_map_with_pygmt.ipynb`.\n", @@ -25,7 +25,7 @@ { "cell_type": "code", "execution_count": null, - "id": "26ee75f4", + "id": "52059e55", "metadata": {}, "outputs": [], "source": [