154 test for 20 release

README.md

@@ -42,7 +42,7 @@ eprint = {https://rmets.onlinelibrary.wiley.com/doi/pdf/10.1002/gdj3.185},
 ## Dependencies
 
 - [pyGPlates](https://www.gplates.org/docs/pygplates/pygplates_getting_started.html#installation)
-- [plate-model-manager](https://pypi.org/project/plate-model-manager/)
+- [plate-model-manager](https://pypi.org/project/plate-model-manager/) >= 1.2.0
 - [Shapely](https://shapely.readthedocs.io/en/stable/project.html#installing-shapely)
 - [NumPy](https://numpy.org/install/) > 1.16
 - [SciPy](https://scipy.org/install/) > 1.0

gplately/reconstruction.py

@@ -59,8 +59,10 @@ def __init__(
         self.static_polygons = _load_FeatureCollection(static_polygons)
 
     def __getstate__(self):
-        filenames = {"rotation_model": self.rotation_model.filenames,
-                     "anchor_plate_id": self.anchor_plate_id}
+        filenames = {
+            "rotation_model": self.rotation_model.filenames,
+            "anchor_plate_id": self.anchor_plate_id,
+        }
 
         if self.topology_features:
             filenames["topology_features"] = self.topology_features.filenames
@@ -79,7 +81,9 @@ def __getstate__(self):
 
     def __setstate__(self, state):
         # reinstate unpicklable items
-        self.rotation_model = _RotationModel(state["rotation_model"], default_anchor_plate_id=state["anchor_plate_id"])
+        self.rotation_model = _RotationModel(
+            state["rotation_model"], default_anchor_plate_id=state["anchor_plate_id"]
+        )
 
         self.anchor_plate_id = state["anchor_plate_id"]
         self.topology_features = None
@@ -175,6 +179,7 @@ def tessellate_subduction_zones(
         The delta time interval used for velocity calculations is, by default, assumed to be 1Ma.
         """
         from . import ptt as _ptt
+
         anchor_plate_id = kwargs.pop("anchor_plate_id", self.anchor_plate_id)
 
         if ignore_warnings:
@@ -474,6 +479,7 @@ def tessellate_mid_ocean_ridges(
             * length of arc segment (in degrees) that current point is on
         """
         from . import ptt as _ptt
+
         anchor_plate_id = kwargs.pop("anchor_plate_id", self.anchor_plate_id)
 
         if ignore_warnings:
@@ -614,7 +620,9 @@ def total_ridge_length(self, time, use_ptt=False, ignore_warnings=False):
 
             return total_ridge_length_kms
 
-    def reconstruct(self, feature, to_time, from_time=0, anchor_plate_id=None, **kwargs):
+    def reconstruct(
+        self, feature, to_time, from_time=0, anchor_plate_id=None, **kwargs
+    ):
         """Reconstructs regular geological features, motion paths or flowlines to a specific geological time.
 
         Parameters
@@ -865,10 +873,13 @@ def create_motion_path(
             query_plate_id = False
             plate_ids = np.ones(len(lons), dtype=int) * plate_id
 
-        if not anchor_plate_id:
+        if anchor_plate_id is None:
             anchor_plate_id = self.anchor_plate_id
 
         seed_points = zip(lats, lons)
+        if return_rate_of_motion is True:
+            StepTimes = np.empty(((len(time_array) - 1) * 2, len(lons)))
+            StepRates = np.empty(((len(time_array) - 1) * 2, len(lons)))
         for i, lat_lon in enumerate(seed_points):
             seed_points_at_digitisation_time = pygplates.PointOnSphere(
                 pygplates.LatLonPoint(float(lat_lon[0]), float(lat_lon[1]))
@@ -887,15 +898,15 @@ def create_motion_path(
                 seed_points_at_digitisation_time,
                 time_array,
                 valid_time=(time_array.max(), time_array.min()),
-                relative_plate=int(plate_id),
-                reconstruction_plate_id=int(anchor_plate_id),
+                relative_plate=int(anchor_plate_id),
+                reconstruction_plate_id=int(plate_id),
             )
 
             reconstructed_motion_paths = self.reconstruct(
                 motion_path_feature,
                 to_time=0,
                 reconstruct_type=pygplates.ReconstructType.motion_path,
-                anchor_plate_id=int(plate_id),
+                anchor_plate_id=int(anchor_plate_id),
             )
             # Turn motion paths in to lat-lon coordinates
             for reconstructed_motion_path in reconstructed_motion_paths:
@@ -908,9 +919,6 @@ def create_motion_path(
 
             # Obtain step-plot coordinates for rate of motion
             if return_rate_of_motion is True:
-                StepTimes = np.empty(((len(time_array) - 1) * 2, len(lons)))
-                StepRates = np.empty(((len(time_array) - 1) * 2, len(lons)))
-
                 # Get timestep
                 TimeStep = []
                 for j in range(len(time_array) - 1):
@@ -2573,6 +2581,7 @@ def _activate_deactivate_points(self):
 
     def _find_resolved_topologies_containing_points(self):
         from . import ptt as _ptt
+
         current_time = self.get_current_time()
 
         # Resolve the plate polygons for the current time.

tests-dir/unittest/common.py

@@ -1,7 +1,6 @@
 import sys
 from pathlib import Path
 
-import cartopy.crs as ccrs
 import matplotlib.pyplot as plt
 
 sys.path.insert(0, "../..")
@@ -11,6 +10,8 @@
 
 Path(OUTPUT_DIR).mkdir(parents=True, exist_ok=True)
 
+MODEL_REPO_DIR = "plate-model-repo"
+
 
 def save_fig(filename):
     output_file = f"{OUTPUT_DIR}/{Path(filename).stem}.png"

tests-dir/unittest/run_all.sh

@@ -14,4 +14,6 @@ $BASEDIR/test_plot_with_raster.py save
 
 $BASEDIR/test_reconstruct_points.py save
 
-$BASEDIR/test_tessellate.py save
+$BASEDIR/test_tessellate.py save
+
+$BASEDIR/test_motion_path.py save

tests-dir/unittest/test_anchor_plate_id.py

@@ -7,7 +7,7 @@
 from plate_model_manager import PlateModelManager
 
 sys.path.insert(0, "../")
-from common import save_fig
+from common import MODEL_REPO_DIR, save_fig
 
 import gplately
 
@@ -17,10 +17,9 @@
 
 def main(show=True):
     pm_manger = PlateModelManager()
-    model = pm_manger.get_model(MODEL_NAME)
+    model = pm_manger.get_model(MODEL_NAME, data_dir=MODEL_REPO_DIR)
     if not model:
         raise Exception(f"Unable to get model {MODEL_NAME}!!!")
-    model.set_data_dir("test-plate-model-folder")
 
     C2020_rotation_file = model.get_rotation_model()
     C2020_topology_features = model.get_topologies()

tests-dir/unittest/test_motion_path.py

@@ -0,0 +1,173 @@
+#!/usr/bin/env python3
+import sys
+
+import cartopy.crs as ccrs
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.gridspec import GridSpec
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+from plate_model_manager import PlateModelManager
+
+sys.path.insert(0, "../..")
+from common import MODEL_REPO_DIR, save_fig
+
+import gplately
+
+
+def main(show=True):
+    pm_manager = PlateModelManager()
+    muller2019_model = pm_manager.get_model("Matthews2016", data_dir=MODEL_REPO_DIR)
+
+    rotation_model = muller2019_model.get_rotation_model()
+    topology_features = muller2019_model.get_topologies()
+    static_polygons = muller2019_model.get_static_polygons()
+
+    coastlines = muller2019_model.get_layer("Coastlines")
+
+    # Call the PlateReconstruction object to create a plate motion model
+    model = gplately.reconstruction.PlateReconstruction(
+        rotation_model, topology_features, static_polygons
+    )
+
+    # Call the PlotTopologies object
+    time = 0
+    gplot = gplately.PlotTopologies(
+        model, coastlines=coastlines, continents=None, COBs=None, time=time
+    )
+
+    # *****************************************do the first plot***************************
+
+    lats = np.array([19])
+    lons = np.array([-155])
+
+    # Create the time array for the motion path - must be float
+    start_reconstruction_time = 0.0
+    time_step = 2.0
+    max_reconstruction_time = 100.0
+    time_array = np.arange(
+        start_reconstruction_time, max_reconstruction_time + time_step, time_step
+    )
+
+    # Creating the motion path and obtaining rate plot arrays using the gplately Points alternative
+    gpts = gplately.reconstruction.Points(model, lons, lats, time=0.0)
+    lon, lat, _, _ = gpts.motion_path(
+        time_array, anchor_plate_id=2, return_rate_of_motion=True
+    )
+
+    fig = plt.figure(figsize=(12, 6))
+    gs = GridSpec(2, 2, figure=fig)
+    ax = fig.add_subplot(gs[0, 0], projection=ccrs.PlateCarree(central_longitude=180))
+    ax.set_title("Motion path of the Hawaiian-Emperor Chain")
+
+    # --- Limit map extent
+    lon_min = 130.0
+    lon_max = 220
+    lat_min = 15.0
+    lat_max = 60.0
+    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
+    ax.gridlines(
+        draw_labels=True, xlocs=np.arange(-180, 180, 20), ylocs=np.arange(-90, 90, 10)
+    )
+
+    # --- Make sure motion path longitudes are wrapped correctly to the dateline
+    lon360 = gplately.tools.correct_longitudes_for_dateline(lon)
+
+    # --- plot motion path
+    ax.scatter(
+        lon360,
+        lat,
+        100,
+        marker=".",
+        c=time_array,
+        cmap=plt.cm.inferno,
+        edgecolor="k",
+        transform=ccrs.PlateCarree(),
+        vmin=time_array[0],
+        vmax=time_array[-1],
+        zorder=4,
+    )
+
+    gplot.plot_coastlines(ax, color="DarkGrey")
+
+    # ***************************** do the second plot **********************************
+
+    lats = np.array([30, 30, 30])
+    lons = np.array([73, 78, 83])
+
+    # Plate ID attributes
+    # motion of the Indian craton (501) with respect to the Northern European craton (301)
+    plate_ID = 501
+    anchor_plate_ID = 301
+
+    # Create the time array for the motion path
+    start_reconstruction_time = 0.0
+    time_step = 5.0
+    max_reconstruction_time = 105.0
+    time_array = np.arange(
+        start_reconstruction_time, max_reconstruction_time + 1, time_step
+    )
+
+    # Get the latitudes and longitudes of all points along the motion path
+    lon, lat, step_times, step_rates_of_motion = model.create_motion_path(
+        lons, lats, time_array, plate_ID, anchor_plate_ID, return_rate_of_motion=True
+    )
+
+    ax2 = fig.add_subplot(gs[0, 1], projection=ccrs.PlateCarree(central_longitude=180))
+    ax2.set_title("Motion path of the Indian craton to the North European Craton")
+
+    gplot.plot_coastlines(ax2, color="DarkGrey")
+
+    # --- Limit map extent
+    lon_min = 0.0
+    lon_max = 100
+    lat_min = -40.0
+    lat_max = 40.0
+    ax2.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
+
+    ax2.gridlines(
+        draw_labels=True, xlocs=np.arange(-180, 180, 20), ylocs=np.arange(-90, 90, 10)
+    )
+    # --- Make sure negtive motion path longitudes are wrapped correctly to the dateline
+    lons = gplately.tools.correct_longitudes_for_dateline(lons)
+
+    # --- plot motion path
+    for i in np.arange(0, len(lons)):
+        ax2.scatter(
+            lon[:, i],
+            lat[:, i],
+            200,
+            marker="*",
+            c=time_array,
+            cmap=plt.cm.inferno,
+            edgecolor="C{}".format(i),
+            linewidth=0.5,
+            transform=ccrs.PlateCarree(),
+            vmin=time_array[0],
+            vmax=time_array[-1],
+            zorder=4,
+        )
+
+    # ********************************do the third plot**********************************
+
+    ax3 = fig.add_subplot(gs[1, :])
+    ax3.set_title("Rate of motion between seed points and the North European Craton")
+    ax3.plot(step_times, step_rates_of_motion)
+    ax3.set_xlabel("Time (Ma)", fontsize=12)
+    ax3.set_ylabel("Rate of motion (km/Myr)", fontsize=12)
+    ax3.invert_xaxis()
+    ax3.set_xlim([max_reconstruction_time, 0])
+    ax3.grid(alpha=0.3)
+
+    plt.subplots_adjust(wspace=0.25)
+
+    if show:
+        plt.show()
+    else:
+        save_fig(__file__)
+
+
+if __name__ == "__main__":
+    if len(sys.argv) == 2 and sys.argv[1] == "save":
+        main(show=False)
+    else:
+        main(show=True)

tests-dir/unittest/test_plate_model.py

@@ -1,12 +1,12 @@
 #!/usr/bin/env python
 
+from common import MODEL_REPO_DIR
 from plate_model_manager import PlateModelManager
 
 
 def main():
     pm_manger = PlateModelManager()
-    model = pm_manger.get_model("Muller2019")
-    model.set_data_dir("test-plate-model-folder")
+    model = pm_manger.get_model("Muller2019", data_dir=MODEL_REPO_DIR)
 
     print(model.get_avail_layers())
 

tests-dir/unittest/test_plot.py

@@ -7,7 +7,7 @@
 from plate_model_manager import PlateModelManager
 
 sys.path.insert(0, "../..")
-from common import save_fig
+from common import MODEL_REPO_DIR, save_fig
 
 from gplately import PlateReconstruction, PlotTopologies
 
@@ -21,8 +21,7 @@ def main(show=True):
     pm_manager = PlateModelManager()
 
     age = 55
-    model = pm_manager.get_model(MODEL_NAME)
-    model.set_data_dir("plate-model-repo")
+    model = pm_manager.get_model(MODEL_NAME, data_dir=MODEL_REPO_DIR)
 
     test_model = PlateReconstruction(
         model.get_rotation_model(),

tests-dir/unittest/test_plot_with_raster.py

@@ -4,7 +4,7 @@
 import cartopy.crs as ccrs
 import matplotlib.pyplot as plt
 import numpy as np
-from common import save_fig
+from common import MODEL_REPO_DIR, save_fig
 from plate_model_manager import PlateModelManager
 
 sys.path.insert(0, "../..")
@@ -14,7 +14,7 @@
 def main(show=True):
     # Call GPlately's PlateModelManager object and request data from the Müller et al. 2019 study
     pm_manager = PlateModelManager()
-    muller2019_model = pm_manager.get_model("Muller2019", data_dir="plate-model-repo")
+    muller2019_model = pm_manager.get_model("Muller2019", data_dir=MODEL_REPO_DIR)
     rotation_model = muller2019_model.get_rotation_model()
     topology_features = muller2019_model.get_topologies()
     static_polygons = muller2019_model.get_static_polygons()