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 using GPlately.
Data¶
📌 IMPORTANT!!!
The notebook needs three datasets in the folder WorkflowData/12-GPlately-Regional-Plots/source_data/:
- Mutschler_WorldPorphyryCopperDeposits_GPlates.xlsx
- CarbonateThickness grids (0-170 Ma)
- SedimentThickness grids (0-170 Ma)
You can download the datasets at https://repo.gplates.org/webdav/gplately/12-GPlately-Regional-Plots-data.zip (2.3G).
Alternatively, the datasets are also available at https://zenodo.org/records/13777155.
This notebook requires several dependencies. The easiest way is to create a Conda environment from this YAML file. Otherwise, make sure all required dependencies are installed before running the notebook.
In [1]:
import gplately
import warnings, os
from pathlib import Path
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import pandas as pd
from IPython.display import clear_output
# Need joblib to run the script on multiple cores, and moviepy to make the movies
from joblib import Parallel, delayed
try:
from moviepy.editor import ImageSequenceClip # moviepy 1.x
except ImportError:
from moviepy import ImageSequenceClip # moviepy 2.x
# Plotting dependencies
from cmcrameri import cm
warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=FutureWarning)
data_dir = Path("WorkflowData/12-GPlately-Regional-Plots")
output_dir = data_dir / "output"
output_dir.mkdir(parents=True, exist_ok=True)
source_data_dir = "12-GPlately-Regional-Plots-Source-Data"
# Set quick_run to False for a full run, or True for a quick run
quick_run = False
#quick_run = True
if not quick_run:
time_steps = np.arange(170, -1, -1)
else:
# This notebook takes a long time to run. For a quick run, we can plot every 34 million years (170, 136, 102, 68, 34, 0)
time_steps = np.arange(170, -1, -34)
Matplotlib is building the font cache; this may take a moment.
Check whether the required datasets are available in the expected directory before continuing.
In [2]:
deposits_file = f"{source_data_dir}/Mutschler_WorldPorphyryCopperDeposits_GPlates.xlsx"
# Template paths to the sediment thickness and carbonate sediment thickness grid
# files, with the reconstruction time filled in later via .format()
total_sed_grid_filename = f"{source_data_dir}/SedimentThickness/sed_thick_0.1d_{{}}.nc"
carbonate_sed_grid_filename = f"{source_data_dir}/CarbonateThickness/uncompacted_carbonate_thickness_{{}}Ma.nc"
try:
if not os.path.exists(deposits_file):
raise FileNotFoundError(f"Deposits file not found: {deposits_file}")
for time in time_steps:
if not os.path.exists(total_sed_grid_filename.format(time)):
raise FileNotFoundError(f"Total sediment grid file not found: {total_sed_grid_filename.format(time)}")
if not os.path.exists(carbonate_sed_grid_filename.format(time)):
raise FileNotFoundError(f"Carbonate sediment grid file not found: {carbonate_sed_grid_filename.format(time)}")
except FileNotFoundError as e:
print(e)
print("Please ensure that the required data files are present in the source_data directory.")
print("You can download the data at https://repo.gplates.org/webdav/gplately/12-GPlately-Regional-Plots-data.zip (2.3G).")
raise
In [3]:
pmm_model = gplately.PlateModelManager().get_model(
"Alfonso2024", # model name
data_dir="plate-model-repo", # the folder to save the model files
)
assert pmm_model is not None, "Failed to load the plate model"
model = gplately.PlateReconstruction(
pmm_model.get_rotation_model(),
topology_features=pmm_model.get_layer("Topologies"),
static_polygons=pmm_model.get_layer("StaticPolygons"),
)
gplot = gplately.PlotTopologies(
model,
coastlines=pmm_model.get_layer("Coastlines"),
time=170,
)
Load deposits¶
Point to the location of your deposits (.xlsx) file — running this cell will generate a pandas DataFrame containing the deposit coordinates.
In [4]:
# Load the Excel file
df = pd.read_excel(deposits_file, 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 <= 171 Ma (keeps ages up to and just beyond 170 Ma)
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
In [5]:
def plot_deposits(time, extents, proj, grid_type, save_fig=None):
"""Plot categorized copper deposits on a reconstructed map, with Cu grade shown by colour and deposit size shown by symbol size."""
warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=FutureWarning)
gplot.time = time
fig = plt.figure(figsize=(14, 10))
ax1 = fig.add_subplot(111, projection=proj)
cmap = cm.batlow_r # type: ignore
# 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(f"Grid type '{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()) # type: ignore
# Cu grade bins (in %), reversed color order
cu_bins = [0, 0.3, 0.6, 1.0, float("inf")]
cu_labels = [
"Low (<0.3%)",
"Moderate (0.3-0.6%)",
"High (0.6-1.0%)",
"Very High (>1.0%)",
]
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"].astype("string").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) # type: ignore
ax1.scatter(
lons, # type: ignore
lats, # type: ignore
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( # type: ignore
draw_labels=True, linewidth=0.5, color="gray", alpha=0.7, linestyle="--"
)
gl.top_labels = False # don't show top labels (longitude)
gl.right_labels = False # don't show right labels (latitude)
gl.bottom_labels = True # show bottom labels (longitude)
gl.left_labels = True # show left labels (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¶
In [6]:
plot_deposits(
0,
extents=(-190, -30, -10, 65),
proj=ccrs.PlateCarree(float(np.mean((-190, -30)))),
grid_type="carbonate",
)
In [7]:
plot_deposits(
0,
extents=(-190, -30, -10, 65),
proj=ccrs.PlateCarree(float(np.mean((-190, -30)))),
grid_type="total",
)
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.
In [8]:
use_parallel = True
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(float(np.mean((-190, -30)))),
grid_type=grid_type,
save_fig=f"{output_dir}/NorthAmerica_{grid_type}_{reconstruction_time}.png"
)
for reconstruction_time in time_steps
)
else:
for grid_type in grid_types:
for reconstruction_time in time_steps:
plot_deposits(
reconstruction_time,
extents=(-190, -30, -10, 65),
proj=ccrs.PlateCarree(float(np.mean((-190, -30)))),
grid_type=grid_type,
save_fig=f"{output_dir}/NorthAmerica_{grid_type}_{reconstruction_time}.png"
)
print(f"Finished plotting {reconstruction_time} Ma")
clear_output()
Generating Movies (mp4) from Plots (png)¶
In [9]:
for grid_type in grid_types:
frame_list = []
for time in time_steps:
frame_list.append(
f"{output_dir}/NorthAmerica_{grid_type}_{time}.png"
)
clip = ImageSequenceClip(frame_list, fps=25)
clip.write_videofile(
f"{output_dir}/NorthAmerica_{grid_type}.mp4",
fps=25,
codec="libx264",
bitrate="5000k",
audio=False,
logger=None,
ffmpeg_params=[
"-vf",
"pad=ceil(iw/2)*2:ceil(ih/2)*2",
"-pix_fmt",
"yuv420p",
],
)
print(f"The video file has been saved to {output_dir}/NorthAmerica_{grid_type}.mp4")
The video file has been saved to WorkflowData/12-GPlately-Regional-Plots/output/NorthAmerica_carbonate.mp4
The video file has been saved to WorkflowData/12-GPlately-Regional-Plots/output/NorthAmerica_total.mp4
East Asia¶
In [10]:
plot_deposits(
0,
extents=(110, 170, -20, 40),
proj=ccrs.PlateCarree(float(np.mean((120, 40)))),
grid_type="carbonate",
)
In [11]:
plot_deposits(
0,
extents=(110, 170, -20, 40),
proj=ccrs.PlateCarree(float(np.mean((120, 40)))),
grid_type="total",
)
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.
In [12]:
use_parallel = True
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(float(np.mean((120, 40)))),
grid_type=grid_type,
save_fig=f"{output_dir}/EastAsia_{grid_type}_{reconstruction_time}.png",
)
for reconstruction_time in time_steps
)
else:
for grid_type in grid_types:
for reconstruction_time in time_steps:
plot_deposits(
reconstruction_time,
extents=(110, 170, -20, 40),
proj=ccrs.PlateCarree(float(np.mean((120, 40)))),
grid_type=grid_type,
save_fig=f"{output_dir}/EastAsia_{grid_type}_{reconstruction_time}.png",
)
print(f"Finished plotting {reconstruction_time} Ma")
clear_output()
Generating Movies (mp4) from Plots (png)¶
In [13]:
for grid_type in grid_types:
frame_list = []
for time in time_steps:
frame_list.append(
f"{output_dir}/EastAsia_{grid_type}_{time}.png"
)
clip = ImageSequenceClip(frame_list, fps=25)
clip.write_videofile(
f"{output_dir}/EastAsia_{grid_type}.mp4",
fps=25,
codec="libx264",
bitrate="5000k",
audio=False,
logger=None,
ffmpeg_params=[
"-vf",
"pad=ceil(iw/2)*2:ceil(ih/2)*2",
"-pix_fmt",
"yuv420p",
],
)
print(f"The video file has been saved to {output_dir}/EastAsia_{grid_type}.mp4")
The video file has been saved to WorkflowData/12-GPlately-Regional-Plots/output/EastAsia_carbonate.mp4
The video file has been saved to WorkflowData/12-GPlately-Regional-Plots/output/EastAsia_total.mp4
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)