from seispy.rfcorrect import RFStation, psrf2depth, Mod3DPerturbation, psrf_1D_raytracing,\
psrf_3D_migration, time2depth, psrf_3D_raytracing
import numpy as np
from seispy.ccppara import ccppara
from seispy.setuplog import setuplog
from seispy.geo import latlon_from, deg2km, rad2deg
from os.path import join, dirname, exists
import argparse
import sys
import glob
[docs]class Station(object):
def __init__(self, sta_lst):
"""
Read station list
"""
dtype = {'names': ('station', 'stla', 'stlo', 'stel'), 'formats': ('U20', 'f4', 'f4', 'f2')}
try:
self.station, self.stla, self.stlo, self.stel = np.loadtxt(sta_lst, dtype=dtype, unpack=True, ndmin=1)
except:
dtype = {'names': ('station', 'stla', 'stlo'), 'formats': ('U20', 'f4', 'f4')}
self.station, self.stla, self.stlo = np.loadtxt(sta_lst, dtype=dtype, unpack=True, ndmin=1)
self.stel = np.zeros(self.stla.size)
self.sta_num = self.stla.shape[0]
[docs]def init_mat(sta_num):
dtype = {'names': ('Station', 'stalat', 'stalon', 'Depthrange', 'events', 'bazi', 'rayp', 'phases', 'moveout_correct',
'Piercelat', 'Piercelon', 'StopIndex'),
'formats': tuple(['O']*12)}
return np.zeros([1, sta_num], dtype=dtype)
def _load_mod(datapath, staname):
"""Load 1D velocity model files with suffix of ".vel". The model file should be including 3 columns with depth, vp and vs.
:param datapath: Folder name with 1D velocity model files.
:type datapath: string
:param staname: The station name as a part of file name of 1D velocity model files.
:type staname: string
"""
expresion = join(datapath, "*"+staname+"*.vel")
modfiles = glob.glob(expresion)
if len(modfiles) == 0:
raise FileNotFoundError("The model file of {} were not found.".format(expresion))
elif len(modfiles) > 1:
raise ValueError('More then 1 file were found as the expresion: {}'.format(expresion))
else:
return modfiles[0]
[docs]def makedata(cpara, velmod3d=None, modfolder1d=None, log=setuplog()):
ismod1d = False
if velmod3d is not None:
if isinstance(velmod3d, str):
velmod = velmod3d
else:
raise ValueError('Path to 3d velocity model should be in str')
elif modfolder1d is not None:
if isinstance(modfolder1d, str):
if exists(modfolder1d):
ismod1d = True
else:
raise FileNotFoundError('No such folder of {}'.format(modfolder1d))
else:
ValueError('Path to 1d velocity model files should be in str')
else:
ismod1d = True
# cpara = ccppara(cfg_file)
sta_info = Station(cpara.stalist)
RFdepth = []
for i in range(sta_info.stla.shape[0]):
rfdep = {}
rfpath = join(cpara.rfpath, sta_info.station[i])
stadatar = RFStation(rfpath, only_r=True)
stadatar.stel = sta_info.stel[i]
stadatar.stla = sta_info.stla[i]
stadatar.stlo = sta_info.stlo[i]
log.RF2depthlog.info('the {}th/{} station with {} events'.format(i + 1, sta_info.stla.shape[0], stadatar.ev_num))
piercelat = np.zeros([stadatar.ev_num, cpara.depth_axis.shape[0]])
piercelon = np.zeros([stadatar.ev_num, cpara.depth_axis.shape[0]])
if stadatar.prime_phase == 'P':
sphere = True
else:
sphere = False
if ismod1d:
if modfolder1d is not None:
velmod = _load_mod(modfolder1d, sta_info.station[i])
else:
velmod = cpara.velmod
PS_RFdepth, end_index, x_s, _ = psrf2depth(stadatar, cpara.depth_axis,
velmod=velmod, srayp=cpara.rayp_lib, sphere=sphere, phase=cpara.phase)
for j in range(stadatar.ev_num):
piercelat[j], piercelon[j] = latlon_from(sta_info.stla[i], sta_info.stlo[i],
stadatar.bazi[j], rad2deg(x_s[j]))
rfdep['station'] = sta_info.station[i]
rfdep['stalat'] = sta_info.stla[i]
rfdep['stalon'] = sta_info.stlo[i]
rfdep['depthrange'] = cpara.depth_axis
# rfdep['events'] = _convert_str_mat(stadatar.event)
rfdep['bazi'] = stadatar.bazi
rfdep['rayp'] = stadatar.rayp
# rfdep['phases'] = stadatar.phase[i]
rfdep['moveout_correct'] = PS_RFdepth
rfdep['piercelat'] = piercelat
rfdep['piercelon'] = piercelon
rfdep['stopindex'] = end_index
RFdepth.append(rfdep)
# savemat(cpara.depthdat, {'RFdepth': RFdepth})
np.save(cpara.depthdat, RFdepth)
[docs]def makedata3d(cpara, velmod3d, log=setuplog(), raytracing3d=True):
mod3d = Mod3DPerturbation(velmod3d, cpara.depth_axis, velmod=cpara.velmod)
sta_info = Station(cpara.stalist)
if cpara.rayp_lib is not None:
srayp = np.load(cpara.rayp_lib)
else:
srayp = None
RFdepth = []
for i in range(sta_info.stla.shape[0]):
rfdep = {}
rfpath = join(cpara.rfpath, sta_info.station[i])
stadatar = RFStation(rfpath, only_r=True)
stadatar.stel = sta_info.stel[i]
stadatar.stla = sta_info.stla[i]
stadatar.stlo = sta_info.stlo[i]
if stadatar.prime_phase == 'P':
sphere = True
else:
sphere = False
log.RF2depthlog.info('the {}th/{} station with {} events'.format(i + 1, sta_info.stla.shape[0], stadatar.ev_num))
if raytracing3d:
pplat_s, pplon_s, pplat_p, pplon_p, newtpds = psrf_3D_raytracing(stadatar, cpara.depth_axis, mod3d, srayp=srayp, sphere=sphere)
else:
pplat_s, pplon_s, pplat_p, pplon_p, raylength_s, raylength_p, tps = psrf_1D_raytracing(
stadatar, cpara.depth_axis, srayp=srayp, sphere=sphere, phase=cpara.phase)
newtpds = psrf_3D_migration(pplat_s, pplon_s, pplat_p, pplon_p, raylength_s, raylength_p,
tps, cpara.depth_axis, mod3d)
amp3d, end_index = time2depth(stadatar, cpara.depth_axis, newtpds)
rfdep['station'] = sta_info.station[i]
rfdep['stalat'] = sta_info.stla[i]
rfdep['stalon'] = sta_info.stlo[i]
rfdep['depthrange'] = cpara.depth_axis
# rfdep['events'] = _convert_str_mat(stadatar.event)
rfdep['bazi'] = stadatar.bazi
rfdep['rayp'] = stadatar.rayp
# rfdep['phases'] = _convert_str_mat(stadatar.phase)
rfdep['moveout_correct'] = amp3d
rfdep['piercelat'] = pplat_s
rfdep['piercelon'] = pplon_s
rfdep['stopindex'] = end_index
RFdepth.append(rfdep)
np.save(cpara.depthdat, RFdepth)
[docs]def rf2depth():
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='')
parser.add_argument('-m', help='Folder path to 1d vel model files with staname.vel as the file name',
metavar='1d_velmodel_folder', type=str, default='')
parser.add_argument('-r', help='Path to 3d vel model in npz file for 3D ray tracing',
metavar='3d_velmodel_path', type=str, default='')
parser.add_argument('cfg_file', type=str, help='Path to configure file')
arg = parser.parse_args()
if len(sys.argv) == 1:
parser.print_help()
sys.exit(1)
cpara = ccppara(arg.cfg_file)
if arg.d != '' and arg.r != '':
raise ValueError('Specify only 1 argument in \'-d\' and \'-r\'')
elif arg.d != '' and arg.r == '' and arg.m == '':
makedata3d(cpara, arg.d, raytracing3d=False)
elif arg.d == '' and arg.r != '' and arg.m == '':
makedata3d(cpara, arg.r, raytracing3d=True)
elif arg.d == '' and arg.r == '' and arg.m != '':
makedata(cpara, modfolder1d=arg.m)
else:
makedata(cpara)
if __name__ == '__main__':
cfg_file = '/Users/xumj/Researches/Tibet_MTZ/process/paraCCP.cfg'
vel3d_file = '/Users/xumj/Researches/Tibet_MTZ/models/GYPSUM.npz'
cpara = ccppara(cfg_file)
makedata3d(cpara, vel3d_file)