Ondřej Šrámek, Charles University, ondrej.sramek@gmail.com, 6 March 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
#%reset -f
from mpl_toolkits.basemap import Basemap, cm
import matplotlib.pyplot as plt
import numpy as np
from math import pi
# 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
# 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
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)
plt.rcParams['figure.figsize'] = (18,9)
z = c1_iscont
fig, ax = fig1(x1, y1, z, 'Continent', '1/0')
plt.show()
z = l1_d1-c1_d1
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in water thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_d2-c1_d2
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in ice thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_d3-c1_d3
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in sediment1 thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_d4-c1_d4
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in sediment2 thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_d5-c1_d5
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in sediment3 thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_d6-c1_d6
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in upper crust thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_d7-c1_d7
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in middle crust thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_d8-c1_d8
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in lower crust thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_rho3-c1_rho3
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in sediment1 density', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_rho4-c1_rho4
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in sediment2 density', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_rho5-c1_rho5
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in sediment3 density', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_rho6-c1_rho6
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in upper crust density', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_rho7-c1_rho7
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in middle crust density', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_rho8-c1_rho8
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in lower crust density', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_vp3-c1_vp3
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in sediment1 Vp', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_vp4-c1_vp4
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in sediment2 Vp', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_vp5-c1_vp5
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in sediment3 Vp', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_vp6-c1_vp6
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in upper crust Vp', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_vp7-c1_vp7
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in middle crust Vp', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_vp8-c1_vp8
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in lower crust Vp', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_vs3-c1_vs3
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in sediment1 Vs', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_vs4-c1_vs4
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in sediment2 Vs', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_vs5-c1_vs5
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in sediment3 Vs', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_vs6-c1_vs6
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in upper crust Vs', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_vs7-c1_vs7
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in middle crust Vs', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = l1_vs8-c1_vs8
fig, ax = fig1(x1, y1, z, 'LITHO1.0–CRUST1.0 difference in lower crust Vs', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())