Gplately Regional Plots

Lauren Ilano, Hojat Shirmard, Dietmar Muller
EarthByte Group, School of Geosciences, University of Sydney, NSW 2006, Australia

This notebook generates regional plots given a certain:

• model_dir - The absolute path to your plate reconstruction model directory
• grid_filename The absolute path to a grid, i.e. crustal, sediment thickness, seafloor age, spreading rate etc. If the grid(s) are unique per timestep, the time in the filename must be replaced by curly brackets {:.0f} with a 0 if the age is an integer, and 1 if the age is to 1 decimal place, and so on.

Instructions

1. Ensure all dependencies are installed to your Python environment (the packages listed below)
   • Note: A dependency is cmcrameri from GitHub. Install by cloning the repository, changing the directory on your command line to the repository's top level, and typing pip install . into your terminal.
2. Change model_dir to your plate model directory, and grid_filename to your required grids for plotting, ensuring to generalise the timesteps as mentioned above.

Data

The datasets used in this ntoebook can be found at https://zenodo.org/records/13777155. Download the source_data.zip and unzip it in the same folder as this notebook file.

import gplately

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np
import pandas as pd
import os
import gplately.tools as tools
import matplotlib.colors as mcolors
from cartopy.io import shapereader as shpreader
import netCDF4
import warnings
from scipy import ndimage
import glob
from rasterio.features import rasterize
from rasterio.transform import from_bounds
import matplotlib.gridspec as gridspec
import matplotlib
import pandas as pd

# Need joblib to run the script on multiple cores, and moviepy to make the movies
from joblib import Parallel, delayed
import joblib
import moviepy.editor as mpy

# Plotting dependencies
from cmcrameri import cm


import cartopy.mpl.ticker as cticker
import matplotlib.ticker as mticker

recon_model = gplately.PlateModelManager().get_model(
    "Alfonso2024",  # model name
    data_dir="plate-model-repo",  # the folder to save the model files
)

model = gplately.PlateReconstruction(
    recon_model.get_rotation_model(),
    topology_features=recon_model.get_layer("Topologies"),
    static_polygons=recon_model.get_layer("StaticPolygons"),
)

gplot = gplately.PlotTopologies(
    model,
    coastlines=recon_model.get_layer("Coastlines"),
    time=170,
)

# Path to source data zip containing sed thickness grids
#source_data_dir = os.path.join('/source_data', 'source_data')
total_sed_grid_filename ="./source_data/SedimentThickness/sed_thick_0.1d_{}.nc"
carbonate_sed_grid_filename = "./source_data/CarbonateThickness/uncompacted_carbonate_thickness_{}Ma.nc"

# Path to save all output plots to
output_dir = "./outputs"
os.makedirs(output_dir+"/Plots_v2", exist_ok=True)

Load deposits

Point to where your deposits (xlsx) file is, and this will generate a pandas dataframe with the deposit coordinates.

# Load the Excel file
file_path = "./inputs/Mutschler_WorldPorphyryCopperDeposits_GPlates.xlsx"
df = pd.read_excel(file_path, sheet_name="Deposits")

# Clean column names
df.columns = df.columns.str.strip().str.replace(' ', '_').str.replace('–', '-')

# Rename odd characters if present
df.rename(columns={
    'Au_grade_gPERt': 'Au_grade',
    'Ag_grade_gPERt': 'Ag_grade',
    'w’W_mt': 'W_mt'
}, inplace=True)

# Keep only relevant columns
columns_to_keep = [
    'Camp', 'Lat_Dec', 'Long_Dec', 'AGE', 'Type',
    'Ore_mt', 'Cu_mt', 'Mo_mt', 'Au_kg', 'Ag_kg',
    'Cu_grade', 'Mo_grade', 'Au_grade', 'Ag_grade'
]
df = df[columns_to_keep]

# Filter to Porphyry deposits (case-insensitive)
df = df[df['Type'].str.contains("Porphyry", case=False, na=False)]

# Rename coordinate columns
df.rename(columns={'Lat_Dec': 'LAT', 'Long_Dec': 'LON'}, inplace=True)

# Convert AGE and Cu_mt to numeric
df['Cu_mt'] = pd.to_numeric(df['Cu_mt'], errors='coerce')
df['AGE'] = pd.to_numeric(df['AGE'], errors='coerce')

# Filter out invalid or missing data
df = df.dropna(subset=['Cu_mt', 'AGE'])
df = df[df['Cu_mt'] > 0]

# Filter deposits with AGE <= 170
df = df[df['AGE'] <= 171]

# Define custom bins and labels
custom_bins = [0, 2e6, 10e6, 30e6, np.inf]
custom_labels = [
    "Minor (<2 Mt)",
    "Moderate (2–10 Mt)",
    "Major (10–30 Mt)",
    "Giant (>30 Mt)"
]

# Categorize SIZE
df['SIZE'] = pd.cut(df['Cu_mt'], bins=custom_bins, labels=custom_labels, include_lowest=True)

# Show distribution
print(df['SIZE'].value_counts())

SIZE
Minor (<2 Mt)         80
Moderate (2–10 Mt)    64
Major (10–30 Mt)      13
Giant (>30 Mt)         7
Name: count, dtype: int64

def plot_deposits(time, extents, proj, grid_type, save_fig=None):
    """Plots categorized copper deposits on a reconstructed map (Cu grade by color, size by symbol)."""

    gplot.time = time
    
    fig = plt.figure(figsize=(14, 10))
    ax1 = fig.add_subplot(111, projection=proj)
    cmap = cm.batlow_r

    # Load grid
    if grid_type.lower() == "total":
        grid_filename = total_sed_grid_filename
        vmin, vmax = 50, 400
        grid_label = "Sediment thickness (m)"
    elif grid_type.lower() == "carbonate":
        grid_filename = carbonate_sed_grid_filename
        vmin, vmax = 0, 500
        grid_label = "Carbonate sediment thickness (m)"
    else:
        raise ValueError("This grid type is not supported.")
        
    im = gplot.plot_grid_from_netCDF(
        ax1, grid_filename.format(int(gplot.time)), 
        cmap=cmap, vmin=vmin, vmax=vmax
    )

    # Plot tectonic features
    gplot.plot_coastlines(ax1, facecolor='silver', edgecolor='None')
    gplot.plot_all_topological_sections(ax1, color='grey', tessellate_degrees=1)
    gplot.plot_ridges(ax1, color='darkred', linewidth=2, tessellate_degrees=1)
    gplot.plot_trenches(ax1, color='white', linewidth=5, tessellate_degrees=1)
    gplot.plot_trenches(ax1, color='k', tessellate_degrees=1)
    gplot.plot_subduction_teeth(ax1, color='k', spacing=0.02, zorder=10)
    gplot.plot_plate_motion_vectors(ax1, spacingX=10, spacingY=10, normalise=True, zorder=10, alpha=0.5)

    ax1.set_extent(extents, ccrs.PlateCarree())

    # Cu grade bins (in %), reversed color order
    cu_bins = [0, 0.3, 0.6, 1.0, float("inf")]
    cu_labels = ["Low (<0.4%)", "Moderate (0.4–0.6%)", "High (0.6–0.8%)", "Very High (>0.8%)"]
    cu_colors = ["blue", "green", "orange", "red"]  # reversed order

    # New Size categories (renamed labels)
    size_styles = {
        "Small (<2 Mt)": {"size": 70},
        "Medium (2–10 Mt)": {"size": 100},
        "Large (10–30 Mt)": {"size": 130},
        "Giant (>30 Mt)": {"size": 170},
    }

    # Filter and classify deposits
    cu_df = df[(df['AGE'] != 'Unknown') & (df['AGE'] >= time) & df['Cu_grade'].notna()].copy()
    cu_df['Cu_grade_%'] = cu_df['Cu_grade'] * 100  # Convert to percent
    cu_df['Cu_bin'] = pd.cut(cu_df['Cu_grade_%'], bins=cu_bins, labels=cu_labels, include_lowest=True)

    # Replace size label column if needed (optional, based on your data)
    size_map = {
        "Minor (<2 Mt)": "Small (<2 Mt)",
        "Moderate (2–10 Mt)": "Medium (2–10 Mt)",
        "Major (10–30 Mt)": "Large (10–30 Mt)",
        "Giant (>30 Mt)": "Giant (>30 Mt)"
    }
    cu_df["SIZE"] = cu_df["SIZE"].replace(size_map)

    # Plot deposits: by Cu bin (color) and SIZE (size)
    for cu_label, cu_color in zip(cu_labels, cu_colors):
        for size_label, style in size_styles.items():
            subset = cu_df[(cu_df['Cu_bin'] == cu_label) & (cu_df['SIZE'] == size_label)]
            if subset.empty:
                continue
            points = gplately.Points(model, subset['LON'].to_numpy(), subset['LAT'].to_numpy())
            lons, lats = points.reconstruct(time, return_array=True)
            ax1.scatter(
                lons, lats, marker="o", color=cu_color, edgecolors='black',
                linewidths=0.2, alpha=0.7, s=style['size'],
                transform=ccrs.PlateCarree()
            )

    # Legends
    from matplotlib.lines import Line2D

    # Cu Grade Legend (color)
    cu_legend_elements = [
        Line2D([0], [0], marker='o', color=color, label=label,
               markersize=10, linestyle='', markeredgecolor='black', alpha=0.7)
        for label, color in zip(cu_labels, cu_colors)
    ]
    # legend1 = ax1.legend(handles=cu_legend_elements, title="Cu Grade (%)",
    #                      loc='upper left', bbox_to_anchor=(0, 1), fontsize=10, title_fontsize=11)

    # Deposit Size Legend (size)
    size_legend_elements = [
        Line2D([0], [0], marker='o', color='gray', label=label,
               markersize=style['size'] * 0.05, linestyle='',
               alpha=0.7, markeredgecolor='black')
        for label, style in size_styles.items()
    ]

    # legend2 = ax1.legend(handles=size_legend_elements, title="Deposit Size (Mt)",
    #                      loc='lower left', bbox_to_anchor=(0, 0), fontsize=10, title_fontsize=11)


    legend1 = ax1.legend(
        handles=cu_legend_elements,
        title="Cu Grade (%)",
        loc='upper left',
        bbox_to_anchor=(0, 1),
        fontsize=10,
        title_fontsize=11,
        framealpha=1,  # fully opaque
        facecolor='white',  # background color
        edgecolor='black'   # optional: border color
        )

    legend2 = ax1.legend(
        handles=size_legend_elements,
        title="Deposit Size (Mt)",
        loc='lower left',
        bbox_to_anchor=(0, 0),
        fontsize=10,
        title_fontsize=11,
        framealpha=1,        # fully opaque
        facecolor='white',   # background color
        edgecolor='black'    # optional: border color
        )
    legend1.set_zorder(100)
    legend2.set_zorder(100)

    # Add both legends to the plot
    ax1.add_artist(legend1)
    ax1.add_artist(legend2)

    # Final touches
    #latlonticks(ax1)
    gl = ax1.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.7, linestyle='--')
    gl.top_labels = False        # don't show top (longitude)
    gl.right_labels = False      # don't show right (latitude)
    gl.bottom_labels = True      # ✅ show bottom (longitude)
    gl.left_labels = True        # ✅ show left (latitude)
    gl.xlabel_style = {"size": 10}
    gl.ylabel_style = {"size": 10}

    plt.title(f"{int(gplot.time)} Ma", fontsize=20)
    cbar = plt.colorbar(im, shrink=0.5)
    cbar.set_label(grid_label, fontsize=14)

    if save_fig:
        fig.savefig(save_fig, dpi=300, bbox_inches='tight')
    else:
        plt.show()
    plt.close()

North America

plot_deposits(0, extents=(-190, -30, -10, 65), proj=ccrs.PlateCarree(np.mean((-190, -30))), grid_type="carbonate")

2025-06-23 20:33:52 - py.warnings - WARNING - /var/folders/gz/437jwwyj4ld14qvq5715jv6w0000gr/T/ipykernel_13422/4240219516.py:63: FutureWarning: The behavior of Series.replace (and DataFrame.replace) with CategoricalDtype is deprecated. In a future version, replace will only be used for cases that preserve the categories. To change the categories, use ser.cat.rename_categories instead.

plot_deposits(0, extents=(-190, -30, -10, 65), proj=ccrs.PlateCarree(np.mean((-190, -30))), grid_type="total")

2025-06-23 20:33:55 - py.warnings - WARNING - /var/folders/gz/437jwwyj4ld14qvq5715jv6w0000gr/T/ipykernel_13422/4240219516.py:63: FutureWarning: The behavior of Series.replace (and DataFrame.replace) with CategoricalDtype is deprecated. In a future version, replace will only be used for cases that preserve the categories. To change the categories, use ser.cat.rename_categories instead.

Generating Plots (png)

This code generates and saves geological deposit plots for different times and grid types, either sequentially or in parallel depending on the use_parallel flag.

reconstruction_times = np.arange(170,-1,-1)

use_parallel=False

grid_types = ["carbonate", "total"]
if use_parallel:
    for grid_type in grid_types:
        # Use LokyBackend to protect the netCDF routine
        slabdip_img = Parallel(n_jobs=-1, backend='loky', verbose=1) \
        (delayed(plot_deposits) \
         (reconstruction_time, 
          extents=(-190, -30, -10, 65), 
          proj=ccrs.PlateCarree(np.mean((-190, -30))),
          grid_type=grid_type,
          save_fig=output_dir+"/Plots_v2/NorthAmerica_{}_{}.{}".format(grid_type, reconstruction_time, "png"),
         ) for reconstruction_time in reconstruction_times)
else:
    for grid_type in grid_types:
        for reconstruction_time in reconstruction_times:
            plot_deposits(
                reconstruction_time, 
                extents=(-190, -30, -10, 65), 
                proj=ccrs.PlateCarree(np.mean((-120, -30))),
                grid_type=grid_type,
                save_fig=output_dir+"/Plots_v2/NorthAmerica_{}_{}.{}".format(grid_type, reconstruction_time, "png"),
            )
            print("Time {} done".format(reconstruction_time))

Generating Movies (mp4) from Plots (png)

for grid_type in grid_types:
    frame_list = []
    for time in np.arange(170,-1,-1):
        frame_list.append(
            output_dir+"/Plots_v2/NorthAmerica_{}_{}.png".format(grid_type, time), 

        )

    clip = mpy.ImageSequenceClip(frame_list, fps=25)

    clip.write_videofile(
    output_dir+"/Plots_v2/NorthAmerica_{}.mp4".format(grid_type),
    fps=100,
    codec="libx264",
    bitrate="5000k",
    audio=False,
    logger=None,
    ffmpeg_params=[
        "-vf",
        "pad=ceil(iw/2)*2:ceil(ih/2)*2",
        "-pix_fmt",
        "yuv420p",
    ],
)

East Asia

plot_deposits(0, extents=(110, 170, -20, 40), proj=ccrs.PlateCarree(np.mean((120, 40))), grid_type="carbonate")

2025-06-23 20:58:35 - py.warnings - WARNING - /var/folders/gz/437jwwyj4ld14qvq5715jv6w0000gr/T/ipykernel_13422/4240219516.py:63: FutureWarning: The behavior of Series.replace (and DataFrame.replace) with CategoricalDtype is deprecated. In a future version, replace will only be used for cases that preserve the categories. To change the categories, use ser.cat.rename_categories instead.

plot_deposits(0, extents=(110, 170, -20, 40), proj=ccrs.PlateCarree(np.mean((120, 40))), grid_type="total")

2025-06-23 20:58:38 - py.warnings - WARNING - /var/folders/gz/437jwwyj4ld14qvq5715jv6w0000gr/T/ipykernel_13422/4240219516.py:63: FutureWarning: The behavior of Series.replace (and DataFrame.replace) with CategoricalDtype is deprecated. In a future version, replace will only be used for cases that preserve the categories. To change the categories, use ser.cat.rename_categories instead.

Generating Plots (png)

This code generates and saves geological deposit plots for different times and grid types, either sequentially or in parallel depending on the use_parallel flag.

reconstruction_times = np.arange(170,-1,-1)

use_parallel=False

grid_types = ["carbonate", "total"]
if use_parallel:
    for grid_type in grid_types:
        # Use LokyBackend to protect the netCDF routine
        slabdip_img = Parallel(n_jobs=-1, backend='loky', verbose=1) \
        (delayed(plot_deposits) \
         (reconstruction_time, 
          extents=(110, 170, -20, 40), 
          proj=ccrs.PlateCarree(np.mean((120, 40))),
          grid_type=grid_type,
          save_fig=output_dir+"/Plots_v2/EastAsia_{}_{}.{}".format(grid_type, reconstruction_time, "png"),
         ) for reconstruction_time in reconstruction_times)
else:
    for grid_type in grid_types:
        for reconstruction_time in reconstruction_times:
            plot_deposits(
                reconstruction_time, 
                extents=(110, 170, -20, 40), 
                proj=ccrs.PlateCarree(np.mean((120, 40))),
                grid_type=grid_type,
                save_fig=output_dir+"/Plots_v2/EastAsia_{}_{}.{}".format(grid_type, reconstruction_time, "png"),
            )
            print("Time {} done".format(reconstruction_time))

Generating Movies (mp4) from Plots (png)

for grid_type in grid_types:
    frame_list = []
    for time in np.arange(170,-1,-1):
        frame_list.append(
            output_dir+"/Plots_v2/EastAsia_{}_{}.png".format(grid_type, time), 

        )

    clip = mpy.ImageSequenceClip(frame_list, fps=25)

    clip.write_videofile(
    output_dir+"/Plots_v2/EastAsia_{}.mp4".format(grid_type),
    fps=100,
    codec="libx264",
    bitrate="5000k",
    audio=False,
    logger=None,
    ffmpeg_params=[
        "-vf",
        "pad=ceil(iw/2)*2:ceil(ih/2)*2",
        "-pix_fmt",
        "yuv420p",
    ],
)

Generate plots and movies for other regions by updating the extents values in the Generating Plots (png) and Generating Movies (mp4) from Plots (png) sections—for example:

• South America: extents = (-120, -30, -60, 5)
• Mediterranean Sea: extents = (0, 70, 10, 60)