Commit f9e2c3ef authored by Mijian Xu's avatar Mijian Xu 😷
Browse files

update hk stacking

parent 71e68fd5
import numpy as np
from os.path import join
import re
from obspy.io.sac.sactrace import SACTrace
import matplotlib.pyplot as plt
from os.path import dirname, join
from scipy.io import loadmat
from matplotlib.colors import ListedColormap
def transarray(array, axis=0):
if not isinstance(array, np.ndarray):
......@@ -89,13 +91,74 @@ def hkstack(seis, t0, dt, p, h, kappa, vp=6.3, weight=(0.7, 0.2, 0.1)):
allstackvar = np.var(allstack, axis=2)
allstack = np.mean(allstack, axis=2)
return stack, stackvar, allstack, allstackvar
Normed_stack = allstack - np.min(allstack)
Normed_stack = Normed_stack / np.max(Normed_stack)
return stack, stackvar, Normed_stack, allstackvar
def load_cyan_map():
path = join(dirname(__file__), 'data', 'cyan.mat')
carray = loadmat(path)['cyan']
return ListedColormap(carray)
def plot(stack, allstack, h, kappa, besth, bestk, cvalue, cmap=load_cyan_map()):
f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex='col', sharey='row')
xlim = (h[0], h[-1])
ylim = (kappa[0], kappa[-1])
ax1.imshow(stack[:, :, 0], cmap=cmap, extent=[xlim[0], xlim[1], ylim[0], ylim[1]], aspect='auto', origin='lower')
ax1.set_ylabel('$V_P/V_S$')
ax1.set_title('Ps')
ax2.imshow(stack[:, :, 1], cmap=cmap, extent=[xlim[0], xlim[1], ylim[0], ylim[1]], aspect='auto', origin='lower')
ax2.set_title('PpPs')
ax3.imshow(stack[:, :, 2], cmap=cmap, extent=[xlim[0], xlim[1], ylim[0], ylim[1]], aspect='auto', origin='lower')
ax3.set_title('PsPs+PpSs')
ax3.set_xlabel('Moho depth (km)')
ax3.set_ylabel('$V_P/V_S$')
ax4.imshow(allstack, cmap=cmap, extent=[xlim[0], xlim[1], ylim[0], ylim[1]], aspect='auto', origin='lower')
ax4.plot(besth, bestk, color='red', marker='s', markerfacecolor='none')
ax4.contour(allstack, [cvalue, 1], colors='k', extent=[xlim[0], xlim[1], ylim[0], ylim[1]], origin='lower')
ax4.plot(xlim, [bestk, bestk], color='red', linestyle='--', linewidth=0.6)
ax4.plot([besth, besth], ylim, color='red', linestyle='--', linewidth=0.6)
ax4.set_xlabel('Moho depth (km)')
plt.show(f)
def ci(allstack, h, kappa, ev_num):
"""
Search best H and kappa from stacked matrix.
Calculate error for H and kappa
:param allstack: stacked HK matrix
:param h: 1-D array of H
:param kappa: 1-D array of kappa
:param ev_num: event number
:return:
"""
[i, j] = np.unravel_index(allstack.argmax(), allstack.shape)
bestk = kappa[i]
besth = h[j]
cvalue = 1 - np.std(allstack.reshape(allstack.size)) / np.sqrt(ev_num)
cs = plt.contour(h, kappa, allstack, [cvalue])
cs_path = cs.collections[0].get_paths()[0].vertices
maxhsig = (np.max(cs_path[:, 0]) - np.min(cs_path[:, 0])) / 2
maxksig = (np.max(cs_path[:, 1]) - np.min(cs_path[:, 1])) / 2
return besth, bestk, cvalue, maxhsig, maxksig
def print_result(besth, bestk, maxhsig, maxksig, print_comment=True):
header = 'H\tH_error\tk\tk_error\n'
if print_comment:
msg = '{}{:.1f}\t{:.2f}\t{:.2f}\t{:.2f}'.format(header, besth, maxhsig, bestk, maxksig)
else:
msg = '{:.1f}\t{:.2f}\t{:.2f}\t{:.2f}'.format(besth, maxhsig, bestk, maxksig)
print(msg)
def hktest():
h = np.arange(30, 70, 1)
h = np.arange(40, 80, 0.1)
kappa = np.arange(1.6, 1.9, 0.01)
path = '/Users/xumj/Researches/YN_crust/dip_study/syn_hk/benchmark'
path = '/Users/xumj/Researches/YN_crust/dip_study/syn_hk/A1'
baz, rayp = np.loadtxt(join(path, 'sample.geom'), usecols=(0, 1), unpack=True)
rayp *= 1000
with open(join(path, 'sample.tr')) as f:
......@@ -110,12 +173,11 @@ def hktest():
sac = SACTrace.read(join(path, 'tr_{:d}_{:.3f}.r'.format(int(b), r)))
seis[_i] = sac.data
_, _, allstack, _ = hkstack(seis, shift, dt, rayp, h, kappa, vp=6.54)
[i, j] = np.unravel_index(allstack.argmax(), allstack.shape)
print('Moho depth = {:.0f}\nkappa = {:.2f}'.format(h[j], kappa[i]))
stack, _, allstack, _ = hkstack(seis, shift, dt, rayp, h, kappa, vp=6.3)
besth, bestk, cvalue, maxhsig, maxksig = ci(allstack, h, kappa, ev_num)
print_result(besth, bestk, maxhsig, maxksig, print_comment=True)
plot(stack, allstack, h, kappa, besth, bestk, cvalue)
if __name__ == '__main__':
hktest()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment