Skip to content

Commit

Permalink
debug PSRF2depth
Browse files Browse the repository at this point in the history
  • Loading branch information
@xumi1993
xumi1993 committed Dec 21, 2023
1 parent bdcb488 commit 919f428
Show file tree
Hide file tree
Showing 3 changed files with 44 additions and 36 deletions.
24 changes: 11 additions & 13 deletions seispy/para.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -238,40 +238,40 @@ def read_para(cls, cfg_file):
elif key == 'catalogpath':
pa.catalogpath = value
else:
pa.__dict__[key] = value
exec('pa.{} = {}'.format(key, value))
sections.remove('path')
if 'fetch' in sections:
for key, value in cf.items('fetch'):
if key == 'use_remote_data':
pa.__dict__[key] = cf.getboolean('fetch', 'use_remote_data')
pa.use_remote_data = cf.getboolean('fetch', 'use_remote_data')
elif key == 'data_server_user':
if value == '':
continue
else:
pa.__dict__[key] = value
pa.data_server_user = value
elif key == 'data_server_password':
if value == '':
continue
else:
pa.__dict__[key] = value
pa.data_server_password = value
else:
pa.stainfo.__dict__[key] = value
exec('pa.stainfo.{} = {}'.format(key, value))
sections.remove('fetch')
for sec in sections:
for key, value in cf.items(sec):
if key == 'date_begin':
pa.__dict__[key] = obspy.UTCDateTime(value)
pa.date_begin = obspy.UTCDateTime(value)
elif key == 'date_end':
pa.__dict__[key] = obspy.UTCDateTime(value)
pa.date_end = obspy.UTCDateTime(value)
elif key == 'offset':
try:
pa.__dict__[key] = float(value)
pa.offset = cf.getfloat(sec, 'offset')
except:
pa.__dict__[key] = None
pa.offset = None
elif key == 'itmax':
pa.__dict__[key] = int(value)
pa.itmax = int(value)
elif key == 'only_r':
pa.__dict__[key] = cf.getboolean(sec, 'only_r')
pa.only_r = cf.getboolean(sec, 'only_r')
elif key == 'criterion':
pa.criterion = value
elif key == 'decon_method':
Expand All @@ -289,8 +289,6 @@ def read_para(cls, cfg_file):
else:
try:
exec('pa.{} = {}'.format(key, float(value)))
# pa.__dict__[key] = float(value)
except ValueError:
exec('pa.{} = value'.format(key))
# pa.__dict__[key] = value
return pa
32 changes: 28 additions & 4 deletions seispy/rf2depth_makedata.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,9 +12,12 @@


class Station(object):
def __init__(self, sta_lst):
def __init__(self, sta_lst:str):
"""
Read station list
:param sta_lst: Path to station list
:type sta_lst: string
"""
dtype = {'names': ('station', 'stla', 'stlo', 'stel'), 'formats': ('U20', 'f4', 'f4', 'f2')}
try:
Expand Down Expand Up @@ -46,8 +49,22 @@ def _load_mod(datapath, staname):


class RFDepth():
"""Convert receiver function to depth axis
"""
def __init__(self, cpara:CCPPara, log=setuplog(),
raytracing3d=False, velmod3d=None, modfolder1d=None) -> None:
"""
:param cpara: CCPPara object
:type cpara: CCPPara
:param log: Log object
:type log: setuplog
:param raytracing3d: If True, use 3D ray tracing to calculate the travel time
:type raytracing3d: bool
:param velmod3d: Path to 3D velocity model in npz file
:type velmod3d: str
:param modfolder1d: Folder path to 1D velocity model files with staname.vel as the file name
:type modfolder1d: str
"""
self.ismod1d = False
self.cpara = cpara
self.modfolder1d = modfolder1d
Expand Down Expand Up @@ -91,7 +108,12 @@ def _test_comp(self):
'\'Q\', \'L\', and \'Z\' components')
sys.exit(1)

def makedata(self):
def makedata(self, psphase=1):
"""Convert receiver function to depth axis
:param psphase: 1 for Ps, 2 for PpPs, 3 for PpSs
:type psphase: int
"""
for i in range(self.sta_info.stla.shape[0]):
rfpath = join(self.cpara.rfpath, self.sta_info.station[i])
stadatar = RFStation(rfpath, only_r=True, prime_comp=self.prime_comp)
Expand All @@ -109,15 +131,15 @@ def makedata(self):
else:
velmod = self.cpara.velmod
ps_rfdepth, end_index, x_s, _ = psrf2depth(stadatar, self.cpara.depth_axis,
velmod=velmod, srayp=self.cpara.rayp_lib, sphere=sphere, phase=stadatar.prime_phase)
velmod=velmod, srayp=self.cpara.rayp_lib, sphere=sphere, phase=psphase)
piercelat, piercelon = latlon_from(self.sta_info.stla[i], self.sta_info.stlo[i],
stadatar.bazi, rad2deg(x_s))
else:
if self.raytracing3d:
pplat_s, pplon_s, pplat_p, pplon_p, newtpds = psrf_3D_raytracing(stadatar, self.cpara.depth_axis, self.mod3d, srayp=self.srayp, sphere=sphere)
else:
pplat_s, pplon_s, pplat_p, pplon_p, raylength_s, raylength_p, tps = psrf_1D_raytracing(
stadatar, self.cpara.depth_axis, srayp=self.srayp, sphere=sphere, phase=stadatar.prime_phase)
stadatar, self.cpara.depth_axis, srayp=self.srayp, sphere=sphere, phase=psphase)
newtpds = psrf_3D_migration(pplat_s, pplon_s, pplat_p, pplon_p, raylength_s, raylength_p,
tps, self.cpara.depth_axis, self.mod3d)
if stadatar.prime_phase == 'P':
Expand Down Expand Up @@ -145,6 +167,8 @@ def _write_rfdep(self, stadata, amp, pplat, pplon, end_index):


def rf2depth():
"""Convert receiver function to depth axis
"""
parser = argparse.ArgumentParser(description="Convert Ps RF to depth axis")
parser.add_argument('-d', help='Path to 3d vel model in npz file for moveout correcting',
metavar='3d_velmodel_path', type=str, default='')
Expand Down
24 changes: 5 additions & 19 deletions seispy/rfcorrect.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -278,17 +278,9 @@ def moveoutcorrect(self, ref_rayp=0.06, dep_range=np.arange(0, 150), velmod='ias
t_corr, _ = moveoutcorrect_ref(self, skm2srad(ref_rayp), dep_range, chan='t', velmod=velmod, **kwargs)
else:
t_corr = None
if 'datar' in self.__dict__:
chan = 'r'
elif 'dataz' in self.__dict__:
chan = 'z'
elif 'datal' in self.__dict__:
chan = 'l'
else:
pass
rf_corr, _ = moveoutcorrect_ref(self, skm2srad(ref_rayp), dep_range, chan=chan, velmod=velmod, **kwargs)
rf_corr, _ = moveoutcorrect_ref(self, skm2srad(ref_rayp), dep_range, chan='', velmod=velmod, **kwargs)
if replace:
self.__dict__['data{}'.format(chan)] = rf_corr
self.data_prime= rf_corr
if not self.only_r:
self.datat = t_corr
else:
Expand Down Expand Up @@ -434,7 +426,7 @@ def _imag2nan(arr):


def moveoutcorrect_ref(stadatar, raypref, YAxisRange,
chan='r', velmod='iasp91', sphere=True, phase=1):
chan='', velmod='iasp91', sphere=True, phase=1):
"""Moveout correction refer to a specified ray-parameter
:param stadatar: data class of RFStation
Expand All @@ -447,16 +439,10 @@ def moveoutcorrect_ref(stadatar, raypref, YAxisRange,
"""
sampling = stadatar.sampling
shift = stadatar.shift
if chan == 'r':
data = stadatar.datar
elif chan == 't':
if chan == 't':
data = stadatar.datat
elif chan == 'z':
data = stadatar.dataz
elif chan == 'l':
data = stadatar.datal
else:
raise ValueError('Field \'datar\' or \'datal\' must be in the SACStation')
data = stadatar.data_prime
dep_mod = DepModel(YAxisRange, velmod, stadatar.stel)
# x_s = np.zeros([stadatar.ev_num, YAxisRange.shape[0]])
# x_p = np.zeros([stadatar.ev_num, YAxisRange.shape[0]])
Expand Down

0 comments on commit 919f428

Please sign in to comment.