Moho depth: GEMMA vs LITHO vs CRUST¶
Ondřej Šrámek, Charles University, ondrej.sramek@gmail.com, 25 Nov 2017
- CRUST2.0 http://igppweb.ucsd.edu/~gabi/crust2.html
- CRUST1.0 http://igppweb.ucsd.edu/~gabi/crust1.html
- LITHO1.0 http://igppweb.ucsd.edu/~gabi/litho1.0.html
- GEMMA1.0 http://gocedata.como.polimi.it
The dataset was downloaded from http://gocedata.como.polimi.it, files t6.asc and moho_err.asc packed in moho.zip. No description or README file included. Comparison with Reguzzioni & Sampietro 2015 (doi:10.1016/j.jag.2014.04.002), especially their Figures 8 and 9, suggests the dataset is Moho depth.
Question remains: How is GEMMA1.0 Moho depth is defined? Could be:
- Moho depth relative to zero elevation surface.
- Moho depth relative to actual surface of the solid/liquid Earth (i.e., ocean level in oceans).
- Moho depth relative to actual surface of the solid Earth (i.e., ocean bottom in oceans).
This version assumes 2
In [1]:
#%reset -f
case = 2
from mpl_toolkits.basemap import Basemap, cm
import matplotlib.pyplot as plt
import numpy as np
from math import pi, sqrt
plt.rcParams['figure.figsize'] = (18,9)
In [2]:
gemma = -np.loadtxt('t6.asc', skiprows=6)
gemma_err = np.loadtxt('moho_err.asc', skiprows=6)
crust2_md = np.loadtxt('crust2_md.z')
crust2_th1 = np.loadtxt('crust2_th1.z')
crust2_th3 = np.loadtxt('crust2_th3.z')
crust1_md = np.loadtxt('crust1_md.z')
crust1_th1 = np.loadtxt('crust1_th1.z')
crust1_th3 = np.loadtxt('crust1_th3.z')
litho1_md = np.loadtxt('litho1_md.z')
litho1_th1 = np.loadtxt('litho1_th1.z')
litho1_th3 = np.loadtxt('litho1_th3.z')
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)
y1 = np.arange(89.5, -90, -1)
x1, y1 = np.meshgrid(x1, y1)
x2 = np.arange(-179, 180, 2)
y2 = np.arange(89, -90, -2)
x2, y2 = np.meshgrid(x2, y2)
In [3]:
if case == 1:
# comparing gemma to Moho depths relative to zero elevation
crust2 = crust2_md
crust1 = crust1_md
litho1 = litho1_md
elif case == 2:
# comparing gemma to Moho depths relative to surface of the solid/liquid Earth
crust2 = crust2_th1
crust1 = crust1_th1
litho1 = litho1_th1
elif case == 3:
# comparing gemma to Moho depths relative to surface of the solid Earth
crust2 = crust2_th3
crust1 = crust1_th3
litho1 = litho1_th3
In [4]:
def fig1(x, y, z, title):
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('km')
plt.title(title)
return fig, ax
def fig2(x, y, z, title, clevs):
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)
cs = m.contourf(x, y, z, clevs, cmap='seismic')
cbar = m.colorbar(cs)
cbar.set_label('km')
plt.title(title)
return fig, ax
In [5]:
fig, ax = fig1(x05, y05, gemma, 'GEMMA1.0')
plt.show()
fig, ax = fig1(x05, y05, gemma_err, 'GEMMA1.0 uncertainty')
plt.show()
In [6]:
fig, ax = fig1(x1, y1, litho1, 'LITHO1.0')
plt.show()
fig, ax = fig1(x1, y1, crust1, 'CRUST1.0')
plt.show()
fig, ax = fig1(x2, y2, crust2, 'CRUST2.0')
plt.show()
In [7]:
plt.rcParams['figure.figsize'] = (18,9)
vmin = min(gemma.min(), litho1.min(), crust1.min(), crust2.min())
vmax = max(gemma.max(), litho1.max(), crust1.max(), crust2.max())
clevs = np.linspace(vmin, vmax, 19)
fig = plt.figure()
#
plt.subplot(2, 2, 1)
z = gemma
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(x05, y05)
cs = m.contourf(x, y, z, clevs, cmap='rainbow')
plt.title('GEMMA1.0')
#
plt.subplot(2, 2, 2)
z = litho1
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(x1, y1)
cs = m.contourf(x, y, z, clevs, cmap='rainbow')
cbar = m.colorbar(cs)
cbar.set_label('km')
plt.title('LITHO1.0')
#
plt.subplot(2, 2, 3)
z = crust2
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(x2, y2)
cs = m.contourf(x, y, z, clevs, cmap='rainbow')
plt.title('CRUST2.0')
#
plt.subplot(2, 2, 4)
z = crust1
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(x1, y1)
cs = m.contourf(x, y, z, clevs, cmap='rainbow')
cbar = m.colorbar(cs)
cbar.set_label('km')
plt.title('CRUST1.0')
#
plt.show()
In [8]:
# extrapolate to 0.5 deg mesh
litho1_05 = np.zeros(gemma.shape)
crust1_05 = np.zeros(gemma.shape)
crust2_05 = np.zeros(gemma.shape)
for i in np.arange(litho1.shape[0]):
for j in np.arange(litho1.shape[1]):
litho1_05[2*i:2*i+2, 2*j:2*j+2] = litho1[i][j]
crust1_05[2*i:2*i+2, 2*j:2*j+2] = crust1[i][j]
for i in np.arange(crust2.shape[0]):
for j in np.arange(crust2.shape[1]):
crust2_05[4*i:4*i+4, 4*j:4*j+4] = crust2[i][j]
diff_litho1 = gemma - litho1_05
diff_crust1 = gemma - crust1_05
diff_crust2 = gemma - crust2_05
vmin = min(diff_litho1.min(), diff_crust1.min(), diff_crust2.min())
vmax = max(diff_litho1.max(), diff_crust1.max(), diff_crust2.max())
Note that LITHO1.0 is dumbly extrapolated from 1 degree grid to 0.5 degree grid.
A proper way would be to interpolate using the litho interpolation code originally from Pasyanos.
In [9]:
plt.rcParams['figure.figsize'] = (18,9)
cmap = 'rainbow'
fig, ax = fig1(x05, y05, diff_litho1, 'Difference GEMMA1.0 – LITHO1.0 [km]')
plt.show()
fig, ax = fig1(x05, y05, diff_crust1, 'Difference GEMMA1.0 – CRUST1.0 [km]')
plt.show()
fig, ax = fig1(x05, y05, diff_crust2, 'Difference GEMMA1.0 – CRUST2.0 [km]')
plt.show()
And plotting with the same scale on the three difference plots:
In [10]:
clevs = np.zeros(13)
clevs[0:7] = np.linspace(np.floor(vmin), 0, 7)
clevs[6:13] = np.linspace(0, np.ceil(vmax), 7)
fig, ax = fig2(x05, y05, diff_litho1, 'Difference GEMMA1.0 – LITHO1.0 [km]', clevs)
plt.show()
fig, ax = fig2(x05, y05, diff_crust1, 'Difference GEMMA1.0 – CRUST1.0 [km]', clevs)
plt.show()
fig, ax = fig2(x05, y05, diff_crust2, 'Difference GEMMA1.0 – CRUST2.0 [km]', clevs)
plt.show()
Show min max values
Calculate mean surface area-averaged value of differences
In [11]:
print (' MIN MAX [km]')
print ('gemma1 {:6.1f} {:6.1f}'.format(gemma.min(), gemma.max()))
print ('crust2 {:6.1f} {:6.1f}'.format(crust2.min(), crust2.max()))
print ('crust1 {:6.1f} {:6.1f}'.format(crust1.min(), crust1.max()))
print ('litho1 {:6.1f} {:6.1f}\n'.format(litho1.min(), litho1.max()))
def mean_globe(z):
nlat = z.shape[0]
sintheta = np.sin( np.arange(0.5/nlat, 1, 1/nlat) * pi )
zz = np.zeros(z.shape)
for i in range(0,z.shape[1]):
zz[:,i] = z[:,i] * sintheta[:]
return np.sum(zz) / (zz.shape[1]*np.sum(sintheta))
def rms_globe(z):
nlat = z.shape[0]
zz = z*z
sintheta = np.sin( np.arange(0.5/nlat, 1, 1/nlat) * pi )
for i in range(0,z.shape[1]):
zz[:,i] = zz[:,i] * sintheta[:]
return sqrt(np.sum(zz) / (zz.shape[1]*np.sum(sintheta)))
print ('Differences between models (GEMMA-XXXXX)')
print (' MEAN RMS [km]')
print ('litho1 {:6.2f} {:6.2f}'.format(mean_globe(diff_litho1), rms_globe(diff_litho1)))
print ('crust1 {:6.2f} {:6.2f}'.format(mean_globe(diff_crust1), rms_globe(diff_crust1)))
print ('crust2 {:6.2f} {:6.2f}'.format(mean_globe(diff_crust2), rms_globe(diff_crust2)))
In [12]:
plt.rcParams['figure.figsize'] = (18,12)
nlat = gemma.shape[0]
sintheta = np.sin( np.arange(0.5/nlat, 1, 1/nlat) * pi )
weights = np.zeros(gemma.shape)
for i in range(0,gemma.shape[1]):
weights[:,i] = sintheta[:]
hist, bins_litho1 = np.histogram(diff_litho1.reshape(gemma.size), bins=40, weights=weights.reshape(gemma.size))
hist, bins_crust1 = np.histogram(diff_crust1.reshape(gemma.size), bins=40, weights=weights.reshape(gemma.size))
hist, bins_crust2 = np.histogram(diff_crust2.reshape(gemma.size), bins=40, weights=weights.reshape(gemma.size))
vmin = min(bins_litho1.min(), bins_crust1.min(), bins_crust2.min())
vmax = max(bins_litho1.max(), bins_crust1.max(), bins_crust2.max())
vrange = (vmin/2,vmax/2) # only to half of extreme values in order to see the interesting stuff
bins = 120
fig = plt.figure()
#
plt.subplot(3, 1, 1)
plt.hist(diff_litho1.reshape(gemma.size), bins=bins, range=vrange, weights=weights.reshape(gemma.size))
plt.title('Difference GEMMA1.0 – LITHO1.0')
plt.ylabel('#')
#
plt.subplot(3, 1, 2)
plt.hist(diff_crust1.reshape(gemma.size), bins=bins, range=vrange, weights=weights.reshape(gemma.size))
plt.title('Difference GEMMA1.0 – CRUST1.0')
plt.ylabel('#')
#
plt.subplot(3, 1, 3)
plt.hist(diff_crust2.reshape(gemma.size), bins=bins, range=vrange, weights=weights.reshape(gemma.size))
plt.title('Difference GEMMA1.0 – CRUST2.0')
plt.xlabel('km')
plt.ylabel('#')
#
plt.show()
