Global Plate Motions

Global Plate Motions

Figure 1 from my SIAM News publication in 2015:

The surface of the Earth is continually moving in response to convection currents in the underlying mantle. A map of the surface velocities shows the motions are organised into a strikingly simple pattern of surface “plates” that are almost rigid bodies. The surface strain rate - coloured contours - is concentrated at the plate boundaries. On closer inspection, there are many regions within the continental crust where the motions are diffuse, not at all plate-like. The arrows are colored by their angle to North as a means of distinguishing the different plates. The scale is such that the longest vector represents 10cm/yr.

The image is composed of a background grayscale topography / bathymetry map, plate motion vectors (in a Europe-fixed frame) and contours of the observed strain rate second invariant. These data are available after running Notebook 0 to download / install everything.

%pylab inline

from osgeo import gdal
from osgeo import gdal_array

import as ccrs
import matplotlib.pyplot as plt

import cartopy.feature as cfeature

from import netcdf
from cloudstor import cloudstor
teaching_data = cloudstor(url="L93TxcmtLQzcfbk", password='')

teaching_data.download_file_if_distinct("BlueMarbleNG-TB_2004-12-01_rgb_3600x1800.TIFF", "Resources/BlueMarbleNG-TB_2004-12-01_rgb_3600x1800.TIFF")
teaching_data.download_file_if_distinct("color_etopo1_ice_low.tif", "Resources/color_etopo1_ice_low.tif")
teaching_data.download_file_if_distinct("EMAG2_image_V2_no_compr.tif", "Resources/EMAG2_image_V2_no_compr.tif")
teaching_data.download_file_if_distinct("global_age_data.3.6.z.npz", "Resources/global_age_data.3.6.z.npz")
teaching_data.download_file_if_distinct("etopo1_grayscale_hillshade.tif", "Resources/etopo1_grayscale_hillshade.tif")
teaching_data.download_file_if_distinct("ETOPO1_Ice_c_geotiff.tif", "Resources/ETOPO1_Ice_c_geotiff.tif")

teaching_data.download_file_if_distinct("HYP_50M_SR_W/HYP_50M_SR_W.tif", "Resources/HYP_50M_SR_W/HYP_50M_SR_W.tif")
teaching_data.download_file_if_distinct("OB_50M/OB_50M.tif", "Resources/OB_50M/OB_50M.tif")
teaching_data.download_file_if_distinct("", "Resources/")
# The colormap routine creates enormous arrays in intermediary calculations. This is
# a way to avoid memory errors: process to RGB (int8) in advance

def apply_colormap_to_image(rawimage, colormap, norm):

    greyimage = norm(rawimage)
    rgbimage = np.empty((greyimage.shape[0], greyimage.shape[1] , 4), dtype=uint8)

    for i in range(0, greyimage.shape[0]):
        rgbimage[i,:,:] = colormap(greyimage[i,:]) * 256
    rgbimage2 = rgbimage[:,:,0:3]        
    return rgbimage2
base_projection     = ccrs.PlateCarree() 
global_extent     = [-180.0, 180.0, -90.0, 90.0]

strainrate = numpy.loadtxt("Resources/sec_invariant_strain_0.2.dat")
strainrate_data = strainrate.reshape(741,1800,3)  # I had to look at the data to work this out !
# strainrate_img  = strainrate_data[:,:,2] # Not actually used here !
# Note: we need to manage the memory here
# memory constraints than the native installation. Be sure to del() unused arrays/images

# Etopo Height field as geotiff

etopoH = gdal.Open("Resources/ETOPO1_Ice_c_geotiff.tif") 
etopoH_img = etopoH.ReadAsArray()[::2,::2].astype(numpy.float16)

colormap = plt.get_cmap('Greys')
norm = matplotlib.colors.Normalize(vmin=-10000, vmax=5000)
etopoH_img_grey = apply_colormap_to_image(etopoH_img, colormap, norm)

# These are alternative images 

# etopo1       = gdal.Open("Resources/color_etopo1_ice_low.tif")
# etopo_img    = etopo1.ReadAsArray().transpose(1,2,0)
# del(etopo1)

# global_shrelief = gdal.Open("Resources/etopo1_grayscale_hillshade.tif")
# global_shrelief_img = global_shrelief.ReadAsArray()[::2,::2].astype(numpy.float16)
# del(global_shrelief)
# print "global_shrelief_img - ", global_shrelief_img.shape
# How does that look in the Native Projection ?

fig = plt.figure(figsize=(12, 6), facecolor="none")
ax = plt.axes(projection=base_projection)
ax.imshow(etopoH_img_grey[::1,::1], transform=ccrs.PlateCarree(), origin="upper", 
          alpha=1.0, extent=global_extent,  zorder=0)
# Some features we need for the map 

coastline = cfeature.NaturalEarthFeature('physical', 'coastline', '50m',

rivers = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '50m',
                                        edgecolor='Blue', facecolor="none")

lakes = cfeature.NaturalEarthFeature('physical', 'lakes', '50m',
                                        edgecolor="blue", facecolor="blue")

ocean = cfeature.NaturalEarthFeature('physical', 'ocean', '50m',

graticules_5 = cfeature.NaturalEarthFeature('physical', 'graticules_5', '10m',
                           edgecolor="black", facecolor=None)
## Reading the velocity vector data from the EU-fixed dataset

rootgrp = netcdf.netcdf_file(filename="Resources/", version=2)

ve = rootgrp.variables["ve"]
vn = rootgrp.variables["vn"]

lonv = rootgrp.variables["lon"]
latv = rootgrp.variables["lat"]

lons = lonv[::5]
lats = latv[::5]
llX, llY = np.meshgrid(lons,lats)

#llX = llX.reshape(-1)
#llY = llY.reshape(-1)

Veast = (np.array(ve[::5,::5]).T)
Vnorth = (np.array(vn[::5,::5]).T)

Vorientation = np.arctan2(Vnorth,Veast)
# I used Robinson in the article but it keeps crashing for me now !

#projection = ccrs.PlateCarree() # For debugging, this one is quick
#projection = ccrs.Orthographic(central_longitude=-150.0, central_latitude=60.0, globe=None)
#projection = ccrs.Orthographic(central_longitude=120.0, central_latitude=5.0, globe=None)
projection = ccrs.Mollweide(central_longitude=120)
#projection = ccrs.Robinson(central_longitude=120) # This is the one in the paper

fig = plt.figure(figsize=(12, 6), facecolor="none")

ax = plt.axes(projection=projection)

# Still too large for the memory 
ax.imshow(etopoH_img_grey[::4,::4,:], transform=ccrs.PlateCarree(), origin="upper", 
          alpha=1.0, extent=global_extent, interpolation="spline16", zorder=0)

mappable2 = ax.contourf(strainrate_data[:,:,0], strainrate_data[:,:,1], strainrate_data[:,:,2], 
         levels=[10, 20, 50, 100, 200, 500], linestyle=None, vmin=1.0, vmax=200,
         transform=base_projection,  cmap=cm.OrRd_r, alpha=0.75, 
         extent=strainrate_extent, extend="max", zorder=12)

plt.colorbar(mappable=mappable2, shrink=0.5)

mappable1 = ax.quiver(llX, llY, Veast, Vnorth, Vorientation, scale=2000, transform=ccrs.PlateCarree(),
        cmap=cm.Blues, alpha=0.7, zorder=13, pivot="mid")

ax.add_feature(coastline, linewidth=1.0,  edgecolor="#000000", zorder=2, alpha=0.75)
# ax.add_feature(rivers,    linewidth=1,    edgecolor="Blue", zorder=2)
# ax.add_feature(lakes,     linewidth=1,    edgecolor="Blue", zorder=3, alpha=0.25)
# ax.add_feature(graticules_5, linewidth=0.25, edgecolor="black", linestyle=":", zorder=1, alpha=0.75)
#fig.savefig("PlateMotionsPlateCarree.png", dpi=300) 
#fig.savefig("PlateMotionsRobinson.png", dpi=500)
fig.savefig("PlateMotionsMollweide.png", dpi=500)
## This was an alternative to the background image which I thought might look good:

globalrelief      = gdal.Open("Resources/HYP_50M_SR_W/HYP_50M_SR_W.tif")
globalrelief_img  = globalrelief.ReadAsArray().transpose(1,2,0)

globalbathym      = gdal.Open("Resources/OB_50M/OB_50M.tif")
globalbathym_img  = globalbathym.ReadAsArray().transpose(1,2,0)

print ("etopoH_img - ", etopoH_img.shape)
print ("globalrelief_img - ", globalrelief_img.shape)

blended_img = np.empty_like(globalrelief_img)
blended_img[...,0] = np.where( etopoH_img < 0.0, globalbathym_img[...,0], globalrelief_img[...,0] )
blended_img[...,1] = np.where( etopoH_img < 0.0, globalbathym_img[...,1], globalrelief_img[...,1] )
blended_img[...,2] = np.where( etopoH_img < 0.0, globalbathym_img[...,2], globalrelief_img[...,2] )


fig = plt.figure(figsize=(12, 6), facecolor="none")
ax = plt.axes(projection=base_projection)
ax.imshow(blended_img[::1,::1], transform=ccrs.PlateCarree(), origin="upper", 
          alpha=1.0, extent=global_extent,  zorder=0)
# Do this if the relief / bathym sizes don't match the etopo data (to make the blended image)
# The datasets we downloaded can be manipulated trivially without the need for this and I have
# commented it all out so you can run all cells without reprocessing the data files. 

import scipy.ndimage
import scipy.misc

etopoH = gdal.Open("Resources/ETOPO1_Ice_g_geotiff.tif")
etopoH_img = etopoH.ReadAsArray()


etopoH_transform = etopoH.GetGeoTransform()
globalrelief_transform = globalrelief.GetGeoTransform()

# Resize to match globalrelief ... this resize is int only ??

globaletopoH = scipy.misc.imresize(etopoH_img, globalrelief_img.shape, mode='F')

## How to turn this array back into the appropriate geotiff

from osgeo import gdal
from osgeo import osr

# data exists in 'ary' with values range 0 - 255
# Uncomment the next line if ary[0][0] is upper-left corner
#ary = numpy.flipup(ary)

Ny, Nx = globaletopoH.shape
driver = gdal.GetDriverByName("GTiff")
# Final argument is optional but will produce much smaller output file
ds = driver.Create('output.tif', Nx, Ny, 1, gdal.GDT_Float64, ['COMPRESS=LZW'])

# this assumes the projection is Geographic lat/lon WGS 84
srs = osr.SpatialReference()

ds.SetGeoTransform( globalrelief_transform  ) # define GeoTransform tuple
ds = None