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 sourced from EarthByte's GPlates 2.3 Dataset library: https://www.earthbyte.org/gplates-2-3-software-and-data-sets/ and loaded in with GPlately's DataServer object.)

Let's set up all our packages:

import gplately

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

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/).

# Initialise gplately data server and extract data from Muller et al. 2019
gdownload = gplately.download.DataServer("Muller2019")
rotation_model, topology_features, static_polygons = gdownload.get_plate_reconstruction_files()

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

Checking whether the requested files need to be updated...
Requested files are up-to-date!

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

# Obtain Muller et al. 2019 geometries with gdownload
coastlines, continents, COBs = gdownload.get_topology_geometries()

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

Checking whether the requested files need to be updated...
Requested files are up-to-date!

Plotting shapely geometries using GPlately's PlotTopologies object

We can use get_feature_data on GPlately's DataServer object 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 GPlately's DataServer.

**(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.)

# Set present day
time = 0 # Ma
seafloor_fabric = gdownload.get_feature_data("SeafloorFabric")

Checking whether the requested files need to be updated...
Requested files are up-to-date!

DataServer 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.

# 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])
    

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.

# 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

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...