add two examples (#122)

* add two examples of plate-model-manager

* save unittest output

.gitignore

@@ -3,3 +3,5 @@ gplately.egg-info
 build
 __pycache__
 .ipynb_checkpoints
+
+output

Notebooks/Examples/introducing-plate-model-manager.py

@@ -0,0 +1,39 @@
+#!/usr/bin/env python3
+
+from plate_model_manager import PlateModelManager
+
+
+def main():
+    pm_manager = PlateModelManager()
+
+    print("available models: ")
+    print("*****************")
+    for name in pm_manager.get_available_model_names():
+        print(name)
+    print()
+
+    model = pm_manager.get_model("Muller2019")
+    model.set_data_dir("plate-model-repo")
+
+    print("available layers in model Muller2019:")
+    print("*************************************")
+    for layer in model.get_avail_layers():
+        print(layer)
+    print()
+
+    print("rotation file(s):")
+    print("*****************")
+    print(model.get_rotation_model())
+    print()
+
+    print("StaticPolygons file(s):")
+    print("***********************")
+    print(model.get_layer("StaticPolygons"))
+    print()
+
+    # for layer in model.get_avail_layers():
+    #    print(model.get_layer(layer))
+
+
+if __name__ == "__main__":
+    main()

Notebooks/Examples/working-with-plate-model-manager.py

@@ -0,0 +1,51 @@
+#!/usr/bin/env python3
+
+
+import cartopy.crs as ccrs
+import matplotlib.pyplot as plt
+from plate_model_manager import PlateModelManager
+
+from gplately import PlateReconstruction, PlotTopologies
+
+# This is a simple example of how to use the Plate Model Manager with GPlately
+
+
+def main():
+    # Here is how to use PlateModelManager to create PlateReconstruction and PlotTopologies objects
+    pm_manager = PlateModelManager()
+    model = pm_manager.get_model("Muller2019")
+    model.set_data_dir("plate-model-repo")
+
+    age = 55
+    test_model = PlateReconstruction(
+        model.get_rotation_model(),
+        topology_features=model.get_layer("Topologies"),
+        static_polygons=model.get_layer("StaticPolygons"),
+    )
+    gplot = PlotTopologies(
+        test_model,
+        coastlines=model.get_layer("Coastlines"),
+        COBs=model.get_layer("COBs"),
+        time=age,
+    )
+
+    # Now do some plotting
+    fig = plt.figure(figsize=(12, 6), dpi=72)
+    ax = fig.add_subplot(111, projection=ccrs.Robinson(central_longitude=180))
+
+    gplot.plot_continent_ocean_boundaries(ax, color="cornflowerblue")
+    gplot.plot_coastlines(ax, color="black")
+    gplot.plot_ridges_and_transforms(ax, color="red")
+    gplot.plot_trenches(ax, color="orange")
+    gplot.plot_subduction_teeth(ax, color="orange")
+    ax.set_global()
+
+    ids = set([f.get_reconstruction_plate_id() for f in gplot.topologies])
+    for id in ids:
+        gplot.plot_plate_id(ax, id, facecolor="None", edgecolor="lightgreen")
+    plt.title(f"{age} Ma")
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()

gplately/examples

@@ -0,0 +1 @@
+../Notebooks/Examples

gplately/plot.py

@@ -331,6 +331,7 @@ def _calculate_triangle_vertices(
         length = geometry.length
         tessellated_x = []
         tessellated_y = []
+
         for distance in np.arange(spacing, length, spacing):
             point = Point(geometry.interpolate(distance))
             tessellated_x.append(point.x)
@@ -2361,9 +2362,19 @@ def plot_subduction_teeth(
 
         central_meridian = _meridian_from_ax(ax)
         tessellate_degrees = np.rad2deg(spacing)
-        # michael chin made this change. if spacing is in meters, it is too large to plot the triangles
-        spacing = math.degrees(spacing)
-        # spacing = spacing * EARTH_RADIUS * 1e3
+
+        try:
+            projection = ax.projection
+        except AttributeError:
+            print(
+                "The ax.projection does not exist. You must set project to plot Cartopy maps, such as ax = plt.subplot(211, projection=cartopy.crs.PlateCarree())"
+            )
+            projection = None
+
+        if isinstance(projection, ccrs.PlateCarree):
+            spacing = math.degrees(spacing)
+        else:
+            spacing = spacing * EARTH_RADIUS * 1e3
 
         if aspect is None:
             aspect = 2.0 / 3.0
@@ -2389,6 +2400,7 @@ def plot_subduction_teeth(
             "l",
             height,
             spacing,
+            projection=projection,
             ax=ax,
             color=color,
             **kwargs,
@@ -2399,6 +2411,7 @@ def plot_subduction_teeth(
             "r",
             height,
             spacing,
+            projection=projection,
             ax=ax,
             color=color,
             **kwargs,

tests-dir/unittest/common.py

@@ -0,0 +1,19 @@
+import sys
+from pathlib import Path
+
+import cartopy.crs as ccrs
+import matplotlib.pyplot as plt
+
+sys.path.insert(0, "../..")
+
+
+OUTPUT_DIR = "output"
+
+Path(OUTPUT_DIR).mkdir(parents=True, exist_ok=True)
+
+
+def save_fig(filename):
+    output_file = f"{OUTPUT_DIR}/{Path(filename).stem}.png"
+    plt.gcf().savefig(output_file, dpi=120, bbox_inches="tight")  # transparent=True)
+    print(f"Done! The {output_file} has been saved.")
+    plt.close(plt.gcf())

tests-dir/unittest/run_all.sh

@@ -0,0 +1,5 @@
+#!/bin/bash
+
+./test_subduction_teeth.py save
+
+./test_plot.py save 

tests-dir/unittest/test_plot.py

@@ -1,28 +1,28 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 
-import matplotlib.pyplot as plt
-import cartopy.crs as ccrs
 import sys
 
+import cartopy.crs as ccrs
+import matplotlib.pyplot as plt
 from plate_model_manager import PlateModelManager
 
-sys.path.insert(0, "../")
+sys.path.insert(0, "../..")
+from common import save_fig
+
 from gplately import PlateReconstruction, PlotTopologies
 
 # test the plot function with the new PlateModel class
 
 
-def main():
+def main(show=True):
     pm_manager = PlateModelManager()
 
     age = 55
     model = pm_manager.get_model("Muller2019")
     model.set_data_dir("plate-model-repo")
 
-    rotation_model = model.get_rotation_model()
-
     test_model = PlateReconstruction(
-        rotation_model,
+        model.get_rotation_model(),
         topology_features=model.get_layer("Topologies"),
         static_polygons=model.get_layer("StaticPolygons"),
     )
@@ -33,7 +33,7 @@ def main():
         time=age,
     )
 
-    fig = plt.figure(figsize=(10, 10), dpi=100)
+    fig = plt.figure(figsize=(10, 5), dpi=96)
     ax = fig.add_subplot(111, projection=ccrs.Robinson(central_longitude=180))
 
     gplot.plot_continent_ocean_boundaries(ax, color="cornflowerblue")
@@ -46,8 +46,16 @@ def main():
     ids = set([f.get_reconstruction_plate_id() for f in gplot.topologies])
     for id in ids:
         gplot.plot_plate_id(ax, id, facecolor="None", edgecolor="lightgreen")
-    plt.show()
+    plt.title(f"{age} Ma")
+
+    if show:
+        plt.show()
+    else:
+        save_fig(__file__)
 
 
 if __name__ == "__main__":
-    main()
+    if len(sys.argv) == 2 and sys.argv[1] == "save":
+        main(show=False)
+    else:
+        main(show=True)

tests-dir/unittest/test_subduction_teeth.py

@@ -4,35 +4,54 @@
 
 import cartopy
 import matplotlib.pyplot as plt
+from plate_model_manager import PlateModelManager
+
+sys.path.insert(0, "../..")
+from common import save_fig
 
-sys.path.insert(0, "../")
 import gplately
 
+# test plotting subduction teeth
 
-def main():
-    gdownload = gplately.download.DataServer("Muller2019")
-    (
-        C2020_rotation_file,
-        C2020_topology_features,
-        C2020_static_polygons,
-    ) = gdownload.get_plate_reconstruction_files()
 
-    fig, ax = plt.subplots(
-        subplot_kw={"projection": cartopy.crs.PlateCarree()}, figsize=(12, 8)
-    )
-    ax.set_extent([50, 105, -10, 40], crs=cartopy.crs.PlateCarree())
+def main(show=True):
+    pm_manager = PlateModelManager()
+    model = pm_manager.get_model("Muller2019")
+    model.set_data_dir("plate-model-repo")
+
+    plt.figure(figsize=(12, 8))
+
+    ax1 = plt.subplot(211, projection=cartopy.crs.PlateCarree())
+    ax2 = plt.subplot(212, projection=cartopy.crs.Robinson())
+
     plate_model = gplately.PlateReconstruction(
-        C2020_rotation_file, C2020_topology_features, C2020_static_polygons
+        model.get_rotation_model(),
+        topology_features=model.get_layer("Topologies"),
+        static_polygons=model.get_layer("StaticPolygons"),
     )
     plot_plates = gplately.PlotTopologies(plate_model, time=100)
     plot_plates.time = 100
-    plot_plates.plot_ridges_and_transforms(ax, color="r")
-    plot_plates.plot_trenches(ax, color="b")
-    plot_plates.plot_faults(ax, color="k")
-    plot_plates.plot_subduction_teeth(ax, color="green")
 
-    plt.show()
+    ax1.set_extent([50, 105, -10, 40], crs=cartopy.crs.PlateCarree())
+    plot_plates.plot_ridges_and_transforms(ax1, color="r")
+    plot_plates.plot_trenches(ax1, color="b")
+    plot_plates.plot_faults(ax1, color="k")
+    plot_plates.plot_subduction_teeth(ax1, color="green")
+
+    ax2.set_global()
+    plot_plates.plot_ridges_and_transforms(ax2, color="r")
+    plot_plates.plot_trenches(ax2, color="b")
+    plot_plates.plot_faults(ax2, color="k")
+    plot_plates.plot_subduction_teeth(ax2, color="green")
+
+    if show:
+        plt.show()
+    else:
+        save_fig(__file__)
 
 
 if __name__ == "__main__":
-    main()
+    if len(sys.argv) == 2 and sys.argv[1] == "save":
+        main(show=False)
+    else:
+        main(show=True)