9 - Creating Motion Paths and Flowlines¶

In this notebook, we will show how GPlately's PlateReconstruction object:

  • Creates a motion path to illustrate the trajectory of a tectonic plate through geological time, and
  • Creates a flowlines to track the motion of tectonic plates from spreading features like mid-ocean ridges.

Also, the rate of motion through time, of the motion path and flowlines, are quantified.

In [1]:
import gplately

import numpy as np
import pygplates
import glob, os
import matplotlib.pyplot as plt
import matplotlib.axes as axs
import cartopy.crs as ccrs
import matplotlib
from matplotlib import image
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from mpl_toolkits.axes_grid1 import make_axes_locatable
from shapely.geometry import LineString
import geopandas as gpd
from plate_model_manager import PlateModelManager, PresentDayRasterManager

%matplotlib inline

We use rotations from the WK08-A absolute plate motion model by Wessel and Kroenke (2008), which can be found in the Matthews et al. 2016 plate model. Let's request this model using the PlateModelManager with the model name Matthews2016, and set up a PlateReconstruction object.

Wessel, P., & Kroenke, L. W. (2008). Pacific absolute plate motion since 145 Ma: An assessment of the fixed hot spot hypothesis. Journal of Geophysical Research: Solid Earth, 113(B6). https://doi.org/10.1029/2007JB005499

Set up a PlotTopologies object for plotting.

In [2]:
# Download Matthews2016 model
pm_manager = PlateModelManager()
matthews2016_model = pm_manager.get_model("Matthews2016", data_dir="plate-model-repo")

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

coastlines = matthews2016_model.get_layer('Coastlines')

# Create a PlateReconstruction object
model = gplately.reconstruction.PlateReconstruction(rotation_model, topology_features, static_polygons)

# Create a PlotTopologies object
gplot = gplately.PlotTopologies(model, coastlines=coastlines, continents=None, COBs=None, time=0)

1. The motion path of a Hawaiian-Emperor chain seed point¶

The longitude and latitude coordinates of present-day Hawaii is (-155, 19). The easiest way to understand motion paths is as a hot spot trail, where we want to trace the path that a hot spot in the relatively stationary reference frame (moving_plate_ID = 2) imprints on the overriding Pacific plate (relative_plate_ID = 901). This 901-2 reconstruction ID pair is from Wessel and Kroenke's (2008) WK08-A model.

Step plots of rates of motion¶

By setting return_times and return_rate_of_motion to True in create_motion_path, the function returns two additional arrays: one contains the times array, and another contains the rate of plate motion.

In [3]:
# Latitude and longitude.
# Note: These can be arrays if you have more than one hot spot trail.
lat = 19
lon = -155

# Plate IDs
moving_plate_ID = 2
relative_plate_ID = 901

# Create the time array for the motion path - must be float
start_reconstruction_time = 0.
time_step = 2.
max_reconstruction_time = 100.
time_array = np.arange(start_reconstruction_time, max_reconstruction_time + time_step, time_step)

# Get the latitudes and longitudes of all points along the motion path
rlons, rlats, rtimes, rates_of_motion = model.create_motion_path(
    lon,
    lat,
    time_array,
    plate_id=moving_plate_ID,
    relative_plate_id=relative_plate_ID,
    return_times=True,
    return_rate_of_motion=True)

Let's plot this motion path with the ETOPO1 global topology and bathymetry model. We can download this grid to the GPlately cache using PresentDayRasterManager:

In [ ]:
from matplotlib import image
raster_manager = PresentDayRasterManager()
etopo_tif_img = gplately.Raster(
    data=image.imread(raster_manager.get_raster("ETOPO1_tif")),
    plate_reconstruction=model)
etopo_tif_img.resample(0.25, 0.25, inplace=True)
etopo_tif_img.lats = etopo_tif_img.lats[::-1]
In [5]:
from matplotlib.colors import ListedColormap

def plot_motion_path(rlons, rlats, rtimes, rates_of_motion, time=0):
    # --- Create figure
    fig = plt.figure(figsize=(12,5), dpi=200)
    plt.suptitle("Motion path of the Hawaiian-Emperor Chain (at {} Ma)".format(time), weight="demi")
    gs = fig.add_gridspec(1, 7)
    
    ax = fig.add_subplot(gs[0, :4], projection=ccrs.Mercator(central_longitude=180))
    ax.set_title("Motion paths", fontsize=14)
    
    # --- Plot tick labels
    ax.set_xticks(np.arange(0, 360, 20), crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(-90, 100, 10), crs=ccrs.PlateCarree())
    ax.set_xticks(np.arange(0, 360, 10), minor=True, crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(-90, 90, 5), minor=True, crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=True)
    lat_formatter = LatitudeFormatter()
    ax.xaxis.tick_bottom()  # labelled ticks on top
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    ax.xaxis.set_ticks_position('both')   # ticks on both top/bottom
    ax.yaxis.set_ticks_position('both')   # ticks on both left/right
    ax.grid()
    
    # --- Limit map extent
    lon_min = 130
    lon_max = 220
    lat_min = 0
    lat_max = 60
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())

    gplot.time = time
    
    # --- Plot filled coastlines
    gplot.plot_coastlines(ax, facecolor="white", edgecolor=None, zorder=2)
    
    # ---- Plot bathymetry
    etopo_tif_img.reconstruct(time=time, fill_value='darkblue').imshow(alpha=0.5, zorder=1)
    
    # --- Make sure negtive motion path longitudes are wrapped correctly to the dateline
    rlons = gplately.tools.correct_longitudes_for_dateline(rlons)
    
    # --- plot motion path
    colors = ['yellow']
    # The mapping of 'rtimes' to colours should not change.
    # Ie, each time should map to the same colour (regardless of the time range).
    rtimes_min, rtimes_max = rtimes[0], rtimes[-1]
    cmap = ListedColormap(plt.cm.inferno(np.linspace(rtimes_min/rtimes_max, 1.0, 256)))
    for i in range(len(rlons)):
        ax.plot(rlons[i], rlats[i], colors[i], linewidth=0.75, transform=ccrs.PlateCarree(), zorder=3)
        mp = ax.scatter(rlons[i], rlats[i], 100, marker='.', c=rtimes, cmap=cmap, edgecolor='k',
                        transform=ccrs.PlateCarree(), vmin=rtimes_min, vmax=rtimes_max, zorder=4)
        
    # --- Set a colorbar to help visualise the passage of every Xth Myr.
    divider = make_axes_locatable(ax)
    cax = divider.new_horizontal(size="4%", pad=0.4, axes_class=plt.Axes)
    fig.add_axes(cax)
    cbar = plt.colorbar(mp, cax=cax, orientation="vertical")
    cbar.set_label('Age (Ma)', fontsize=12)
    cbar.ax.minorticks_on()

    # --- Plot times and rates of motion on 2nd subplot
    p = fig.add_subplot(gs[0, 5:], xlabel='Time (Ma)', ylabel='Rate of motion (cm/yr)',
                        xlim=[rtimes.max(), 0],ylim=[0.0,12.0])
    p.set_title("Rate of motion")
    p.tick_params(direction="in", length=5, top=True, right=True)
    p.tick_params(direction="in", which='minor', length=2.5, top=True, right=True)
    p.minorticks_on()
    # For each seed point plot its rate of motion.
    for rate_of_motion in rates_of_motion:
        plt.stairs(rate_of_motion, rtimes)

    plt.show()

    fig.savefig("Hawaii_Emperor_motion_path_{}Ma.pdf".format(time), bbox_inches='tight')

Plot the motion path at present day.

In [6]:
plot_motion_path(rlons, rlats, rtimes, rates_of_motion)
No description has been provided for this image

Now plot the motion path at a past time.

To do this we specify the reconstruction time using the to_time parameter of PlateReconstruction.create_motion_path.

The from_time parameter defaults to present day. This is what we want since the hot spot coordinates are present day positions.

Also, we are using the relative_plate_id parameter (instead of anchor_plate_id). This allows us to view the motion path in the default frame of reference (when reconstructing the motion path to past times).

In [7]:
time = 20

# Reconstruct the motion path to a past time.
rlons, rlats, rtimes, rates_of_motion = model.create_motion_path(
    lon,
    lat,
    time_array,
    # Specify the reconstruction time...
    to_time=time,
    plate_id=moving_plate_ID,
    # Note that we specify 'relative_plate_id' and leave 'anchor_plate_id' as default (of 'model')...
    relative_plate_id=relative_plate_ID,
    return_times=True,
    return_rate_of_motion=True)

plot_motion_path(rlons, rlats, rtimes, rates_of_motion, time=time)
No description has been provided for this image

Now do the same but using a Points object (instead of PlateReconstruction).

To do this we first create a Points object at present day since the hot spot coordinates are present day positions.

Then we create a motion path using Points.motion_path and specify the reconstruction time with its time parameter.

In [8]:
time = 20

# Create a 'Points' object at present day (since point coordinates are present day positions).
gpts = gplately.Points(model, lon, lat, time=0, plate_id=moving_plate_ID)

# Reconstruct the motion path to a past time using our 'Points' object.
rlons, rlats, rtimes, rates_of_motion = gpts.motion_path(
    time_array,
    # Specify the reconstruction time...
    time=time,
    # Note that we specify 'relative_plate_id' and leave 'anchor_plate_id' as default (of 'model')...
    relative_plate_id=relative_plate_ID,
    return_times=True,
    return_rate_of_motion=True)

plot_motion_path(rlons, rlats, rtimes, rates_of_motion, time=time)
No description has been provided for this image

So far, our hot spot's initial location has been at present day. But we can create a motion path using a reconstructed hot spot location.

To demonstrate this we'll take our hot spot initial location and reconstruct it to an arbitrary past time (10 Ma). Then we'll create a motion path using that reconstructed hot spot location. We'll also specify that past time (10 Ma) using the from_time parameter in PlateReconstruction.create_motion_path. Finally we'll specify the reconstruction time (20 Ma) using the to_time parameter.

This should then match the above plot since it's essentially achieving the same result (which is reconstructing the same motion path to 20 Ma).

On a side note, if we don't specify the motion path plate ID (in PlateReconstruction.create_motion_path) then PlateReconstruction will assign it for us. In our case, however, we need to explicitly specify the plate ID since it's a hot spot (plate ID 2) that does not have an associated plate in the PlateReconstruction's static polygons (used to assign plate IDs).

In [9]:
# Use a 'Points' object containing the 'present day' hot spot location.
gpts = gplately.Points(model, lon, lat, time=0, plate_id=moving_plate_ID)

initial_hot_spot_time = 10

# Reconstruct the 'present day' hot spot location to its position at 'initial_hot_spot_time'.
initial_lons, initial_lats = gpts.reconstruct(initial_hot_spot_time, return_array=True)

time = 20

# Create the motion path using our reconstructed hot spot location at 'initial_hot_spot_time'.
# Note, however, that the motion path is reconstructed to 'time'.
rlons, rlats, rtimes, rates_of_motion = model.create_motion_path(
    initial_lons,
    initial_lats,
    time_array,
    # Specify the initial time of the hot spot...
    from_time=initial_hot_spot_time,
    # Specify the reconstruction time...
    to_time=time,
    # Note that we need to specify the hot spot plate ID (it cannot be assigned using static polygons)...
    plate_id=moving_plate_ID,
    # Note that we specify 'relative_plate_id' and leave 'anchor_plate_id' as default (of 'model')...
    relative_plate_id=relative_plate_ID,
    return_times=True,
    return_rate_of_motion=True)

plot_motion_path(rlons, rlats, rtimes, rates_of_motion, time=time)
No description has been provided for this image

2. Flowlines from a mid-ocean ridge in the South Atlantic¶

A flowline is a series of point locations reconstructed through time to track the motion of tectonic plates from spreading features like mid-ocean ridges.

We can generate one or more flowlines with create_flowline from GPlately's PlateReconstruction object. All we need is an array of latitudes and longitudes of points along a mid-ocean ridge, the IDs of the plate to the left and right of the mid-ocean ridge, and an array of reconstruction times.

In [10]:
# Longitudes and latitudes of points along the South Atlantic mid ocean ridge
lons = np.array([-15.4795, -14.9688, -14.0632, -17.7739])
lats = np.array([-1.60584,-11.9764, -25.7209, -35.7228])

# Left and right plate IDs
left_plate_ID = 201
right_plate_ID = 701

# Constructing the time array
time_step = 5.
max_reconstruction_time = 120
time_array = np.arange(0, max_reconstruction_time+time_step, time_step)

# Generate the latitudes and longitude coordinates of the flowlines to the left and right of the MOR!
left_rlons, left_rlats, right_rlons, right_rlats, rtimes, rates_of_motion = model.create_flowline(
    lons,
    lats,
    time_array,
    left_plate_ID,
    right_plate_ID,
    return_times=True,
    return_rate_of_motion=True,
)

Let's illustrate these flowlines with mid-ocean ridge/transform boundaries from Seton et al. (2012) using the PlotTopologies object.

In [11]:
def plot_flowline(left_rlons, left_rlats, right_rlons, right_rlats, rtimes, rates_of_motion, time=0):
    # --- Create figure
    fig = plt.figure(figsize=(12,5), dpi=200)
    plt.suptitle("Flowlines of seed points along a South-Atlantic mid-ocean ridge (at {} Ma)".format(time), weight="demi")
    gs = fig.add_gridspec(1, 7)
    
    ax = fig.add_subplot(gs[0, :4], projection=ccrs.Mercator())
    ax.set_title("Flowlines", fontsize=14)
    
    # --- Plot tick labels. This doesn't work for all projections in cartopy yet...
    ax.set_xticks(np.arange(-180, 190, 10), crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(-90, 95, 5), crs=ccrs.PlateCarree())
    ax.set_xticks(np.arange(-180, 185, 5), minor=True, crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(-90, 90, 2.5), minor=True, crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=True)
    lat_formatter = LatitudeFormatter()
    ax.xaxis.tick_bottom()  # labelled ticks on top
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    ax.xaxis.set_ticks_position('both')   # ticks on both top/bottom
    ax.yaxis.set_ticks_position('both')   # ticks on both left/right
    
    # --- Limit map extent
    lon_min = -60
    lon_max = 30
    lat_min = -40.
    lat_max = 10
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    
    # ---- Plot bathymetry
    etopo_tif_img.reconstruct(time=time, fill_value='darkblue').imshow(alpha=0.3, zorder=1)
    
    gplot.time = time
    
    # --- Plot filled coastlines
    gplot.plot_coastlines(ax, facecolor="white", edgecolor=None, zorder=2)
    
    gplot.plot_ridges(ax, color="r", linewidth=1)
    gplot.plot_transforms(ax, color="r", linewidth=1)
    
    # --- Make sure flowline longitudes are wrapped correctly to the dateline
    #left_rlons = gplately.tools.correct_longitudes_for_dateline(left_rlons)
    #right_rlons = gplately.tools.correct_longitudes_for_dateline(right_rlons)
    
    # The mapping of 'rtimes' to colours should not change.
    # Ie, each time should map to the same colour (regardless of the time range).
    rtimes_min, rtimes_max = rtimes[0], rtimes[-1]
    cmap = ListedColormap(plt.cm.inferno(np.linspace(rtimes_min/rtimes_max, 1.0, 256)))
    
    # --- Iterate over the reconstructed flowlines. Each seed point results in a 'left' and 'right' flowline 
    for i in range(len(lons)):
         
        # plot left flowlines as circles
        ax.plot(left_rlons[i], left_rlats[i], color="C{}".format(i), linewidth=2, transform=ccrs.PlateCarree(), zorder=3)
        fl = ax.scatter(left_rlons[i], left_rlats[i], s=50, marker='D', c=rtimes, cmap=cmap,
                        transform=ccrs.PlateCarree(), edgecolor="C{}".format(i),
                        vmin=rtimes_min, vmax=rtimes_max, zorder=4)
        
        # plot right flowlines as stars
        ax.plot(right_rlons[i], right_rlats[i], color="C{}".format(i), linewidth=2, transform=ccrs.PlateCarree(), zorder=3)
        ax.scatter(right_rlons[i], right_rlats[i], s=50, marker='s', c=rtimes, cmap=cmap,
                   transform=ccrs.PlateCarree(), edgecolor="C{}".format(i),
                   vmin=rtimes_min, vmax=rtimes_max, zorder=4)
        
    # --- Set a colorbar to help visualise the passage of every Xth Myr.
    divider = make_axes_locatable(ax)
    cax = divider.new_horizontal(size="4%", pad=0.4, axes_class=plt.Axes)
    fig.add_axes(cax)
    cbar = plt.colorbar(fl, cax=cax, orientation="vertical")
    cbar.set_label('Age (Ma)', fontsize=12)
    cbar.ax.minorticks_on()

    # --- Plot times and rates of motion on 2nd subplot
    p = fig.add_subplot(gs[0, 5:], xlabel='Time (Ma)', ylabel='Spreading Rate (cm/yr)',
                        xlim=[rtimes.max(), 0],ylim=[0.0,6.0])
    
    p.set_title("Flowline spreading rates")
    p.tick_params(direction="in", length=5, top=True, right=True)
    p.tick_params(direction="in", which='minor', length=2.5, top=True, right=True)
    p.minorticks_on()
    # For each seed point plot its rate of motion.
    for rate_of_motion in rates_of_motion:
        plt.stairs(rate_of_motion, rtimes)

    plt.show()

    fig.savefig("Flowlines_South_Atlantic_{}Ma.pdf".format(time), bbox_inches='tight')

Plot the flowline at present day.

In [12]:
plot_flowline(left_rlons, left_rlats, right_rlons, right_rlats, rtimes, rates_of_motion)
No description has been provided for this image

Now plot the flowline at a past time.

To do this we specify the reconstruction time using the to_time parameter of PlateReconstruction.create_flowline.

The from_time parameter defaults to present day. This is what we want since the seed point coordinates are present day positions.

In [13]:
time = 55

# Reconstruct the flowline to a past time.
left_rlons, left_rlats, right_rlons, right_rlats, rtimes, rates_of_motion = model.create_flowline(
    lons,
    lats,
    time_array,
    left_plate_ID,
    right_plate_ID,
    # Specify the reconstruction time...
    to_time=time,
    return_times=True,
    return_rate_of_motion=True,
)

plot_flowline(left_rlons, left_rlats, right_rlons, right_rlats, rtimes, rates_of_motion, time=time)
No description has been provided for this image

Now do the same but using a Points object (instead of PlateReconstruction).

To do this we first create a Points object at present day since the seed point coordinates are present day positions.

Then we create a flowline using Points.flowline and specify the reconstruction time with its time parameter.

In [14]:
time = 55

# Create a 'Points' object at present day (since point coordinates are present day positions).
gpts = gplately.Points(model, lons, lats, time=0)

# Reconstruct the flowline to a past time using our 'Points' object.
left_rlons, left_rlats, right_rlons, right_rlats, rtimes, rates_of_motion = gpts.flowline(
    time_array,
    left_plate_ID,
    right_plate_ID,
    # Specify the reconstruction time...
    time=time,
    return_times=True,
    return_rate_of_motion=True,
)

plot_flowline(left_rlons, left_rlats, right_rlons, right_rlats, rtimes, rates_of_motion, time=time)
No description has been provided for this image

So far, our seed locations have been at present day. But we can create a flowline using reconstructed seed locations.

To demonstrate this we've taken our present day seed locations and reconstructed them (as a mid-ocean ridge feature using GPlates) to an arbitrary past time (30 Ma) - these are the initial_lons and initial_lats values. Then we create flowlines using those reconstructed seed locations. We'll also specify that past time (30 Ma) using the from_time parameter in PlateReconstruction.create_flowline. Finally we specify the reconstruction time (55 Ma) using the to_time parameter.

This should then match the above plot since it's essentially achieving the same result (which is reconstructing the same flowline to 55 Ma).

In [15]:
initial_seed_time = 30

# Longitudes and latitudes of the same seed points - but at 'initial_seed_time' (not present day).
initial_lons = np.array([-13.886, -13.859, -13.625, -17.951])
initial_lats = np.array([-1.950, -12.333, -26.143, -35.954])

time = 55

# Create the flowlines using our reconstructed seed locations at 'initial_seed_time'.
# Note, however, that the flowlines are reconstructed to 'time'.
left_rlons, left_rlats, right_rlons, right_rlats, rtimes, rates_of_motion = model.create_flowline(
    initial_lons,
    initial_lats,
    time_array,
    left_plate_ID,
    right_plate_ID,
    # Specify the initial time of the seed points...
    from_time=initial_seed_time,
    # Specify the reconstruction time...
    to_time=time,
    return_times=True,
    return_rate_of_motion=True,
)

plot_flowline(left_rlons, left_rlats, right_rlons, right_rlats, rtimes, rates_of_motion, time=time)
No description has been provided for this image