5 - Working With Feature Geometries¶

GPlately can be used to read, process and plot geological features contained in ESRI Shapefiles (with file extension .shp) and GPML files (with file extension .gpml or .gpmlz). In this notebook, we will look at:

  • Plotting shapely geometries using GPlately's PlotTopologies object
  • Plotting extra .gpml features such as:
    • Polylines
    • Polygons
    • Points

(These GPML files are downloaded with PlateModelManager.)

Let's set up all our packages:

In [1]:
import gplately

import numpy as np
import gplately.pygplates as pygplates
import glob, os
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from plate_model_manager import PlateModelManager

GPlately's PlotTopologies object uses the PlateReconstruction object to reconstruct geological features. It then turns the reconstructed features into Shapely MultiPolygon, MultiPoint and/or MultiLine geometries and plots them onto GeoAxis maps. To call PlotTopologies, we need to pass:

  • the PlateReconstruction plate motion model we just created
  • a specific reconstruction time (Ma),
  • a coastline filename or <pygplates.FeatureCollection> object,
  • a continent filename or <pygplates.FeatureCollection> object,
  • and a continent-ocean boundary (COBs) filename or <pygplates.FeatureCollection> object,

We'll first construct the PlateReconstruction object, which needs a rotation model, a topology_feature collection and static_polygons. Let's use GPlately's DataServer to download these files from Müller et al. (2019) (https://www.earthbyte.org/muller-et-al-2019-deforming-plate-reconstruction-and-seafloor-age-grids-tectonics/).

In [ ]:
# Download Muller et al. 2019 model data
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()

# Create the plate motion model!
model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)

We can also download the coastlines, continents and COBs needed for our PlotTopologies object using DataServer.

In [ ]:
# Obtain Muller et al. 2019 geometries 
coastlines = muller2019_model.get_layer('Coastlines')
continents = muller2019_model.get_layer('ContinentalPolygons')
COBs =  muller2019_model.get_layer('COBs')

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

Plotting shapely geometries using GPlately's PlotTopologies object¶

We can use PlateModelManager to download assorted feature data to reconstruct through geological time.

Plotting Shapely Polylines¶

Let's visualise a set of polylines from Matthews et al. 2011** that define the tectonic fabric of global seafloors at present day digitised from vertical gravity gradient (VGG) maps. These polylines represent global fracture zones (FZs), v-shaped structures (VANOMs), discordant zones (DZs), fracture zones (hand-traced with less certainty) (FZLCs), unclassified V-anomalies (UNCVs) and extinct ridges.

Let's load these features in with PlateModelManager.

**(Matthews, K. J., Müller, R. D., Wessel, P., Whittaker, J. M. 2011. The tectonic fabric of the ocean basins, The Journal of Geophysical Research. Doi: 10.1029/2011JB008413.)

In [ ]:
# Set present day
time = 0 # Ma
seafloor_fabric = muller2019_model.get_layer("SeafloorFabric")

PlateModelManager returned the seafloor fabric as a list of pygplates feature collections, one for each type of fabric. To reconstruct and plot them, we loop through each Feature and:

  • Reconstruct them using the PlateReconstruction object
  • Turn the reconstructed feature into a set of Shapely polylines using the shapelify_feature_lines method in the GPlately plot module.

We want to identify each feature by colour, so these steps are done in a for loop. Each seafloor feature is plotted onto a Cartopy GeoAxis map.

In [5]:
# Set up GeoAxis and plot shapefile topologies
ax2 = plt.figure(figsize=(16,12)).add_subplot(111, projection=ccrs.Robinson(central_longitude=10))
gplot.time = time
gplot.plot_continents(ax2, facecolor='0.8')
gplot.plot_continent_ocean_boundaries(ax2, color='0.98')
gplot.plot_coastlines(ax2, color='0.9')
gplot.plot_ridges_and_transforms(ax2, color='r')
gplot.plot_trenches(ax2, color='navy')
gplot.plot_subduction_teeth(ax2, color='navy')
plt.title('Global seafloor fabric at %i Ma' % (time))

# Seafloor fabric topology identification variables
colours = ['powderblue', 'k', 'm', 'g', 'b', 'y']

# Loop through all seafloor fabric filenames, reconstruct each topology and plot onto ax2 using GPlately
for i, fabric in enumerate(seafloor_fabric):
    reconstructed_seafloor_topology = model.reconstruct(fabric, time)
    polylines = gplately.plot.shapelify_feature_lines(reconstructed_seafloor_topology)
    ax2.add_geometries(polylines, crs=ccrs.PlateCarree(), facecolor=colours[i], edgecolor=colours[i])
No description has been provided for this image

If you have moviepy available, you can create a gif that illustrates the generation of tectonic seafloor fabric over many Ma. Let's reconstruct plate movements over 150 Ma in intervals of 10 Ma.

In [6]:
# Seafloor fabric topology identification variables
colours = ['powderblue', 'k', 'm', 'g', 'b', 'y']
feat = ['Fracture Zones', 'V-Shaped Structures', 'Discordant Zones', 'Fracture Zones (less certainty)', 
        'Unclassified V-Anomalies', 'Extinct Ridges']

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

# Create a plot for each 10 Ma interval
for time in np.arange(oldest_seed_time,0.,-time_step):
    
    # Set up GeoAxis and plot shapefile topologies
    ax5 = plt.figure(figsize=(16,12), dpi=200).add_subplot(111, projection=ccrs.Robinson(central_longitude=10))
    gplot.time = time
    gplot.plot_continents(ax5, facecolor='0.8')
    gplot.plot_continent_ocean_boundaries(ax5, color='0.98')
    gplot.plot_coastlines(ax5, color='0.9')
    gplot.plot_ridges_and_transforms(ax5, color='r')
    gplot.plot_trenches(ax5, color='navy')
    gplot.plot_subduction_teeth(ax5, color='navy')
    plt.title('Global seafloor fabric at %i Ma' % (time))
    
    # Loop through all seafloor fabric filenames, reconstruct each topology and plot onto ax2 using GPlately
    for i, fabric in enumerate(seafloor_fabric):
        reconstructed_seafloor_topology = model.reconstruct(fabric, time)
        polylines = gplately.plot.shapelify_feature_lines(reconstructed_seafloor_topology)
        ax5.add_geometries(polylines, crs=ccrs.PlateCarree(), facecolor=colours[i], edgecolor=colours[i])
        ax5.set_global()
        
    plt.savefig('/tmp/seafloor_fabric_%d_Ma.png' % time)
    plt.close()
    print('Image for %d Ma saved' % time)
    

Image for 150 Ma saved
Image for 140 Ma saved
Image for 130 Ma saved
Image for 120 Ma saved
Image for 110 Ma saved
Image for 100 Ma saved
Image for 90 Ma saved
Image for 80 Ma saved
Image for 70 Ma saved
Image for 60 Ma saved
Image for 50 Ma saved
Image for 40 Ma saved
Image for 30 Ma saved
Image for 20 Ma saved
Image for 10 Ma saved
In [7]:
import moviepy.editor as mpy

frame_list = []

for time in np.arange(oldest_seed_time,0.,-time_step):
    frame_list.append('/tmp/seafloor_fabric_%d_Ma.png' % time)

clip = mpy.ImageSequenceClip(frame_list, fps=5)
clip.write_gif('/tmp/seafloor_fabric_movie.gif')

from IPython.display import Image
print('The movie will show up in a few seconds. Please be patient...')
with open('/tmp/seafloor_fabric_movie.gif','rb') as f:
    display(Image(data=f.read(), format='png', width = 3000, height = 1000))

MoviePy - Building file /tmp/seafloor_fabric_movie.gif with imageio.

                                                                                                                                                                                                                      

The movie will show up in a few seconds. Please be patient...
No description has been provided for this image

Plotting Shapely Polygons¶

Let's visualise two polygon-feature data sets:

  1. Global volcanic provinces at present day from "The interplay between the eruption and weathering of Large Igneous Provinces and the deep-time carbon cycle" by Johansson et al. (2018)

  2. Global large igneous provinces (LIPs) from the "Long-term interaction between mid-ocean ridges and mantle plumes" by Whittaker et al. (2015).

We'll obtain these files using GPlately's PlateModelManager object:

In [ ]:
# Locate the Johansson et al. (2018) and Whittaker et al. (2015) gpmlz files containing volcanic provinces & LIPs
lip_volcanic_provinces = [muller2019_model.get_layer("Johansson2018LIPs"),
                          muller2019_model.get_layer("Whittaker2015LIPs")]

We now have a list containing the LIP and Volcanic Province topology files. Let's use a for loop to loop through each <pygplates.FeatureCollection> object in the list and:

  • Reconstruct topologies to a specific geological time,
  • Turn topologies into shapely polygons
  • Add shapely geometries onto a GeoAxis map with formatting keyword arguments
In [9]:
# Set time to present day
time = 0

# Set up GeoAxis and plot shapefile topologies to present day
ax3 = plt.figure(figsize=(16,12)).add_subplot(111, projection=ccrs.Robinson(central_longitude=10))
gplot.time = time
gplot.plot_continents(ax3, facecolor='0.8')
gplot.plot_continent_ocean_boundaries(ax3, color='0.98')
gplot.plot_coastlines(ax3, color='0.9')
gplot.plot_ridges_and_transforms(ax3, color='r')
gplot.plot_trenches(ax3, color='navy')
gplot.plot_subduction_teeth(ax3, color='navy')
plt.title('Global volcanic & large igneous provinces at %i Ma' % (time))

# Loop through all seafloor fabric filenames, reconstruct each topology and plot onto ax2 using GPlately
feat = ["Johansson et al. 2018", "Whittaker et al. 2015"]
colours = ['cyan', 'maroon']
for i, topology in enumerate(lip_volcanic_provinces):
    reconstructed_topology = model.reconstruct(topology, time)
    polygons = gplately.plot.shapelify_feature_polygons(reconstructed_topology)
    ax3.add_geometries(polygons, crs=ccrs.PlateCarree(), facecolor=colours[i], edgecolor=colours[i], label=feat[i])
    ax3.set_global()
No description has been provided for this image

Plotting Shapely Points¶

Let's visualise present day surface hotspot/plume locations from the "Long-term interaction between mid-ocean ridges and mantle plumes" by Whittaker et al. (2015). These locations are point data split into Pacific and Indo/Atlantic domains. They were compiled from studies by Montelli et al. (2004), Courtillot et al. 2003, Anderson and Schramm (2005) and Steinberger et al. (2000). Any plume points are separated by 500 km or less have been combined into an average point location.

Let's obtain these topology files with GPlately's PlateModelManager object and reconstruct these point features to present-day.

In [10]:
# Set reconstruction time to present day
time = 0 # Ma
hotspot_plumes = muller2019_model.get_layer("Hotspots")
        
# Reconstruct hotspot and plume point locations to present day.
reconstructed_hotspot_plumes = model.reconstruct(hotspot_plumes, time)

downloading https://repo.gplates.org/webdav/pmm/muller2019/Hotspots.zip
The local file(s) is/are still good. Will not download again at this moment.

Plotting point data requires a matplotlib scatterplot. You'll need two arrays: 1 for the latitudes of the point locations, and 1 for the longitudes of the point locations. We can use pygplates to extract these lat-lon coordinates. Once these arrays are created, we can plot them using "scatter" onto a GeoAxis map!

In [11]:
# Extract lat-lon coordinates of point locations
plon = np.zeros(len(reconstructed_hotspot_plumes))
plat = np.zeros(len(reconstructed_hotspot_plumes))
for i, feature in enumerate(reconstructed_hotspot_plumes):
    geometry = feature.get_reconstructed_geometry()
    plat[i], plon[i] = geometry.to_lat_lon()
    
# Set up GeoAxis, plot shapefile topologies and hotspot/plume point features to present day
ax4 = plt.figure(figsize=(18,10)).add_subplot(111, projection=ccrs.Mollweide(central_longitude=0))
gplot.time = time
gplot.plot_continents(ax4, facecolor='0.8')
gplot.plot_continent_ocean_boundaries(ax4, color='0.98')
gplot.plot_coastlines(ax4, color='0.9')
gplot.plot_ridges_and_transforms(ax4, color='r')
gplot.plot_trenches(ax4, color='navy')
gplot.plot_subduction_teeth(ax4, color='navy')
plt.title('Global surface hotspot & plume locations at %i Ma' % (time))
ax4.scatter(plon, plat, transform=ccrs.PlateCarree(), marker='o', color='greenyellow', edgecolor='k', s=30, zorder=2)
ax4.set_global()
No description has been provided for this image

Dataset sources/citations¶

GPML feature topology data used for this notebook have been sourced from EarthByte's GPlates 2.3 software and dataset database: https://www.earthbyte.org/gplates-2-3-software-and-data-sets/

Global seafloor fabric:

  • Matthews, K.J., Müller, R.D., Wessel, P. and Whittaker, J.M., 2011. The tectonic fabric of the ocean basins. Journal of Geophysical Research, 116(B12): B12109, DOI: 10.1029/2011JB008413.

Volcanic provinces:

  • Johansson, L., Zahirovic, S., and Müller, R. D., 2018, The interplay between the eruption and weathering of Large Igneous Provinces and the deep-time carbon cycle: Geophysical Research Letters, doi:10.1029/2017GL076691.
  • Bryan, S. E., and Ernst, R. E., 2008, Revised definition of large igneous provinces (LIPs): Earth-Science Reviews, v. 86, no. 1, p. 175-202, doi:10.1016/j.earscirev.2007.08.008.
  • Courtillot, V. E., and Renne, P. R., 2003, On the ages of flood basalt events: Comptes Rendus Geoscience, v. 335, no. 1, p. 113-140, doi:10.1016/S1631-0713(03)00006-3.

Large igneous provinces:

  • Whittaker, J., Afonso, J., Masterton, S., Müller, R., Wessel, P., Williams, S., and Seton, M., 2015, Long-term interaction between mid-ocean ridges and mantle plumes: Nature Geoscience, v. 8, no. 6, p. 479-483, doi: 10.1038/ngeo2437.

Surface hotspot + plume locations:

  • Whittaker, J., Afonso, J., Masterton, S., Müller, R., Wessel, P., Williams, S., and Seton, M., 2015, Long-term interaction between mid-ocean ridges and mantle plumes: Nature Geoscience, v. 8, no. 6, p. 479-483, doi: 10.1038/ngeo2437.
In [ ]: