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
import cartopy.feature as cfeature
from mpl_toolkits.axes_grid1 import make_axes_locatable
from shapely.geometry import LineString
import geopandas as gpd
%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
gdownload = gplately.download.DataServer("Matthews2016")
rotation_model, topology_features, static_polygons = gdownload.get_plate_reconstruction_files()
# Call the PlateReconstruction object to create a plate motion model
model = gplately.reconstruction.PlateReconstruction(rotation_model, topology_features, static_polygons)
# Attempt to fetch Matthews et al. 2016 coastlines, continents and COBs from EarthByte servers
coastlines, continents, COBs = gdownload.get_topology_geometries()
# Call the PlotTopologies object
time = 0
gplot = gplately.PlotTopologies(model, coastlines=coastlines, continents=continents, COBs=COBs, time=time)
Checking whether the requested files need to be updated... Requested files are up-to-date! Checking whether the requested files need to be updated... Requested files are up-to-date! No continent-ocean boundaries in Matthews2016.
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_rate_of_motion
to True
in create_motion_path
, the function returns two additional arrays: one is a timestep array, and another is a rate of plate motion array.
# Latitude and longitude array (can have >1 lats and lons)
lats = np.array([19])
lons = np.array([-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
lon, lat, step_times, rates_of_motion = model.create_motion_path(
lons, lats, time_array, anchor_plate_ID, moving_plate_ID, 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:
etopo_tif_img = gdownload.get_raster("ETOPO1_tif")
Checking whether the requested files need to be updated... Requested files are up-to-date!
# Creating the motion path and obtaining rate plot arrays using the gplately Points alternative
gpts = gplately.reconstruction.Points(model, lons, lats, time=0.)
lon, lat, step_times, rates_of_motion = gpts.motion_path(
time_array,
anchor_plate_ID,
return_rate_of_motion=True
)
# --- Create figure
fig = plt.figure(figsize=(12, 3), dpi=300)
plt.rcParams['font.family'] = 'Arial'
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 = 130.
lon_max = 220 # if I specify as -150, it will plot as the 'other' way around the world
lat_min = 15.
lat_max = 60.
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
ax.add_feature(cfeature.LAND, color='white', zorder=2)
# --- Make sure motion path longitudes are wrapped correctly to the dateline
lon360 = gplately.tools.correct_longitudes_for_dateline(lon)
# --- plot motion path
ax.plot(lon360, lat, color='yellow', linewidth=0.75, transform=ccrs.PlateCarree(), zorder=3)
mp = ax.scatter(lon360, lat, 100, marker='.', c=time_array, cmap=plt.cm.inferno, edgecolor='k',
transform=ccrs.PlateCarree(), vmin=time_array[0], vmax=time_array[-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 step times and step rates of motion on 2nd subplot
p = fig.add_subplot(122, xlabel='Time (Ma)', ylabel='Rate of motion (cm/yr)',
xlim=[time_array.max(), 0],ylim=[3.5,12])
p.set_title("Rate of Pacific Plate motion (cm/yr)")
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.plot(step_times, rates_of_motion, color='k')
plt.subplots_adjust(wspace=0.25)
fig.savefig("Hawaii_Emperor_motion_path.pdf", bbox_inches='tight')