fix slant stack

seispy/rfcorrect.py

@@ -103,7 +103,6 @@ class SACStation(object):
             for _i, evt, ph in zip(range(self.ev_num), self.event, self.phase):
                 sac = SACTrace.read(join(data_path, evt + '_' + ph + '_{}.sac'.format(self.comp)))
                 self.datar[_i] = sac.data
-        raise FutureWarning('This "SACStation" class will be removed in the new version and replaced by "RFStation"')
 
     def normalize(self):
         maxamp = np.nanmax(np.abs(self.datar), axis=1)

seispy/slantstack.py

 import numpy as np
 from scipy.interpolate import interp1d
 import matplotlib.pyplot as plt
+from obspy.taup import TauPyModel
+from seispy.geo import srad2skm, skm2sdeg
 
 class SlantStack():
     def __init__(self, seis, timeaxis, dis) -> None:
@@ -19,30 +21,52 @@ class SlantStack():
         self.ref_dis = 65
         self.rayp_range = np.arange(-0.35, 0.35, 0.01)
         self.tau_range = np.arange(0, 100)
+        self.tau = []
+        self.drayp = []
 
     def stack(self, ref_dis=None, rayp_range=None, tau_range=None):
         if ref_dis is not None and isinstance(ref_dis, (int, float)):
             self.ref_dis = ref_dis
+        elif ref_dis is None:
+            pass
         else:
             raise TypeError('{} should be in int or float type.'.format(ref_dis))
         if rayp_range is not None and isinstance(rayp_range, np.ndarray):
             self.rayp_range = rayp_range
+        elif rayp_range is None:
+            pass
         else:
             raise TypeError('{} should be in numpy.ndarray type.'.format(rayp_range))
         if tau_range is not None and isinstance(tau_range, np.ndarray):
             self.tau_range = tau_range
+        elif tau_range is None:
+            pass
         else:
             raise TypeError('{} should be in numpy.ndarray type.'.format(tau_range))
         ev_num = self.datar.shape[0]
-        tmp = np.zeros([ev_num, tau_range.shape[0]])
-        self.stack_amp = np.zeros([rayp_range.shape[0], tau_range.shape[0]])
-        for j in range(rayp_range.shape[0]):
-            for i in range(ev_num):
-                self.datar[i, :] = self.datar[i, :] / np.max(np.abs(self.datar[i, :]))
-                tps = tau_range - rayp_range[j] * (self.dis[i] - ref_dis)
-                tmp[i, :] = interp1d(self.time_axis, self.datar[i, :], fill_value='extrapolate')(tps)
-            self.stack_amp[j, :] = np.mean(tmp, axis=0)
-        
+        taus, rayps = np.meshgrid(self.tau_range, self.rayp_range)
+        self.stack_amp = np.zeros([self.rayp_range.shape[0], self.tau_range.shape[0]])
+        for i in range(ev_num):
+            tps = taus - rayps * (self.dis[i] - self.ref_dis)
+            self.stack_amp += interp1d(self.time_axis, self.datar[i, :], fill_value='extrapolate')(tps)
+        self.stack_amp /= ev_num
+
+        # for j in range(self.rayp_range.shape[0]):
+        #     for i in range(ev_num):
+        #         self.datar[i, :] = self.datar[i, :] / np.max(np.abs(self.datar[i, :]))
+        #         tps = self.tau_range - self.rayp_range[j] * (self.dis[i] - self.ref_dis)
+        #         tmp[i, :] = interp1d(self.time_axis, self.datar[i, :], fill_value='extrapolate')(tps)
+        #     self.stack_amp[j, :] = np.mean(tmp, axis=0)
+
+    def syn_tps(self, phase_list, model='iasp91'):
+        model = TauPyModel(model=model)
+        phase_list.insert(0, 'P')
+        arrs = model.get_travel_times(10, self.ref_dis, phase_list=phase_list)
+        p_arr = arrs[0].time
+        p_rayp = skm2sdeg(srad2skm(arrs[0].ray_param))
+        self.tau = [arr.time - p_arr for arr in arrs[1:]]
+        self.drayp = [p_rayp - skm2sdeg(srad2skm(arr.ray_param))for arr in arrs[1:]]
+
     def plot(self, cmap='jet', xlim=None, vmin=None, vmax=None, figpath=None):
         self.fig = plt.figure()
         self.ax = self.fig.add_subplot()
@@ -57,13 +81,14 @@ class SlantStack():
             pass
         else:
             raise TypeError('vmin and vmax must be in int or in float type')
-        cs = self.ax.pcolor(self.tau_range, self.rayp_range, self.stack_amp,
+        im = self.ax.pcolor(self.tau_range, self.rayp_range, self.stack_amp,
                             vmax=vmax, vmin=vmin, cmap=cmap)
-        self.ax.colorbar(cs)
+        self.fig.colorbar(im, ax=self.ax)
+        self.ax.scatter(self.tau, self.drayp, color='k', marker='x')
         if xlim is not None and isinstance(xlim, (list, np.ndarray)):
             self.ax.set_xlim(xlim)
         self.ax.set_xlabel('Time (s)')
-        self.ax.set_ylabel('Slowness (s/km)')
+        self.ax.set_ylabel('Slowness (s/$^\circ$)')
         if figpath is not None and isinstance(figpath, str):
             self.fig.savefig(figpath, format='png', dpi=400, bbox_inches='tight')
         else:

setup.py

@@ -6,7 +6,7 @@ with open("README.md", "r") as fh:
     long_description = fh.read()
 
 
-VERSION = "1.2.8"
+VERSION = "1.2.9"
 setup(name='python-seispy',
       version=VERSION,
       author='Mijian Xu',
@@ -19,7 +19,6 @@ setup(name='python-seispy',
       package_data={'': ['data/*']},
       install_requires=[
                 'netcdf4>=1.5.2',
-                'pyerf>=1.0.1',
                 'obspy>=1.2.0',
                 'pandas>=1.0.0',
                 'numpy>=1.19.0',