6 - Working with Rasters¶

Any numpy array can be turned into a gplately.Raster:

raster = gplately.Raster(
    plate_reconstruction=model,
    data=array,
    extent="global",  # equivalent to (-180, 180, -90, 90)
    origin="lower",  # or set extent to (-180, 180, -90, 90)
)

In this notebook, we will:

  • Download rasters from EarthByte's webDAV server
  • Plot rasters
  • Resize and respace rasters
  • Reconstruct rasters
  • Linearly interpolate point data on rasters.

Import all needed packages, and create PlateReconstruction and Plot objects for the Muller et al. (2019) plate tectonic model.

In [1]:
import os

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

Using GPlately's download module, we can easily download a variety of plate tectonic models!

In [ ]:
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()

coastlines = muller2019_model.get_layer('Coastlines')
continents = muller2019_model.get_layer('ContinentalPolygons')
COBs =  muller2019_model.get_layer('COBs')

model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)
gplot = gplately.PlotTopologies(model, coastlines, continents, COBs)

PlateModelManager for rasters¶

Let's use PlateModelManager to download Muller et al. 2019 netCDF age grids.

There is a unique age grid for each millionth year - let's access the 0 Ma age grid by passing time to get_age_grid. It is returned as a gplately Raster object which we call muller_2019_age_grid.

In [ ]:
time = 0  # Ma
muller_2019_age_grid = gplately.Raster(
    data=muller2019_model.get_raster("AgeGrids", time),
    plate_reconstruction=model,
    extent=[-180, 180, -90, 90],
    )

Plotting rasters¶

The muller_2019_age_grid is a Raster object - this object allows us to work with age grids and other rasters. Let's visualise the data with imshow.

In [4]:
fig = plt.figure(figsize=(16, 12))
muller_2019_age_grid.imshow(cmap="YlGnBu")
plt.show()
No description has been provided for this image

Let's plot this netCDF grid along with coastlines, mid-ocean ridges and subduction zones (with teeth) resolved from the Muller et al. 2019 plate model.

In [5]:
fig = plt.figure(figsize=(16, 12))
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))

gplot.time = time

muller_2019_age_grid.imshow(ax1, cmap='YlGnBu', vmin=0, vmax=200)
gplot.plot_coastlines(ax1, edgecolor='k', facecolor='1', alpha=0.1)
gplot.plot_trenches(ax1, color='magenta', zorder=5)
gplot.plot_ridges(ax1, color='b', zorder=5)
gplot.plot_transforms(ax1, color='lightblue', linewidth=0.75)
gplot.plot_subduction_teeth(ax1, color='magenta')
ax1.set_title("{} Ma".format(time))
plt.show()
No description has been provided for this image

Resizing and resampling rasters¶

Let's resize and resample the the present-day Muller et al. 2019 agegrid. We can do this using resize.

  • resize - provide the number of points in each direction to resize the raster (e.g. 100 cols by 200 rows)
  • resample - provide the grid spacing in each direction (e.g. 0.1 degrees by 0.2 degrees)
In [6]:
model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)

muller_2019_age_grid = gplately.Raster(
    data=muller2019_model.get_raster("AgeGrids", time),
    plate_reconstruction=model,
    extent=[-180, 180, -90, 90],
    )

# Set grid size in x and y directions
muller_2019_age_grid.resize(20, 10, inplace=True)

# Plot resampled and resized age grid
fig = plt.figure(figsize=(16,12))
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude = 20))
muller_2019_age_grid.imshow(ax1, cmap='YlGnBu', vmin=0, vmax=200)
gplot.plot_coastlines(ax1, edgecolor='k', facecolor='1', alpha=0.1)
gplot.plot_trenches(ax1, color='magenta', zorder=5)
gplot.plot_ridges(ax1, color='b', zorder=5)
gplot.plot_subduction_teeth(ax1, color='magenta')
ax1.set_title("{} Ma".format(time))
plt.show()

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-0.nc
The local file(s) is/are still good. Will not download again at this moment.
No description has been provided for this image

Downloading general rasters¶

Let's visualise an ETOPO1 relief raster. GPlately also makes it easy to import, using get_raster

Alternatively, you can import a local netcdf file using by passing the filename to gplately.grids.Raster or gplately.grids.read_netcdf_grid

In [7]:
# This returns ETOPO1 as a `Raster` object.
from matplotlib import image
raster_manager = PresentDayRasterManager()
etopo = gplately.Raster(data=image.imread(raster_manager.get_raster("ETOPO1_tif")))
etopo.lats = etopo.lats[::-1]

print(np.shape(etopo))

downloading https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/image/color_etopo1_ice_low.tif.gz
The local file(s) is/are still good. Will not download again at this moment.
(2700, 5400, 3)

etopo is a large (5400 x 2700 pixels) RGB image, with a total of 14,580,000 grid points. We can visualise it using a number of methods, including Raster.imshow.

In [8]:
fig = plt.figure(figsize=(16, 12))
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))
etopo.imshow(ax=ax1, interpolation="none")
plt.show()
No description has been provided for this image

We can also resize this RGB raster using the Raster.resize and Raster.resample methods. Here we resample it to a 0.5 by 0.5 degree resolution, then plot it using Raster.imshow.

In [9]:
etopo_downscaled = etopo.resample(0.5,0.5)
print(etopo_downscaled.shape)

fig = plt.figure(figsize=(16, 12))
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))
etopo_downscaled.imshow(ax1, interpolation="none")

(361, 721, 3)
Out[9]:
<matplotlib.image.AxesImage at 0x1672c0890>
No description has been provided for this image

Reconstructing rasters¶

The ETOPO1 raster can be reconstructed back in time by assigning a plate_reconstruction to the object. We use the plate reconstruction we defined earlier (model). Note that a Raster object will also accept a plate_reconstruction at initialisation.

After this, we use reconstruct to reconstruct the raster to a given time.

In [10]:
# assign a plate reconstruction in order to reconstruct the raster
etopo.plate_reconstruction = model

white_rgb = (255, 255, 255)  # RGB code for white, to fill gaps in output
etopo_reconstructed = etopo.reconstruct(50, threads=4, fill_value=white_rgb)

fig = plt.figure(figsize=(16, 12))
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=180))
etopo_reconstructed.imshow(ax1)
Out[10]:
<matplotlib.image.AxesImage at 0x28b22c850>
No description has been provided for this image

Raster can be reconstructed in-place (inplace=True), and fill_value can be set to any valid matplotlib colour when reconstructing RGB images.

In [11]:
# create a duplicate of the obj
etopo_dup = etopo.copy()

etopo_dup.reconstruct(time=75, threads=4, fill_value="darkblue", inplace=True)
etopo_dup.imshow(projection=ccrs.Robinson())
Out[11]:
<matplotlib.image.AxesImage at 0x28b236050>
No description has been provided for this image

By default, Raster.reconstruct uses self.plate_reconstruction.static_polygons to assign plate IDs to grid points. To override this behaviour, pass any collection of pygplates.Feature (e.g. list, pygplates.FeatureCollection, etc) to the partitioning_features argument.

In [12]:
etopo_reconstructed = etopo.reconstruct(140, partitioning_features=continents, threads=4, fill_value="grey")
etopo_reconstructed.imshow(projection=ccrs.Orthographic(0, -80))
plt.gca().set_title("Reconstructed to 140 Ma")
Out[12]:
Text(0.5, 1.0, 'Reconstructed to 140 Ma')
No description has been provided for this image

Reverse reconstructions¶

Rasters can be also be reverse reconstructed forward in time!

In [13]:
reconstructed = etopo.reconstruct(50, fill_value="white", threads=4)

reverse_reconstructed = reconstructed.reconstruct(0, fill_value="white", threads=4)

# plot
fig, axs = plt.subplots(3, 1, figsize=(8, 12), subplot_kw={"projection": ccrs.Mollweide(central_longitude=0)})
etopo.imshow(ax=axs[0])
reconstructed.imshow(ax=axs[1])
reverse_reconstructed.imshow(ax=axs[2])

axs[0].set_title("Original image (present day)")
axs[1].set_title("Reconstructed to 50 Ma")
axs[2].set_title("Reverse reconstructed back to present day")
plt.show()
No description has been provided for this image

Reconstructing netCDF rasters¶

Similar to above, we can also reconstruct numeric netcdf grids! GPlately also includes ETOPO1 as a netcdf for download.

In [14]:
# ETOPO1 downloaded as a Raster object
raster_manager = PresentDayRasterManager()
etopo_nc = gplately.Raster(data=raster_manager.get_raster("ETOPO1_grd"))
etopo_nc._data = etopo_nc._data.astype(float)

# resample to a more manageable size
etopo_nc.resample(0.5, 0.5, inplace=True)
print(etopo_nc.shape)

# Assign plate reconstruction
etopo_nc.plate_reconstruction = model

# Reconstruct raster to 50 Ma
etopo_nc_reconstructed = etopo_nc.reconstruct(50, threads=4)

# plot
fig = plt.figure(figsize=(16, 12))
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))
im1 = etopo_nc_reconstructed.imshow(ax1)
plt.colorbar(im1)
plt.show()

downloading https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/grid_registered/netcdf/ETOPO1_Ice_g_gmt4.grd.gz
The local file(s) is/are still good. Will not download again at this moment.

RuntimeWarning: Invalid values encountered in lat; clipping to [-90, 90]

(361, 721)
No description has been provided for this image

The reconstructed netcdf can also be saved.

In [15]:
# Save the reconstructed ETOPO grid to netCDF 
# (note: this ETOPO netCDF grid is high-res, ~190+ MB)
save_filename = os.path.join(
    "NotebookFiles",
    "reconstructed_etopo.nc",
)

etopo_nc_reconstructed.save_to_netcdf4(save_filename)

Linear interpolation¶

As an example of linear interpolation, let's find which subduction zones form continental arcs. For this we use continental netcdf4 files and sample these at points projected 250 km from the trench (in the direction of the subducting plate) to determine if they are inside a continent.

In [16]:
# Define the continental grid
continental_grid_filename = os.path.join(
    "NotebookFiles",
    "continental_grid_0.nc",
)
continental_raster = gplately.Raster(
    continental_grid_filename,
    model,
    extent="global",
)

# tessellate trenches and extract subduction polarity angle, and the lat-lon coordinates
trench_data = model.tessellate_subduction_zones(time, ignore_warnings=True)
trench_normal_azimuthal_angle = trench_data[:, 7]
trench_pt_lon = trench_data[:, 0]
trench_pt_lat = trench_data[:, 1]
    
# calculate 250 km arc distance
arc_distance = 250 / (gplately.tools.geocentric_radius(trench_pt_lat) / 1e3)
    
# Lat and lon coordinates of all trench points after being projected out 250 km in the direction of subduction.
dlon = arc_distance * np.sin(np.radians(trench_normal_azimuthal_angle))
dlat = arc_distance * np.cos(np.radians(trench_normal_azimuthal_angle))
ilon = trench_pt_lon + np.degrees(dlon)
ilat = trench_pt_lat + np.degrees(dlat)

Let's use GPlately to linearly interpolate these projected trench points onto the continental grids using Raster.interpolate

In [17]:
# Use GPlately to interpolate these new points on the defined grid
sampled_points = continental_raster.interpolate(ilon, ilat, method='linear')

# sampled_point[0] is a list of points in the grid (ascribed the integer 1). Collect their indices. 
in_raster_indices = sampled_points > 0

# Get the lat-lon coordinates of the in_raster points
lat_in = ilat[in_raster_indices]
lon_in = ilon[in_raster_indices]

Plot the in-raster points along with the raster, trenches, and coastlines.

In [18]:
fig = plt.figure(figsize=(16, 12), dpi=100)
ax1 = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))
gplot.time = time
gplot.plot_grid_from_netCDF(ax1, continental_grid_filename, cmap="twilight", alpha=0.5, vmin=0, vmax=200)
gplot.plot_coastlines(ax1, edgecolor='k', facecolor='1', alpha=0.1)
gplot.plot_trenches(ax1, color='r', zorder=5)

# Plot the trench points in-continent
ax1.plot(
    lon_in, lat_in,
    linestyle="none",
    marker="o",
    markersize=0.25,
    markerfacecolor="blue",
    markeredgecolor="blue",
    transform=ccrs.PlateCarree(),
)
plt.title("%i Ma" %time)

# Custom legend
from matplotlib.lines import Line2D

handles = [
    Line2D([0], [0], color="red"),
    Line2D([0], [0], color="blue"),
]
labels = [
    "Trenches",
    "Arc segments in continental grids",
]

plt.legend(handles, labels, loc="lower left")
plt.show()
No description has been provided for this image
In [ ]: