LITHO1.0 relative to CRUST1.0

Ondřej Šrámek, Charles University, ondrej.sramek@gmail.com, 13 April 2018

LITHO1.0 model: http://igppweb.ucsd.edu/~gabi/litho1.0.html

CRUST1.0 model: http://igppweb.ucsd.edu/~gabi/crust1.html

Pasyanos, Masters, Laske, Ma (2014): "LITHO1.0: An updated crust and lithospheric model of the Earth." J. Geophys. Res., 119(3): 2153–2173, doi:10.1002/2013JB010626

In [1]:
#%reset -f
from mpl_toolkits.basemap import Basemap, cm
import matplotlib.pyplot as plt
import numpy as np
from math import pi
In [2]:
# read-in CRUST1.0 data
#c1_top1 = np.loadtxt('crust1_1deg_top1.txt')
#c1_top2 = np.loadtxt('crust1_1deg_top2.txt')
#c1_top3 = np.loadtxt('crust1_1deg_top3.txt')
#c1_top4 = np.loadtxt('crust1_1deg_top4.txt')
#c1_top5 = np.loadtxt('crust1_1deg_top5.txt')
#c1_top6 = np.loadtxt('crust1_1deg_top6.txt')
#c1_top7 = np.loadtxt('crust1_1deg_top7.txt')
#c1_top8 = np.loadtxt('crust1_1deg_top8.txt')
#c1_top9 = np.loadtxt('crust1_1deg_top9.txt')

#c1_rho1 = np.loadtxt('crust1_1deg_rho1.txt')
#c1_rho2 = np.loadtxt('crust1_1deg_rho2.txt')
#c1_rho3 = np.loadtxt('crust1_1deg_rho3.txt')
#c1_rho4 = np.loadtxt('crust1_1deg_rho4.txt')
#c1_rho5 = np.loadtxt('crust1_1deg_rho5.txt')
c1_rho6 = np.loadtxt('crust1_1deg_rho6.txt')
c1_rho7 = np.loadtxt('crust1_1deg_rho7.txt')
c1_rho8 = np.loadtxt('crust1_1deg_rho8.txt')

#c1_vp1 = np.loadtxt('crust1_1deg_vp1.txt')
#c1_vp2 = np.loadtxt('crust1_1deg_vp2.txt')
#c1_vp3 = np.loadtxt('crust1_1deg_vp3.txt')
#c1_vp4 = np.loadtxt('crust1_1deg_vp4.txt')
#c1_vp5 = np.loadtxt('crust1_1deg_vp5.txt')
c1_vp6 = np.loadtxt('crust1_1deg_vp6.txt')
c1_vp7 = np.loadtxt('crust1_1deg_vp7.txt')
c1_vp8 = np.loadtxt('crust1_1deg_vp8.txt')

#c1_vs1 = np.loadtxt('crust1_1deg_vs1.txt')
#c1_vs2 = np.loadtxt('crust1_1deg_vs2.txt')
#c1_vs3 = np.loadtxt('crust1_1deg_vs3.txt')
#c1_vs4 = np.loadtxt('crust1_1deg_vs4.txt')
#c1_vs5 = np.loadtxt('crust1_1deg_vs5.txt')
#c1_vs6 = np.loadtxt('crust1_1deg_vs6.txt')
#c1_vs7 = np.loadtxt('crust1_1deg_vs7.txt')
#c1_vs8 = np.loadtxt('crust1_1deg_vs8.txt')

c1_iscont = np.loadtxt('crust1_1deg_iscont.txt')

#c1_d1 = c1_top1 - c1_top2
#c1_d2 = c1_top2 - c1_top3
#c1_d3 = c1_top3 - c1_top4
#c1_d4 = c1_top4 - c1_top5
#c1_d5 = c1_top5 - c1_top6
#c1_d6 = c1_top6 - c1_top7
#c1_d7 = c1_top7 - c1_top8
#c1_d8 = c1_top8 - c1_top9

#c1_d678 = c1_top6 - c1_top9
In [3]:
# read-in LITHO1.0 data
#l1_top1 = np.loadtxt('litho1_1deg_top1.txt')
#l1_top2 = np.loadtxt('litho1_1deg_top2.txt')
#l1_top3 = np.loadtxt('litho1_1deg_top3.txt')
#l1_top4 = np.loadtxt('litho1_1deg_top4.txt')
#l1_top5 = np.loadtxt('litho1_1deg_top5.txt')
#l1_top6 = np.loadtxt('litho1_1deg_top6.txt')
#l1_top7 = np.loadtxt('litho1_1deg_top7.txt')
#l1_top8 = np.loadtxt('litho1_1deg_top8.txt')
#l1_top9 = np.loadtxt('litho1_1deg_top9.txt')
#l1_top10 = np.loadtxt('litho1_1deg_top10.txt')

#l1_rho1 = np.loadtxt('litho1_1deg_rho1.txt')
#l1_rho2 = np.loadtxt('litho1_1deg_rho2.txt')
#l1_rho3 = np.loadtxt('litho1_1deg_rho3.txt')
#l1_rho4 = np.loadtxt('litho1_1deg_rho4.txt')
#l1_rho5 = np.loadtxt('litho1_1deg_rho5.txt')
l1_rho6 = np.loadtxt('litho1_1deg_rho6.txt')
l1_rho7 = np.loadtxt('litho1_1deg_rho7.txt')
l1_rho8 = np.loadtxt('litho1_1deg_rho8.txt')
#l1_rho9 = np.loadtxt('litho1_1deg_rho9.txt')

#l1_vp1 = np.loadtxt('litho1_1deg_vp1.txt') * 1e-3 # in km/s
#l1_vp2 = np.loadtxt('litho1_1deg_vp2.txt') * 1e-3
#l1_vp3 = np.loadtxt('litho1_1deg_vp3.txt') * 1e-3
#l1_vp4 = np.loadtxt('litho1_1deg_vp4.txt') * 1e-3
#l1_vp5 = np.loadtxt('litho1_1deg_vp5.txt') * 1e-3
l1_vp6 = np.loadtxt('litho1_1deg_vp6.txt') * 1e-3
l1_vp7 = np.loadtxt('litho1_1deg_vp7.txt') * 1e-3
l1_vp8 = np.loadtxt('litho1_1deg_vp8.txt') * 1e-3
#l1_vp9 = np.loadtxt('litho1_1deg_vp9.txt') * 1e-3

#l1_vs1 = np.loadtxt('litho1_1deg_vs1.txt') * 1e-3
#l1_vs2 = np.loadtxt('litho1_1deg_vs2.txt') * 1e-3
#l1_vs3 = np.loadtxt('litho1_1deg_vs3.txt') * 1e-3
#l1_vs4 = np.loadtxt('litho1_1deg_vs4.txt') * 1e-3
#l1_vs5 = np.loadtxt('litho1_1deg_vs5.txt') * 1e-3
#l1_vs6 = np.loadtxt('litho1_1deg_vs6.txt') * 1e-3
#l1_vs7 = np.loadtxt('litho1_1deg_vs7.txt') * 1e-3
#l1_vs8 = np.loadtxt('litho1_1deg_vs8.txt') * 1e-3
#l1_vs9 = np.loadtxt('litho1_1deg_vs9.txt') * 1e-3

#l1_d1 = l1_top1 - l1_top2
#l1_d2 = l1_top2 - l1_top3
#l1_d3 = l1_top3 - l1_top4
#l1_d4 = l1_top4 - l1_top5
#l1_d5 = l1_top5 - l1_top6
#l1_d6 = l1_top6 - l1_top7
#l1_d7 = l1_top7 - l1_top8
#l1_d8 = l1_top8 - l1_top9
#l1_d9 = l1_top9 - l1_top10

#l1_d678 = l1_top6 - l1_top9
#l1_d6789 = l1_top6 - l1_top10
In [4]:
def fig1(x, y, z, title, unit):
    fig, ax = plt.subplots()
    m = Basemap(projection='cyl', llcrnrlat=-90, urcrnrlat=90, llcrnrlon=-180, urcrnrlon=180, resolution='c')
    m.drawcoastlines()
    m.drawparallels(np.arange(-90.,120.,30.))
    m.drawmeridians(np.arange(0.,360.,60.))
    x, y = m(x, y)
    clevs = np.linspace(z.min(), z.max(), 19)
    cs = m.contourf(x, y, z, clevs, cmap='rainbow')
    #cs = m.pcolormesh(x, y, z, vmin=z.min(), vmax=z.max(), cmap='rainbow')
    cbar = m.colorbar(cs)
    cbar.set_label(unit)
    plt.title(title)
    return fig, ax

x1 = np.arange(-179.5, 180, 1.0)
y1 = np.arange(89.5, -90, -1.0)
x1, y1 = np.meshgrid(x1, y1)

Continental mask from CRUST1.0

In [5]:
plt.rcParams['figure.figsize'] = (18,9)
In [6]:
#z = c1_iscont
#fig, ax = fig1(x1, y1, z, 'Continent', '1/0')
#plt.show()

Plots of LITHO1.0–CRUST1.0 Vp difference in %

In [7]:
z = (l1_vp6-c1_vp6)/c1_vp6*100
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in upper crust Vp [percent]', '%')
plt.show()
print('min =', z.min())
print('max =', z.max())
/Users/ondrej/anaconda/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3608: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.
  b = ax.ishold()
/Users/ondrej/anaconda/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3677: MatplotlibDeprecationWarning: axes.hold is deprecated.
    See the API Changes document (http://matplotlib.org/api/api_changes.html)
    for more details.
  ax.hold(b)
min = -21.216
max = 24.5646
In [8]:
z = (l1_vp7-c1_vp7)/c1_vp7*100
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in middle crust Vp in [percent]', '%')
plt.show()
print('min =', z.min())
print('max =', z.max())
/Users/ondrej/anaconda/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3608: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.
  b = ax.ishold()
/Users/ondrej/anaconda/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3677: MatplotlibDeprecationWarning: axes.hold is deprecated.
    See the API Changes document (http://matplotlib.org/api/api_changes.html)
    for more details.
  ax.hold(b)
min = -14.4634848485
max = 12.3920338983
In [9]:
z = (l1_vp8-c1_vp8)/c1_vp8*100
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in lower crust Vp in [percent]', '%')
plt.show()
print('min =', z.min())
print('max =', z.max())
/Users/ondrej/anaconda/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3608: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.
  b = ax.ishold()
/Users/ondrej/anaconda/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3677: MatplotlibDeprecationWarning: axes.hold is deprecated.
    See the API Changes document (http://matplotlib.org/api/api_changes.html)
    for more details.
  ax.hold(b)
min = -10.5630985915
max = 12.5692424242

Plots of LITHO1.0–CRUST1.0 density difference in percent %

In [10]:
z = (l1_rho6-c1_rho6)/c1_rho6*100
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in upper crust density [percent]', '%')
plt.show()
print('min =', z.min())
print('max =', z.max())
/Users/ondrej/anaconda/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3608: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.
  b = ax.ishold()
/Users/ondrej/anaconda/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3677: MatplotlibDeprecationWarning: axes.hold is deprecated.
    See the API Changes document (http://matplotlib.org/api/api_changes.html)
    for more details.
  ax.hold(b)
min = -12.3136029412
max = 11.5617977528
In [11]:
z = (l1_rho7-c1_rho7)/c1_rho7*100
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in middle crust density [percent]', '%')
plt.show()
print('min =', z.min())
print('max =', z.max())
/Users/ondrej/anaconda/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3608: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.
  b = ax.ishold()
/Users/ondrej/anaconda/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3677: MatplotlibDeprecationWarning: axes.hold is deprecated.
    See the API Changes document (http://matplotlib.org/api/api_changes.html)
    for more details.
  ax.hold(b)
min = -9.79929824561
max = 11.0912408759
In [12]:
z = (l1_rho8-c1_rho8)/c1_rho8*100
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in lower crust density [percent]', '%')
plt.show()
print('min =', z.min())
print('max =', z.max())
/Users/ondrej/anaconda/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3608: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.
  b = ax.ishold()
/Users/ondrej/anaconda/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3677: MatplotlibDeprecationWarning: axes.hold is deprecated.
    See the API Changes document (http://matplotlib.org/api/api_changes.html)
    for more details.
  ax.hold(b)
min = -9.90426229508
max = 9.87517482517