remove 'continue' statement in reconstruction.py. see Issue #174

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()

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"],

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,28 +1891,28 @@ 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
         ----------
         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,35 +1939,33 @@ 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()
 
         # reconstruct points to reconstruction_time
         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,19 +2115,19 @@ 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)
             )
             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
+                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()
+            )