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 warnings

#warnings.filterwarnings("ignore", category=UserWarning, module="gplately.utils.dev_warning")

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

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

Matplotlib is building the font cache; this may take a moment.

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=gplately.plot.get_age_grid_cmap(), 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=gplately.plot.get_age_grid_cmap(), 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=gplately.plot.get_age_grid_cmap(), vmin=0, vmax=200)
else:
    graster.plot(ax=ax, cmap=gplately.plot.get_age_grid_cmap(), 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.