Commit 9d321342 authored by Mijian Xu's avatar Mijian Xu 😷
Browse files

update hk stacking

parent f9e2c3ef
...@@ -97,4 +97,4 @@ def ccppara(cfg_file): ...@@ -97,4 +97,4 @@ def ccppara(cfg_file):
stack_end = cf.getfloat('stack', 'stack_end') stack_end = cf.getfloat('stack', 'stack_end')
cpara.stack_val = cf.getfloat('stack', 'stack_val') cpara.stack_val = cf.getfloat('stack', 'stack_val')
cpara.stack_range = np.append(np.arange(stack_start, stack_end, cpara.stack_val), stack_end) cpara.stack_range = np.append(np.arange(stack_start, stack_end, cpara.stack_val), stack_end)
return cpara return cpara
\ No newline at end of file
...@@ -2,9 +2,13 @@ import numpy as np ...@@ -2,9 +2,13 @@ import numpy as np
import re import re
from obspy.io.sac.sactrace import SACTrace from obspy.io.sac.sactrace import SACTrace
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from os.path import dirname, join from os.path import dirname, join, basename
from scipy.io import loadmat from scipy.io import loadmat
from matplotlib.colors import ListedColormap from matplotlib.colors import ListedColormap
from seispy.rfcorrect import SACStation
from seispy.hkpara import hkpara
import argparse
def transarray(array, axis=0): def transarray(array, axis=0):
if not isinstance(array, np.ndarray): if not isinstance(array, np.ndarray):
...@@ -102,10 +106,12 @@ def load_cyan_map(): ...@@ -102,10 +106,12 @@ def load_cyan_map():
return ListedColormap(carray) return ListedColormap(carray)
def plot(stack, allstack, h, kappa, besth, bestk, cvalue, cmap=load_cyan_map()): def plot(stack, allstack, h, kappa, besth, bestk, cvalue, cmap=load_cyan_map(), title=None, path=None):
f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex='col', sharey='row') f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex='col', sharey='row')
xlim = (h[0], h[-1]) xlim = (h[0], h[-1])
ylim = (kappa[0], kappa[-1]) ylim = (kappa[0], kappa[-1])
if title is not None:
f.suptitle(title, fontsize='large')
ax1.imshow(stack[:, :, 0], cmap=cmap, extent=[xlim[0], xlim[1], ylim[0], ylim[1]], aspect='auto', origin='lower') 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_ylabel('$V_P/V_S$')
ax1.set_title('Ps') ax1.set_title('Ps')
...@@ -121,7 +127,10 @@ def plot(stack, allstack, h, kappa, besth, bestk, cvalue, cmap=load_cyan_map()): ...@@ -121,7 +127,10 @@ def plot(stack, allstack, h, kappa, besth, bestk, cvalue, cmap=load_cyan_map()):
ax4.plot(xlim, [bestk, bestk], color='red', linestyle='--', linewidth=0.6) ax4.plot(xlim, [bestk, bestk], color='red', linestyle='--', linewidth=0.6)
ax4.plot([besth, besth], ylim, color='red', linestyle='--', linewidth=0.6) ax4.plot([besth, besth], ylim, color='red', linestyle='--', linewidth=0.6)
ax4.set_xlabel('Moho depth (km)') ax4.set_xlabel('Moho depth (km)')
plt.show(f) if path is None:
f.show()
else:
f.savefig(path, format='pdf')
def ci(allstack, h, kappa, ev_num): def ci(allstack, h, kappa, ev_num):
...@@ -155,6 +164,45 @@ def print_result(besth, bestk, maxhsig, maxksig, print_comment=True): ...@@ -155,6 +164,45 @@ def print_result(besth, bestk, maxhsig, maxksig, print_comment=True):
print(msg) print(msg)
def hksta(hpara, isplot=False):
station = basename(hpara.rfpath)
stadata = SACStation(join(hpara.rfpath, station+'finallist.dat'), only_r=True)
stack, _, allstack, _ = hkstack(stadata.datar, stadata.shift, stadata.sampling, stadata.rayp,
hpara.hrange, hpara.krange, vp=hpara.vp)
besth, bestk, cvalue, maxhsig, maxksig = ci(allstack, hpara.hrange, hpara.krange, stadata.ev_num)
with open(hpara.hklist, 'a') as f:
f.write('{}\t{:.3f}\t{:.3f}\t{:.1f}\t{:.2f}\t{:.2f}\t{:.3f}\n'.format(station, stadata.stla, stadata.stlo,
besth, maxhsig, bestk, maxksig))
title = '{}\nMoho depth = ${:.1f}\pm{:.2f}$\nV_P/V_S = ${:.2f}\pm{:.3f}$'.format(station, besth,
maxhsig, bestk, maxksig)
if isplot:
img_path = join(hpara.hkpath, station+'.pdf')
plot(stack, allstack, hpara.hrange, hpara.krange, besth, bestk, cvalue, title=title, path=img_path)
else:
plot(stack, allstack, hpara.hrange, hpara.krange, besth, bestk, cvalue, title=title)
def hk():
parser = argparse.ArgumentParser(description="HK stacking for single station")
parser.add_argument('cfg_file', type=str, help='Path to HK configure file')
parser.add_argument('-s', help='Path to PRFs with folder name of the station name',
dest='station', type=str, default=None)
parser.add_argument('-H', help='Range for searching best H: <hmin>/<hmax>', type=str, default='')
parser.add_argument('-K', help='Range for searching best K: <kmin>/<kmax>', type=str, default='')
parser.add_argument('-p', help='If save the figure', dest='isplot', action='store_false')
arg = parser.parse_args()
hpara = hkpara(arg.cfg_file)
if arg.station is not None:
hpara.hkpath = arg.station
if arg.H != '':
hmin, hmax = [float(val) for val in arg.H.split('/')]
hpara.hrange = np.arange(hmin, hmax, 0.1)
if arg.K != '':
kmin, kmax = [float(val) for val in arg.K.split('/')]
hpara.krange = np.arange(kmin, kmax, 0.1)
hksta(hpara, isplot=arg.isplot)
def hktest(): def hktest():
h = np.arange(40, 80, 0.1) h = np.arange(40, 80, 0.1)
kappa = np.arange(1.6, 1.9, 0.01) kappa = np.arange(1.6, 1.9, 0.01)
......
from os.path import expanduser
import configparser
import numpy as np
class HKPara(object):
def __init__(self):
self.rfpath = expanduser('~')
self.hkpath = expanduser('~')
self.hklist = 'hk.dat'
self.hrange = np.arange(20, 80, 0.1)
self.krange = np.arange(1.6, 1.9, 0.01)
self.vp = 6.3
@property
def hrange(self):
return self._hrange
@hrange.setter
def hrange(self, value):
if not (isinstance(value, np.ndarray) or value is None):
raise TypeError('Error type of hrange')
else:
self._hrange = value
@property
def krange(self):
return self._krange
@krange.setter
def krange(self, value):
if not (isinstance(value, np.ndarray) or value is None):
raise TypeError('Error type of krange')
else:
self._krange = value
def hkpara(cfg_file):
hpara = HKPara()
cf = configparser.ConfigParser()
try:
cf.read(cfg_file)
except Exception:
raise FileNotFoundError('Cannot open configure file %s' % cfg_file)
# para for FileIO section
hpara.rfpath = cf.get('FileIO', 'rfpath')
hpara.hkpath = cf.get('FileIO', 'hkpath')
hpara.hklst = cf.get('FileIO', 'hklst')
hmin = cf.getfloat('hk', 'hmin')
hmax = cf.getfloat('hk', 'hmax')
kmin = cf.getfloat('hk', 'kmin')
kmax = cf.getfloat('hk', 'kmax')
hpara.hrange = np.arange(hmin, hmax, 0.1)
hpara.krange = np.arange(kmin, kmax, 0.01)
vp = cf.get('hk', 'vp')
if vp != '':
hpara.vp = float(vp)
return hpara
...@@ -26,6 +26,8 @@ class SACStation(object): ...@@ -26,6 +26,8 @@ class SACStation(object):
self.rayp = skm2srad(self.rayp) self.rayp = skm2srad(self.rayp)
self.ev_num = self.evla.shape[0] self.ev_num = self.evla.shape[0]
sample_sac = SACTrace.read(join(data_path, self.event[0] + '_' + self.phase[0] + '_R.sac')) sample_sac = SACTrace.read(join(data_path, self.event[0] + '_' + self.phase[0] + '_R.sac'))
self.stla = sample_sac.stla
self.stlo = sample_sac.stlo
self.RFlength = sample_sac.npts self.RFlength = sample_sac.npts
self.shift = -sample_sac.b self.shift = -sample_sac.b
self.sampling = sample_sac.delta self.sampling = sample_sac.delta
......
...@@ -19,7 +19,8 @@ setup(name='seispy', ...@@ -19,7 +19,8 @@ setup(name='seispy',
'updatecatalog=seispy.updatecatalog:main', 'updatecatalog=seispy.updatecatalog:main',
'ndk2dat=seispy.updatecatalog:ndk2dat', 'ndk2dat=seispy.updatecatalog:ndk2dat',
'lsnc=seispy.io:lsnc', 'lsnc=seispy.io:lsnc',
'ccp_profile=seispy.ccp:ccp_profile']}, 'ccp_profile=seispy.ccp:ccp_profile',
'hk=seispy.hk:hk']},
include_package_data=True, include_package_data=True,
zip_safe=False zip_safe=False
) )
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