1 - Getting Started

Welcome to GPlately!

GPlately uses object-oriented programming to make life simple. In this notebook we will explore some of the main objects you will use:

• PlateReconstruction - reconstruct features, tessellate mid ocean ridges, subduction zones
• Points - partition points onto plates, rotate back through time
• Raster - read in NetCDF grids, interpolation, resampling.
• PlotTopologies - plotting topologies e.g. ridges, trenches, subduction teeth on maps
• PlateModelManager - downloading plate models, features, and rasters e.g. .rot, .gpml, .shp and .nc files

import cartopy.crs as ccrs
import gplately
import matplotlib.pyplot as plt
import numpy as np
from plate_model_manager import PlateModelManager

Tectonic plate reconstructions

We simply supply a rotation model, plate topologies, and static polygons to initialise a plate reconstruction model. You can download these files into your machine's cache using GPlately's PlateModelManager object.

# Call GPlately's PlateModelManager object and request data from the Müller et al. 2019 study
pm_manager = PlateModelManager()
muller2019_model = pm_manager.get_model("Muller2019", data_dir="plate-model-repo")
rotation_model = muller2019_model.get_rotation_model()
topology_features = muller2019_model.get_topologies()
static_polygons = muller2019_model.get_static_polygons()

# Tessellate the subduction zones to 0.5 degrees.
tessellation_threshold_radians = np.radians(0.05)

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

downloading https://repo.gplates.org/webdav/pmm/muller2019/Rotations.zip
downloading https://repo.gplates.org/webdav/pmm/muller2019/Topologies.zip
downloading https://repo.gplates.org/webdav/pmm/muller2019/StaticPolygons.zip

Now let's find the subduction zones and mid-ocean ridges at 10 Ma.

time = 10

# these bundle a lot of information - check PTT docs for more info
subduction_data = model.tessellate_subduction_zones(time, ignore_warnings=True)
ridge_data = model.tessellate_mid_ocean_ridges(time, ignore_warnings=True)

Mapping

The PlotTopologies function injests the plate model we have defined as well as the coastlines, continents, and COB. It computes all of the plate topologies for a given reconstruction time.

This object has been designed to work specifically with cartopy. Define your figure and supply your axes to these plotting routines. Some common favourites include:

• coastlines
• continents
• ridges and transforms
• trenches
• subduction teeth (!!)
• netCDF grids
• plate motion vectors

You can still supply optional keywords as you normally would.

# Obtain features for the PlotTopologies object with PlateModelManager
coastlines = muller2019_model.get_layer('Coastlines')
continents = muller2019_model.get_layer('ContinentalPolygons')
COBs =  muller2019_model.get_layer('COBs')

# Call the PlotTopologies object
gplot = gplately.plot.PlotTopologies(model, coastlines=coastlines, continents=continents, COBs=COBs)

# Download all Muller et al. 2019 netCDF age grids with PlateModelManager. This is returned as a Raster object.
agegrid = gplately.Raster(data=muller2019_model.get_raster("AgeGrids",time))

downloading https://repo.gplates.org/webdav/pmm/muller2019/Coastlines.zip
downloading https://repo.gplates.org/webdav/pmm/muller2019/ContinentalPolygons.zip
downloading https://repo.gplates.org/webdav/pmm/muller2019/COBs.zip
downloading https://www.earthbyte.org/webdav/ftp/Data_Collections/Muller_etal_2019_Tectonics/Muller_etal_2019_Agegrids/Muller_etal_2019_Tectonics_v2.0_netCDF/Muller_etal_2019_Tectonics_v2.0_AgeGrid-10.nc

Set the time attribute to reconstruct all topologies to the specified time.

IMPORTANT: You must set gplot.time or provide a time at initialisation before plotting anything.

gplot.time = 10 # Ma

Create a map with some useful geological information

fig = plt.figure(figsize=(16,12))

ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(190))

gplot.plot_continents(ax1, facecolor='0.8')
gplot.plot_coastlines(ax1, color='0.5')
gplot.plot_ridges_and_transforms(ax1, color='red')
gplot.plot_trenches(ax1, color='k')
gplot.plot_subduction_teeth(ax1, color='k')
im = gplot.plot_grid(ax1, agegrid.data, cmap='YlGnBu', vmin=0, vmax=200)
gplot.plot_plate_motion_vectors(ax1, spacingX=10, spacingY=10, normalise=True, zorder=10, alpha=0.5)

fig.colorbar(im, orientation='horizontal', shrink=0.4, pad=0.05, label='Age (Ma)')

<matplotlib.colorbar.Colorbar at 0x13d99e150>

# update the time to regenerate topologies
time = 100
gplot.time = time
agegrid = gplately.Raster(data=muller2019_model.get_raster("AgeGrids",time))

fig = plt.figure(figsize=(16,12))

ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(190))

gplot.plot_continents(ax1, facecolor='0.8')
gplot.plot_coastlines(ax1, color='0.5')
gplot.plot_ridges_and_transforms(ax1, color='red')
gplot.plot_trenches(ax1, color='k')
gplot.plot_subduction_teeth(ax1, color='k')
im = gplot.plot_grid(ax1, agegrid.data, cmap='YlGnBu', vmin=0, vmax=200)
gplot.plot_plate_motion_vectors(ax1, spacingX=10, spacingY=10, normalise=True, zorder=10, alpha=0.5)

fig.colorbar(im, orientation='horizontal', shrink=0.4, pad=0.05, label='Age (Ma)')

downloading https://www.earthbyte.org/webdav/ftp/Data_Collections/Muller_etal_2019_Tectonics/Muller_etal_2019_Agegrids/Muller_etal_2019_Tectonics_v2.0_netCDF/Muller_etal_2019_Tectonics_v2.0_AgeGrid-100.nc

<matplotlib.colorbar.Colorbar at 0x1485ee750>

Working with points

Now that we have defined our reconstruction object, we can reconstruct point data.

pt_lons = np.array([140., 150., 160.])
pt_lats = np.array([-30., -40., -50.])

gpts = gplately.Points(model, pt_lons, pt_lats)

vel_x, vel_y = gpts.plate_velocity(0)
vel_mag = np.hypot(vel_x, vel_y)

print("average point velocity (cm/yr)", vel_mag)

average point velocity (cm/yr) [6.33037749 5.61546482 4.93440741]

Plot their position from time=0 to time=20

rlons = np.empty((21, pt_lons.size))
rlats = np.empty((21, pt_lons.size))


for time in range(0, 21):
    rlons[time], rlats[time] = gpts.reconstruct(time, return_array=True)

gplot.time = 0 # present day

fig = plt.figure(figsize=(16,12))
ax1 = fig.add_subplot(111, projection=ccrs.Mercator(190)) 
ax1.set_extent([130,180,-60,-10])

gplot.plot_coastlines(ax1, color='0.8')

for i in range(0, len(pt_lons)):
    ax1.plot(rlons[:,i], rlats[:,i], 'o', transform=ccrs.PlateCarree())

Rasters

You can initialise a Raster object by providing a raster (either a NetCDF file path or a regular numpy array), and optionally passing a PlateReconstruction model. The data attribute stores the raster data as a 2D numpy array.

In this example, we will pass a NetCDF file path to the Raster object.

time = 100

# download a netcdf age grid Raster and assign a plate reconstruction `model'
graster = gplately.Raster(data=muller2019_model.get_raster("AgeGrids",time))
graster.plate_reconstruction = model

# the underlying numpy (masked) array can be accessed using the data attribute
print(type(graster.data))

# alternatively initialise a Raster from a numpy array
graster = gplately.Raster(data=graster.data, # 2D numpy array
                          plate_reconstruction=model, # PlateReconstruction object
                          extent='global', # equivalent to [-180, 180, -90, 90]
                          time=100, # time in Ma
                         )

downloading https://www.earthbyte.org/webdav/ftp/Data_Collections/Muller_etal_2019_Tectonics/Muller_etal_2019_Agegrids/Muller_etal_2019_Tectonics_v2.0_netCDF/Muller_etal_2019_Tectonics_v2.0_AgeGrid-100.nc
The local file(s) is/are still good. Will not download again at this moment.
<class 'numpy.ma.core.MaskedArray'>

gplot.time = time

fig = plt.figure(figsize=(16,12))

ax1 = fig.add_subplot(111, projection=ccrs.Robinson())
gplot.plot_grid(ax1, graster, cmap='YlGnBu', vmin=0, vmax=100)
# Alternatively:
# graster.imshow(ax=ax1, cmap="YlGnBu", vmin=0, vmax=100)
gplot.plot_coastlines(ax1, edgecolor='k', facecolor='none')

<GeoAxes: >

There are a bunch of routines such as,

• filling masked (NaN) regions
• interpolation
• resampling
• reconstructing rasters

In-place operations can be achieved using inplace=True which will update the internal data structures.