CRUST1.0¶
Ondřej Šrámek, Charles University, ondrej.sramek@gmail.com, 6 March 2018
CRUST1.0 model: http://igppweb.ucsd.edu/~gabi/crust1.html
CRUST1.0 describes the lithosphere as an 8-layer sandwich.
At each (latitude,longitude) point, the density is constant within each layer.
Elevations refer to a reference Earth shape, positive elevations above, negative below.
Density is set to zero where layer thickness is zero.
In the following, I just plot those files as maps.
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]:
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
x05 = np.arange(-179.75, 180, 0.5)
y05 = np.arange(89.75, -90, -0.5)
x05, y05 = np.meshgrid(x05, y05)
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 [4]:
plt.rcParams['figure.figsize'] = (18,9)
In [5]:
z = c1_iscont
fig, ax = fig1(x1, y1, z, 'Continent', '1/0')
plt.show()
Plots of CRUST1.0 elevation (km) of tops of layers¶
In [6]:
z = c1_top1
fig, ax = fig1(x1, y1, z, 'CRUST1.0 top1 (water)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [7]:
z = c1_top2
fig, ax = fig1(x1, y1, z, 'CRUST1.0 top2 (ice)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [8]:
z = c1_top3
fig, ax = fig1(x1, y1, z, 'CRUST1.0 top3 (sediment1)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [9]:
z = c1_top4
fig, ax = fig1(x1, y1, z, 'CRUST1.0 top4 (sediment2)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [10]:
z = c1_top5
fig, ax = fig1(x1, y1, z, 'CRUST1.0 top5 (sediment3)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [11]:
z = c1_top6
fig, ax = fig1(x1, y1, z, 'CRUST1.0 top6 (upper crust)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [12]:
z = c1_top7
fig, ax = fig1(x1, y1, z, 'CRUST1.0 top7 (middle crust)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [13]:
z = c1_top8
fig, ax = fig1(x1, y1, z, 'CRUST1.0 top8 (lower crust)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [14]:
z = c1_top9
fig, ax = fig1(x1, y1, z, 'CRUST1.0 top9 (mantle)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
Plots of CRUST1.0 layer thicknesses (km)¶
In [15]:
z = c1_d1
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Water thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [16]:
z = c1_d2
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Ice thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [17]:
z = c1_d3
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Sediment1 thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [18]:
z = c1_d4
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Sediment2 thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [19]:
z = c1_d5
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Sediment3 thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [20]:
z = c1_d6
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Upper crust thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [21]:
z = c1_d7
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Middle crust thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [22]:
z = c1_d8
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Lower crust thickness', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
CRUST1.0 crystalline crust thicknesses, layers 6–8 (km)¶
In [23]:
z = c1_d678
fig, ax = fig1(x1, y1, z, 'CRUST1.0 crustal thickness (layers 6–8)', 'km')
plt.show()
print('min =', z.min())
print('max =', z.max())
Plots of CRUST1.0 density (kg/m3) in layers¶
In [24]:
z = c1_rho1
print('min =', z.min())
print('max =', z.max())
In [25]:
z = c1_rho2
print('min =', z.min())
print('max =', z.max())
In [26]:
z = c1_rho3
fig, ax = fig1(x1, y1, z, 'CRUST1.0 rho3 (sediment1)', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [27]:
z = c1_rho4
fig, ax = fig1(x1, y1, z, 'CRUST1.0 rho4 (sediment2)', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [28]:
z = c1_rho5
fig, ax = fig1(x1, y1, z, 'CRUST1.0 rho5 (sediment3)', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [29]:
z = c1_rho6
fig, ax = fig1(x1, y1, z, 'CRUST1.0 rho6 (upper crust)', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [30]:
z = c1_rho7
fig, ax = fig1(x1, y1, z, 'CRUST1.0 rho7 (middle crust)', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [31]:
z = c1_rho8
fig, ax = fig1(x1, y1, z, 'CRUST1.0 rho8 (lower crust)', 'kg/m3')
plt.show()
print('min =', z.min())
print('max =', z.max())
Plots of CRUST1.0 Vp (km/s)¶
In [32]:
z = c1_vp1
print('min =', z.min())
print('max =', z.max())
In [33]:
z = c1_vp2
print('min =', z.min())
print('max =', z.max())
In [34]:
z = c1_vp3
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Vp3 (sediment1)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [35]:
z = c1_vp4
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Vp4 (sediment2)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [36]:
z = c1_vp5
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Vp5 (sediment3)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [37]:
z = c1_vp6
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Vp6 (upper crust)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [38]:
z = c1_vp7
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Vp7 (middle crust)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [39]:
z = c1_vp8
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Vp8 (lower crust)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
Plots of CRUST1.0 Vs (km/s)¶
In [40]:
z = c1_vs1
print('min =', z.min())
print('max =', z.max())
In [41]:
z = c1_vs2
print('min =', z.min())
print('max =', z.max())
In [42]:
z = c1_vs3
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Vs3 (sediment1)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [43]:
z = c1_vs4
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Vs4 (sediment2)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [44]:
z = c1_vs5
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Vs5 (sediment3)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [45]:
z = c1_vs6
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Vs6 (upper crust)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [46]:
z = c1_vs7
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Vs7 (middle crust)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
In [47]:
z = c1_vs8
fig, ax = fig1(x1, y1, z, 'CRUST1.0 Vs8 (lower crust)', 'km/s')
plt.show()
print('min =', z.min())
print('max =', z.max())
Histograms Vp of continental crust¶
In [48]:
plt.rcParams['figure.figsize'] = (18,12)
In [49]:
# weights to account for latitude effect on surface area
nlat = c1_vp6.shape[0]
sintheta = np.sin( np.arange(0.5/nlat, 1, 1/nlat) * pi )
c1_wt = np.zeros(c1_vp6.shape)
for i in range(0,c1_vp6.shape[1]):
c1_wt[:,i] = sintheta[:]
In [50]:
dz = 0.025
zmin = min(c1_vp6.min(), c1_vp7.min(), c1_vp8.min())
zmax = max(c1_vp6.max(), c1_vp7.max(), c1_vp8.max())
zrange = ( (int(zmin/dz)+1)*dz, int(zmax/dz)*dz )
nbins = int((zrange[1]-zrange[0])/dz)
fig = plt.figure()
wt1 = c1_wt.flatten()*c1_iscont.flatten()
#
plt.subplot(3, 1, 1)
plt.hist(c1_vp6.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('CRUST1.0 Vp6 (upper crust) – CONTINENTS ONLY')
plt.ylabel('#')
#
plt.subplot(3, 1, 2)
plt.hist(c1_vp7.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('CRUST1.0 Vp7 (middle crust) – CONTINENTS ONLY')
plt.ylabel('#')
#
plt.subplot(3, 1, 3)
plt.hist(c1_vp8.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('CRUST1.0 Vp8 (lower crust) – CONTINENTS ONLY')
plt.xlabel('km/s')
plt.ylabel('#')
#
plt.show()
Histograms fractional thickness of continental crustal layers¶
In [51]:
# thickness of crystalline crustal layers
c1_d678 = c1_top6 - c1_top9
c1_fd6 = c1_d6/c1_d678
c1_fd7 = c1_d7/c1_d678
c1_fd8 = c1_d8/c1_d678
dz = 0.01
zmin = min(c1_fd6.min(), c1_fd7.min(), c1_fd8.min())
zmax = max(c1_fd6.max(), c1_fd7.max(), c1_fd8.max())
zrange = ( (int(zmin/dz)+1)*dz, int(zmax/dz)*dz )
nbins = int((zrange[1]-zrange[0])/dz)
fig = plt.figure()
wt1 = c1_wt.flatten()*c1_iscont.flatten()
#
plt.subplot(3, 1, 1)
plt.hist(c1_fd6.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('CRUST1.0 6/678 thickness (upper crust / crust) – CONTINENTS ONLY')
plt.ylabel('#')
#
plt.subplot(3, 1, 2)
plt.hist(c1_fd7.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('CRUST1.0 7/678 thickness (middle crust / crust) – CONTINENTS ONLY')
plt.ylabel('#')
#
plt.subplot(3, 1, 3)
plt.hist(c1_fd8.flatten(), bins=nbins, range=zrange, weights=wt1)
plt.title('CRUST1.0 8/678 thickness (lower crust / crust) – CONTINENTS ONLY')
plt.xlabel('m/s')
plt.ylabel('#')
#
plt.show()
In [52]:
z = c1_fd6
fig, ax = fig1(x1, y1, z, 'CRUST1.0 lay6/lays678 thickness (upper crust / crust)', 'fraction')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = c1_fd7
fig, ax = fig1(x1, y1, z, 'CRUST1.0 lay7/lays678 thickness (upper crust / crust)', 'fraction')
plt.show()
print('min =', z.min())
print('max =', z.max())
z = c1_fd8
fig, ax = fig1(x1, y1, z, 'CRUST1.0 lay8/lays678 thickness (upper crust / crust)', 'fraction')
plt.show()
print('min =', z.min())
print('max =', z.max())
