9 - Creating Motion Paths and Flowlines¶
In this notebook, we will show how GPlately's PlateReconstruction
object:
- Creates a motion path of points to illustrate the trajectory of a tectonic plate through geological time, and
- Quantifies the rate of motion of a tectonic plate through time.
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 DataServer object with the string 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.
# Request Matthews2016 from GPlately
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')
# Call the PlateReconstruction object to create a plate motion model
model = gplately.reconstruction.PlateReconstruction(rotation_model, topology_features, static_polygons)
# Call the 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 find the trail of points on the overriding Pacific plate (moving_plate_ID = 901
) relative to a relatively stationary reference frame (anchor_plate_ID = 2
). 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.
# Latitude and longitude
lat = 19
lon = -155
# Plate IDs
moving_plate_ID = 901
anchor_plate_ID = 2
# 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,
anchor_plate_id=anchor_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 DataServer:
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]
# Creating the motion path and obtaining rate plot arrays using the gplately Points alternative
gpts = gplately.Points(model, lon, lat, time=0, plate_id=moving_plate_ID)
rlons, rlats, rtimes, rates_of_motion = gpts.motion_path(
time_array,
anchor_plate_id=anchor_plate_ID,
return_times=True,
return_rate_of_motion=True)
# --- Create figure
fig = plt.figure(figsize=(12, 3), dpi=300)
ax = fig.add_subplot(121, projection=ccrs.Mercator(central_longitude=180))
ax.set_title("Motion path of the Hawaiian-Emperor Chain")
# --- Limit map extent
lon_min = 180+10
lon_max = 180+80
lat_min = -20
lat_max = 30
ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
ax.gridlines(draw_labels=True, xlocs=np.arange(-180,180,20), ylocs=np.arange(-90,90,10))
# ---- Plot bathymetry
etopo_tif_img.imshow(alpha=0.5, zorder=1)
# --- Plot filled coastlines
gplot.plot_coastlines(ax, facecolor="white", edgecolor=None, zorder=2)
# --- Make sure motion path longitudes are wrapped correctly to the dateline
rlons360 = gplately.tools.correct_longitudes_for_dateline(rlons)
# --- plot motion path
ax.plot(rlons360, rlats, color='yellow', linewidth=0.75, transform=ccrs.PlateCarree(), zorder=3)
mp = ax.scatter(rlons360, rlats, 100, marker='.', c=rtimes, cmap=plt.cm.inferno, edgecolor='k',
transform=ccrs.PlateCarree(), vmin=rtimes[0], vmax=rtimes[-1], 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)')
cbar.ax.minorticks_on()
# --- Plot times and rates of motion on 2nd subplot
p = fig.add_subplot(122, xlabel='Time (Ma)', ylabel='Rate of motion (cm/yr)',
xlim=[rtimes.max(), 0],ylim=[rates_of_motion.min()-0.5,rates_of_motion.max()+0.5])
p.set_title("Rate of Pacific Plate 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()
p.stairs(rates_of_motion[0], rtimes, color='k')
plt.subplots_adjust(wspace=0.25)
fig.savefig("Hawaii_Emperor_motion_path.pdf", bbox_inches='tight')
The motion path of several seed points¶
Below, we model the motion of three seed points on the Indian plate (allocated plate ID 501) with respect to the Eurasian plate (allocated anchor plate ID 301). The 501-301 rotation is derived through a plate circuit with the Matthews et al. 2016 model.
Also, we can use 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).
# Latitude and longitude array (can have >1 lats and lons)
lats = np.array([30, 30, 30])
lons = np.array([73, 78, 83])
# Plate ID attributes
# motion of the Indian plate (501) with respect to the Eurasian plate (301)
plate_ID = 501
relative_plate_ID = 301
# Create the time array for the motion path
start_reconstruction_time = 0.
time_step = 5.
max_reconstruction_time = 105.
time_array = np.arange(start_reconstruction_time, max_reconstruction_time + 1, time_step)
# Get the latitudes and longitudes of all points along the motion path
rlons, rlats, rtimes, rates_of_motion = model.create_motion_path(
lons,
lats,
time_array,
plate_id=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)
def plot_motion_path(rlons, rlats, rtimes, rates_of_motion, time=0):
# --- Create figure
fig = plt.figure(figsize=(12,5), dpi=100)
plt.suptitle("Motion of the Indian plate relative to the Eurasian plate (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
# --- Limit map extent
lon_min = 0.
lon_max = 100
lat_min = -40.
lat_max = 40.
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 = ['c0', 'c1', 'c2']
for i in range(len(lons)):
ax.plot(rlons[i], rlats[i], color="C{}".format(i), linewidth=2, transform=ccrs.PlateCarree(), zorder=3)
mp = ax.scatter(rlons[i], rlats[i], 200, marker='*', c=rtimes, cmap=plt.cm.inferno, edgecolor="C{}".format(i),
linewidth=0.5, transform=ccrs.PlateCarree(), vmin=rtimes[0], vmax=rtimes[-1], 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=[rates_of_motion.min()-0.5,rates_of_motion.max()+0.5])
p.set_title("Rates 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()
Plot the motion path at present day.
plot_motion_path(rlons, rlats, rtimes, rates_of_motion)
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 seed point coordinates are present day positions.
time = 55
# Reconstruct the motion path to a past time.
rlons, rlats, rtimes, rates_of_motion = model.create_motion_path(
lons,
lats,
time_array,
# Specify the reconstruction time...
to_time=time,
plate_id=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)
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 motion path using Points.motion_path
and specifying the reconstruction time with its time
parameter.
time = 55
# Create a 'Points' object at present day (since point coordinates are present day positions).
gpts = gplately.Points(model, lons, lats, time=0, plate_id=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)
So far, our seed locations have been at present day. But we can create a motion path using reconstructed seed locations.
To demonstrate this we'll take our present day seed locations and reconstruct them to an arbitrary past time (30 Ma). Then we'll create motion paths using those reconstructed seed locations. We'll also specify that past time (30 Ma) using the from_time
parameter in PlateReconstruction.create_motion_path
. Finally we'll 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 motion path to 55 Ma).
On a side note, if we don't specify the motion path plate ID(s) (in PlateReconstruction.create_motion_path
) then PlateReconstruction
will assign them for us. In our case it will assign plate 501
since the (reconstructed) seed points are on the (reconstructed) Indian plate.
# Use a 'Points' object containing the 'present day' seed locations.
gpts = gplately.Points(model, lons, lats, time=0, plate_id=plate_ID)
initial_seed_time = 30
# Reconstruct the 'present day' seed locations to their position at 'initial_seed_time'.
initial_lons, initial_lats = gpts.reconstruct(initial_seed_time, return_array=True)
time = 55
# Create the motion paths with our reconstructed seed locations at 'initial_seed_time'.
#
# Note: We didn't specify the 'plate_id' parameter, so 'model' will assign
# plate IDs for us by reconstructing its static polygons to 'initial_seed_time' and
# finding which reconstructed static polygons contain the reconstructed seed locations.
# In our case, all three seed points are on the Indian plate (501).
rlons, rlats, rtimes, rates_of_motion = model.create_motion_path(
initial_lons,
initial_lats,
time_array,
# Specify the initial time of the seed points...
from_time=initial_seed_time,
# Specify the reconstruction time...
to_time=time,
assign_plate_id_using_static_polygons=True,
# 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)
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.
# 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.
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=100)
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)
# --- 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,
transform=ccrs.PlateCarree(), cmap=plt.cm.inferno, edgecolor="C{}".format(i),
vmin=rtimes[0], vmax=rtimes[-1], 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,
transform=ccrs.PlateCarree(), cmap=plt.cm.inferno, edgecolor="C{}".format(i),
vmin=rtimes[0], vmax=rtimes[-1], 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=[rates_of_motion.min()-0.5,rates_of_motion.max()+0.5])
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()
Plot the flowline at present day.
plot_flowline(left_rlons, left_rlats, right_rlons, right_rlats, rtimes, rates_of_motion)
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.
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)
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 specifying the reconstruction time with its time
parameter.
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)
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).
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 with our reconstructed seed locations at 'initial_seed_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)