Merge branch 'dev' of https://github.com/xumi1993/seispy into dev

seispy/catalog.py

@@ -5,12 +5,32 @@
 import pandas as pd
 
 def download_catalog(fname, server='IRIS', format='seispy', **kwargs):
+    """
+    Download catalog to local file
+    :param fname: file name
+    :type fname: str
+    :param server: service provider, defaults to IRIS
+    :type server: str
+    :param format: file format, defaults to seispy
+    :type format: str
+    :param kwargs: parameters for Query.get_events
+    :type kwargs: dict
+    """
     query = Query(server=server)
     query.get_events(**kwargs)
     write_catalog(query, fname, format)
 
 
 def write_catalog(query, fname, format):
+    """
+    Write catalog to file
+    :param query: Query object
+    :type query: seispy.io.Query
+    :param fname: file name
+    :type fname: str
+    :param format: file format
+    :type format: str
+    """
     if format == 'seispy':
         with open(fname, 'w') as f:
             for i, row in query.events.iterrows():
@@ -23,6 +43,12 @@ def write_catalog(query, fname, format):
 
 
 def read_catalog_file(fname):
+    """
+    Read catalog file
+    :param fname: catalog file name
+    :type fname: str
+    :return: pandas.DataFrame
+    """
     col = ['date', 'evla', 'evlo', 'evdp', 'mag']
     colcata = ['year', 'month', 'day', 'hour', 'minute', 'second']
     data_cata = pd.read_table(fname, header=None, sep='\s+')

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/ccppara.py

@@ -61,6 +61,12 @@ def shape(self, value):
 
 
 def ccppara(cfg_file):
+    """ Read configure file for CCP stacking
+    :param cfg_file: Path to configure file
+    :type cfg_file: str
+    :return: CCPPara object
+    :rtype: CCPPara
+    """
     cpara = CCPPara()
     cf = configparser.ConfigParser()
     try:

seispy/ccpprofile.py

@@ -15,29 +15,38 @@
 import sys
 
 
-def prof_range(lat, lon):
+def _prof_range(lat, lon):
     dis = [0]
     for i in range(lat.size-1):
         dis.append(distaz(lat[i], lon[i], lat[i+1], lon[i+1]).degreesToKilometers())
     return np.cumsum(dis)
 
 
 def create_center_bin_profile(stations, val=5, method='linear'):
+    """Create bins along a profile with given stations
+    :param stations: Stations along a profile
+    :type stations: :class:`seispy.rf2depth_makedata.Station`
+    :param val: The interval between two points in km
+    :type val: float
+    :param method: Method for interpolation
+    :type method: str
+    :return: The location of bins (bin_loca), and length between each bin and the start point (profile_range)
+    :rtype: (numpy.array, numpy.array)
+    """
     if not isinstance(stations, Station):
         raise TypeError('Stations should be seispy.rf2depth_makedata.Station')
-    dis_sta = prof_range(stations.stla, stations.stlo)
+    dis_sta = _prof_range(stations.stla, stations.stlo)
     dis_inter = np.append(np.arange(0, dis_sta[-1], val), dis_sta[-1])
     r, theta, phi = geo2sph(np.zeros(stations.stla.size), stations.stla, stations.stlo)
     # t_po = np.arange(stations.stla.size)
     # ip_t_po = np.linspace(0, stations.stla.size, bin_num)
     theta_i = interp1d(dis_sta, theta, kind=method, bounds_error=False, fill_value='extrapolate')(dis_inter)
     phi_i = interp1d(dis_sta, phi, kind=method, bounds_error=False, fill_value='extrapolate')(dis_inter)
     _, lat, lon = sph2geo(r, theta_i, phi_i)
-    # dis = prof_range(lat, lon)
     return lat, lon, dis_inter
 
 
-def line_proj(lat1, lon1, lat2, lon2):
+def _line_proj(lat1, lon1, lat2, lon2):
     daz = distaz(lat1, lon1, lat2, lon2)
     az1_begin = (daz.baz - 90) % 360
     az2_begin = (daz.baz + 90) % 360
@@ -46,7 +55,7 @@ def line_proj(lat1, lon1, lat2, lon2):
     return az1_begin, az2_begin, az1_end, az2_end, daz
 
 
-def fix_filename(filename, typ='dat'):
+def _fix_filename(filename, typ='dat'):
     dname = dirname(filename)
     if not exists(dname) and dname != '':
         raise FileExistsError('internal error')
@@ -61,7 +70,7 @@ def fix_filename(filename, typ='dat'):
         return filename
 
 
-def bin_shape(cpara):
+def _bin_shape(cpara):
     if cpara.bin_radius is None:
         depmod = DepModel(cpara.stack_range)
         fzone = km2deg(np.sqrt(0.5*cpara.domperiod*depmod.vs*cpara.stack_range))
@@ -98,11 +107,14 @@ def init_profile(lat1, lon1, lat2, lon2, val):
 
 
 class CCPProfile():
-    def __init__(self, cfg_file=None, log=None):
-        if log is None:
-            self.logger = SetupLog()
-        else:
-            self.logger = log
+    def __init__(self, cfg_file=None, log:SetupLog=SetupLog()):
+        """CCPProfile class for CCP stacking along a profile
+        :param cfg_file: Configure file for CCP stacking
+        :type cfg_file: str
+        :param log: Log class for logging
+        :type log: :class:`seispy.setuplog.SetupLog`
+        """
+        self.logger = log
         if cfg_file is None:
             self.cpara = CCPPara()
         elif isinstance(cfg_file, str):
@@ -149,7 +161,7 @@ def initial_profile(self):
             self.bin_loca = np.vstack((lat, lon)).T
         else:
             self.bin_loca, self.profile_range = init_profile(*self.cpara.line, self.cpara.slide_val)
-        self.fzone = bin_shape(self.cpara)
+        self.fzone = _bin_shape(self.cpara)
         self._get_sta()
         self._select_sta()
 
@@ -199,7 +211,7 @@ def _write_sta(self):
                 self.stalst[idx,0], self.stalst[idx, 1]))
 
     def _proj_sta(self, width):
-        az1_begin, az2_begin, az1_end, az2_end, daz = line_proj(*self.cpara.line)
+        az1_begin, az2_begin, az1_end, az2_end, daz = _line_proj(*self.cpara.line)
         az_sta_begin = distaz(self.stalst[:, 0], self.stalst[:, 1], self.cpara.line[0], self.cpara.line[1]).az
         az_sta_end = distaz(self.stalst[:, 0], self.stalst[:, 1], self.cpara.line[2], self.cpara.line[3]).az
         if 0 <= daz.baz < 90 or 270 <= daz.baz < 360:
@@ -275,7 +287,7 @@ def save_stack_data(self, format='npz'):
         :param format: Format for stacked data
         :type format: str
         """
-        self.cpara.stackfile = fix_filename(self.cpara.stackfile, format)
+        self.cpara.stackfile = _fix_filename(self.cpara.stackfile, format)
         self.logger.CCPlog.info('Saving stacked data to {}'.format(self.cpara.stackfile))
         if not isinstance(self.cpara.stackfile, str):
             self.logger.CCPlog.error('fname should be in \'str\'')

seispy/decon.py

@@ -274,6 +274,20 @@ def deconwater(uin, win, dt, tshift=10., wlevel=0.05, f0=2.0, normalize=False, p
 
 
 def deconvolute(uin, win, dt, method='iter', **kwargs):
+    """ Deconvolute receiver function from waveforms.
+    :param uin: R or Q component for the response function
+    :type uin: np.ndarray
+    :param win: Z or L component for the source function
+    :type win: np.ndarray
+    :param dt: sample interval in second
+    :type dt: float
+    :param method: Method for deconvolution, defaults to 'iter'
+    :type method: str, optional
+    :param kwargs: Parameters for deconvolution
+    :type kwargs: dict
+    :return: RF, rms, [iter]
+    :rtype: np.ndarray, float, int
+    """
     if method.lower() == 'iter':
         return deconit(uin, win, dt, **kwargs)
     elif method.lower() == 'water':

seispy/distaz.py

@@ -72,7 +72,7 @@ class distaz:
     and I vaguely remember a perl port. Long live distaz! 
 
     Mijian Xu, Jan 01, 2016
-    Add np.ndarray to availible input varible.
+    Add np.ndarray to available input.
     """
 
     def __init__(self, lat1, lon1, lat2, lon2):