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
#%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_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())
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())
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())
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())
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())
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())