diff --git a/seispy/ccp/rayplib.py b/seispy/ccp/rayplib.py new file mode 100644 index 00000000..f8f40b16 --- /dev/null +++ b/seispy/ccp/rayplib.py @@ -0,0 +1,95 @@ +""" +use obspy.taup instead of taup in java +""" + +import numpy as np + +from obspy.taup import TauPyModel + +def rayp_input(phase_main:str, + source_dist_in_degree:str, + conver_layers:str, + source_depth_in_km:str, + out_put_path:str): + """ + anaylize input informations + only allow P, SKS, S, ScS + """ + if phase_main not in ['P', 'S', 'SKS', 'ScS']: + raise ValueError('not a valid phase') + + try: + # set default step value for 2 input + source_dist_in_degree +='/1';source_dist_in_degree = [float(_i) for _i in source_dist_in_degree.split('/')[:4]] + conver_layers +='/1';conver_layers = [float(_i) for _i in conver_layers.split('/')[:4]] + source_depth_in_km +='/2';source_depth_in_km = [float(_i) for _i in source_depth_in_km.split('/')] + if conver_layers[0] < 0: + conver_layers[0] = 0 + conver_layers = np.arange(conver_layers[0], conver_layers[1], conver_layers[2]) + source_dist_in_degree = np.arange(source_dist_in_degree[0], source_dist_in_degree[1], source_dist_in_degree[2]) + source_depth_in_km = np.arange(source_depth_in_km[0],source_depth_in_km[1], source_depth_in_km[2]) + except: + raise ValueError('invalid inputs in dist and layers') + # here we limit layer_max to 30km + if conver_layers[-1] < 30: + raise ValueError('theres no need to do such a thing') + + + rayp = cal_taup('iasp91', phase_main, + conver_layers, source_dist_in_degree, source_depth_in_km) + + np.savez(out_put_path, dis = source_dist_in_degree, dep = source_depth_in_km, layers = conver_layers, rayp = rayp) + + + + +def cal_taup(model:str, phase:str, + layer:np.ndarray, dist:np.ndarray, depth:np.ndarray): + """ + cal rayp from obspy.taup + wont cal layers above 11km, change min_con + + phase_list: list of phase names, like ["P30s", "P50s"] + layer : conversion layers + dist : numpy.ndarray, dist in deg + depth : numpy.ndarray, source_depth in km + + """ + min_con = 11 + if phase[-1] == 'P': + conv = 's' + else: + conv = 'p' + ## init taup and matrix + rayp_lib = np.zeros((dist.shape[0], depth.shape[0], layer.shape[0])) + model = TauPyModel(model=model) + cal_layers_slice = np.where(layer > min_con); head = cal_layers_slice[0][0] + phase_list = ["{}{}{}".format(phase,_d, conv) for _d in layer[cal_layers_slice]] + + ### main + for _i, _deg in enumerate(dist): + for _j, _dep in enumerate(depth): + #print("now at degree {}, dep{}".format(_deg, _dep)) + for _k, _phase in enumerate(phase_list): + arrivals=model.get_travel_times(source_depth_in_km=_dep, + distance_in_degree=_deg, + phase_list=[_phase]) + if len(arrivals) == 0: + continue + rayp_lib[_i,_j,head+_k] = arrivals[0].ray_param_sec_degree + if head > 0: + for _layer in range(head): + rayp_lib[:,:,_layer] = rayp_lib[:,:,head] + return rayp_lib + + +if __name__ == "__main__": + # paras + phase = "S" + dist = "30/90/5" + layers = "0/100/2" + depth = '0/100/5' + out_put = "G:/Ps_rayp.npz" + + # run + rayp_input(phase, dist, layers, depth, out_put) diff --git a/seispy/rf2depth_makedata.py b/seispy/rf2depth_makedata.py index 2a45b4e5..e56b5ef0 100644 --- a/seispy/rf2depth_makedata.py +++ b/seispy/rf2depth_makedata.py @@ -1,4 +1,18 @@ -from seispy.core.depmodel import _load_mod +""" +class RF2depth : process cal RF2depth +class sta_part, sta_full, _RFInd + +""" +from os.path import join, exists +import argparse +import sys +import glob +from collections import namedtuple +from logging import Logger + +import numpy as np + +from seispy.core.depmodel import _load_mod, DepModel from seispy.rfcorrect import RFStation, psrf2depth, psrf_1D_raytracing,\ psrf_3D_migration, time2depth, psrf_3D_raytracing from seispy.core.pertmod import Mod3DPerturbation @@ -6,10 +20,30 @@ from seispy.ccppara import ccppara, CCPPara from seispy.setuplog import SetupLog from seispy.geo import latlon_from, rad2deg -from os.path import join, exists -import argparse -import sys -import glob + +# sta_part is not used +sta_part = namedtuple('sta_part', ['station', 'stla', 'stlo']) +sta_full = namedtuple('sta_full', ['station', 'stla', 'stlo', 'stel']) +_RfInDepth = namedtuple('RfInDepth', + ['station', 'stalat', 'stalon', + 'depthrange', 'bazi', 'rayp', + 'moveout_correct', 'piercelat', 'piercelon', 'stopindex']) +class RFInDepth(_RfInDepth): + """ + this class is inheried from a namedtuple, and can be called as a dict + mainly to restrict behaviour while output and + provide a better speed while + + warnings: this class is still in develpment + """ + def __getitem__(self, item): + """ + make RFInDepth can be called like dict + >>> ll = RFInDepth(station='sta', stalat='stla', stalon='stlo', depthrange='dep', bazi='baz',rayp='rayp', moveout_correct=0, piercelat=0, piercelon=0 ) + >>> ll['stalat'], ll.stalat + ('stla', 'stla') + """ + return getattr(self, item) class Station(object): @@ -28,20 +62,43 @@ def __init__(self, sta_lst:str): 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] + def __getitem__(self,index): + """ + allows [] do indexing + """ + if hasattr(self,'stel'): + return sta_full(self.station[index], + self.stla[index], self.stlo[index], self.stel[index]) + else: + return sta_full(self.station[index], + self.stla[index], self.stlo[index], 0) + def __iter__(self): + """ + allow for...in to iter + """ + for index in range(self.stla.shape[0]): + if hasattr(self, 'stel'): + yield sta_full(self.station[index], + self.stla[index], self.stlo[index], self.stel[index]) + else: + yield sta_full(self.station[index], self.stla[index], self.stlo[index], 0) + def __len__(self): + return self.station.shape[0] + class RFDepth(): """Convert receiver function to depth axis """ - def __init__(self, cpara:CCPPara, log=SetupLog(), - raytracing3d=False, velmod3d=None, modfolder1d=None) -> None: + def __init__(self, cpara:CCPPara, log:Logger=SetupLog().RF2depthlog, + 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 + :type log: logging.Logger + :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 @@ -51,22 +108,22 @@ def __init__(self, cpara:CCPPara, log=SetupLog(), self.cpara = cpara self.modfolder1d = modfolder1d self.log = log - self.raytracing3d = raytracing3d + self.RAYTRACING3D = RAYTRACING3D if velmod3d is not None: if isinstance(velmod3d, str): self.mod3d = Mod3DPerturbation(velmod3d, cpara.depth_axis, velmod=cpara.velmod) else: - log.RF2depthlog.error('Path to 3d velocity model should be in str') + log.error('Path to 3d velocity model should be in str') sys.exit(1) elif modfolder1d is not None: if isinstance(modfolder1d, str): if exists(modfolder1d): self.ismod1d = True else: - log.RF2depthlog.error('No such folder of {}'.format(modfolder1d)) + log.error('No such folder of {}'.format(modfolder1d)) sys.exit(1) else: - log.RF2depthlog.error('Folder to 1d velocity model files should be in str') + log.error('Folder to 1d velocity model files should be in str') sys.exit(1) else: self.ismod1d = True @@ -74,6 +131,7 @@ def __init__(self, cpara:CCPPara, log=SetupLog(), self.srayp = np.load(cpara.rayp_lib) else: self.srayp = None + self.sta_info = Station(cpara.stalist) self.rfdepth = [] self._test_comp() @@ -86,7 +144,7 @@ def _test_comp(self): self.prime_comp = comp break if not self.prime_comp: - self.log.RF2depthlog.error('No such any RF files in \'R\',' + self.log.error('No such any RF files in \'R\',' '\'Q\', \'L\', and \'Z\' components') sys.exit(1) @@ -96,29 +154,38 @@ def makedata(self, psphase=1): :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]) + for _i, _sta in enumerate(self.sta_info): + rfpath = join(self.cpara.rfpath, _sta.station) stadatar = RFStation(rfpath, only_r=True, prime_comp=self.prime_comp) - stadatar.stel = self.sta_info.stel[i] - stadatar.stla = self.sta_info.stla[i] - stadatar.stlo = self.sta_info.stlo[i] + stadatar.stel = _sta.stel + stadatar.stla = _sta.stla + stadatar.stlo = _sta.stlo if stadatar.prime_phase == 'P': sphere = True else: sphere = False - self.log.RF2depthlog.info('the {}th/{} station with {} events'.format(i + 1, self.sta_info.stla.shape[0], stadatar.ev_num)) + self.log.info('the {}th/{} station with {} events'.format(_i + 1, self.sta_info.shape[0], stadatar.ev_num)) + + + #### 1d model for each station if self.ismod1d: if self.modfolder1d is not None: - velmod = _load_mod(self.modfolder1d, self.sta_info.station[i]) + velmod = _load_mod(self.modfolder1d, _sta.station) + 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=psphase) + velmod=velmod, srayp=self.cpara.rayp_lib, + sphere=sphere, phase=psphase) + piercelat, piercelon = np.zeros_like(x_s, dtype=np.float64), np.zeros_like(x_s, dtype=np.float64) + for j in range(stadatar.ev_num): piercelat[j], piercelon[j] = latlon_from(self.sta_info.stla[i], self.sta_info.stlo[i], stadatar.bazi[j], rad2deg(x_s[j])) else: + ### 3d model interp 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: @@ -130,59 +197,81 @@ def makedata(self, psphase=1): piercelat, piercelon = pplat_s, pplon_s else: piercelat, piercelon = pplat_p, pplon_p + ps_rfdepth, end_index = time2depth(stadatar, self.cpara.depth_axis, newtpds) - rfdep = self._write_rfdep(stadatar, ps_rfdepth, piercelat, piercelon, end_index) - self.rfdepth.append(rfdep) + + rfdep = self._append_rfdep(self.cpara,stadatar, ps_rfdepth, piercelat, piercelon, end_index) np.save(self.cpara.depthdat, self.rfdepth) - def _write_rfdep(self, stadata, amp, pplat, pplon, end_index): + def _write_rfdep(self,cpara, stadata, amp, pplat, pplon, end_index): rfdep = {} rfdep['station'] = stadata.staname rfdep['stalat'] = stadata.stla rfdep['stalon'] = stadata.stlo - rfdep['depthrange'] = self.cpara.depth_axis + rfdep['depthrange'] = cpara.depth_axis rfdep['bazi'] = stadata.bazi rfdep['rayp'] = stadata.rayp rfdep['moveout_correct'] = amp rfdep['piercelat'] = pplat rfdep['piercelon'] = pplon rfdep['stopindex'] = end_index - return rfdep + self.rfdepth.append(rfdep) + def rf2depth(): - """Convert receiver function to depth axis + """ + CLI for Convert receiver function to depth axis + There's 4 branch provided to do RF 2 depth conversion + + 1. only -d :do moveout correction + 2. only -r : do raytracing but no moveout correction + 3. -d and -r : do moveout correction and raytracing + 4. -m : use {staname}.vel file for RF2depth conversion + + """ 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', + 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', + 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', + 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') + parser.add_argument('cfg_file', + help='Path to configure file', + metavar='ccp.cfg', type=str) arg = parser.parse_args() cpara = ccppara(arg.cfg_file) if arg.d != '' and arg.r != '': + #### print help issue raise ValueError('Specify only 1 argument in \'-d\' and \'-r\'') elif arg.d != '' and arg.r == '' and arg.m == '': + #### do 3d moveout correction but 1d rf2depth raytracing3d = False velmod3d = arg.d modfolder1d = None elif arg.d == '' and arg.r != '' and arg.m == '': + #### do 3d raytraying raytracing3d = True velmod3d = arg.d modfolder1d = None elif arg.d == '' and arg.r == '' and arg.m != '': + #### use multiple 1d velmod for time2depth convertion raytracing3d = False velmod3d = None modfolder1d = arg.m else: + ### all last, use default 1D vel model do time2depth conversion raytracing3d = False velmod3d = None modfolder1d = None rfd = RFDepth( - cpara, raytracing3d=raytracing3d, + cpara, + RAYTRACING3D=raytracing3d, velmod3d=velmod3d, modfolder1d=modfolder1d, )