2 - Plate Reconstructions¶

In this notebook, we will set up and use a plate reconstruction model with GPlately's PlateReconstruction object.

model = gplately.PlateReconstruction(
    rotation_model, # required
    topology_features, # optional
    static_polygons # optional
)

The PlateReconstruction object provides methods for reconstructing topology features to a specific geological time. To use this object, you simply need to supply a rotation model, a set of topology features (or feature collection), and a set of static polygons.

In [1]:
import gplately
from pathlib import Path
import numpy as np
import glob, os, warnings
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

from matplotlib.lines import Line2D

from plate_model_manager import PlateModelManager
from gplately.auxiliary import get_gplot

model_repo_dir = "plate-model-repo"
model_name = "zahirovic2022"
data_dir = Path("WorkflowData/02-Plate-Reconstructions")
data_dir.mkdir(parents=True, exist_ok=True)

To create a PlateReconstruction object, we need three file sets:

  1. A set of rotation files — files with the ".rot" extension
  2. A set of topology feature files — typically of type ".gpml" or ".gpmlz"
  3. A set of static polygons

We will demonstrate two ways to load plate model files below.

In [2]:
# Set to True to use Method 2 (load files from the local hard drive)
use_local_files = False

Method 1: Loading Files with PlateModelManager¶

You can use PlateModelManager to download the files required for a supported plate reconstruction model. Run the command pmm ls to list all supported plate models.

Use PlateModelManager's get_model() method to obtain a PlateModel object.

Now let's retrieve the following files from the PlateModel object:

  • Rotation files
  • Topology files
  • Static polygon files
In [3]:
# Use the model name to create a PlateModel object
pmm_model = PlateModelManager().get_model(model_name, data_dir=model_repo_dir)
assert pmm_model is not None, f"Failed to load the {model_name} plate model."
if not use_local_files:
    # Obtain all rotation files, topology features, and static polygons
    rotation_model = pmm_model.get_rotation_model()
    topology_features = pmm_model.get_topologies()
    static_polygons = pmm_model.get_static_polygons()

Method 2: Loading Local Files¶

The cell below shows how the glob and os libraries can be used to locate rotation files, topology feature files, and static polygon files from one or more directories on your computer.

In [4]:
if use_local_files:
    # For demonstration purposes only, we download the plate model files once.
    # This call will not re-download the files unless they have been updated on the server.
    PlateModelManager().get_model(model_name, data_dir=model_repo_dir).download_all_layers() # type: ignore

    model_path = os.path.join(model_repo_dir, model_name)
    rotation_model = glob.glob(os.path.join(model_path, "Rotations", "*.rot"))
    topology_features = []
    for topology_filename in glob.glob(
        os.path.join(model_path, "Topologies", "*.gpml")
    ):
        # Skip files whose name contains "Inactive"
        if "Inactive" not in topology_filename:
            topology_features.append(topology_filename)
    static_polygons = glob.glob(os.path.join(model_path, "StaticPolygons", "*.shp"))

Constructing a Plate Reconstruction Model Using the PlateReconstruction Object¶

Once we have our rotation model, topology features, and static polygons, we can supply them to the PlateReconstruction object to construct the plate motion model.

In [5]:
model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)

Reconstructing Feature Geometries¶

The plate motion model we created can be used to generate plate reconstructions through geological time. Let's reconstruct subduction zones and mid-ocean ridges to 50 Ma.

tessellate_subduction_zones() samples points along subduction zone trenches and returns subduction data at a given geological time as a 10-column, vertically-stacked tuple.

tessellate_mid_ocean_ridges() samples points along resolved spreading features (e.g. mid-ocean ridges) and returns spreading rates and ridge segment lengths at a given geological time as a 4-column, vertically-stacked tuple.

You may notice that the tessellated subduction zones and MORs in the plot below are missing some sections compared to the maps plotted by PlotTopologies. This is because the tessellation algorithm applies stricter rules and excludes some ineligible sections. See the online documentation for details.

In [6]:
time = 50  # Ma
subduction_data = model.tessellate_subduction_zones(time)
ridge_data = model.tessellate_mid_ocean_ridges(time)

print(subduction_data.shape)
print(ridge_data.shape)

fig = plt.figure(figsize=(8, 6), dpi=72)
ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=0))
ax.gridlines(  # type: ignore
    color="0.7",
    linestyle="--",
    xlocs=np.arange(-180, 180, 15),
    ylocs=np.arange(-90, 90, 15),
)
ax.set_global()  # type: ignore
ax.scatter(
    subduction_data[:, 0],  # longitude
    subduction_data[:, 1],  # latitude
    color="blue",
    s=1,
    transform=ccrs.PlateCarree(),
)
ax.scatter(
    ridge_data[:, 0],  # longitude
    ridge_data[:, 1],  # latitude
    color="red",
    s=1,
    transform=ccrs.PlateCarree(),
)
plt.title(f"Tessellated Subduction Zones and MORs at {time} Ma")
plt.show()

(11093, 10)
(11745, 4)
No description has been provided for this image

Plotting Plate Reconstructions with the PlotTopologies Object¶

Let's visualise this reconstruction on a GeoAxes plot using GPlately's PlotTopologies object. To create the object, we need to supply:

  • the PlateReconstruction plate motion model we just created,
  • a coastline filename or <pygplates.FeatureCollection> object,
  • a continent filename or <pygplates.FeatureCollection> object,
  • a continent-ocean boundary (COBs) filename or <pygplates.FeatureCollection> object, and
  • a specific reconstruction time (Ma).
gplot = gplately.PlotTopologies(
    plate_reconstruction,
    coastlines=None,
    continents=None,
    COBs=None,
    time=None,
    anchor_plate_id=None, # by default uses anchor plate of 'plate_reconstruction'
)

We demonstrate the same methods used above to locate the coastline, continent, and COB files.

Method 1: Loading Files with PlateModelManager¶

We already defined a PlateModel object above (pmm_model) to obtain a rotation model, topology features, and static polygons. Let's reuse this object to locate the coastlines, continents, and COBs downloaded to the local folder.

In [7]:
if not use_local_files:
    assert pmm_model is not None, f"Failed to load the {model_name} plate model."
    coastlines = pmm_model.get_coastlines()
    continents = pmm_model.get_continental_polygons()
    COBs = pmm_model.get_COBs()

Method 2: Loading Local Files¶

Load the local coastlines, continents, and COBs files.

In [8]:
if use_local_files:
    coastlines = glob.glob(os.path.join(model_path, "Coastlines", "*.shp"))
    continents = glob.glob(os.path.join(model_path, "ContinentalPolygons", "*.shp"))
    COBs = glob.glob(os.path.join(model_path, "COBs", "*.shp"))

Define the PlotTopologies Object¶

Let's create a PlotTopologies object named gplot and set it up to visualise geologic features at 50 Ma:

In [9]:
# Create the PlotTopologies object
time = 50  # Ma
gplot = gplately.PlotTopologies(model, coastlines=coastlines, continents=continents, COBs=COBs, time=time)

To plot using GPlately's PlotTopologies object, first create a GeoAxes plot (here we call it ax) and select a projection using Cartopy. This is the plot we then supply to our gplot object.

In [10]:
with warnings.catch_warnings():
    warnings.filterwarnings("ignore", category=UserWarning)
    # Set up a GeoAxes plot
    fig = plt.figure(figsize=(8, 6), dpi=100)
    ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=0))
    ax.gridlines(  # type: ignore
        color="0.7",
        linestyle="--",
        xlocs=np.arange(-180, 180, 15),
        ylocs=np.arange(-90, 90, 15),
    )
    plt.title(f"Subduction Zones and MORs at {time} Ma")

    # Plot shapefile features, subduction zones, and MOR boundaries at the given time
    gplot.time = time  # Ma
    gplot.plot_continent_ocean_boundaries(ax, color="b", alpha=0.05)
    gplot.plot_continents(ax, facecolor="palegoldenrod", alpha=0.2)
    gplot.plot_coastlines(ax, color="DarkKhaki")
    gplot.plot_all_topological_sections(
        ax,
        plot_subduction_teeth=True,
        other_kwargs={"color": "grey", "linewidth": 0.8},
        ridge_kwargs={"color": "red", "linewidth": 1.0},
        transform_kwargs={"color": "green", "linewidth": 1.0},
        trench_kwargs={"color": "blue", "linewidth": 1.0},
    )
    ax.set_global()  # type: ignore
No description has been provided for this image

If you have moviepy installed, you can create a GIF that illustrates plate motions through geological time. Let's reconstruct plate movements up to 100 Ma in intervals of 10 Ma!

In [11]:
def generate_frame(output_filename, gplot_m):
    time = gplot_m.time
    # Set up a GeoAxes plot
    fig = plt.figure(figsize=(8, 6), dpi=72)
    ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=0))
    ax.gridlines(  # type: ignore
        color="0.7",
        linestyle="--",
        xlocs=np.arange(-180, 180, 15),
        ylocs=np.arange(-90, 90, 15),
    )
    plt.title(f"Subduction Zones and MORs at {int(time)} Ma")

    # Plot shapefile features, subduction zones, and MOR boundaries at the given time
    gplot_m.plot_continent_ocean_boundaries(ax, color="b", alpha=0.05)
    gplot_m.plot_continents(ax, facecolor="palegoldenrod", alpha=0.2)
    gplot_m.plot_coastlines(ax, color="DarkKhaki", alpha=0.4)
    gplot_m.plot_all_topological_sections(
        ax,
        plot_subduction_teeth=True,
        other_kwargs={"color": "grey", "linewidth": 0.8},
        ridge_kwargs={"color": "red", "linewidth": 1.0},
        transform_kwargs={"color": "green", "linewidth": 1.0},
        trench_kwargs={"color": "blue", "linewidth": 1.0},
    )
    ax.set_global()  # type: ignore
    plt.savefig(output_filename, bbox_inches="tight")
    plt.close()
In [12]:
import tempfile
from IPython.display import Image

try:
    from moviepy.editor import ImageSequenceClip  # moviepy 1.x
except ImportError:
    from moviepy import ImageSequenceClip  # moviepy 2.x #type: ignore

# Time variables
oldest_seed_time = 100  # Ma
time_step = 10  # Ma

with tempfile.TemporaryDirectory() as tmpdir, warnings.catch_warnings():
    warnings.filterwarnings("ignore", category=UserWarning)

    frame_list = []

    # Create a plot for each time step
    gplot_m = get_gplot("Zahirovic2022", time=oldest_seed_time, model_repo_dir=model_repo_dir)
    for time in np.arange(oldest_seed_time, 0.0, -time_step):
        gplot_m.time = time
        print("Generating %d Ma frame..." % time)
        frame_filename = os.path.join(tmpdir, "frame_%d_Ma.png" % time)
        generate_frame(frame_filename, gplot_m)
        frame_list.append(frame_filename)

    video_filename = data_dir / "subd_mor_boundary_features.gif"

    clip = ImageSequenceClip(frame_list, fps=5)
    clip.write_gif(video_filename)

    print("Displaying the animation below...")
    with open(video_filename, "rb") as f:
        display(Image(data=f.read(), format="png", width=500, height=250))

Generating 100 Ma frame...

Generating 90 Ma frame...

Generating 80 Ma frame...

Generating 70 Ma frame...

Generating 60 Ma frame...

Generating 50 Ma frame...

Generating 40 Ma frame...

Generating 30 Ma frame...

Generating 20 Ma frame...

Generating 10 Ma frame...

MoviePy - Building file WorkflowData/02-Plate-Reconstructions/subd_mor_boundary_features.gif with imageio.

frame_index:   0%|          | 0/10 [00:00<?, ?it/s, now=None]
                                                             












Displaying the animation below...












No description has been provided for this image























Comparing Two Different Plate Models¶
Let's create two PlotTopologies objects using two different plate models and plot the resulting maps side by side for comparison.


















In [13]:






reconstruction_time = 100  # Ma

# You can change the model names below to compare different plate models.
# Make sure both models contain the required data files.
model_name_1 = "Zahirovic2022"
model_name_2 = "Muller2025"

gplot_1 = get_gplot(model_name_1, time=reconstruction_time, model_repo_dir=model_repo_dir)
gplot_2 = get_gplot(model_name_2, time=reconstruction_time, model_repo_dir=model_repo_dir)

pmm_1 = gplot_1.plate_reconstruction.plate_model
pmm_2 = gplot_2.plate_reconstruction.plate_model

agegrid_1 = gplately.Raster(data=pmm_1.get_raster("AgeGrids", reconstruction_time))
agegrid_2 = gplately.Raster(data=pmm_2.get_raster("AgeGrids", reconstruction_time))

warnings.filterwarnings("ignore", category=UserWarning)

# Set up a GeoAxes plot
fig = plt.figure(figsize=(20, 10), dpi=72)

# --- First subplot ---
ax1 = fig.add_subplot(121, projection=ccrs.Mollweide(central_longitude=0))
plt.title(f"{model_name_1} {reconstruction_time} Ma", fontsize=15)
ax1.set_title("a", loc="left", fontsize="16", weight="demi")
ax1.set_global()  # type: ignore

# Plot seafloor age grid, coastlines, subduction zones, and MOR/transform boundaries at the reconstruction time
gplot_1.plot_coastlines(ax1, color="0.5")
im = gplot_1.plot_grid(ax1, agegrid_1.data, cmap="YlGnBu", vmin=0, vmax=250, alpha=0.4)
gplot_1.plot_all_topological_sections(
    ax1,
    plot_subduction_teeth=True,
    other_kwargs={"color": "grey", "linewidth": 0.8},
    ridge_kwargs={"color": "red", "linewidth": 1.0},
    transform_kwargs={"color": "green", "linewidth": 1.0},
    trench_kwargs={"color": "blue", "linewidth": 1.0},
)
gplot_1.plot_plate_motion_vectors(
    ax1, spacingX=10, spacingY=10, normalise=False, zorder=4, alpha=0.4
)

# --- Second subplot ---
ax2 = fig.add_subplot(122, projection=ccrs.Mollweide(central_longitude=0))
plt.title(f"{model_name_2} {reconstruction_time} Ma", fontsize=15)
ax2.set_title("b", loc="left", fontsize="16", weight="demi")
ax2.set_global()  # type: ignore

# Plot seafloor age grid, coastlines, subduction zones, and MOR/transform boundaries at the reconstruction time
gplot_2.plot_coastlines(ax2, color="0.5")

# Use the age grid for the current time step
im = gplot_2.plot_grid(ax2, agegrid_2.data, cmap="YlGnBu", vmin=0, vmax=250, alpha=0.4)
gplot_2.plot_all_topological_sections(
    ax2,
    plot_subduction_teeth=True,
    other_kwargs={"color": "grey", "linewidth": 0.8},
    ridge_kwargs={"color": "red", "linewidth": 1.0},
    transform_kwargs={"color": "green", "linewidth": 1.0},
    trench_kwargs={"color": "blue", "linewidth": 1.0},
)
gplot_2.plot_plate_motion_vectors(
    ax2, spacingX=10, spacingY=10, normalise=False, zorder=10, alpha=0.4
)


# --- Plot properties ---
plt.subplots_adjust(wspace=0.025)  # spacing between subplots

# Colorbar settings
cb_ax = fig.add_axes((0.17, 0.25, 0.15, 0.02))
assert im is not None
cb = fig.colorbar(im, cax=cb_ax, orientation="horizontal", shrink=0.4, pad=0.05)
cb.set_label(label="Age (Ma)", fontsize=15)
tick_values = np.arange(0, 251, 50)
cb.set_ticks(
    tick_values.tolist(), labels=[str(int(v)) for v in tick_values], fontsize=15
)

# Legend settings
legend_elements = [
    Line2D(
        [0],
        [0],
        linestyle="-",
        color="r",
        label="Mid-ocean Ridges and Transform Boundaries",
    ),
    Line2D(
        [0],
        [0],
        marker="^",
        linestyle="-",
        color="blue",
        label="Subduction Zones with Polarity Teeth",
        markerfacecolor="blue",
        markersize=5,
    ),
]

lg = fig.legend(
    handles=legend_elements,
    bbox_to_anchor=(0.65, 0.3),
    ncol=1,
    fontsize=15,
    frameon=False,
)
plt.show()

























2026-08-15 10:29:19 - matplotlib.font_manager - WARNING - findfont: Failed to find font weight demi, now using 700.












No description has been provided for this image