4 - Velocity Basics

In this notebook, we will show how to calculate and visualise plate velocity data with GPlately's PlateReconstruction and Points objects.

Let's import all needed packages:

import os, warnings
from pathlib import Path
import cartopy.crs as ccrs
import gplately
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import Normalize
from plate_model_manager import PlateModelManager
import pygplates

warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning)

data_dir = Path("WorkflowData/04-Velocity-Basics")
output_dir = data_dir / "output"
output_dir.mkdir(parents=True, exist_ok=True)

Let's first create a plate motion model using GPlately's PlateReconstruction object. To create the object, we need to pass a rotation_model, a set of pygplates topology_features and a path to a static_polygons file.

model_name = "Zahirovic2022"
pmm_model = PlateModelManager().get_model(model_name, data_dir="plate-model-repo")
assert pmm_model is not None

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

model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)

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

gplot = gplately.PlotTopologies(model, coastlines, continents, COBs)

The two approaches to calculating velocities at points

Velocities can be calculated by either:

• Using topologies:

  This is useful when you have a grid of static points that you want to calculate velocities at. The points do not reconstruct, so they remain unmoved regardless of the reconstruction time.

  This is done by calling the get_point_velocities method of a PlateReconstruction object.

  It is implemented to find which topology (topological rigid plate or deforming network) each point is contained within (at the requested reconstruction time) and calculate the plate (or network) velocity at that point. This can be thought of as dynamically assigning a plate ID to each point (at the reconstruction time). In other words, each point is essentially assigned a plate ID every time velocities are calculated (in contrast to static polygons below). This concept applies well to topological plates. For deforming networks it's a little more complicated, since a network is essentially a deforming triangulation of different plate IDs — however, the same principle applies.

• Using static polygons:

  This is useful when you have points that you want attached to plates so that they reconstruct through time.

  This is done by calling the plate_velocity method of a Points object.

  It is implemented as a once-only step of assigning a plate ID to each point using the static polygons (and storing those assigned plate IDs in the Points object) and then calculating velocities (using that Points object) at any reconstruction time (where the static plate IDs are used to calculate velocities). This differs from calculating velocities using topologies, in that the plate IDs are essentially static (unchanging through time), whereas the plate IDs are dynamic with topologies. Also, topologies typically have global coverage over the model's entire time range, whereas the static polygons only have global coverage at present day (and progressively disappear further back in time).

Velocities using topologies

Calculating velocities using topologies is done by calling the get_point_velocities method of a PlateReconstruction object.

This is useful when you have a grid of static points that you want to calculate velocities at. The points do not reconstruct, so they remain unmoved regardless of the reconstruction time.

# The distribution of points in the velocity domain: set global extent with 5 degree intervals
Xnodes = np.arange(-180, 180, 5)
Ynodes = np.arange(-90, 90, 5)

# Create a lat-lon mesh and convert to 1d lat-lon arrays
x, y = np.meshgrid(Xnodes, Ynodes)
x = x.flatten()
y = y.flatten()

# Obtain plate velocities at 249 Ma.
time = 249
#
# 'PlateReconstruction.get_point_velocities()' defaults to kms/myr (not cms/yr, the default for 'Points.plate_velocity()').
# So we explicitly specify cms/yr.
#
# Also, we return two separate velocity arrays (the east and north components) instead of a single 2D array of (north, east).
vel_x, vel_y = model.get_point_velocities(
    x,
    y,
    time,
    velocity_units=pygplates.VelocityUnits.cms_per_yr,
    return_east_north_arrays=True,
)
vel_mag = np.hypot(vel_x, vel_y)

print(f"Number of points in our velocity domain = {len(vel_x)}")
print(f"Average velocity at {time} Ma = {vel_mag.mean():.2f} cm/yr")

Number of points in our velocity domain = 2592
Average velocity at 249 Ma = 5.20 cm/yr

Calculating velocity data using the PlateReconstruction object

Let's calculate plate velocity data using get_point_velocities, a method of the PlateReconstruction object. It returns the east and north velocity components for each point we give it, at a specific reconstruction time.

We need the following parameters:

• the PlateReconstruction model
• 2 1D flattened arrays of mesh-node longitudes and latitudes representing the velocity domain
• the reconstruction time (Ma)

By default it returns a single array of (north, east) pairs for each point; passing return_east_north_arrays=True (as we did above) instead returns two separate arrays for the east and north components.

Calculating average global plate velocity through time

Global average plate velocities (cm/yr) are obtained by averaging point velocities in the velocity domain and looping over a time range.

time_range = np.arange(0, 251)

vel_av = np.zeros(time_range.size)
vel_std = np.zeros(time_range.size)

for t, time in enumerate(time_range):
    vel_x, vel_y = model.get_point_velocities(
        x,
        y,
        time,
        velocity_units=pygplates.VelocityUnits.cms_per_yr,
        return_east_north_arrays=True,
    )
    vel_mag = np.hypot(vel_x, vel_y)

    # an optional setting: if there are points in the velocity domain with a large plate velocity,
    # we can ignore these outliers. This should not be used when debugging plate models.
    ignore_outliers = True

    # Set the outlier velocity to be 50 cm/yr
    outlier_velocity = 50.0

    if ignore_outliers is True:
        vel_mag_new = [v for v in vel_mag if v < outlier_velocity]
        vel_av[t] = np.mean(vel_mag_new)
        vel_std[t] = np.std(vel_mag_new)
    else:
        vel_av[t] = np.mean(vel_mag_new)
        vel_std[t] = np.std(vel_mag)

# save to a CSV file

output_data = np.column_stack([time_range, vel_av, vel_std])

header = "Time (Ma),Mean plate velocities (cm/yr),Standard deviation (cm/yr)"

np.savetxt(
    os.path.join(
        output_dir,
        "GlobalAveragePlateVelocities.csv",
    ),
    output_data,
    delimiter=",",
    header=header,
    comments="",
)

fig = plt.figure(figsize=(8, 4), dpi=100)
ax1 = fig.add_subplot(
    111,
    xlim=(250, 0),
    ylim=(0, 10),
    xlabel="Age (Ma)",
    ylabel="Velocity (cm/yr)",
    title="Global average plate velocity (cm/yr)",
)

ax1.fill_between(
    time_range,
    vel_av - vel_std,
    vel_av + vel_std,
    color="0.8",
    label="Standard deviation (cm/yr)",
)
ax1.plot(time_range, vel_av, c="k", label="Mean plate velocity (cm/yr)")

ax1.legend(loc="upper right", frameon=False)
fig.savefig(
    os.path.join(
        output_dir,
        "average_plate_velocity.pdf",
    ),
    bbox_inches="tight",
)

Visualising PlateReconstruction velocity data

As a first example, let's reconstruct all topological plates and boundaries to 50 Ma and illustrate the velocity of each moving plate! One way to do this is by plotting a velocity vector field using the plot_plate_motion_vectors method on the PlotTopologies object (which internally uses PlateReconstruction.get_point_velocities to calculate velocities on a regular longitude-latitude grid of points).

Since plot_plate_motion_vectors uses Cartopy's quiver function, it accepts quiver keyword arguments like regrid_shape. This is useful if you'd like your vectors interpolated onto a regular grid in projection space.

time = 50

# Set up a GeoAxis plot
fig = plt.figure(figsize=(16, 12))
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=0))
ax1.gridlines(
    color="0.7",
    linestyle="--",
    xlocs=np.arange(-180, 180, 15),
    ylocs=np.arange(-90, 90, 15),
)
plt.title("Global plate motion velocity field at %i Ma" % (time))

# Plot all topologies
gplot.time = time
gplot.plot_continents(ax1, facecolor="navajowhite")
gplot.plot_coastlines(ax1, color="orange")
gplot.plot_all_topological_sections(
    ax1,
    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},
)
ax1.set_global()

# Plot a velocity vector field
#
# Use a 10 degree longitude-latitude spacing between points.
gplot.plot_plate_motion_vectors(
    ax1, spacingX=10, spacingY=10, regrid_shape=20, alpha=0.5, color="green", zorder=2
)
plt.show()

Plotting a velocity streamplot

We can visualise the same data with streamplot from matplotlib.

# Set up a GeoAxis plot
fig = plt.figure(figsize=(16, 12))
ax2 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=0))
ax2.gridlines(
    color="0.7",
    linestyle="--",
    xlocs=np.arange(-180, 180, 15),
    ylocs=np.arange(-90, 90, 15),
)
plt.title("Global plate motion velocity streamplot at %i Ma" % (time))

# Plot all topologies
gplot.time = time  # Ma
gplot.plot_continents(ax2, facecolor="0.95")
gplot.plot_coastlines(ax2, color="0.9")
gplot.plot_all_topological_sections(
    ax2,
    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},
)
ax2.set_global()

vel_x, vel_y = model.get_point_velocities(x, y, time, return_east_north_arrays=True)
vel_mag = np.hypot(vel_x, vel_y)

ax2.streamplot(
    x,
    y,
    vel_x,
    vel_y,
    color=vel_mag,
    transform=ccrs.PlateCarree(),
    linewidth=0.02 * vel_mag,
    cmap=plt.cm.turbo,
    density=2,
)
plt.show()

# Set up a GeoAxis plot
fig = plt.figure(figsize=(12, 4))
norm = Normalize(0, 10)

for i, time in enumerate([80, 60, 40, 20]):
    ax2 = fig.add_subplot(
        1, 4, i + 1, projection=ccrs.Orthographic(70, 0), title="{:.0f} Ma".format(time)
    )
    ax2.set_global()
    ax2.gridlines(
        color="0.7",
        linestyle=":",
        xlocs=np.arange(-180, 180, 15),
        ylocs=np.arange(-90, 90, 15),
    )
    # plt.title('Global plate motion velocity streamplot at %i Ma' % (time))

    #     gplot.plot_grid(ax2, rgb)

    # Plot  topologies
    gplot.time = time  # Ma
    gplot.plot_continents(ax2, facecolor="0.9")
    gplot.plot_coastlines(ax2, color="0.7")
    gplot.plot_all_topological_sections(
        ax2,
        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},
    )

    vel_x, vel_y = model.get_point_velocities(x, y, time, return_east_north_arrays=True)
    vel_mag = np.hypot(vel_x, vel_y)

    sp = ax2.streamplot(
        x,
        y,
        vel_x,
        vel_y,
        color=vel_mag * 0.1,
        transform=ccrs.PlateCarree(),
        norm=norm,
        linewidth=0.01 * vel_mag,
        cmap="plasma",
        density=1,
    )

fig.subplots_adjust(wspace=0.05)

cax = plt.axes((0.36, 0.1, 0.3, 0.04))
fig.colorbar(
    sp.lines, # type: ignore
    cax=cax,
    orientation="horizontal",
    label="Plate velocity (cm/yr)",
    extend="max",
)

fig.savefig(
    os.path.join(
        output_dir,
        "India_collision.pdf",
    ),
    bbox_inches="tight",
)

Create a velocity animation

If you have moviepy installed, you can animate the motion of topological plates through geological time with a streamplot of domain point velocities (in cm/yr) overlying the plates. Let's reconstruct plate movements over 100 Ma in intervals of 10 Ma.

def generate_frame(output_filename, time):
    vel_x, vel_y = model.get_point_velocities(x, y, time, return_east_north_arrays=True)
    vel_mag = np.hypot(vel_x, vel_y)

    # Set up a GeoAxis plot
    fig = plt.figure(figsize=(16, 12))
    ax2 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=0))
    ax2.gridlines(
        color="0.7",
        linestyle="--",
        xlocs=np.arange(-180, 180, 15),
        ylocs=np.arange(-90, 90, 15),
    )
    plt.title("Global plate motion velocity streamplot at %i Ma" % (time))

    # Reconstruct topological plates and boundaries with PlotTopologies
    gplot.time = time  # Ma
    gplot.plot_continents(ax2, facecolor="0.95")
    gplot.plot_coastlines(ax2, color="0.9")
    gplot.plot_all_topological_sections(
        ax2,
        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},
    )
    ax2.set_global()

    # Create the streamplot, using speed as a colormap.
    ax2.streamplot(
        x,
        y,
        vel_x,
        vel_y,
        color=vel_mag,
        transform=ccrs.PlateCarree(),
        linewidth=0.02 * vel_mag,
        cmap=plt.cm.turbo,
        density=2,
    )

    fig.savefig(output_filename, bbox_inches="tight")
    plt.close(fig)

import tempfile
from IPython.display import Image

try:
    from moviepy.editor import ImageSequenceClip  # moviepy 1.x
except ImportError:
    from moviepy import ImageSequenceClip  # moviepy 2.x

# Time variables
oldest_seed_time = 100  # Ma
time_step = 10  # Ma

frame_list = []
# Create a plot for each 10 Ma interval
for time in np.arange(oldest_seed_time, 0.0, -time_step):
    print(f"Generating {time} Ma frame...")
    frame_filename = os.path.join(
        output_dir, f"plate_velocity_stream_plot_frame_{time}_Ma.png"
    )
    generate_frame(frame_filename, time)
    frame_list.append(frame_filename)

video_filename = os.path.join(output_dir, "plate_velocity_stream_plot.gif")

clip = ImageSequenceClip(frame_list, fps=5)
clip.write_gif(video_filename)

print("The movie will show up in a few seconds...")
with open(video_filename, "rb") as f:
    display(Image(data=f.read(), format="png", width=1000, height=500))

Generating 100.0 Ma frame...

Generating 90.0 Ma frame...

Generating 80.0 Ma frame...

Generating 70.0 Ma frame...

Generating 60.0 Ma frame...

Generating 50.0 Ma frame...

Generating 40.0 Ma frame...

Generating 30.0 Ma frame...

Generating 20.0 Ma frame...

Generating 10.0 Ma frame...

MoviePy - Building file WorkflowData/04-Velocity-Basics/output/plate_velocity_stream_plot.gif with imageio.

frame_index:   0%|          | 0/10 [00:00<?, ?it/s, now=None]

frame_index:  70%|███████   | 7/10 [00:00<00:00, 62.65it/s, now=None]

The movie will show up in a few seconds...

Velocities using static polygons

Calculating velocities using static polygons is done by calling the plate_velocity method of a Points object.

This is useful when you have points that you want attached to plates so that they reconstruct through time.

Calculating velocities of reconstructed points using the Points object

PlateReconstruction.get_point_velocities calculates velocities of topological plates (and topological deforming networks) on a static grid of points (as shown above). What is meant by static grid is that the point positions don't change over time. This is different to the static polygons which don't change shape over time (but do change position, or get reconstructed, over time).

Another way to calculate velocities at point locations is to specify points at an initial time in a Points object. Then we can reconstruct them to any reconstruction time and calculate velocities at the reconstructed locations. Hence these points are not static. So, instead of using topological features to calculate velocities (with PlateReconstruction.get_point_velocities), a Points object optionally assigns a plate ID to each point (using a model's static polygons dataset) and uses that (along with the model's rotation model) to both reconstruct each point and calculate each point's velocity at any reconstruction time (using Points.plate_velocity).

A Points object requires a PlateReconstruction object for its rotation model and static polygons. A Points object also requires initial point locations and an initial time. By default, the initial time is present day (0 Ma), in which case the initial locations represent the present-day locations of the points. However, if the initial time is not present day then the initial locations represent locations at the initial time.

Let's calculate plate velocity data for a global grid of longitude-latitude points (x and y). We'll use the plate_velocity method in the Points object. It returns the east and north components of velocities for each point, at the specified reconstruction time.

# The distribution of points in the velocity domain: set global extent with 5 degree intervals
Xnodes = np.arange(-180, 180, 5)
Ynodes = np.arange(-90, 90, 5)

# Create a lat-lon mesh and convert to 1d lat-lon arrays
x, y = np.meshgrid(Xnodes, Ynodes)
x = x.flatten()
y = y.flatten()

# Create a 'gplately.Points' object using the longitude-latitude grid of points (x, y).
#
# These points represent locations at 0 Ma (the default initial time).
#
# Note: These positions are not static throughout time. They get reconstructed using their plate IDs
#       (assigned using 'model's static polygons).
gpts = gplately.Points(model, x, y)

# Obtain plate velocities at 0 Ma.
#
# 'Point.plate_velocity()' defaults to cms/yr (not kms/myr, the default for 'PlateReconstruction.get_point_velocities()').
# So we don't need to explicitly specify cms/yr.
#
# Also, by default it returns two separate velocity arrays (the east and north components) instead of a single 2D array of (north, east),
# unlike 'PlateReconstruction.get_point_velocities()' for which you must specify `return_east_north_arrays=True` for separate arrays.
vel_x, vel_y = gpts.plate_velocity(time=0)
vel_mag = np.hypot(vel_x, vel_y)

print(f"Number of points in our velocity domain = {len(vel_x)}")
print(f"Average velocity at 0 Ma = {vel_mag.mean():.2f} cm/yr")

Number of points in our velocity domain = 2590
Average velocity at 0 Ma = 2.88 cm/yr

Visualising Points object velocity data

We can visualise all the points on our velocity domain on a scatterplot - they will be colour mapped with their velocity magnitudes. We use Matplotlib's scatter function for the scatterplot.

In addition to the scatterplot, we can optionally also visualise a plot of the velocity arrows using Matplotlib's quiver function, where we also colour map their velocity magnitudes.

For these plots we use an array of latitudes and an array of longitudes to plot point data on. These will be the initial point locations that are subsequently reconstructed to the time at which velocities are calculated. Note that, when plotting velocities at the reconstructed locations, we cannot use the initial point locations (unless the initial time matches the reconstruction time).

def plot_point_velocities(
    vel_time, vel_x, vel_y, rlons, rlats, draw_arrows=False, plot_time=None
):
    # It's possible to calculate velocities at one time but plot them at another.
    if plot_time is None:
        plot_time = vel_time

    # Velocity magnitudes.
    vel_mag = np.hypot(vel_x, vel_y)

    # Set up a GeoAxis plot
    fig = plt.figure(figsize=(18, 14))
    ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=0))

    # Plot continents, coastlines and topologies.
    gplot.time = plot_time  # Ma
    gplot.plot_continents(ax, facecolor="0.95")
    gplot.plot_coastlines(ax, color="0.9")
    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},
    )

    title = f"Global point velocities at {vel_time} Ma"
    if plot_time != vel_time:
        title += f" (plotted at {plot_time} Ma)"
    plt.title(title)

    # Plot the velocity domain points with their velocity magnitudes as a colour scale.
    im = ax.scatter(
        rlons,
        rlats,
        transform=ccrs.PlateCarree(),
        c=vel_mag,
        s=10,
        cmap=plt.cm.afmhot_r,
        vmin=0,
        vmax=10,
        zorder=2,
    )

    if draw_arrows:
        # Plot the velocity arrows (with their velocity magnitudes as a colour scale).
        im = ax.quiver(
            rlons,
            rlats,
            vel_x,
            vel_y,
            vel_mag,
            transform=ccrs.PlateCarree(),
            cmap=plt.cm.afmhot_r,
            clim=(0, 10),
            regrid_shape=None,
            zorder=2,
        )

    # Add colorbar and set global extent.
    fig.colorbar(im, ax=ax, shrink=0.5).set_label(
        "Velocity magnitude (cm/yr)", fontsize=12
    )
    ax.set_global()

    plt.show()

Present day Points

First we calculate velocities at present day (0 Ma).

Since the gpts object is being created at present day (the default initial time), we can just use the present day locations of its points (gpts.lons and gpts.lats) when plotting velocities at present day. However, this is just to demonstrate that it's possible in this particular scenario. Ideally you would always use the reconstructed locations, especially when the initial time is not present day (as we'll show later).

# Create points with initial positions at present day (0 Ma).
#
# These points represent locations at 0 Ma (the default initial time).
gpts = gplately.Points(model, x, y)

# Present day (0 Ma).
time = 0

# Calculate the velocities at present day (0 Ma), returning locations aligned with velocity arrays.
vel_x, vel_y, rlons, rlats = gpts.plate_velocity(time, return_reconstructed_points=True)

# Plot the velocities at present day (0 Ma).
plot_point_velocities(time, vel_x, vel_y, rlons, rlats, draw_arrows=True)

Velocities at a past time

However, if we calculate velocities at a past time (more specifically at a time that's different to the initial time of the gpts object) then we can no longer use the initial point locations (in our case present day locations) when plotting. This is because the points have been reconstructed to a new location.

So we ask gpts.plate_velocity() to also return the reconstructed point locations (using the return_reconstructed_points argument). And then we plot velocities at those reconstructed locations.

# A time in the past.
time = 100

# Calculate the velocities at the past time.
vel_x, vel_y, rlons, rlats = gpts.plate_velocity(time, return_reconstructed_points=True)

# Plot the velocities at the past time.
#
# Note: We need the locations of the initial points (provided to the 'gpts' object) reconstructed to 'time'.
#       We can no longer use the initial point locations ('gpts.lons' and 'gpts.lats').
plot_point_velocities(time, vel_x, vel_y, rlons, rlats)

Dealing with oceanic points

By default, all points exist for all time. However this is not true for oceanic points. And this is why the above plot is drawing points that it should not.

To fix this we can specify age=None to gpts.plate_velocity() (instead of relying on the default age=numpy.inf). This will use the static polygons of the model to assign a time of appearance to each point. This way, if we calculate velocities at a time prior to the age of appearance of some oceanic points, then those points will now be missing from the output (as desired).

# Create points with initial positions at present day (0 Ma).
#
# With 'age=None': the age of each point (its time of appearance) is determined from the static polygons of 'model'.
gpts = gplately.Points(model, x, y, age=None)

# Calculate the velocities at the past time.
vel_x, vel_y, rlons, rlats = gpts.plate_velocity(time, return_reconstructed_points=True)

# Plot the velocities at the past time.
plot_point_velocities(time, vel_x, vel_y, rlons, rlats)

Plotting past velocities at present day

Normally we plot the velocities at the reconstructed point locations (as above). However we can plot velocities at the initial point locations (present day locations in our case). This is sometimes called the plate reference frame.

To do this we ask gpts.plate_velocity() to also return the indices of the reconstructed points (using the return_point_indices argument). These are indices into the initial points passed into the gpts object. Then we can index into the initial points and plot at those locations (present day locations in our case).

NOTE: These velocities still represent velocities at the reconstruction time, not velocities at present day. This also means the velocity directions are not un-rotated back to present day (only the point locations are essentially un-rotated back to present day, because we're plotting the initial locations). So it might be confusing to plot velocity arrows. Instead we just map velocity magnitude using a colour scale.

# Calculate the velocities at 'time'.
#
# But also return the indices into the initial points (passed into 'gpts').
vel_x, vel_y, point_indices = gpts.plate_velocity(time, return_point_indices=True)

# Use the point indices to find the corresponding initial point locations (at the initial time).
#
# Note that the length of 'point_indices' is the same length as 'vel_x' and 'vel_y',
# but is smaller than the length of 'gpts.lons' and 'gpts.lats'.
initial_lons = gpts.lons[point_indices]
initial_lats = gpts.lats[point_indices]

# Plot the velocities calculated at 'time'.
# But at the locations of the points at the initial time 'gpts.time' (instead of the reconstructed point locations at 'time').
# In our case 'gpts.time' is zero (ie, present day).
plot_point_velocities(
    time, vel_x, vel_y, initial_lons, initial_lats, plot_time=gpts.time
)

Initialising Points at a past time

So far we've only created a Points object using present day point locations. Because the static polygons all exist at present day, they also have global coverage at present day, and so we could assign plate IDs (and ages) to all our points (passed into the Points object).

However, if we create a Points object using point locations in the past then the static polygons are also reconstructed to that same past time (to assign plate IDs). But the reconstructed static polygons will no longer have global coverage (ie, some oceanic polygons will disappear) and hence some of the points will not be assigned plate IDs (and/or ages).

# A time in the past.
time = 100

# Create points with initial positions at 100 Ma (ie, not present day).
#
# Note that we're now specifying the 'time' argument.
# Previously we did not specify it, and so it would default to 0 Ma.
#
# With 'age=None': the age of each point (its time of appearance) is determined from the static polygons of 'model'.
#
# With 'remove_unreconstructable_points=True': we're removing any points that were not assigned plate IDs and ages.
gpts = gplately.Points(
    model, x, y, time, age=None, remove_unreconstructable_points=True
)

Note that the number of points stored in the gpts object is less than the number of points we passed into it (x and y).

This is because we requested that gplately.Points only store points that are reconstructable (ie, that were assigned a plate ID and age using the static polygons reconstructed to the initial time 100 Ma).

print(
    f'Number of points passed into "gpts" at initial time ({gpts.time} Ma) = {len(x)}'
)
print(
    f'Number of points actually in "gpts" at initial time ({gpts.time} Ma) = {gpts.size}'
)

Number of points passed into "gpts" at initial time (100 Ma) = 2592
Number of points actually in "gpts" at initial time (100 Ma) = 1515

From now on, let's just keep all points, including those that are not reconstructable.

Even though not all points will be reconstructable, we can still access the initial locations (and associated plate IDs and ages) of all points (using gpts.lons, gpts.lats, gpts.plate_id and gpts.age; which all have a length of gpts.size).

# Create points with initial positions at 100 Ma (ie, not present day).
#
# If we do not specify 'remove_unreconstructable_points' then it defaults to False.
gpts = gplately.Points(model, x, y, time, age=None)

Note that the number of points stored in the gpts object now equals the number of points we passed into it (x and y). But that doesn't mean they are all reconstructable.

print(
    f'Number of points passed into "gpts" at initial time ({gpts.time} Ma) = {len(x)}'
)
print(
    f'Number of points actually in "gpts" at initial time ({gpts.time} Ma) = {gpts.size}'
)

Number of points passed into "gpts" at initial time (100 Ma) = 2592
Number of points actually in "gpts" at initial time (100 Ma) = 2592

Velocities at the initial time

Next we calculate and plot the velocities at the initial time (100 Ma).

# Calculate the velocities at the initial time.
vel_x, vel_y, rlons, rlats = gpts.plate_velocity(time, return_reconstructed_points=True)

# Plot the velocities at the initial time.
plot_point_velocities(time, vel_x, vel_y, rlons, rlats, draw_arrows=True)

Note that the number of points in the above plot is now less than the number of points stored in the gpts object, even though we calculated velocities at the initial time (ie, the same time the gpts object was created).

This is because not all the initial points could be assigned ages (and plate IDs) at the initial time (100 Ma). And hence not all of the initial points are reconstructable, which also means we cannot calculate velocities for them.

print(f'Number of points in "gpts" at initial time ({gpts.time} Ma) = {gpts.size}')
print(f"Number of reconstructed points at initial time ({time} Ma) = {len(rlons)}")

Number of points in "gpts" at initial time (100 Ma) = 2592
Number of reconstructed points at initial time (100 Ma) = 1515

Velocities at present day

Now let's calculate and plot velocities at present day. These are the initial points at 100 Ma reconstructed to present day (and velocities calculated at present day).

# Present day (0 Ma).
time = 0

# Calculate the velocities at present day.
vel_x, vel_y, rlons, rlats = gpts.plate_velocity(time, return_reconstructed_points=True)

# Plot the velocities at present day.
plot_point_velocities(time, vel_x, vel_y, rlons, rlats, draw_arrows=True)

Note that the number of points in the above plot at present day is equal to the number of points in the plot at the initial time (100 Ma).

Again, this is because not all of the initial points are reconstructable. And it's the same number of points that are not reconstructable.

print(f'Number of points in "gpts" at initial time ({gpts.time} Ma) = {gpts.size}')
print(f'Number of points in "gpts" reconstructed to present day (0 Ma) = {len(rlons)}')

Number of points in "gpts" at initial time (100 Ma) = 2592
Number of points in "gpts" reconstructed to present day (0 Ma) = 1515

Velocities at a past time

Let's calculate and plot at one more time. This time at 200 Ma, which is 100 Myr older than the initial time (100 Ma).

# A past time that is 100 Myr older than the initial time of 'gpts'.
time = gpts.time + 100

# Calculate the velocities at the past time.
vel_x, vel_y, rlons, rlats = gpts.plate_velocity(time, return_reconstructed_points=True)

# Plot the velocities at the past time.
plot_point_velocities(time, vel_x, vel_y, rlons, rlats, draw_arrows=True)

Note that the number of points in the above plot (at a time older time than the initial time) is less than the number of points in the plot at the initial time (100 Ma).

This is because fewer of the reconstructable points exist at the older time. In other words, some of the oceanic points that existed at the initial time (100 Ma) did not exist at the older time (200 Ma).

print(f'Number of points in "gpts" at initial time ({gpts.time} Ma) = {gpts.size}')
print(f'Number of points in "gpts" reconstructed to {time} Ma = {len(rlons)}')

Number of points in "gpts" at initial time (100 Ma) = 2592
Number of points in "gpts" reconstructed to 200 Ma = 1197

Create a velocity animation

If you have moviepy installed, you can animate the motion of topological plates through geological time with a scatterplot of domain point velocities (in cm/yr) overlying the plates. Let's reconstruct plate movements from 0-100 Ma in intervals of 10 Ma. With each iteration of the time loop we re-calculate velocity data.

def generate_frame(output_filename, time):
    # Get all point velocities and their magnitudes
    vel_x, vel_y, rlons, rlats = gpts.plate_velocity(
        time, return_reconstructed_points=True
    )
    vel_mag = np.hypot(vel_x, vel_y)

    # Set up a GeoAxis plot
    fig = plt.figure(figsize=(18, 14))
    ax3 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=0))
    plt.title(f"Global point velocity scatterplot at {time} Ma")

    # Plot all topologies reconstructed to the current Ma
    gplot.time = time  # Ma
    gplot.plot_continents(ax3, facecolor="0.95")
    gplot.plot_coastlines(ax3, color="0.9")
    gplot.plot_all_topological_sections(
        ax3,
        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},
    )

    # Plot the velocity domain points with their velocity magnitudes as a colour scale.
    im = ax3.scatter(
        rlons,
        rlats,
        transform=ccrs.PlateCarree(),
        c=vel_mag,
        s=30,
        cmap=plt.cm.afmhot_r,
        vmin=0,
        vmax=10,
        zorder=2,
    )

    # Add colorbar, set global extent and show plot
    cbar = fig.colorbar(im, ax=ax3, shrink=0.5)
    cbar.set_label("Velocity magnitude (cm/yr)", fontsize=12)

    ax3.set_global()
    fig.savefig(output_filename, dpi=300, bbox_inches="tight")
    plt.close(fig)

from IPython.display import Image

try:
    from moviepy.editor import ImageSequenceClip  # moviepy 1.x
except ImportError:
    from moviepy import ImageSequenceClip  # moviepy 2.x

# Time variables
oldest_seed_time = 100  # Ma
time_step = 10  # Ma

frame_list = []
# Create a plot for each 10 Ma interval
for time in np.arange(oldest_seed_time, 0.0, -time_step):
    print(f"Generating {time} Ma frame...")
    frame_filename = os.path.join(
        output_dir, f"plate_velocity_scatter_plot_frame_{time}_Ma.png"
    )
    generate_frame(frame_filename, time)
    frame_list.append(frame_filename)

video_filename = os.path.join(output_dir, "plate_velocity_scatter_plot.gif")

clip = ImageSequenceClip(frame_list, fps=5)
clip.write_gif(video_filename)

print("The movie will show up in a few seconds...")
with open(video_filename, "rb") as f:
    display(Image(data=f.read(), format="png", width=1000, height=500))

Generating 100.0 Ma frame...

Generating 90.0 Ma frame...

Generating 80.0 Ma frame...

Generating 70.0 Ma frame...

Generating 60.0 Ma frame...

Generating 50.0 Ma frame...

Generating 40.0 Ma frame...

Generating 30.0 Ma frame...

Generating 20.0 Ma frame...

Generating 10.0 Ma frame...

MoviePy - Building file WorkflowData/04-Velocity-Basics/output/plate_velocity_scatter_plot.gif with imageio.

frame_index:   0%|          | 0/10 [00:00<?, ?it/s, now=None]

frame_index:  30%|███       | 3/10 [00:00<00:00, 13.92it/s, now=None]

frame_index:  50%|█████     | 5/10 [00:00<00:00,  8.70it/s, now=None]

frame_index:  60%|██████    | 6/10 [00:00<00:00,  7.95it/s, now=None]

frame_index:  70%|███████   | 7/10 [00:00<00:00,  7.38it/s, now=None]

frame_index:  80%|████████  | 8/10 [00:01<00:00,  7.04it/s, now=None]

frame_index:  90%|█████████ | 9/10 [00:01<00:00,  6.81it/s, now=None]

frame_index: 100%|██████████| 10/10 [00:01<00:00,  6.67it/s, now=None]

The movie will show up in a few seconds...

You may have noticed the banding problem of the color bar in the .gif above. GIF's 256-color limit means smooth gradients get reduced to visible steps. To avoid this banding problem, we can create a .mp4 file instead.

mp4_filename = os.path.join(output_dir, "plate_velocity_scatter_plot.mp4")
clip = ImageSequenceClip(frame_list, fps=5)
clip.write_videofile(
    mp4_filename,
    fps=5,
    codec="libx264",
    audio=False,
    ffmpeg_params=[
        "-pix_fmt",
        "yuv420p",
        "-vf",
        "scale=trunc(iw/2)*2:trunc(ih/2)*2",
    ],
)

from IPython.display import Video

display(Video(mp4_filename, embed=True, mimetype="video/mp4", width=600))

MoviePy - Building video WorkflowData/04-Velocity-Basics/output/plate_velocity_scatter_plot.mp4.
MoviePy - Writing video WorkflowData/04-Velocity-Basics/output/plate_velocity_scatter_plot.mp4

frame_index:   0%|          | 0/10 [00:00<?, ?it/s, now=None]

frame_index:  30%|███       | 3/10 [00:00<00:00,  7.98it/s, now=None]

frame_index:  40%|████      | 4/10 [00:00<00:01,  5.59it/s, now=None]

frame_index:  50%|█████     | 5/10 [00:00<00:01,  4.62it/s, now=None]

frame_index:  60%|██████    | 6/10 [00:01<00:00,  4.19it/s, now=None]

frame_index:  70%|███████   | 7/10 [00:01<00:00,  3.90it/s, now=None]

frame_index:  80%|████████  | 8/10 [00:01<00:00,  3.75it/s, now=None]

frame_index:  90%|█████████ | 9/10 [00:02<00:00,  3.65it/s, now=None]

frame_index: 100%|██████████| 10/10 [00:02<00:00,  3.56it/s, now=None]

MoviePy - Done !
MoviePy - video ready WorkflowData/04-Velocity-Basics/output/plate_velocity_scatter_plot.mp4