Merge pull request #60 from xumi1993/dev

1.3.4

.github/workflows/doc_build.yml

@@ -8,6 +8,14 @@ on:
     branches: [ docs ]
   pull_request:
     branches: [ master ]
+  workflow_call:
+    inputs:
+      config-path:
+        required: true
+        type: string
+    secrets:
+      token:
+        required: true
 
 
 # A workflow run is made up of one or more jobs that can run sequentially or in parallel
@@ -49,4 +57,4 @@ jobs:
           force_orphan: true
           cname: seispy.xumijian.me
           user_name: 'github-actions[bot]'
-          user_email: 'github-actions[bot]@users.noreply.github.com'
+          user_email: 'github-actions[bot]@users.noreply.github.com'

README.md

@@ -28,16 +28,20 @@
 
 Seispy is a Python module for processing seismological data and calculating Receiver Functions. The advanced functions are available to improve the Obspy.
 
+## Acknowledgements
+
+For the use of the Seispy package, please cite as:
+
+- Xu, M. & He, J. (2023). Seispy: Python Module for Batch Calculation and Postprocessing of Receiver Functions. Seismological Research Letters, 94 (2A): 935–943. doi: https://doi.org/10.1785/0220220288
+
+For 3D time-difference correction, please also consider citing:
+
+- Xu, M., Huang, H., Huang, Z., Wang, P., Wang, L., Xu, M., ... & Yuan, X. (2018). Insight into the subducted Indian slab and origin of the Tengchong volcano in SE Tibet from receiver function analysis. Earth and Planetary Science Letters, 482, 567-579. doi: https://doi.org/10.1016/j.epsl.2017.11.048
+
+For 2D and 3D CCP stacking, please also consider citing:
+
+- Xu, M., Huang, Z., Wang, L., Xu, M., Zhang, Y., Mi, N., ... & Yuan, X. (2020). Sharp lateral Moho variations across the SE Tibetan margin and their implications for plateau growth. Journal of Geophysical Research: Solid Earth, 125(5), e2019JB018117. doi: https://doi.org/10.1029/2019JB018117
 
-## Dependencies
-  * [Python]() >= 3.8
-  * [ObsPy](http://docs.obspy.org) >= 1.2.0
-  * [NumPy](http://www.numpy.org/) >= 1.16
-  * [SciPy](http://www.scipy.org/) >= 1.2.0
-  * [matplotlib](https://matplotlib.org/) >= 3.6.0
-  * [PySide6](https://www.riverbankcomputing.com/software/pyqt/) >= 6.3.0
-  * [scikits.bootstrap](https://github.com/cgevans/scikits-bootstrap) >= 1.0.0
-
 ## Installation
 
 See [Seispy documentation](https://seispy.xumijian.me/installation.html) in detail.

seispy/catalog.py

@@ -25,10 +25,10 @@ def write_catalog(query, fname, format):
 def read_catalog_file(fname):
     col = ['date', 'evla', 'evlo', 'evdp', 'mag']
     colcata = ['year', 'month', 'day', 'hour', 'minute', 'second']
-    data_cata = pd.read_table(fname, header=None, sep=' |\t', engine='python')
+    data_cata = pd.read_table(fname, header=None, sep='\s+')
     date_cols = data_cata.loc[:, [0,1,2,4,5,6]]
     date_cols.columns = colcata
-    dates = pd.to_datetime(date_cols)
+    dates = pd.DataFrame(pd.to_datetime(date_cols))[0].apply(UTCDateTime)
     eq_lst = pd.concat([dates, data_cata.loc[:, 7:]], axis=1)
     eq_lst.columns = col
     return eq_lst

seispy/ccp3d.py

@@ -1,5 +1,5 @@
 import numpy as np
-from seispy.geo import km2deg, latlon_from, cosd, extrema, skm2srad, rad2deg
+from seispy.geo import km2deg, extrema, skm2srad, rad2deg
 from seispy import distaz
 from seispy.rfcorrect import DepModel
 from seispy.setuplog import setuplog
@@ -8,50 +8,48 @@
 from seispy.signal import smooth
 from seispy.utils import check_stack_val, read_rfdep
 from scipy.interpolate import interp1d
+import pyproj
 import warnings
 import sys
 
 
 def gen_center_bin(center_lat, center_lon, len_lat, len_lon, val):
     """
-Create spaced grid point with coordinates of the center point in the area in spherical coordinates.
+    Create spaced grid point with coordinates of the center point in the area in spherical coordinates.
 
-:param center_lat: Latitude of the center point.
-:type center_lat: float
-:param center_lon: Longitude of the center point.
-:type center_lon: float
-:param len_lat: Half length in degree along latitude axis.
-:type len_lat: float
-:param len_lon: Half length in degree along longitude axis.
-:type len_lon: float
-:param val: Interval in degree between adjacent grid point.
-:type val: float
-:return: Coordinates of Grid points.
-:rtype: 2-D ndarray of floats with shape (n, 2), where n is the number of grid points.
+    :param center_lat: Latitude of the center point.
+    :type center_lat: float
+    :param center_lon: Longitude of the center point.
+    :type center_lon: float
+    :param len_lat: Half length in degree along latitude axis.
+    :type len_lat: float
+    :param len_lon: Half length in degree along longitude axis.
+    :type len_lon: float
+    :param val: Interval in degree between adjacent grid point.
+    :type val: float
+    :return: Coordinates of Grid points.
+    :rtype: 2-D ndarray of floats with shape (n, 2), where n is the number of grid points.
     """
-    lats = np.arange(0, 2*len_lat, val)
-    lons = np.arange(0, 2*len_lon, val)
-    plat, plon = latlon_from(center_lat, center_lon, 0, 90)
-    da = distaz(plat, plon, center_lat, center_lon)
-    begx = -len_lon 
-    begy = -len_lat
-    bin_loca = []
-    bin_mat = np.zeros([lats.size, lons.size, 2])
-    bin_map = np.zeros([lats.size, lons.size]).astype(int)
+    ellps="WGS84"
+    deg2km = 111.19
+    val_m = val*deg2km*1000
+    len_lat_m = len_lat*deg2km*1000
+    len_lon_m = len_lon*deg2km*1000
+    proj = pyproj.Proj(proj='aeqd', ellps=ellps, datum=ellps, lat_0=center_lat, lon_0=center_lon)
+    dx = np.arange(-len_lon_m, len_lon_m + val_m, val_m)
+    dy = np.arange(-len_lat_m, len_lat_m + val_m, val_m)
+    bin_mat = np.zeros([dy.size, dx.size, 2])
+    bin_map = np.zeros([dy.size, dx.size]).astype(int)
     n = 0
-    for j in range(lats.size):
-        delyinc = j * val + begy
-        delt = da.delta + delyinc
-        for i in range(lons.size):
-            azim = da.az + (begx + i * val) / cosd(delyinc)
-            glat, glon = latlon_from(plat, plon, azim, delt)
-            if glon > 180:
-                glon -= 360
+    bin_loca = []
+    for j, dyy in enumerate(dy):
+        for i, dxx in enumerate(dx):
+            glon, glat = proj(dxx, dyy, inverse=True)
             bin_loca.append([glat, glon])
             bin_mat[j, i, 0] = glat
             bin_mat[j, i, 1] = glon
             bin_map[j, i] = n
-            n += 1
+            n += 1 
     return np.array(bin_loca), bin_mat, bin_map
 
 

seispy/ccppara.py

@@ -2,6 +2,7 @@
 from os.path import expanduser, join, dirname, exists
 import configparser
 from seispy.geo import km2deg
+from seispy.utils import scalar_instance
 import warnings
 
 
@@ -40,7 +41,7 @@ def bin_radius(self):
 
     @bin_radius.setter
     def bin_radius(self, value):
-        if not (isinstance(value, (int, float)) or value is None):
+        if not (scalar_instance(value) or value is None):
             raise TypeError('Error type of bin_radius')
         else:
             self._bin_radius = value
@@ -118,7 +119,8 @@ def ccppara(cfg_file):
         lon2 = cf.getfloat('line', 'profile_lon2')
         cpara.line = np.array([lat1, lon1, lat2, lon2])
     except:
-        warnings.warn('line section not found. Setup it for ccp_profile')
+        # warnings.warn('line section not found. Setup it for ccp_profile')
+        pass
     try:
         # para for center bins
         cla = cf.getfloat('spacedbins', 'center_lat')
@@ -127,7 +129,10 @@ def ccppara(cfg_file):
         hllo = cf.getfloat('spacedbins', 'half_len_lon')
         cpara.center_bin = [cla, clo, hlla, hllo, km2deg(cpara.slide_val)]
     except:
-        warnings.warn('No such section of spacedbins and line. Setup them for CCP stacking')
+        # warnings.warn('No such section of spaced bins. Setup them for ccp3d')
+        pass
+    if cpara.line.size == 0 and len(cpara.center_bin) == 0:
+        raise ValueError('Please setup line or spacedbins section')
     # para for depth section
     dep_end = cf.getfloat('depth', 'dep_end')
     cpara.dep_val = cf.getfloat('depth', 'dep_val')

seispy/ccpprofile.py

@@ -1,15 +1,40 @@
 import numpy as np
-import seispy
-from seispy.geo import km2deg, deg2km, latlon_from, geoproject, sind, rad2deg, skm2srad
+from seispy.geo import km2deg, deg2km, latlon_from, \
+                       geoproject, sind, rad2deg, skm2srad, \
+                       geo2sph, sph2geo
 from seispy.setuplog import setuplog
 from seispy.distaz import distaz
 from seispy.rfcorrect import DepModel
 from seispy.rf2depth_makedata import Station
 from seispy.ccppara import ccppara, CCPPara
 from scikits.bootstrap import ci
 from seispy.ccp3d import boot_bin_stack
-from seispy.utils import check_stack_val, read_rfdep, create_center_bin_profile
+from seispy.utils import check_stack_val, read_rfdep
+from scipy.interpolate import interp1d
 from os.path import exists, dirname, basename, join
+import sys
+
+
+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'):
+    if not isinstance(stations, Station):
+        raise TypeError('Stations should be seispy.rf2depth_makedata.Station')
+    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):
@@ -146,7 +171,9 @@ def _select_sta(self):
                     self.logger.CCPlog.warning('{} does not in RFdepth structure'.format(sta))
                 if self.cpara.shape == 'rect' and self.cpara.adaptive == False:
                     self._pierce_project(self.rfdep[idx])
-        elif self.cpara.width is None and self.cpara.shape == 'circle':
+        elif self.cpara.shape == 'circle':
+            if self.cpara.width is not None:
+                self.logger.CCPlog.warning('The \'width\' will be ignored, when circle bin was set')
             dep_mod = DepModel(self.cpara.stack_range)
             x_s = np.cumsum((dep_mod.dz / dep_mod.R) / np.sqrt((1. / (skm2srad(0.085) ** 2. * (dep_mod.R / dep_mod.vs) ** -2)) - 1))
             dis = self.fzone[-1] + rad2deg(x_s[-1]) + 0.3
@@ -157,21 +184,19 @@ def _select_sta(self):
         elif self.cpara.width is not None and self.cpara.shape == 'rect':
             self.logger.CCPlog.info('Select stations within {} km perpendicular to the profile'.format(self.cpara.width))
             self.idxs = self._proj_sta(self.cpara.width)
-            self._write_sta()  
+            for idx in self.idxs:
+                self._pierce_project(self.rfdep[idx])      
+            if self.cpara.stack_sta_list != '':
+                self._write_sta()  
         else:
-            if self.cpara.width is not None:
-                self.logger.CCPlog.error('Width of profile was set, the bin shape of {} is invalid'.format(self.cpara.shape))
-                raise ValueError('Width of profile was set, the bin shape of {} is invalid'.format(self.cpara.shape))
-            else:
-                self.logger.CCPlog.error('Width of profile was not set, the bin shape of {} is invalid'.format(self.cpara.shape))
-                raise ValueError('Width of profile was not set, the bin shape of {} is invalid'.format(self.cpara.shape))
+            self.logger.CCPlog.error('Width of profile was not set, the bin shape of {} is invalid'.format(self.cpara.shape))
+            sys.exit(1)
 
     def _write_sta(self):
         with open(self.cpara.stack_sta_list, 'w') as f:
             for idx in self.idxs:
                 f.write('{}\t{:.3f}\t{:.3f}\n'.format(self.staname[idx],
                 self.stalst[idx,0], self.stalst[idx, 1]))
-                self._pierce_project(self.rfdep[idx])      
 
     def _proj_sta(self, width):
         az1_begin, az2_begin, az1_end, az2_end, daz = line_proj(*self.cpara.line)
@@ -192,7 +217,7 @@ def _proj_sta(self, width):
         final_idx = np.intersect1d(tmp_idx, proj_idx).astype(int)
         if not final_idx.any():
             self.logger.CCPlog.error('No stations within the profile belt with width of {}'.format(width))
-            raise ValueError('Satisfied stations not found')
+            sys.exit(1)
         return final_idx
 
     def _pierce_project(self, rfsta):
@@ -218,7 +243,7 @@ def stack(self):
             bin_ci = np.zeros([self.cpara.stack_range.size, 2])
             bin_count = np.zeros(self.cpara.stack_range.size)
             self.logger.CCPlog.info('{}/{} bin from {:.2f} km at lat: {:.3f} lon: {:.3f}'.format(i + 1, self.bin_loca.shape[0], self.profile_range[i], bin_info[0], bin_info[1]))
-            if self.cpara.width is None and self.cpara.shape == 'circle':
+            if self.cpara.shape == 'circle' and not exists(self.cpara.stack_sta_list):
                 idxs = self.idxs[i]
             else:
                 idxs = self.idxs
@@ -262,13 +287,21 @@ def save_stack_data(self, format='npz'):
                 for i, bin in enumerate(self.stack_data):
                     for j, dep in enumerate(self.cpara.stack_range):
                         if dep == self.cpara.stack_range[-1]:
-                            f.write('{:.4f}\t{:.4f}\t{:.4f}\t{:.2f} nan nan nan nan\n'.format(bin['bin_lat'], bin['bin_lon'],
-                                                                                              self.profile_range[i], dep))
+                            f.write(
+                                '{:.4f}\t{:.4f}\t{:.4f}\t{:.2f} nan nan nan nan\n'.format(
+                                    bin['bin_lat'], bin['bin_lon'],
+                                    self.profile_range[i], dep
+                                )
+                            )
                         else:
-                            f.write('{:.4f}\t{:.4f}\t{:.4f}\t{:.2f}\t{:.4f}\t{:.4f}\t{:.4f}\t{:d}\n'.format(bin['bin_lat'], bin['bin_lon'],
-                                                                                    self.profile_range[i],
-                                                                                    dep, bin['mu'][j], bin['ci'][j, 0],
-                                                                                    bin['ci'][j, 1], int(bin['count'][j])))
+                            f.write(
+                                '{:.4f}\t{:.4f}\t{:.4f}\t{:.2f}\t{:.4f}\t{:.4f}\t{:.4f}\t{:d}\n'.format(
+                                    bin['bin_lat'], bin['bin_lon'],
+                                    self.profile_range[i],
+                                    dep, bin['mu'][j], bin['ci'][j, 0],
+                                    bin['ci'][j, 1], int(bin['count'][j])
+                                )
+                            )
 
 
 if __name__ == '__main__':

seispy/decon.py

@@ -1,19 +1,16 @@
 # -*- coding: utf-8 -*-
 import numpy as np
 import obspy
-from obspy.io.sac import SACTrace
 from obspy.signal.util import next_pow_2
-from math import pi
-from numpy.fft import fft, ifft, ifftshift
+from numpy.fft import fft, ifft
 # from scipy.linalg import solve_toeplitz
-import matplotlib.pyplot as plt
 
 
 def gaussFilter(dt, nft, f0):
     df = 1.0 / (nft * dt)
     nft21 = 0.5 * nft + 1
     f = df * np.arange(0, nft21)
-    w = 2 * pi * f
+    w = 2 * np.pi * f
 
     gauss = np.zeros([nft, 1])
     gauss1 = np.exp(-0.25 * (w / f0) ** 2) / dt
@@ -41,9 +38,9 @@ def correl(R, W, nfft):
 def phaseshift(x, nfft, dt, tshift):
     Xf = fft(x, nfft)
     shift_i = int(tshift / dt)
-    p = 2 * pi * np.arange(1, nfft + 1) * shift_i / nfft
+    p = 2 * np.pi * np.arange(1, nfft + 1) * shift_i / nfft
     Xf = Xf * np.vectorize(complex)(np.cos(p), -np.sin(p))
-    x = ifft(Xf, nfft) / np.cos(2 * pi * shift_i / nfft)
+    x = ifft(Xf, nfft) / np.cos(2 * np.pi * shift_i / nfft)
     x = x.real
     return x
 
@@ -164,7 +161,7 @@ def deconwater(uin, win, dt, tshift=10., wlevel=0.05, f0=2.0, normalize=False, p
     fny = 1. / (2.* dt);     # nyquist
     delf = fny / (0.5 * nft)
     freq = delf * np.arange(nfpts)
-    w = 2 * pi * freq
+    w = 2 * np.pi * freq
 
     # containers
     # rff = np.zeros(nft); # Rfn in freq domain

seispy/distaz.py

@@ -153,7 +153,7 @@ def __init__(self, lat1, lon1, lat2, lon2):
 
         dbaz_idx = np.where(dbaz < 0.0)[0]
         if len(dbaz_idx) != 0:
-            if isinstance(dbaz, (int, float)):
+            if isinstance(dbaz, (int, float, np.integer, np.floating)):
                 dbaz += 2 * math.pi
             else:
                 dbaz[dbaz_idx] += 2 * math.pi

seispy/eq.py

@@ -7,6 +7,7 @@
 from seispy.decon import RFTrace
 from seispy.geo import snr, srad2skm, rotateSeisENtoTR, \
                        rssq, extrema
+from seispy.utils import scalar_instance
 from obspy.signal.trigger import recursive_sta_lta
 import glob
 
@@ -323,7 +324,7 @@ def deconvolute(self, shift, time_after, f0=2.0, method='iter', only_r=False,
             Time delta for resampling, by default None
         """
         self.method = method
-        if isinstance(f0, (int, float)):
+        if scalar_instance(f0):
             f0 = [f0]
         for ff in f0:
             if method == 'iter':