From e7112cac760e3317e312d6c89a30392f61c8fee7 Mon Sep 17 00:00:00 2001 From: michaelchin Date: Mon, 6 May 2024 16:22:33 +1000 Subject: [PATCH] remove 'continue' statement in reconstruction.py. see Issue #174 --- gplately/commands/create_age_grids.py | 18 ++++-- gplately/oceans.py | 8 +++ gplately/reconstruction.py | 80 ++++++++++++++------------- 3 files changed, 63 insertions(+), 43 deletions(-) diff --git a/gplately/commands/create_age_grids.py b/gplately/commands/create_age_grids.py index 2cb82f4d..ff1c782e 100644 --- a/gplately/commands/create_age_grids.py +++ b/gplately/commands/create_age_grids.py @@ -6,7 +6,7 @@ from typing import Optional, Sequence, Union import pygplates -from plate_model_manager import PlateModelManager +from plate_model_manager import PlateModel, PlateModelManager from gplately import PlateReconstruction, PlotTopologies, SeafloorGrid @@ -133,9 +133,9 @@ def add_parser(parser: argparse.ArgumentParser): __description__ = """Create age grids for a plate model. -example usage: - - gplately ag output -m muller2019 -e 0 -s 10 - - gplately ag plate-model-repo/muller2019/Rotations/Muller_etal_2019_CombinedRotations.rot plate-model-repo/muller2019/Topologies/Muller_etal_2019_PlateBoundaries_DeformingNetworks.gpmlz output -c plate-model-repo/muller2019/ContinentalPolygons/Global_EarthByte_GPlates_PresentDay_ContinentalPolygons_2019_v1.shp -e 0 -s 10 +Example usage: + - gplately ag output -m muller2019 -e 0 -s 10 + - gplately ag plate-model-repo/muller2019/Rotations/Muller_etal_2019_CombinedRotations.rot plate-model-repo/muller2019/Topologies/Muller_etal_2019_PlateBoundaries_DeformingNetworks.gpmlz output -c plate-model-repo/muller2019/ContinentalPolygons/Global_EarthByte_GPlates_PresentDay_ContinentalPolygons_2019_v1.shp -e 0 -s 10 """ @@ -158,8 +158,14 @@ def create_agegrids( """Create age grids for a plate model.""" if model_name: - pm_manager = PlateModelManager() - plate_model = pm_manager.get_model(model_name, data_dir="plate-model-repo") + try: + plate_model = PlateModelManager().get_model( + model_name, data_dir="plate-model-repo" + ) + except: + plate_model = PlateModel( + model_name, data_dir="plate-model-repo", readonly=True + ) rotations = pygplates.FeaturesFunctionArgument( plate_model.get_rotation_model() diff --git a/gplately/oceans.py b/gplately/oceans.py index 07a353fd..d469f6dd 100644 --- a/gplately/oceans.py +++ b/gplately/oceans.py @@ -1335,6 +1335,10 @@ def reconstruct_by_topologies(self): curr_points_including_inactive = ( topology_reconstruction.get_all_current_points() ) + logger.debug(f"the number of current active points is :{len(curr_points)}") + logger.debug( + f"the number of all current points is :{len(curr_points_including_inactive)}" + ) # Collect latitudes and longitudes of currently ACTIVE points in the ocean basin curr_lat_lon_points = [point.to_lat_lon() for point in curr_points] @@ -1428,6 +1432,10 @@ def reconstruct_by_topologies(self): "GPLATELY_DEBUG" in os.environ and os.environ["GPLATELY_DEBUG"].lower() == "true" ): + seafloor_ages = gridding_input_dictionary["SEAFLOOR_AGE"] + logger.debug( + f"The max and min values of seafloor age are: {np.max(seafloor_ages)} - {np.min(seafloor_ages)} ({topology_reconstruction.get_current_time()}Ma)" + ) save_age_grid_sample_points_to_gpml( gridding_input_dictionary["CURRENT_LONGITUDES"], gridding_input_dictionary["CURRENT_LATITUDES"], diff --git a/gplately/reconstruction.py b/gplately/reconstruction.py index db32f14f..4e917e24 100644 --- a/gplately/reconstruction.py +++ b/gplately/reconstruction.py @@ -1,6 +1,7 @@ """Tools that wrap up pyGplates and Plate Tectonic Tools functionalities for reconstructing features, working with point data, and calculating plate velocities at specific geological times. """ + import math import os import warnings @@ -52,7 +53,9 @@ def __init__( self.name = None self._anchor_plate_id = self._check_anchor_plate_id(anchor_plate_id) - self.rotation_model = _RotationModel(rotation_model, default_anchor_plate_id=anchor_plate_id) + self.rotation_model = _RotationModel( + rotation_model, default_anchor_plate_id=anchor_plate_id + ) self.topology_features = _load_FeatureCollection(topology_features) self.static_polygons = _load_FeatureCollection(static_polygons) @@ -97,7 +100,6 @@ def __setstate__(self, state): for polygon in state["static_polygons"]: self.static_polygons.add(_FeatureCollection(polygon)) - @property def anchor_plate_id(self): """Anchor plate ID for reconstruction. Must be an integer >= 0.""" @@ -106,7 +108,9 @@ def anchor_plate_id(self): @anchor_plate_id.setter def anchor_plate_id(self, anchor_plate): self._anchor_plate_id = self._check_anchor_plate_id(anchor_plate) - self.rotation_model = _RotationModel(self.rotation_model, default_anchor_plate_id=self._anchor_plate_id) + self.rotation_model = _RotationModel( + self.rotation_model, default_anchor_plate_id=self._anchor_plate_id + ) @staticmethod def _check_anchor_plate_id(id): @@ -115,7 +119,6 @@ def _check_anchor_plate_id(id): raise ValueError("Invalid anchor plate ID: {}".format(id)) return id - def tessellate_subduction_zones( self, time, @@ -1888,15 +1891,15 @@ def rotate_reference_frames( self, reconstruction_time, from_rotation_features_or_model=None, # filename(s), or pyGPlates feature(s)/collection(s) or a RotationModel - to_rotation_features_or_model=None, # filename(s), or pyGPlates feature(s)/collection(s) or a RotationModel + to_rotation_features_or_model=None, # filename(s), or pyGPlates feature(s)/collection(s) or a RotationModel from_rotation_reference_plate=0, to_rotation_reference_plate=0, non_reference_plate=701, output_name=None, return_array=False, ): - """Rotate a grid defined in one plate model reference frame - within a gplately.Raster object to another plate + """Rotate a grid defined in one plate model reference frame + within a gplately.Raster object to another plate reconstruction model reference frame. Parameters @@ -1904,12 +1907,12 @@ def rotate_reference_frames( reconstruction_time : float The time at which to rotate the reconstructed points. from_rotation_features_or_model : str, list of str, or instance of `pygplates.RotationModel` - A filename, or a list of filenames, or a pyGPlates + A filename, or a list of filenames, or a pyGPlates RotationModel object that defines the rotation model that the input grid is currently associated with. `self.plate_reconstruction.rotation_model` is default. to_rotation_features_or_model : str, list of str, or instance of `pygplates.RotationModel` - A filename, or a list of filenames, or a pyGPlates + A filename, or a list of filenames, or a pyGPlates RotationModel object that defines the rotation model that the input grid shall be rotated with. `self.plate_reconstruction.rotation_model` is default. @@ -1920,7 +1923,7 @@ def rotate_reference_frames( The desired reference plate for the plate model the points to be rotated to. Defaults to the anchor plate 0. non_reference_plate : int, default = 701 - An arbitrary placeholder reference frame with which + An arbitrary placeholder reference frame with which to define the "from" and "to" reference frames. output_name : str, default None If passed, the rotated points are saved as a gpml to this filename. @@ -1936,24 +1939,19 @@ def rotate_reference_frames( if to_rotation_features_or_model is None: to_rotation_features_or_model = self.plate_reconstruction.rotation_model - - # Create the pygplates.FiniteRotation that rotates + # Create the pygplates.FiniteRotation that rotates # between the two reference frames. - from_rotation_model = pygplates.RotationModel( - from_rotation_features_or_model - ) - to_rotation_model = pygplates.RotationModel( - to_rotation_features_or_model - ) + from_rotation_model = pygplates.RotationModel(from_rotation_features_or_model) + to_rotation_model = pygplates.RotationModel(to_rotation_features_or_model) from_rotation = from_rotation_model.get_rotation( - reconstruction_time, - non_reference_plate, - anchor_plate_id=from_rotation_reference_plate + reconstruction_time, + non_reference_plate, + anchor_plate_id=from_rotation_reference_plate, ) to_rotation = to_rotation_model.get_rotation( - reconstruction_time, - non_reference_plate, - anchor_plate_id=to_rotation_reference_plate + reconstruction_time, + non_reference_plate, + anchor_plate_id=to_rotation_reference_plate, ) reference_frame_conversion_rotation = to_rotation * from_rotation.get_inverse() @@ -1961,10 +1959,13 @@ def rotate_reference_frames( lons, lats = self.reconstruct( reconstruction_time, anchor_plate_id=from_rotation_reference_plate, - return_array=True) + return_array=True, + ) # convert FeatureCollection to MultiPointOnSphere - input_points = pygplates.MultiPointOnSphere((lat, lon) for lon, lat in zip(lons, lats)) + input_points = pygplates.MultiPointOnSphere( + (lat, lon) for lon, lat in zip(lons, lats) + ) # Rotate grid nodes to the other reference frame output_points = reference_frame_conversion_rotation * input_points @@ -1978,7 +1979,13 @@ def rotate_reference_frames( if return_array: return out_lon, out_lat else: - return Points(self.plate_reconstruction, out_lon, out_lat, time=reconstruction_time, plate_id=self.plate_id) + return Points( + self.plate_reconstruction, + out_lon, + out_lat, + time=reconstruction_time, + plate_id=self.plate_id, + ) # FROM RECONSTRUCT_BY_TOPOLOGIES.PY @@ -2108,9 +2115,9 @@ def __call__( 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 + 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) ) @@ -2118,9 +2125,9 @@ def __call__( 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 + 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). @@ -2567,7 +2574,6 @@ def reconstruct_to_next_time(self): # 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] - continue # 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. @@ -2760,6 +2766,6 @@ def _find_resolved_topologies_containing_points(self): 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() + self.curr_topology_plate_ids[point_index] = ( + curr_polygon.get_feature().get_reconstruction_plate_id() + )