为rayplib添加了命令行工具，

rf2depth的一点优化，使用namedtuple改善了Station类的使用
depmodel.py 中的ccp_model太繁琐，未来会拆分重写。

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)

seispy/rf2depth_makedata.py

@@ -1,15 +1,49 @@
-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
 import numpy as np
 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,29 +108,30 @@ 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
         if cpara.rayp_lib is not None:
             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,
     )