Convert grid reference frame

Plate-tectonic models reconstruct the positions of continents and ocean basins through deep time. They always store those positions relative to a chosen reference frame — and the choice matters.

Both of the reference frames you'll meet in this project — the palaeomagnetic reference frame (pmag) and the mantle reference frame — orient the reconstructed Earth against its spin axis. They differ in what data fix the orientation, and therefore in which motions appear in the reconstruction:

• The paleomagnetic reference frame (pmag) is derived from palaeomagnetic measurements — the orientation of magnetic minerals locked into rocks when they cooled or were deposited — and restores plate positions relative to the geographic poles. Palaeomagnetic data record the full latitudinal signal, so the reconstructed positions reflect changes in plate position relative to the spin axis, including relative plate motion, absolute plate motion, and True Polar Wander (TPW — solid-body re-orientations of the whole Earth relative to its spin axis). However, because Earth’s long-term magnetic field is approximately radially symmetric about the rotation axis, palaeomagnetic data constrain palaeolatitude and plate rotation, but not palaeolongitude. As a result, palaeomagnetic reference frames require an additional constraint to position continents in longitude, as well as a True Polar Wander correction, if more realistic plate motions relative to the mantle are the objective. However, paleomagnetic reference frames are most appropriate for comparisons with climate proxies, palaeolatitude data, or other plate models that use the pmag convention.

• The mantle reference frame is designed to isolate plate motion relative to the mantle itself, which is what geodynamic plate–mantle models need. TPW rotates the whole solid Earth (mantle included) as one body relative to the spin axis, so a mantle frame has to exclude TPW to give a clean plate-versus-mantle signal. Several methods can be used to achieve that separation — hotspot tracks, no-net rotation (NNR) models and mantle reference frame optimisation which combines NNR with other constraints such as minimising subduction zone migration. True Polar Wander can also be approximated from paleomagnetic data and then subtracted out. The choice of method is part of how a particular mantle frame is defined.

A grid (paleobathymetry, paleo-SST, an age grid, anything in NetCDF) is "in" whatever reference frame it was generated in. If you want to use it alongside another dataset that lives in a different frame, you first have to rotate the grid so the two line up at every age. This notebook does that rotation.

For paleoclimate you'll want a pmag frame. The default configuration below therefore rotates the Alfonso2024 mantle frame (anchor = 0) into the Alfonso2024 pmag frame (anchor = 701701), the convention used in the Müller group: in Alfonso 2024, anchor 701701 matches the reference frame Merdith v6-2 is published in.

---

How to use this notebook

1. Edit the USER CONFIGURATION cell below — that is the only place you need to change anything.
2. Run the cells top to bottom.
3. Inspect the comparison plot at the bottom to confirm the output looks right.

You can also rotate between two different plate models (e.g. an Alfonso2024 grid into the Müller 2022 frame), or use your own unpublished rotation files — see the next section.

---

Choosing rotations

The two big choices are which plate model and which anchor plate.

Plate model — pick one of two ways

• Use a plate model from gplately's online registry by setting ..._model_name to its short name (e.g. "Alfonso2024", "Muller2022", "Merdith2021") and leaving ..._rotation_files = None. gplately downloads and caches the rotation files locally, so the first run takes a moment and subsequent runs are instant.
• Use rotation files you already have on disk by setting ..._rotation_files to a list of .rot / .gpml paths and leaving ..._model_name = None. Useful for in-house or unpublished plate models that aren't in the registry.

You set this independently for the source (from_*) and target (to_*) — they don't have to match.

Anchor plate — which reference frame inside the model

A single plate model usually defines several reference frames; you choose between them by picking the anchor plate ID. Common values in the Müller-group conventions:

• 0 — mantle reference frame
• 701701 — Alfonso2024 pmag reference frame

Rotating within a single plate model means both ..._model_name strings are the same and only the two anchors differ. Rotating between models means they differ.

---

Quick examples

Default — Alfonso2024 mantle → Alfonso2024 pmag (already filled in below):

from_model_name      = "Alfonso2024";  from_rotation_files = None;  from_anchor = 0
to_model_name        = "Alfonso2024";  to_rotation_files   = None;  to_anchor   = 701701

Use unpublished rotation files (e.g. an in-house iteration):

from_model_name      = None
from_rotation_files  = ["data/XYZ/CombinedRotations.rot"]
from_anchor          = 0

Rotate between two different models (Cao 2024 → Müller 2025):

from_model_name      = "Cao2024";  from_rotation_files = None;  from_anchor = 0
to_model_name        = "Muller2025";   to_rotation_files   = None;  to_anchor   = 0

import os
import io
import contextlib
import warnings
import numpy as np
import gplately
import pygplates
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from plate_model_manager import PlateModelManager

USER CONFIGURATION — the only cell you need to edit

Set the source and target rotations, the anchor plates, the time range, the output resolution, and the input/output paths in the cell below. Each option has a brief comment next to it.

The cells after this one don't need any user input. Run them top to bottom.

# ============================================================================
# USER CONFIGURATION
# ============================================================================

# --- Source ('from') reference frame -----------------------------------------
# Provide EITHER a named plate model OR a list of rotation file paths.
# Set the unused one to None.
from_model_name      = "Alfonso2024"     # e.g. "Alfonso2024", "Muller2022", "Merdith2021"; None to use files
from_rotation_files  = None              # e.g. ["./CombinedRotations.rot"]; None to use a named model
from_anchor          = 0                 # plate ID of the source reference frame's anchor

# --- Target ('to') reference frame -------------------------------------------
to_model_name        = "Alfonso2024"
to_rotation_files    = None
to_anchor            = 701701            # plate ID of the target reference frame's anchor

# --- Reconstruction times (Ma) -----------------------------------------------
min_time             = 105
max_time             = 110
timestep_size        = 1                 # 1-Myr cadence

# --- Output grid resolution (degrees) ----------------------------------------
# 1.0° is the project standard for paleoclimate work in this repo: it matches
# the Merdith v6-2 paleogeography (1° native) and the climate-model output,
# so the rotated grids drop straight in alongside them with no resampling.
# Drill-site spacing (~1° equivalent at low–mid latitudes) means 1° also
# resolves the lithology-classifier depth feature at the granularity the
# training data can support.  Bump to 0.5° or 0.2° only if a downstream task
# specifically needs higher resolution.
grid_spacing_degrees = 0.1

# --- Input / output paths (add actual file names below) ----------------------
# Templates accept the reconstruction time via str.format(time), e.g. {:.2f}.
notebook_data_dir = "15-ConvertGridReferenceFrame-notebook-data"
input_grid_template  = f"{notebook_data_dir}/alfonso2024/Rasters/AgeGrids/Alfonso2024_SEAFLOOR_AGE_grid_"+"{:.2f}Ma.nc"
output_grid_template = f"{notebook_data_dir}/alfonso2024/Rasters/AgeGrids-PMAG/Alfonso2024_SEAFLOOR_AGE_grid_PMAG_"+"{:.2f}Ma.nc"

# Download the AgeGrids for this notebook
pm = PlateModelManager().get_model("Alfonso2024", data_dir=notebook_data_dir)
pm.get_rasters("AgeGrids", times=list(range(min_time,max_time+1,timestep_size)))
    
# Flag controlling execution mode: True uses Parallel (multi-process), False uses a single-process loop.
use_parallel = True
# ============================================================================

/home/runner/micromamba/envs/doc-env/lib/python3.14/site-packages/plate_model_manager/plate_model.py:501: RuntimeWarning: coroutine 'PlateModel.download_time_dependent_rasters.<locals>.f' was never awaited

Resolve and validate the configuration

This cell loads the rotation models, lists the time steps, makes the output directory, and prints a one-screen summary of what is about to run. If anything is misconfigured (a missing rotation file, both model_name and rotation_files set, etc.) it stops here with a clear error rather than later in the rotation loop.

It also has a quiet but important second job: when you use a model name from gplately's registry, this cell triggers the download into gplately's local cache. The parallel workers below then read from that cache instead of all racing to download the same files at once.

def build_rotation_model(model_name=None, rotation_files=None):
    """
    Construct a pygplates.RotationModel from EITHER a named plate model
    (downloaded via gplately.download.DataServer) OR a list of local
    rotation file paths.  Exactly one of `model_name` and `rotation_files`
    must be provided (the other should be None).

    Returning the RotationModel here is fine for in-process inspection.
    For parallel rotation (joblib loky workers) prefer to call this helper
    INSIDE the worker rather than passing the model across process
    boundaries.
    """
    has_name  = bool(model_name)
    has_files = bool(rotation_files)
    if has_name and has_files:
        raise ValueError("Provide either model_name OR rotation_files, not both.")
    if not has_name and not has_files:
        raise ValueError("Provide either model_name or rotation_files (one must be set).")

    if has_name:
        ds = gplately.download.DataServer(model_name)
        rotation_model, _, _ = ds.get_plate_reconstruction_files()
        return rotation_model

    # rotation_files path — accept a single string or any iterable of paths.
    if isinstance(rotation_files, (str, os.PathLike)):
        rotation_files = [rotation_files]
    rotation_files = [str(p) for p in rotation_files]
    missing = [p for p in rotation_files if not os.path.exists(p)]
    if missing:
        raise FileNotFoundError(f"Rotation file(s) not found: {missing}")
    return pygplates.RotationModel(rotation_files)


# Build / prime each side once in the main process for cache warm-up + sanity.
_from_model_main = build_rotation_model(from_model_name, from_rotation_files)
_to_model_main   = build_rotation_model(to_model_name,   to_rotation_files)

reconstruction_times = np.arange(min_time, max_time + timestep_size, timestep_size)

# Make sure the output directory exists.
os.makedirs(os.path.dirname(os.path.abspath(output_grid_template.format(0.0))), exist_ok=True)

# Print a one-screen summary of what will run.
def _summary(label, name, files, anchor):
    src = f"model='{name}'" if name else f"files={files}"
    return f"  {label:<6}: {src}, anchor={anchor}"

print("Rotation configuration")
print(_summary("source", from_model_name, from_rotation_files, from_anchor))
print(_summary("target", to_model_name,   to_rotation_files,   to_anchor))
print(f"  times : {min_time}-{max_time} Ma at {timestep_size} Myr cadence "
      f"({len(reconstruction_times)} grids)")
print(f"  res   : {grid_spacing_degrees}°")
print(f"  input : {input_grid_template}")
print(f"  output: {output_grid_template}")

Rotation configuration
  source: model='Alfonso2024', anchor=0
  target: model='Alfonso2024', anchor=701701
  times : 105-110 Ma at 1 Myr cadence (6 grids)
  res   : 0.1°
  input : 15-ConvertGridReferenceFrame-notebook-data/alfonso2024/Rasters/AgeGrids/Alfonso2024_SEAFLOOR_AGE_grid_{:.2f}Ma.nc
  output: 15-ConvertGridReferenceFrame-notebook-data/alfonso2024/Rasters/AgeGrids-PMAG/Alfonso2024_SEAFLOOR_AGE_grid_PMAG_{:.2f}Ma.nc

Rotation function (parallel-safe)

rotate_grid is the worker function that does one age slice. It rebuilds its own copy of the rotation models inside the worker process — that's how parallel rotation stays robust, because the underlying pygplates objects don't always survive being pickled and shipped to another process.

You don't need to call this directly; the cells below do.

def rotate_grid(input_grid_filename, reconstruction_time, output_name):
    """
    Rotate a single NetCDF grid from the source reference frame to the
    target reference frame.

    Parameters
    ----------
    input_grid_filename : str
        Path to the source grid (already formatted with the reconstruction
        time, i.e. no remaining `{}` placeholders).
    reconstruction_time : float
        The age (Ma) the grid represents.
    output_name : str
        Path to write the rotated grid to (already formatted).

    Notes
    -----
    Reads notebook globals: `from_model_name`, `from_rotation_files`,
    `from_anchor`, `to_model_name`, `to_rotation_files`, `to_anchor`,
    `grid_spacing_degrees`.

    Worker-safe: rebuilds rotation models inside the worker process from
    primitive (string / list-of-string / int) inputs.  The main process
    has already primed gplately's local cache (see the validation cell
    above), so worker model construction is a fast cache read.

    Output silencer: per-worker stdout/stderr are redirected and
    `warnings` are suppressed for the duration of this call.  Without
    this, gplately's cache messages and pygplates' NaN-on-land
    RuntimeWarnings accumulate in the parallel cell's notebook output
    across 171 ages, inflating the notebook payload past JupyterLab's
    autosave limit ('Failed to fetch').  Errors that raise exceptions
    are NOT swallowed — they propagate up so the parallel cell still
    surfaces real failures.
    """
    with contextlib.redirect_stdout(io.StringIO()), \
         contextlib.redirect_stderr(io.StringIO()), \
         warnings.catch_warnings():
        warnings.simplefilter("ignore")

        src_model = build_rotation_model(from_model_name, from_rotation_files)
        tgt_model = build_rotation_model(to_model_name,   to_rotation_files)

        rasterobj = gplately.grids.Raster(input_grid_filename, time=reconstruction_time)
        rasterobj.rotate_reference_frames(
            grid_spacing_degrees=grid_spacing_degrees,
            reconstruction_time=reconstruction_time,
            from_rotation_features_or_model=src_model,
            to_rotation_features_or_model=tgt_model,
            from_rotation_reference_plate=from_anchor,
            to_rotation_reference_plate=to_anchor,
            output_name=output_name,
        )

Rotate the grids — parallel (recommended)

This is the fast path. By default it uses every CPU minus one (n_jobs=-2). If anything goes wrong, fall back to the single-core cell below to see the error message clearly without the parallel machinery in the way.

from joblib import Parallel, delayed

# set use_parallel to True if you want to run this code cell
if use_parallel:
    converted_frames = Parallel(n_jobs=-2, backend='loky', verbose=1)(
        delayed(rotate_grid)(
            input_grid_filename = input_grid_template.format(reconstruction_time),
            reconstruction_time = reconstruction_time,
            output_name         = output_grid_template.format(reconstruction_time),
        )
        for reconstruction_time in reconstruction_times
    )

[Parallel(n_jobs=-2)]: Using backend LokyBackend with 3 concurrent workers.

[Parallel(n_jobs=-2)]: Done   6 out of   6 | elapsed:  1.4min finished

Rotate the grids — single core (fallback / debugging)

Equivalent to the parallel cell above, just slower and easier to debug. Run this instead of the parallel cell, not as well as it.

# set use_parallel to False if you want to run this code cell
if not use_parallel:
    for reconstruction_time in reconstruction_times:
        rotate_grid(
            input_grid_filename = input_grid_template.format(reconstruction_time),
            reconstruction_time = reconstruction_time,
            output_name         = output_grid_template.format(reconstruction_time),
        )
        print(f"Finished rotating {reconstruction_time} Ma")

Sample comparison plot

A side-by-side look at the source and target grids at one age, to confirm the rotation did something sensible. Change compare_time to inspect a different slice. Continents should be in different positions in the two panels — that's the rotation working.

compare_time = max_time

proj = ccrs.Mollweide(central_longitude=60)
fig, (ax1, ax2) = plt.subplots(
    1, 2, subplot_kw={'projection': proj}, figsize=(12, 6), dpi=150,
)
input_max = np.ma.max(gplately.Raster(input_grid_template.format(compare_time)).data)
input_min = np.ma.min(gplately.Raster(input_grid_template.format(compare_time)).data)
ax1.imshow(
    np.flipud(gplately.Raster(input_grid_template.format(compare_time)).data),
    transform=ccrs.PlateCarree(),
)
ax1.set_title(f"From (anchor={from_anchor}) @ {compare_time} Ma")

ax2.imshow(
    np.flipud(gplately.Raster(output_grid_template.format(compare_time)).data),
    transform=ccrs.PlateCarree(), vmin=input_min, vmax=input_max
)
ax2.set_title(f"To (anchor={to_anchor}) @ {compare_time} Ma")
plt.show()