3 - Working With Points¶

In this notebook, we use GPlately to manipulate and reconstruct point data.

gpts = gplately.Points(
    plate_reconstruction, # plate reconstruction model
    lons, # list or numpy array of longitudinal coordinates
    lats, # list or numpy array of latitudinal coordinates
    time=0, # time is set to the present day by default
    plate_id=None, # optionally pass an array (or single integer) of pre-determined plate IDs
    age=numpy.inf, # optionally pass an array (or single float) of pre-determined appearance ages (defaults to: appearing for all time)
)

In this example, we will reconstruct data from the Paleobiology Database (PBDB). This data is in CSV format.

In [1]:
import os, warnings
from pathlib import Path
import pandas as pd
import numpy as np
import gplately
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from plate_model_manager import PlateModelManager

data_dir = Path("WorkflowData/03-Working-With-Points")
output_dir = data_dir / "output"
output_dir.mkdir(parents=True, exist_ok=True)

We first download the Zahirovic2022 plate reconstruction model files to use in this notebook, and set up the PlateReconstruction and PlotTopologies objects (calling them model and gplot) from these model files.

In [2]:
model_name = "Zahirovic2022" # You can change this to another model name available in PlateModelManager
# print(PlateModelManager().get_available_model_names()) # uncomment this line to see all available model names
pmm_model = PlateModelManager().get_model(model_name, data_dir="plate-model-repo")
assert pmm_model is not None, f"Model {model_name} not found."

rotation_model = pmm_model.get_rotation_model()
topology_features = pmm_model.get_topologies()
static_polygons = pmm_model.get_static_polygons()

coastlines = pmm_model.get_layer("Coastlines")
continents = pmm_model.get_layer("ContinentalPolygons")
COBs = pmm_model.get_layer("COBs")

model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)
gplot = gplately.PlotTopologies(
    model, coastlines=coastlines, continents=continents, COBs=COBs
)

Download and import PBDB data¶

We can import data from the PBDB directly into pandas using the data URL. Alternatively, we can download the CSV file from their website and import that.

When importing CSV files, it is often easier if the first row is the column names, although pandas does allow you to skip header rows if needed. Conveniently, the PBDB provides an option when downloading data to exclude the metadata at the beginning of the file.

In [3]:
# First, check if the data file already exists locally. If it does, we can read it directly from the file. If not, we will download it from the PBDB website.
data_file_path = "NotebookFiles/pbdb-data.csv"

if os.path.exists(data_file_path):
    pbdb_data = pd.read_csv(data_file_path)
    print(f"Data loaded from {data_file_path}.")
else:
    data_file_path = data_dir / "pbdb-data.csv"
    if os.path.exists(data_file_path):
        pbdb_data = pd.read_csv(data_file_path)
        print(f"Data loaded from {data_file_path}.")
    else:
        # Download data for the Jurassic period, and include the paleoenvironment column.
        # You can use the download page to play with the options and get the download link and/or CSV.
        pbdb_data_url = "https://paleobiodb.org/data1.2/occs/list.csv?datainfo&rowcount&base_name=Foraminifera&interval=Jurassic&show=coords,env"
        pbdb_data = pd.read_csv(pbdb_data_url, sep=",", skiprows=18)
        pbdb_data.to_csv(data_file_path, index=False)
        print(f"Data downloaded from {pbdb_data_url}.")

Data downloaded from https://paleobiodb.org/data1.2/occs/list.csv?datainfo&rowcount&base_name=Foraminifera&interval=Jurassic&show=coords,env.
In [4]:
print(pbdb_data.columns)
print(pbdb_data.shape)

Index(['occurrence_no', 'record_type', 'reid_no', 'flags', 'collection_no',
       'identified_name', 'identified_rank', 'identified_no', 'difference',
       'accepted_name', 'accepted_rank', 'accepted_no', 'early_interval',
       'late_interval', 'max_ma', 'min_ma', 'reference_no', 'lng', 'lat',
       'environment'],
      dtype='str')
(5273, 20)
In [5]:
# Set up a GeoAxis plot
fig = plt.figure(figsize=(16, 8), dpi=300)
ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=0))
ax.set_global()  # type: ignore
ax.gridlines(  # type: ignore
    color="0.7",
    linestyle="--",
    xlocs=np.arange(-180, 180, 15),
    ylocs=np.arange(-90, 90, 15),
)
ax.set_title("Present Day Distribution of Jurassic Foraminifera")

warnings.filterwarnings("ignore", category=UserWarning)
# Plot shapefile features, subduction zones and MOR boundaries at 0 Ma
gplot.time = 0  # Ma
gplot.plot_continent_ocean_boundaries(ax, color="b", alpha=0.05)
gplot.plot_continents(ax, facecolor="palegoldenrod", alpha=0.2)
gplot.plot_coastlines(ax, color="DarkGrey")
gplot.plot_all_topological_sections(
    ax,
    plot_subduction_teeth=True,
    other_kwargs={"color": "grey", "linewidth": 0.8},
    ridge_kwargs={"color": "red", "linewidth": 1.0},
    transform_kwargs={"color": "green", "linewidth": 1.0},
    trench_kwargs={"color": "blue", "linewidth": 1.0},
)

sc = ax.scatter(
    pbdb_data["lng"],
    pbdb_data["lat"],
    color="orange",
    transform=ccrs.PlateCarree(),
    label="Jurassic Foraminifera",
)
ax.legend(frameon=False)
plt.show()
No description has been provided for this image

Reconstruct PBDB data with GPlately¶

We use the lon/lat coordinates of the PBDB data, keeping only occurrences whose age range (min_ma to max_ma) contains a chosen reconstruction_time.

The Points object needs a PlateReconstruction object as a parameter — we already created one above (model), built from a rotation_model, topology_features, and static_polygons.

In [6]:
reconstruction_time = 185 # Ma
# Filter the data to only include occurrences whose age range (min_ma to max_ma) contains reconstruction_time
filtered_data = pbdb_data[
    (pbdb_data["min_ma"] <= reconstruction_time) & (reconstruction_time <= pbdb_data["max_ma"])
]
print(filtered_data["min_ma"].max(), filtered_data["max_ma"].min())
print(reconstruction_time,pbdb_data.shape, filtered_data.shape)
gpts = gplately.Points(model, filtered_data["lng"], filtered_data["lat"])
rlons, rlats = gpts.reconstruct(reconstruction_time, return_array=True) # type: ignore

184.2 192.9
185 (5273, 20) (549, 20)
In [7]:
# Set up a GeoAxis plot
fig = plt.figure(figsize=(16, 8), dpi=300)
ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=0))
ax.set_global()  # type: ignore
ax.gridlines(  # type: ignore
    color="0.7",
    linestyle="--",
    xlocs=np.arange(-180, 180, 15),
    ylocs=np.arange(-90, 90, 15),
)

# Plot shapefile features, subduction zones and MOR boundaries at 0 Ma
gplot.time = reconstruction_time  # Ma
gplot.plot_continent_ocean_boundaries(ax, color="b", alpha=0.05)
gplot.plot_continents(ax, facecolor="palegoldenrod", alpha=0.2)
gplot.plot_coastlines(ax, color="DarkGrey")
gplot.plot_all_topological_sections(
    ax,
    plot_subduction_teeth=True,
    other_kwargs={"color": "grey", "linewidth": 0.8},
    ridge_kwargs={"color": "red", "linewidth": 1.0},
    transform_kwargs={"color": "green", "linewidth": 1.0},
    trench_kwargs={"color": "blue", "linewidth": 1.0},
)

sc = ax.scatter(
    rlons,
    rlats,
    color="orange",
    transform=ccrs.PlateCarree(),
    label="Jurassic Foraminifera",
)
ax.legend(frameon=False)
ax.set_title(f"Jurassic Foraminifera Locations at {reconstruction_time} Ma")
plt.show()
No description has been provided for this image

We can make this map look a little nicer by condensing data points that are close to each other. One way is to bin the data on a longitude/latitude grid. Alternatively, we can use an icosahedral mesh, which has relatively uniform point spacing.

In [8]:
from gplately.lib.icosahedron import get_mesh, xyz2lonlat

mesh_vertices, mesh_faces = get_mesh(level=5)

mesh_lons, mesh_lats = xyz2lonlat(
    mesh_vertices[:, 0], mesh_vertices[:, 1], mesh_vertices[:, 2]
)
In [9]:
# Convert reconstructed lon/lat points to unit Cartesian coordinates.
rlon_rad = np.deg2rad(rlons) # type: ignore
rlat_rad = np.deg2rad(rlats) # type: ignore
rx = np.cos(rlat_rad) * np.cos(rlon_rad)
ry = np.cos(rlat_rad) * np.sin(rlon_rad)
rz = np.sin(rlat_rad)
reconstructed_xyz = np.column_stack((rx, ry, rz))

# On the unit sphere, nearest vertex is the one with the largest dot product.
dot_products = reconstructed_xyz @ mesh_vertices.T
indices = np.argmax(dot_products, axis=1)
distance = np.arccos(np.clip(np.max(dot_products, axis=1), -1.0, 1.0))

uindices, ucount = np.unique(indices, return_counts=True)

ulons = mesh_lons[uindices]
ulats = mesh_lats[uindices]
In [10]:
# Set up a GeoAxis plot
fig = plt.figure(figsize=(16, 20), dpi=300)
ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=0))
ax.set_global() # type: ignore
ax.gridlines( # type: ignore
    color="0.7",
    linestyle="--",
    xlocs=np.arange(-180, 180, 15),
    ylocs=np.arange(-90, 90, 15),
)
ax.set_title(f"Jurassic Foraminifera Locations at {reconstruction_time} Ma")

# Plot shapefile features, subduction zones and MOR boundaries at 0 Ma
gplot.time = reconstruction_time  # Ma
# gplot.plot_continent_ocean_boundaries(ax, color='b', alpha=0.05)
# gplot.plot_continents(ax, facecolor='palegoldenrod', alpha=0.2)
gplot.plot_coastlines(ax, color="DarkGrey")
gplot.plot_all_topological_sections(
    ax,
    plot_subduction_teeth=True,
    other_kwargs={"color": "grey", "linewidth": 0.8},
    ridge_kwargs={"color": "red", "linewidth": 1.0},
    transform_kwargs={"color": "green", "linewidth": 1.0},
    trench_kwargs={"color": "blue", "linewidth": 1.0},
)

mask_interval = np.ones_like(ucount, dtype=bool)

sc = ax.scatter(
    ulons,
    ulats,
    s=50 + ucount, # type: ignore
    color="DarkOrange",
    edgecolor="k",
    alpha=0.5,
    transform=ccrs.PlateCarree(),
    label="Jurassic Foraminifera",
    zorder=10,
)

handles, labels = sc.legend_elements(
    prop="sizes", num=7, color="DarkOrange", markeredgecolor="k"
)
ax.legend(
    handles,
    labels,
    markerscale=2.4,
    loc="upper right",
    title="Number of\nForaminifera",
    labelspacing=4.2,
    handletextpad=2,
    bbox_to_anchor=(1.1, 1.0),
    frameon=False,
)

fig.savefig(f"{output_dir}/Reconstructed_Jurassic_Foraminifera.pdf", bbox_inches="tight")
plt.show()
No description has been provided for this image

Add feature attributes¶

Attributes can be added to each point using the add_attributes method by supplying keyword-value pairs. Some key attributes that can easily be read by GPlates include:

  • FROMAGE: the 'from' age specifies the oldest limit at which the data was active
  • TOAGE: the 'to' age specifies the youngest limit at which the data was active
  • PLATEID: the plate ID

Below, we add the FROMAGE and TOAGE attributes to the Points object and save it to a GPML file, which can be read directly by GPlates.

In [11]:
filtered_data = filtered_data.reset_index(drop=True)
gpts.add_attributes(FROMAGE=filtered_data["max_ma"], TOAGE=filtered_data["min_ma"])
# save to files
gpts.save(f"{output_dir}/output_pbdb_data.csv")
gpts.save(f"{output_dir}/output_pbdb_data.gpmlz")
print(f"Saved reconstructed data to {output_dir}/output_pbdb_data.csv and {output_dir}/output_pbdb_data.gpmlz")

Saved reconstructed data to WorkflowData/03-Working-With-Points/output/output_pbdb_data.csv and WorkflowData/03-Working-With-Points/output/output_pbdb_data.gpmlz