1 - Getting Started

Welcome to GPlately.

GPlately provides an object-oriented interface for plate reconstruction workflows. In this notebook, we introduce several core objects you will use frequently:

• PlateReconstruction - reconstruct features and tessellate mid-ocean ridges and subduction zones
• Points - partition points onto plates and rotate them back through time
• Raster - read NetCDF grids, interpolate, and resample data
• PlotTopologies - plot topologies (for example, ridges, trenches, and subduction teeth) on maps
• PlateModelManager - download plate models and data files (for example, .rot, .gpml, .shp, and .nc)

import os, warnings
from pathlib import Path
import cartopy.crs as ccrs
import gplately
import matplotlib.pyplot as plt
import numpy as np
from plate_model_manager import PlateModelManager
from gplately.plot.gmt_cpt import get_cmap_from_gmt_cpt

# Download the age-grid CPT file if it is not present. This CPT is used for map plotting below.
data_dir = Path("WorkflowData")
data_dir.mkdir(parents=True, exist_ok=True)
cpt_file = data_dir / "agegrid.cpt"
if not os.path.isfile(cpt_file):
    import urllib.request

    urllib.request.urlretrieve(
        "https://raw.githubusercontent.com/GPlates/gplately/refs/heads/master/tests-dir/unittest/create-age-grids-video/agegrid.cpt",
        cpt_file,
    )

Tectonic Plate Reconstructions

To initialize a plate reconstruction model, provide a rotation model, plate topologies, and static polygons. You can download these files into your local folder using GPlately's PlateModelManager.

# Use PlateModelManager to request data from the Zahirovic2022 plate model
pm_manager = PlateModelManager()
z22_model = pm_manager.get_model("Zahirovic2022", data_dir="plate-model-repo")
assert z22_model is not None
rotation_model = z22_model.get_rotation_model()
topology_features = z22_model.get_topologies()
static_polygons = z22_model.get_static_polygons()

# Tessellate subduction zones at 0.05 degrees.
tessellation_threshold_radians = np.radians(0.05)

model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)

Now let's retrieve subduction zones and mid-ocean ridges at 10 Ma.

time = 10

# These methods return rich topology data; see the documentation for details.
subduction_data = model.tessellate_subduction_zones(
    time, tessellation_threshold_radians=tessellation_threshold_radians
)
ridge_data = model.tessellate_mid_ocean_ridges(
    time, tessellation_threshold_radians=tessellation_threshold_radians
)

Mapping

The PlotTopologies object uses the plate model we defined, along with coastlines, continents, and COBs, to compute topologies at a given reconstruction time.

PlotTopologies supports both cartopy and PyGMT as plotting engines.

In this example, we use cartopy. Define a figure and pass axes to the plotting routines. Common layers include:

• coastlines
• continents
• ridges and transforms
• trenches
• subduction teeth
• NetCDF grids
• plate motion vectors

You can also pass optional keyword arguments as usual.

# Retrieve feature layers for the PlotTopologies object
assert z22_model is not None
coastlines = z22_model.get_layer("Coastlines")
continents = z22_model.get_layer("ContinentalPolygons")
COBs = z22_model.get_layer("COBs")

# Create the PlotTopologies object
gplot = gplately.plot.PlotTopologies(
    model, coastlines=coastlines, continents=continents, COBs=COBs
)

# Download an age-grid NetCDF raster from the Zahirovic2022 model and create a Raster object.
agegrid = gplately.Raster(data=z22_model.get_raster("AgeGrids", time))

Set the time attribute to reconstruct all topologies at a specified reconstruction time.

Important: Before plotting, set gplot.time or provide time during initialization.

gplot.time = 10 # Ma

Create a map with key geological information.

with warnings.catch_warnings():
    warnings.filterwarnings("ignore", category=UserWarning)
    fig = plt.figure(figsize=(8, 6))

    ax = fig.add_subplot(111, projection=ccrs.Mollweide(190))

    gplot.plot_continents(ax, facecolor="0.8")
    gplot.plot_coastlines(ax, color="0.5", linewidth=0.5)
    gplot.plot_all_topological_sections(
        ax,
        plot_subduction_teeth=True,
        other_kwargs={"color": "grey", "linewidth": 0.8},
        ridge_kwargs={"color": "black", "linewidth": 1.0},
        transform_kwargs={"color": "green", "linewidth": 1.0},
        trench_kwargs={"color": "blue", "linewidth": 1.0},
    )
    im = gplot.plot_grid(
        ax, agegrid.data, cmap=get_cmap_from_gmt_cpt(cpt_file), vmin=0, vmax=200
    )
    gplot.plot_plate_motion_vectors(
        ax, spacingX=10, spacingY=10, normalise=True, zorder=10, alpha=0.5
    )
    if im:
        fig.colorbar(
            im, orientation="horizontal", shrink=0.4, pad=0.05, label="Age (Ma)"
        )
    time_value = gplot.time
    if time_value is None:
        raise ValueError("gplot.time is not set")
    ax.set_title(f"{int(time_value)} Ma")
    plt.show()

with warnings.catch_warnings():
    warnings.filterwarnings('ignore', category=UserWarning)
    # Update time and recompute topologies.
    time = 100
    gplot.time = time
    assert z22_model is not None
    agegrid = gplately.Raster(data=z22_model.get_raster("AgeGrids",time))
    
    fig = plt.figure(figsize=(8,6))
    
    ax = fig.add_subplot(111, projection=ccrs.Mollweide(190))
    
    gplot.plot_continents(ax, facecolor='0.8')
    gplot.plot_coastlines(ax, color='0.5')
    gplot.plot_all_topological_sections(
            ax,
            plot_subduction_teeth=True,
            other_kwargs={"color": "grey", "linewidth": 0.8},
            ridge_kwargs={"color": "black", "linewidth": 1.0},
            transform_kwargs={"color": "green", "linewidth": 1.0},
            trench_kwargs={"color": "blue", "linewidth": 1.0},
    )
    im = gplot.plot_grid(ax, agegrid.data, cmap=get_cmap_from_gmt_cpt(cpt_file), vmin=0, vmax=200)
    gplot.plot_plate_motion_vectors(ax, spacingX=10, spacingY=10, normalise=True, zorder=10, alpha=0.5)
    assert im is not None
    fig.colorbar(im, orientation='horizontal', shrink=0.4, pad=0.05, label='Age (Ma)')
    time_value = gplot.time
    assert time_value is not None
    ax.set_title(f"{int(time_value)} Ma")
    plt.show()

Working with Points

Now that we have defined our reconstruction object, we can reconstruct point data.

pt_lons = np.array([140.0, 150.0, 160.0])
pt_lats = np.array([-30.0, -40.0, -50.0])

gpts = gplately.Points(model, pt_lons, pt_lats)

# Point velocities at 0 Ma.
vel_x, vel_y = gpts.plate_velocity(0)  # type: ignore
vel_mag = np.hypot(vel_x, vel_y)

print("point velocity (cm/yr)", vel_mag)

point velocity (cm/yr) [7.51712692 6.89469761 6.2383027 ]

Plot point positions from time=0 to time=20.

rlons = np.empty((21, pt_lons.size))
rlats = np.empty((21, pt_lons.size))

for time in range(0, 21):
    rlons[time], rlats[time] = gpts.reconstruct(time, return_array=True)  # type: ignore

gplot.time = 0  # present day

fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111, projection=ccrs.Mercator(190))
ax.set_extent((130, 180, -60, -10)) # type: ignore

gplot.plot_coastlines(ax, color="0.8")

for i in range(0, len(pt_lons)):
    ax.plot(rlons[:, i], rlats[:, i], "o", transform=ccrs.PlateCarree())
plt.show()

Rasters

You can initialize a Raster object by providing raster data (either a NetCDF file path or a NumPy array), and optionally a PlateReconstruction model. The data attribute stores raster data as a 2D NumPy array.

In this example, we pass a NetCDF file path to the Raster object.

time = 100

# Download a NetCDF age-grid raster and assign a plate reconstruction `model`.
assert z22_model is not None
graster = gplately.Raster(data=z22_model.get_raster("AgeGrids", time))
graster.plate_reconstruction = model

# Access the underlying NumPy masked array using the data attribute.
print(type(graster.data))

# Alternatively, initialize a Raster from a NumPy array.
graster = gplately.Raster(
    data=graster.data,  # 2D NumPy array
    plate_reconstruction=model,  # PlateReconstruction object
    extent="global",  # equivalent to [-180, 180, -90, 90]
    time=100,  # time in Ma
)

<class 'numpy.ndarray'>

gplot.time = time

fig = plt.figure(figsize=(8, 6))

ax = fig.add_subplot(111, projection=ccrs.Robinson())

use_raster_plot = (
    False  # Set this flag to True to use Raster.plot() for raster plotting.
)

if not use_raster_plot:
    gplot.plot_grid(ax, graster, cmap=get_cmap_from_gmt_cpt(cpt_file), vmin=0, vmax=200)
else:
    graster.plot(ax=ax, cmap=get_cmap_from_gmt_cpt(cpt_file), vmin=0, vmax=200)

gplot.plot_coastlines(ax, edgecolor="0.5", facecolor="none", linewidth=0.5)
time_value = gplot.time
assert isinstance(time_value, (int, float)), "gplot.time must be a number"
ax.set_title(f"{int(time_value)} Ma")
plt.show()

Additional routines include:

• filling masked (NaN) regions
• interpolation
• resampling
• raster reconstruction

Use inplace=True for in-place operations that update internal data structures.