From d9dd82ef746c70be94f15cb81cc3819f7f467012 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Fri, 29 May 2026 11:36:58 +1000 Subject: [PATCH 1/5] Agegridding now generates MOR seed points at topological lines. Previously mid-ocean ridge features needed spreading properties as outlined in #190. Now we use plate boundary statistics provided by pyGPlates 1.0 to calculate spreading rates from left/right plate velocities as covered in #364. --- gplately/oceans.py | 189 +++++++++++++++++++++++---------------------- 1 file changed, 95 insertions(+), 94 deletions(-) diff --git a/gplately/oceans.py b/gplately/oceans.py index d137623a..095f95b4 100644 --- a/gplately/oceans.py +++ b/gplately/oceans.py @@ -955,114 +955,115 @@ def _build_mid_ocean_ridge_seed_points( self, time: float, ): - """Resolve mid-ocean ridges at ``time`` and divide them into points that make up their shared sub-segments. - Rotate these points to the left and right of the ridge using their stage rotation so that they spread from the ridge. + """Resolve mid-ocean ridges at ``time`` and divide them into uniformly spaced points. + Rotate these points to the left and right of the ridge using the left and right plate velocities so that they spread from the ridge. .. note:: - This assumes that points spread from ridges symmetrically, with the exception of - large ridge jumps at successive timesteps. Therefore, spreading rates of ridge-emerging - points will appear symmetrical until changes in spreading ridge geometries create asymmetries. + This supports spreading asymmetry across the ridge. Each point is shifted by a small angle (0.01 degrees) to the left and right of the ridge. + The left and right spreading rates are calculated as the difference between the left and right plate velocities and the ridge velocity at the point. + """ - .. seealso:: + # Generate statistics at uniformly spaced points along mid-ocean ridges. + mor_boundary_statistics = self.plate_reconstruction.topological_snapshot( + time, + # Ignore topological slab boundaries since they are not *plate* boundaries... + include_topological_slab_boundaries=False, + ).calculate_plate_boundary_statistics( + uniform_point_spacing_radians=np.radians(self.ridge_sampling), + velocity_units=pygplates.VelocityUnits.kms_per_my, # km/Myr is the same as mm/yr # type: ignore + boundary_section_filter=pygplates.FeatureType.gpml_mid_ocean_ridge, # type: ignore + ) - `Get tessellated points along a mid ocean ridge `__. - """ + # How much to rotate each ridge point off the ridge to get a point definitely inside the plate for later topological reconstruction. + boundary_rotation_angle_radians = np.radians(0.01) + + def _calc_shifted_mor_point_and_spreading_rate(stat, plate, plate_velocity): + """Shift the mid-ocean ridge point at 'stat'to the left (or right) of the ridge and calculate the spreading rate for the left (or right) plate.""" + + # If there's no left (or right) plate then don't return a left (or right) shifted point. + if plate.not_located_in_resolved_topology(): + return + # Note: The plate velocity should not be None since there IS a left (or right) plate. + + # Get the pole of rotation to shift the point to the left (or right) of the ridge. + # The pole is perpendicular to the boundary point and the left (or right) plate velocity. + boundary_rotation_pole = pygplates.Vector3D.cross( # type: ignore + stat.boundary_point.to_xyz(), plate_velocity + ) + # If we can't shift the point then don't return a shifted point. + # This can happen when the plate velocity is zero. + # This would result in the point not shifting off the ridge and thus + # not definitely inside the plate for later topological reconstruction. + if boundary_rotation_pole.is_zero_magnitude(): + return + + # Shift the boundary point into the left (or right) plate. + boundary_rotation_pole = boundary_rotation_pole.to_normalised() + boundary_rotation = pygplates.FiniteRotation( # type: ignore + boundary_rotation_pole.to_xyz(), + boundary_rotation_angle_radians, + ) + shifted_boundary_point = boundary_rotation * stat.boundary_point + + # If the left (or right) shifted point is OUTSIDE the left (or right) plate then don't return it. + # We want the shifted point to be inside the left (or right) plate so that it can later be reconstructed + # by that same plate to the next time step. + # + # Note: We give preference to networks since 'stat' (pygplates.PlateBoundaryStatistic) does that. + # This is because it's possible that a network overlays a rigid plate. + plate_resolved_topology = ( + plate.located_in_resolved_network() + or plate.located_in_resolved_boundary() + ) + if ( + not plate_resolved_topology # should always be false since we checked plate.not_located_in_resolved_topology() above + or plate_resolved_topology.get_point_location(shifted_boundary_point) + != plate + ): + return + + # The spreading rate is the left (or right) plate velocity relative to the ridge velocity. + spreading_rate = (plate_velocity - stat.boundary_velocity).get_magnitude() + + return shifted_boundary_point, spreading_rate # Points and their spreading rates that emerge from MORs at this time. shifted_mor_points = [] point_spreading_rates = [] - # Resolve topologies to the current time. - topological_snapshot = self.plate_reconstruction.topological_snapshot(time) - shared_boundary_sections = ( - topological_snapshot.get_resolved_topological_sections() - ) + # Iterate over MOR boundary points and shift each point to the left and right plates. + for stat in mor_boundary_statistics: - # pygplates.ResolvedTopologicalSection objects. - for shared_boundary_section in shared_boundary_sections: - if ( - shared_boundary_section.get_feature().get_feature_type() - == pygplates.FeatureType.gpml_mid_ocean_ridge - ): - spreading_feature = shared_boundary_section.get_feature() - - # Find the stage rotation of the spreading feature in the - # frame of reference of its geometry at the current - # reconstruction time (the MOR is currently actively spreading). - # The stage pole can then be directly geometrically compared - # to the *reconstructed* spreading geometry. - stage_rotation = separate_ridge_transform_segments.get_stage_rotation_for_reconstructed_geometry( - spreading_feature, self.plate_reconstruction.rotation_model, time + # Shift boundary point to the left and caculate spreading rate for the left plate. + left_shifted_mor_point_and_spreading_rate = ( + _calc_shifted_mor_point_and_spreading_rate( + stat, + stat.left_plate, + stat.left_plate_velocity, ) - if not stage_rotation: - # Skip current feature - it's not a spreading feature. - continue - - # Get the stage pole of the stage rotation. - # Note that the stage rotation is already in frame of - # reference of the *reconstructed* geometry at the spreading time. - stage_pole, _ = stage_rotation.get_euler_pole_and_angle() - - # One way rotates left and the other right, but don't know - # which - doesn't matter in our example though. - rotate_slightly_off_mor_one_way = pygplates.FiniteRotation( - stage_pole, np.radians(0.01) + ) + if left_shifted_mor_point_and_spreading_rate: + left_shifted_mor_point, left_spreading_rate = ( + left_shifted_mor_point_and_spreading_rate ) - rotate_slightly_off_mor_opposite_way = ( - rotate_slightly_off_mor_one_way.get_inverse() + shifted_mor_points.append(left_shifted_mor_point) + point_spreading_rates.append(left_spreading_rate) + + # Shift boundary point to the right and caculate spreading rate for the right plate. + right_shifted_mor_point_and_spreading_rate = ( + _calc_shifted_mor_point_and_spreading_rate( + stat, + stat.right_plate, + stat.right_plate_velocity, ) - - # Iterate over the shared sub-segments. - for ( - shared_sub_segment - ) in shared_boundary_section.get_shared_sub_segments(): - # Tessellate MOR section. - mor_points = pygplates.MultiPointOnSphere( - shared_sub_segment.get_resolved_geometry().to_tessellated( - np.radians(self.ridge_sampling) - ) - ) - - coords = mor_points.to_lat_lon_list() - lats = [i[0] for i in coords] - lons = [i[1] for i in coords] - boundary_feature = shared_boundary_section.get_feature() - left_plate = boundary_feature.get_left_plate(None) - right_plate = boundary_feature.get_right_plate(None) - if left_plate is not None and right_plate is not None: - # Get the spreading rates for all points in this sub segment - ( - spreading_rates, - _, - ) = tools.calculate_spreading_rates( - time=time, - lons=lons, - lats=lats, - left_plates=[left_plate] * len(lons), - right_plates=[right_plate] * len(lons), - rotation_model=self.plate_reconstruction.rotation_model, - delta_time=self._ridge_time_step, - ) - - else: - spreading_rates = [np.nan] * len(lons) - - # Loop through all but the 1st and last points in the current sub segment - for point, rate in zip( - mor_points.get_points()[1:-1], - spreading_rates[1:-1], - ): - # Add the point "twice" to the main shifted_mor_points list; once for a L-side - # spread, another for a R-side spread. Then add the same spreading rate twice - # to the list - this therefore assumes spreading rate is symmetric. - shifted_mor_points.append( - rotate_slightly_off_mor_one_way * point - ) - shifted_mor_points.append( - rotate_slightly_off_mor_opposite_way * point - ) - point_spreading_rates.extend([rate] * 2) + ) + if right_shifted_mor_point_and_spreading_rate: + right_shifted_mor_point, right_spreading_rate = ( + right_shifted_mor_point_and_spreading_rate + ) + shifted_mor_points.append(right_shifted_mor_point) + point_spreading_rates.append(right_spreading_rate) logger.info(f"Finished building MOR seedpoints at {time} Ma!") From df8d83560c801a5ebeffb898b827e7d87c58212a Mon Sep 17 00:00:00 2001 From: John Cannon Date: Fri, 29 May 2026 14:03:18 +1000 Subject: [PATCH 2/5] Shift MOR seed points using boundary normal (not plate velocities). Anchors shifted points in their respective left/right plates, and no longer need to test if shifted point within left/right plate. --- gplately/oceans.py | 31 +++++++++++-------------------- 1 file changed, 11 insertions(+), 20 deletions(-) diff --git a/gplately/oceans.py b/gplately/oceans.py index 095f95b4..9bf9b31c 100644 --- a/gplately/oceans.py +++ b/gplately/oceans.py @@ -978,7 +978,9 @@ def _build_mid_ocean_ridge_seed_points( # How much to rotate each ridge point off the ridge to get a point definitely inside the plate for later topological reconstruction. boundary_rotation_angle_radians = np.radians(0.01) - def _calc_shifted_mor_point_and_spreading_rate(stat, plate, plate_velocity): + def _calc_shifted_mor_point_and_spreading_rate( + stat, is_left_plate, plate, plate_velocity + ): """Shift the mid-ocean ridge point at 'stat'to the left (or right) of the ridge and calculate the spreading rate for the left (or right) plate.""" # If there's no left (or right) plate then don't return a left (or right) shifted point. @@ -987,10 +989,14 @@ def _calc_shifted_mor_point_and_spreading_rate(stat, plate, plate_velocity): # Note: The plate velocity should not be None since there IS a left (or right) plate. # Get the pole of rotation to shift the point to the left (or right) of the ridge. - # The pole is perpendicular to the boundary point and the left (or right) plate velocity. + # The pole is perpendicular to both the boundary point (vector from Earth centre) and + # the boundary normal to the left (or right) plate. boundary_rotation_pole = pygplates.Vector3D.cross( # type: ignore - stat.boundary_point.to_xyz(), plate_velocity + stat.boundary_point.to_xyz(), stat.boundary_normal ) + # The boundary normal 'stat.boundary_normal' is defined to point to the left, so reverse it if the plate is to the right. + if not is_left_plate: + boundary_rotation_pole = -boundary_rotation_pole # If we can't shift the point then don't return a shifted point. # This can happen when the plate velocity is zero. # This would result in the point not shifting off the ridge and thus @@ -1006,23 +1012,6 @@ def _calc_shifted_mor_point_and_spreading_rate(stat, plate, plate_velocity): ) shifted_boundary_point = boundary_rotation * stat.boundary_point - # If the left (or right) shifted point is OUTSIDE the left (or right) plate then don't return it. - # We want the shifted point to be inside the left (or right) plate so that it can later be reconstructed - # by that same plate to the next time step. - # - # Note: We give preference to networks since 'stat' (pygplates.PlateBoundaryStatistic) does that. - # This is because it's possible that a network overlays a rigid plate. - plate_resolved_topology = ( - plate.located_in_resolved_network() - or plate.located_in_resolved_boundary() - ) - if ( - not plate_resolved_topology # should always be false since we checked plate.not_located_in_resolved_topology() above - or plate_resolved_topology.get_point_location(shifted_boundary_point) - != plate - ): - return - # The spreading rate is the left (or right) plate velocity relative to the ridge velocity. spreading_rate = (plate_velocity - stat.boundary_velocity).get_magnitude() @@ -1039,6 +1028,7 @@ def _calc_shifted_mor_point_and_spreading_rate(stat, plate, plate_velocity): left_shifted_mor_point_and_spreading_rate = ( _calc_shifted_mor_point_and_spreading_rate( stat, + True, # is_left_plate stat.left_plate, stat.left_plate_velocity, ) @@ -1054,6 +1044,7 @@ def _calc_shifted_mor_point_and_spreading_rate(stat, plate, plate_velocity): right_shifted_mor_point_and_spreading_rate = ( _calc_shifted_mor_point_and_spreading_rate( stat, + False, # is_left_plate stat.right_plate, stat.right_plate_velocity, ) From 06a3d8a76beddd6cdb4c747c123536cff50f86c3 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Fri, 29 May 2026 14:54:42 +1000 Subject: [PATCH 3/5] Set spreading rate to NaN for MOR seed point if no left/right plate. --- gplately/oceans.py | 65 +++++++++++++++++++++------------------------- 1 file changed, 29 insertions(+), 36 deletions(-) diff --git a/gplately/oceans.py b/gplately/oceans.py index 9bf9b31c..7e5ae9cf 100644 --- a/gplately/oceans.py +++ b/gplately/oceans.py @@ -978,42 +978,47 @@ def _build_mid_ocean_ridge_seed_points( # How much to rotate each ridge point off the ridge to get a point definitely inside the plate for later topological reconstruction. boundary_rotation_angle_radians = np.radians(0.01) - def _calc_shifted_mor_point_and_spreading_rate( - stat, is_left_plate, plate, plate_velocity - ): - """Shift the mid-ocean ridge point at 'stat'to the left (or right) of the ridge and calculate the spreading rate for the left (or right) plate.""" + def _calc_shifted_mor_point_and_spreading_rate(stat, is_left_plate): + """Shift the mid-ocean ridge point at 'stat' to the left (or right) of the ridge and calculate the spreading rate for the left (or right) plate. - # If there's no left (or right) plate then don't return a left (or right) shifted point. - if plate.not_located_in_resolved_topology(): - return - # Note: The plate velocity should not be None since there IS a left (or right) plate. + If there's no left (or right) plate then return NaN for the spreading rate. + """ # Get the pole of rotation to shift the point to the left (or right) of the ridge. # The pole is perpendicular to both the boundary point (vector from Earth centre) and # the boundary normal to the left (or right) plate. + # + # Note: The boundary normal should always be perpendicular to the boundary point. + # So the cross product should be a non-zero vector and hence can be normalised. boundary_rotation_pole = pygplates.Vector3D.cross( # type: ignore stat.boundary_point.to_xyz(), stat.boundary_normal - ) + ).to_normalised() # The boundary normal 'stat.boundary_normal' is defined to point to the left, so reverse it if the plate is to the right. if not is_left_plate: boundary_rotation_pole = -boundary_rotation_pole - # If we can't shift the point then don't return a shifted point. - # This can happen when the plate velocity is zero. - # This would result in the point not shifting off the ridge and thus - # not definitely inside the plate for later topological reconstruction. - if boundary_rotation_pole.is_zero_magnitude(): - return # Shift the boundary point into the left (or right) plate. - boundary_rotation_pole = boundary_rotation_pole.to_normalised() + # This should shift it off the ridge and thus definitely inside the plate for later topological reconstruction. boundary_rotation = pygplates.FiniteRotation( # type: ignore boundary_rotation_pole.to_xyz(), boundary_rotation_angle_radians, ) shifted_boundary_point = boundary_rotation * stat.boundary_point - # The spreading rate is the left (or right) plate velocity relative to the ridge velocity. - spreading_rate = (plate_velocity - stat.boundary_velocity).get_magnitude() + # Get the left (or right) plate velocity. + if is_left_plate: + plate_velocity = stat.left_plate_velocity + else: + plate_velocity = stat.right_plate_velocity + + if plate_velocity: + # The spreading rate is the left (or right) plate velocity relative to the ridge velocity. + spreading_rate = ( + plate_velocity - stat.boundary_velocity + ).get_magnitude() + else: + # There's no left (or right) plate so the spreading rate is NaN. + spreading_rate = np.nan return shifted_boundary_point, spreading_rate @@ -1025,36 +1030,24 @@ def _calc_shifted_mor_point_and_spreading_rate( for stat in mor_boundary_statistics: # Shift boundary point to the left and caculate spreading rate for the left plate. - left_shifted_mor_point_and_spreading_rate = ( + left_shifted_mor_point, left_spreading_rate = ( _calc_shifted_mor_point_and_spreading_rate( stat, True, # is_left_plate - stat.left_plate, - stat.left_plate_velocity, ) ) - if left_shifted_mor_point_and_spreading_rate: - left_shifted_mor_point, left_spreading_rate = ( - left_shifted_mor_point_and_spreading_rate - ) - shifted_mor_points.append(left_shifted_mor_point) - point_spreading_rates.append(left_spreading_rate) + shifted_mor_points.append(left_shifted_mor_point) + point_spreading_rates.append(left_spreading_rate) # Shift boundary point to the right and caculate spreading rate for the right plate. - right_shifted_mor_point_and_spreading_rate = ( + right_shifted_mor_point, right_spreading_rate = ( _calc_shifted_mor_point_and_spreading_rate( stat, False, # is_left_plate - stat.right_plate, - stat.right_plate_velocity, ) ) - if right_shifted_mor_point_and_spreading_rate: - right_shifted_mor_point, right_spreading_rate = ( - right_shifted_mor_point_and_spreading_rate - ) - shifted_mor_points.append(right_shifted_mor_point) - point_spreading_rates.append(right_spreading_rate) + shifted_mor_points.append(right_shifted_mor_point) + point_spreading_rates.append(right_spreading_rate) logger.info(f"Finished building MOR seedpoints at {time} Ma!") From 282106a363bd3e7acf06c0e80aebedd31033a6e8 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Mon, 1 Jun 2026 15:40:10 +1000 Subject: [PATCH 4/5] Avoid anomalously high spreading rates due to spurious MOR velocities. This is due to plate ID(s) not existing in rotation tree at 'time + 1' which can happen when MOR only has rotations for the times that it exists (however 'time-of-appearance + 1' used for its velocity at 'time=time-of-appearance' then has no rotation). --- gplately/oceans.py | 136 +++++++++++++++++++++++++++++++++++++-------- 1 file changed, 114 insertions(+), 22 deletions(-) diff --git a/gplately/oceans.py b/gplately/oceans.py index 7e5ae9cf..91cbb317 100644 --- a/gplately/oceans.py +++ b/gplately/oceans.py @@ -37,7 +37,7 @@ ) from .lib.reconstruct_continents import ReconstructContinents from .parallel import get_num_cpus -from .ptt import continent_contours, separate_ridge_transform_segments +from .ptt import continent_contours from .ptt.utils import points_in_polygons, points_spatial_tree from .tools import _deg2pixels from .utils import seafloor_grid_utils @@ -964,24 +964,90 @@ def _build_mid_ocean_ridge_seed_points( The left and right spreading rates are calculated as the difference between the left and right plate velocities and the ridge velocity at the point. """ + velocity_delta_time = 1.0 # Myr - the time interval over which to calculate velocity changes for the spreading rate calculation. + # Generate statistics at uniformly spaced points along mid-ocean ridges. - mor_boundary_statistics = self.plate_reconstruction.topological_snapshot( + mor_boundary_statistics_dict = self.plate_reconstruction.topological_snapshot( time, # Ignore topological slab boundaries since they are not *plate* boundaries... include_topological_slab_boundaries=False, ).calculate_plate_boundary_statistics( uniform_point_spacing_radians=np.radians(self.ridge_sampling), + velocity_delta_time=velocity_delta_time, + velocity_delta_time_type=pygplates.VelocityDeltaTimeType.t_plus_delta_t_to_t, # [t+1, t] # type: ignore velocity_units=pygplates.VelocityUnits.kms_per_my, # km/Myr is the same as mm/yr # type: ignore boundary_section_filter=pygplates.FeatureType.gpml_mid_ocean_ridge, # type: ignore + return_shared_sub_segment_dict=True, ) + def _could_mor_have_anomalous_velocity(mor_feature): + """Returns True if the MOR feature could have an anomalous velocity. + + The MOR velocity is calculated from [time + velocity_delta_time, time] because we specified + 'pygplates.VelocityDeltaTimeType.t_plus_delta_t_to_t' above. + The current 'time' should have valid rotations, but 'time + velocity_delta_time' might not. + If it doesn't then the MOR velocity can be all wrong (eg, an anomalously high value) and we return True. + """ + + # The MOR is reconstructing either by plate ID or by half stage rotation. + mor_reconstruction_method = mor_feature.get_reconstruction_method(None) + if ( + mor_reconstruction_method is not None + and mor_reconstruction_method.startswith("HalfStageRotation") + ): # reconstruction is by half stage rotation... + + # See if MOR feature has left and right plate ids. + left_plate_id = mor_feature.get_left_plate(None) + right_plate_id = mor_feature.get_right_plate(None) + if left_plate_id is not None and right_plate_id is not None: + # Check that valid rotations exist for the left and right plate IDs at 'time' and 'time + velocity_delta_time'. + # We can do this in one call to 'get_rotation' by using 'use_identity_for_missing_plate_ids=False' and checking if the result is None. + if ( + self.plate_reconstruction.rotation_model.get_rotation( + time, + right_plate_id, + time + velocity_delta_time, + left_plate_id, + use_identity_for_missing_plate_ids=False, + ) + is None + ): + return True + # else the MOR feature doesn't have left and right plate IDs so it'll end up using zero plate IDs + # which, while inaccurate, are not anomalous. + + else: # reconstruction is by plate ID... + + # See if MOR feature has a reconstruction plate ID. + reconstruction_plate_id = mor_feature.get_reconstruction_plate_id(None) + if reconstruction_plate_id is not None: + # Check that valid rotations exist for the reconstruction plate ID at 'time' and 'time + velocity_delta_time'. + # We can do this in one call to 'get_rotation' by using 'use_identity_for_missing_plate_ids=False' and checking if the result is None. + if ( + self.plate_reconstruction.rotation_model.get_rotation( + time, + reconstruction_plate_id, + time + velocity_delta_time, + use_identity_for_missing_plate_ids=False, + ) + is None + ): + return True + # else the MOR feature doesn't have a reconstruction plate ID so it'll end up using a zero plate ID + # which, while inaccurate, is not anomalous. + + return False + # How much to rotate each ridge point off the ridge to get a point definitely inside the plate for later topological reconstruction. boundary_rotation_angle_radians = np.radians(0.01) - def _calc_shifted_mor_point_and_spreading_rate(stat, is_left_plate): + def _calc_shifted_mor_point_and_spreading_rate( + stat, is_left_plate, mor_has_anomalous_velocity + ): """Shift the mid-ocean ridge point at 'stat' to the left (or right) of the ridge and calculate the spreading rate for the left (or right) plate. If there's no left (or right) plate then return NaN for the spreading rate. + Also if 'mor_has_anomalous_velocity' is True then return NaN for the spreading rate. """ # Get the pole of rotation to shift the point to the left (or right) of the ridge. @@ -1005,6 +1071,12 @@ def _calc_shifted_mor_point_and_spreading_rate(stat, is_left_plate): ) shifted_boundary_point = boundary_rotation * stat.boundary_point + # If the MOR has an anomalous velocity then use NaN for the spreading rate. + # This is because 'stat.boundary_velocity' could return an anomalously high velocity. + if mor_has_anomalous_velocity: + spreading_rate = np.nan + return shifted_boundary_point, spreading_rate + # Get the left (or right) plate velocity. if is_left_plate: plate_velocity = stat.left_plate_velocity @@ -1018,6 +1090,7 @@ def _calc_shifted_mor_point_and_spreading_rate(stat, is_left_plate): ).get_magnitude() else: # There's no left (or right) plate so the spreading rate is NaN. + # This can happen if there are gaps in plate/network coverage in a global topological model. spreading_rate = np.nan return shifted_boundary_point, spreading_rate @@ -1026,28 +1099,47 @@ def _calc_shifted_mor_point_and_spreading_rate(stat, is_left_plate): shifted_mor_points = [] point_spreading_rates = [] - # Iterate over MOR boundary points and shift each point to the left and right plates. - for stat in mor_boundary_statistics: + # Iterate over the MOR shared sub-segments. + for ( + mor_shared_sub_segment, + mor_boundary_statistics, + ) in mor_boundary_statistics_dict.items(): - # Shift boundary point to the left and caculate spreading rate for the left plate. - left_shifted_mor_point, left_spreading_rate = ( - _calc_shifted_mor_point_and_spreading_rate( - stat, - True, # is_left_plate - ) + # If the MOR feature could have an anomalously high velocity then we should use NaN for its spreading rate. + # + # The MOR velocity is calculated as a stage rotation from [time + velocity_delta_time, time]. + # The current 'time' should have valid rotations, but 'time + velocity_delta_time' might not because + # it might be prior to the time of appearance of the MOR feature (and hence might have no rotations in rotation file). + # That can cause the MOR velocity to be all wrong (eg, an anomalously high value). + mor_has_anomalous_velocity = _could_mor_have_anomalous_velocity( + mor_shared_sub_segment.get_feature() ) - shifted_mor_points.append(left_shifted_mor_point) - point_spreading_rates.append(left_spreading_rate) - - # Shift boundary point to the right and caculate spreading rate for the right plate. - right_shifted_mor_point, right_spreading_rate = ( - _calc_shifted_mor_point_and_spreading_rate( - stat, - False, # is_left_plate + + # Iterate over the MOR boundary points in the current MOR shared sub-segment + # and shift each point to the left and right plates. + for stat in mor_boundary_statistics: + + # Shift boundary point to the left and caculate spreading rate for the left plate. + left_shifted_mor_point, left_spreading_rate = ( + _calc_shifted_mor_point_and_spreading_rate( + stat, + True, # is_left_plate + mor_has_anomalous_velocity, + ) ) - ) - shifted_mor_points.append(right_shifted_mor_point) - point_spreading_rates.append(right_spreading_rate) + shifted_mor_points.append(left_shifted_mor_point) + point_spreading_rates.append(left_spreading_rate) + + # Shift boundary point to the right and caculate spreading rate for the right plate. + right_shifted_mor_point, right_spreading_rate = ( + _calc_shifted_mor_point_and_spreading_rate( + stat, + False, # is_left_plate + mor_has_anomalous_velocity, + ) + ) + shifted_mor_points.append(right_shifted_mor_point) + point_spreading_rates.append(right_spreading_rate) logger.info(f"Finished building MOR seedpoints at {time} Ma!") From 62b59b7bd79205c0d4a097d3c4542f064ca90f7e Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Thu, 25 Jun 2026 05:16:13 +0000 Subject: [PATCH 5/5] 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 eef54ec6..ff30d24a 100644 --- a/Notebooks/14-RuleBasedGPMLProcessingPipeline.ipynb +++ b/Notebooks/14-RuleBasedGPMLProcessingPipeline.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "596b9073", + "id": "3d561d78", "metadata": {}, "source": [ "\n", @@ -19,7 +19,7 @@ }, { "cell_type": "markdown", - "id": "079076fc", + "id": "4ca19dbe", "metadata": {}, "source": [ "\n", @@ -30,7 +30,7 @@ }, { "cell_type": "markdown", - "id": "d01e3ef6", + "id": "743d5977", "metadata": { "lines_to_next_cell": 0 }, @@ -41,7 +41,7 @@ { "cell_type": "code", "execution_count": null, - "id": "fc7bb6ec", + "id": "680ba013", "metadata": {}, "outputs": [], "source": [ @@ -118,7 +118,7 @@ }, { "cell_type": "markdown", - "id": "4fb25240", + "id": "b6bff13f", "metadata": {}, "source": [ "### Search and filter feature collection with pre-defined filters\n", @@ -155,7 +155,7 @@ { "cell_type": "code", "execution_count": null, - "id": "3d2f5be9", + "id": "7b9d374f", "metadata": {}, "outputs": [], "source": [ @@ -363,7 +363,7 @@ }, { "cell_type": "markdown", - "id": "d0d7417e", + "id": "45bd7b67", "metadata": { "lines_to_next_cell": 0 }, @@ -382,7 +382,7 @@ { "cell_type": "code", "execution_count": null, - "id": "66e0286f", + "id": "d87cad2f", "metadata": {}, "outputs": [], "source": [ @@ -418,7 +418,7 @@ }, { "cell_type": "markdown", - "id": "378fb426", + "id": "d8fabe5d", "metadata": { "lines_to_next_cell": 2 }, @@ -435,7 +435,7 @@ { "cell_type": "code", "execution_count": null, - "id": "cd898dc8", + "id": "16182778", "metadata": { "lines_to_next_cell": 0 }, @@ -468,7 +468,7 @@ }, { "cell_type": "markdown", - "id": "c002c664", + "id": "cbb8ee58", "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": "b2e7b237", + "id": "54860ec7", "metadata": {}, "outputs": [], "source": [ @@ -515,7 +515,7 @@ }, { "cell_type": "markdown", - "id": "f10da5f6", + "id": "b4fb459c", "metadata": { "lines_to_next_cell": 0 }, @@ -530,7 +530,7 @@ { "cell_type": "code", "execution_count": null, - "id": "b8f5c679", + "id": "b49c55b4", "metadata": { "lines_to_next_cell": 0 }, @@ -561,7 +561,7 @@ }, { "cell_type": "markdown", - "id": "fb750a60", + "id": "d5125a0d", "metadata": {}, "source": [ "#### Topological reference map\n", @@ -577,7 +577,7 @@ { "cell_type": "code", "execution_count": null, - "id": "71664c67", + "id": "fd649426", "metadata": {}, "outputs": [], "source": [ @@ -595,7 +595,7 @@ }, { "cell_type": "markdown", - "id": "bb4e9803", + "id": "ac0440dc", "metadata": { "lines_to_next_cell": 2 }, @@ -610,7 +610,7 @@ { "cell_type": "code", "execution_count": null, - "id": "1bc1436c", + "id": "6835cda5", "metadata": {}, "outputs": [], "source": [ @@ -641,7 +641,7 @@ }, { "cell_type": "markdown", - "id": "f58ba5c8", + "id": "eca479ff", "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": "1babb6c7", + "id": "a512cac1", "metadata": { "lines_to_next_cell": 2 }, @@ -666,7 +666,7 @@ { "cell_type": "code", "execution_count": null, - "id": "037550fa", + "id": "22b6779a", "metadata": {}, "outputs": [], "source": [ @@ -694,7 +694,7 @@ }, { "cell_type": "markdown", - "id": "91cc50b7", + "id": "4937648e", "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": "5ba1bf9e", + "id": "469147d9", "metadata": { "lines_to_next_cell": 2 }, @@ -720,7 +720,7 @@ { "cell_type": "code", "execution_count": null, - "id": "5aa09b8c", + "id": "194b9e5e", "metadata": { "lines_to_next_cell": 0 }, @@ -766,7 +766,7 @@ }, { "cell_type": "markdown", - "id": "f1ce7142", + "id": "5f3e7dec", "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 ea01accc..bf1c073a 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": "54d4e461", + "id": "a1c4908b", "metadata": {}, "source": [ "This notebook demonstrates how to use the PyGMT plotting maps in GPlately.\n", @@ -14,7 +14,7 @@ }, { "cell_type": "markdown", - "id": "b7374f9e", + "id": "fe001ac0", "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": "d9c53bf2", + "id": "26ee75f4", "metadata": {}, "outputs": [], "source": [