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

.github/workflows/test.yml

@@ -14,7 +14,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, macos-latest]
-        python-version: [ '3.10', '3.11', '3.12']
+        python-version: ['3.9', '3.10', '3.11']
     runs-on: ${{ matrix.os }}
 
     steps:
@@ -38,7 +38,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v2
         with:
-          python-version: '3.12'
+          python-version: '3.11'
       - name: Install dependencies
         run: |
           python -m pip install --upgrade pip

seispy/__init__.py

@@ -1 +1,4 @@
 from seispy.distaz import distaz
+from seispy.core.depmodel import DepModel
+from seispy.rfcorrect import RFStation
+from seispy.seisfwd import SynSeis

seispy/core/depmodel.py

@@ -30,8 +30,8 @@ def _search_vel_file(mode_name):
     if exists(mode_name):
         filename = mode_name
     ## if not found, search in default vel model fold
-    elif exists(join(dirname(__file__), '../data', mode_name.lower() + '.vel')):
-        filename = join(dirname(__file__), '../data', mode_name.lower() + '.vel')
+    elif exists(join(dirname(dirname(__file__)), 'data', mode_name.lower() + '.vel')):
+        filename = join(dirname(dirname(__file__)), 'data', mode_name.lower() + '.vel')
     else:
         raise ValueError('No such file of velocity model')
 
@@ -75,6 +75,7 @@ def _layer2grid(dep_range, model):
         neo_model[_i, :] = model[_j, :]
     return neo_model[:, 1], neo_model[:, 2], neo_model[:, 3]
 
+
 def _intep_mod(model, depths_elev):
     vp = interp1d(model[:,0], model[:,1], bounds_error=False,
                        fill_value=model[0,1])(depths_elev)
@@ -84,6 +85,28 @@ def _intep_mod(model, depths_elev):
                         fill_value=model[0,3])(depths_elev)
     return vp, vs, rho
 
+def _from_layer_model(dep_range, h, vp, vs, rho=None):
+    dep = 0
+    vp_dep = np.zeros_like(dep_range).astype(float)
+    vs_dep = np.zeros_like(dep_range).astype(float)
+    rho_dep = np.zeros_like(dep_range).astype(float)
+    for i, layer in enumerate(h):
+        if (layer == 0 and i == len(h)-1) or \
+           (dep+layer < dep_range[-1] and i == len(h)-1):
+           idx = np.where(dep_range>=dep)[0]
+        else:
+            idx = np.where((dep_range >= dep) & (dep_range < dep+layer))[0]
+        if idx.size == 0:
+            raise ValueError('The thickness of layer {} less than the depth interval'.format(i+1))
+        vp_dep[idx] = vp[i]
+        vs_dep[idx] = vs[i]
+        if rho is not None:
+            rho_dep[idx] = rho[i]
+        dep += layer
+        if dep > dep_range[-1]:
+            break
+    return np.vstack((dep_range, vp_dep, vs_dep, rho_dep)).T
+
 
 class DepModel(object):
     """
@@ -123,8 +146,8 @@ def __init__(self, dep_range, velmod='iasp91', elevation=0., layer_mod=False):
 
         try:
             self.model_array = _search_vel_file(velmod)
-        except (IOError, ValueError):
-            raise IOError(" failed while loading vel model {}".format(velmod))
+        except (IOError, TypeError):
+            return
         else:
             self._elevation()
             if layer_mod:
@@ -137,6 +160,11 @@ def __init__(self, dep_range, velmod='iasp91', elevation=0., layer_mod=False):
     @classmethod
     def read_layer_model(cls, dep_range, h, vp, vs, rho=None, elevation=0):
         mod = cls(dep_range, velmod=None, layer_mod=True, elevation=elevation)
+        if rho is not None:
+            mod.isrho = True
+        mod._elevation()
+        mod.model_array = _from_layer_model(mod.depths, h, vp, vs, rho=rho)
+        mod.vp, mod.vs, mod.rho = _layer2grid(mod.depths, mod.model_array)
         return mod
 
     def _elevation(self):

seispy/io.py

@@ -17,6 +17,15 @@ def __init__(self, server='IRIS', **kwargs):
     def get_events(self, starttime=None,
                    endtime=UTCDateTime.now(), 
                    **kwargs):
+        """Get events from IRIS
+
+        :param starttime: Start time of events, defaults to None
+        :type starttime: :class:`obspy.UTCDateTime`, optional
+        :param endtime: End time of events, defaults to UTCDateTime.now()
+        :type endtime: :class:`obspy.UTCDateTime`, optional
+        :return: Events
+        :rtype: :class:`obspy.Catalog`
+        """
         if endtime > UTCDateTime.now():
             endtime = UTCDateTime.now()
         events = Catalog()

seispy/scripts.py

@@ -5,6 +5,7 @@
 from seispy.utils import read_rfdep
 from seispy.rf import RF
 from seispy.recalrf import ReRF
+import sys
 
 
 def rfharmo():
@@ -267,3 +268,163 @@ def plot_r():
     rfsta = RFStation(arg.rfpath, prime_comp=arg.c)
     rfsta.plotr(rfsta, arg.output, enf=arg.enf, key=arg.k, xlim=[-2, arg.x], outformat=arg.format)
 
+
+def get_events():
+    from obspy import UTCDateTime
+    from seispy.io import Query
+    parser = argparse.ArgumentParser(description="Get seismic events from IRIS Web-Service")
+    parser.add_argument('-b', help='Start time, e.g., 20210101, 20210101020304',
+                        metavar='datetime', default=None)
+    parser.add_argument('-c', help='Catalog type',
+                        metavar='GCMT|NEIC PDE|ISC', default=None)
+    parser.add_argument('-d', help='Radial geographic constraints with center point and distance range',
+                        metavar='<lat>/<lon>/<minradius>/<maxradius>', default=None)
+    parser.add_argument('-e', help='End time, e.g., 20210101, 20210101020304',
+                        metavar='datetime', default=None)
+    parser.add_argument('-m', help='Magnitude range, optional for max magnitude',
+                        metavar='<minmagnitude>[/<maxmagnitude>]', default=None)
+    parser.add_argument('-r', help='Box range with min and max latitude and logitude, omited when specify \'-d\' ',
+                        metavar='<lat>/<lon>/<minradius>/<maxradius>', default=None)
+    parser.add_argument('-p', help='Focal depth, optional for max depth', 
+                        metavar='<mindepth>[/<maxdepth>]', default=None)
+    arg = parser.parse_args()
+    args = {}
+    if arg.c is not None:
+        args['catalog'] = arg.c
+    if arg.p is not None:
+        try:
+            values = [float(value) for value in arg.p.split('/')]
+        except:
+            raise ValueError('Error format with focal depth')
+        if len(values) == 1:
+            args['mindepth'] = values[0]
+        elif len(values) == 2:
+            args['mindepth'] = values[0]
+            args['maxdepth'] = values[1]
+        else:
+            raise ValueError('Error format with focal depth')
+    if arg.m is not None:
+        try:
+            values = [float(value) for value in arg.m.split('/')]
+        except:
+            raise ValueError('Error format with magnitude')
+        if len(values) == 1:
+            args['minmagnitude'] = values[0]
+        elif len(values) == 2:
+            args['minmagnitude'] = values[0]
+            args['maxmagnitude'] = values[1]
+        else:
+            raise ValueError('Error format with focal depth')
+    if arg.r is not None:
+        try:
+            values = [float(value) for value in arg.r.split('/')]
+        except:
+            raise ValueError('Error format with box range')
+        if len(values) == 4:
+            args['minlongitude'] = values[0]
+            args['maxlongitude'] = values[1]
+            args['minlatitude'] = values[2]
+            args['maxlatitude'] = values[3]
+        else:
+            raise ValueError('Error format with box range')
+    elif arg.d is not None:
+        try:
+            values = [float(value) for value in arg.d.split('/')]
+        except:
+            raise ValueError('Error format with Radial geographic constraints')
+        if len(values) == 4:
+            args['latitude'] = values[0]
+            args['longitude'] = values[1]
+            args['minradius'] = values[2]
+            args['maxradius'] = values[3]
+        else:
+            raise ValueError('Error format with radial geographic constraints')
+    else:
+        pass
+    if arg.b is not None:
+        try:
+            args['starttime'] = UTCDateTime(arg.b)
+        except:
+            raise ValueError('-b: Error format with time string')
+    if arg.e is not None:
+        try:
+            args['endtime'] = UTCDateTime(arg.e)
+        except:
+            raise ValueError('-e: Error format with time string')
+    if args == {}:
+        parser.print_usage()
+        sys.exit(1)
+    query = Query()
+    query.get_events(**args)
+    for _, row in query.events.iterrows():
+        print('{} {:.2f} {:.2f} {:.2f} {:.1f} {}'.format(
+              row.date.isoformat(), row.evla, row.evlo,
+              row.evdp, row.mag, row.magtype)
+        )
+
+def get_stations():
+    from obspy import UTCDateTime
+    from seispy.io import Query
+    parser = argparse.ArgumentParser(description="Get stations from IRIS Web-Service")
+    parser.add_argument('-S', help='Server name, defaults to IRIS', metavar='server', default='IRIS')
+    parser.add_argument('-b', help='Start time, e.g., 20210101, 20210101020304', metavar='datetime', default=None)
+    parser.add_argument('-c', help='Channel, wildcard is available like *,?,[EB]...', metavar='channel', default=None)
+    parser.add_argument('-d', help='Radial geographic constraints with center point and distance range',
+                        metavar='<lat>/<lon>/<minradius>/<maxradius>', default=None)
+    parser.add_argument('-e', help='End time, e.g., 20210101, 20210101020304',
+                        metavar='datetime', default=None)
+    parser.add_argument('-n', help='Network name', metavar='network', default=None)
+    parser.add_argument('-r', help='Box range with min and max latitude and logitude, omited when specify \'-d\'',
+                        metavar='<lon1>/<lon2>/<lat1>/<lat2>', default=None)
+    parser.add_argument('-s', help='Station name', metavar='station', default=None)
+    arg = parser.parse_args()
+    args = {}
+    if arg.n is not None:
+        args['network'] = arg.n
+    if arg.s is not None:
+        args['station'] = arg.s
+    if arg.r is not None:
+        try:
+            values = [float(value) for value in arg.r.split('/')]
+        except:
+            raise ValueError('Error format with box range')
+        if len(values) == 4:
+            args['minlongitude'] = values[0]
+            args['maxlongitude'] = values[1]
+            args['minlatitude'] = values[2]
+            args['maxlatitude'] = values[3]
+        else:
+            raise ValueError('Error format with box range')
+    elif arg.d is not None:
+        try:
+            values = [float(value) for value in arg.d.split('/')]
+        except:
+            raise ValueError('Error format with Radial geographic constraints')
+        args['latitude'] = values[0]
+        args['longitude'] = values[1]
+        args['minradius'] = values[2]
+        args['maxradius'] = values[3]
+    else:
+        pass
+    if arg.b is not None:
+        try:
+            args['starttime'] = UTCDateTime(arg.b)
+        except:
+            raise ValueError('-b: Error format with time string')
+    if arg.e is not None:
+        try:
+            args['endtime'] = UTCDateTime(arg.e)
+        except:
+            raise ValueError('-e: Error format with time string')
+    if arg.c is not None:
+        args['channel'] = arg.c
+    if args == {}:
+        parser.print_usage()
+        sys.exit(1)
+    query = Query(arg.S)
+    query.get_stations(**args)
+    for net in query.stations:
+        for sta in net:
+            print('{} {} {:.6f} {:.6f} {:.4f} {} {}'.format(net.code, sta.code,
+                  sta.latitude, sta.longitude, sta.elevation, sta.start_date, sta.end_date,
+                  sta.restricted_status))

seispy/seisfwd.py

@@ -4,9 +4,12 @@
 from seispy.utils import scalar_instance, array_instance
 from obspy import Trace, Stream
 from seispy.decon import RFTrace
+from numba import njit
 
 ei = 0+1j
 
+
+@njit(fastmath=True,cache=True)
 def e_inverse(omega, rho, alpha, beta, p):
     """ E_inverse (Aki & Richards, pp. 161, Eq. (5.71))
 
@@ -23,7 +26,7 @@ def e_inverse(omega, rho, alpha, beta, p):
     p : _type_
         _description_
     """
-    e_inv = np.zeros([4,4], dtype=complex)
+    e_inv = np.zeros((4,4), dtype=np.complex128)
     eta = np.sqrt(1.0/(beta*beta) - p*p)
     xi  = np.sqrt(1.0/(alpha*alpha) - p*p)
     bp = 1.0 - 2.0*beta*beta*p*p
@@ -46,7 +49,7 @@ def e_inverse(omega, rho, alpha, beta, p):
     e_inv[3,3] = e_inv[1,3]
     return e_inv
 
-
+@njit(fastmath=True,cache=True)
 def propagator_sol(omega, rho, alpha, beta, p, z):
     """ 
     propagator (Aki & Richards, pp. 398, Eq. (3) in Box 9.1)
@@ -72,7 +75,7 @@ def propagator_sol(omega, rho, alpha, beta, p, z):
         _description_
     """
 
-    p_mat = np.zeros([4,4], dtype=complex)
+    p_mat = np.zeros((4,4), dtype=np.complex128)
     beta2 = beta*beta
     p2 = p*p
     bp = 1.0 -2.0*beta2*p2
@@ -102,15 +105,16 @@ def propagator_sol(omega, rho, alpha, beta, p, z):
 
     return p_mat
 
+@njit(fastmath=True,cache=True)
 def haskell(omega, p, nl, ipha, alpha, beta, rho, h):
     i0 = 0
     e_inv = e_inverse(omega, rho[-1], alpha[-1], beta[-1], p)
     p_mat = propagator_sol(omega, rho[i0], alpha[i0], beta[i0], p, h[i0] )
     for i in range(i0+1, nl):
         p_mat2 = propagator_sol(omega, rho[i], alpha[i], beta[i], p, h[i])
-        p_mat = np.matmul(p_mat2, p_mat)
+        p_mat = p_mat2 @ p_mat
     if nl > i0+1:
-        sl = np.matmul(e_inv, p_mat)
+        sl = e_inv @ p_mat
     else:
         sl = e_inv
     denom = sl[2,0] * sl[3,1] - sl[2,1] * sl[3,0]
@@ -122,21 +126,17 @@ def haskell(omega, p, nl, ipha, alpha, beta, rho, h):
         uz = sl[2,0] / denom
     return ur, uz
 
-
+@njit(fastmath=True,cache=True)
 def fwd_seis(rayp, dt, npts, ipha, alpha, beta, rho, h):
     nlay = h.size
-    npts_max = next_pow_2(npts)
-    ur_freq = np.zeros(npts_max, dtype=complex)
-    uz_freq = np.zeros(npts_max, dtype=complex)
-    nhalf = int(npts_max / 2 + 1)
+    ur_freq = np.zeros(npts, dtype=np.complex128)
+    uz_freq = np.zeros(npts, dtype=np.complex128)
+    nhalf = int(npts / 2 + 1)
     for i in range(1, nhalf):
-        omg = 2*np.pi * i / (npts_max * dt)
+        omg = 2*np.pi * i / (npts * dt)
         ur_freq[i], uz_freq[i] = haskell(omg, rayp, nlay, ipha, 
                                          alpha, beta, rho, h)
-    ur = ifft(ur_freq).real[::-1]/npts_max
-    uz = -ifft(uz_freq).real[::-1]/npts_max
-
-    return ur[0:npts], uz[0:npts]
+    return ur_freq, uz_freq
 
 
 class SynSeis():
@@ -175,10 +175,13 @@ def run_fwd(self):
         """
         self.rstream = Stream()
         self.zstream = Stream()
+        npts_max = next_pow_2(self.npts)
         for _, rayp in enumerate(self.rayp):
-            ur, uz = fwd_seis(rayp, self.dt, self.npts, self.ipha,
+            ur_freq, uz_freq = fwd_seis(rayp, self.dt, npts_max, self.ipha,
                             self.depmod.vp, self.depmod.vs, self.depmod.rho,
                             self.depmod.thickness)
+            ur = ifft(ur_freq).real[::-1]/npts_max
+            uz = -ifft(uz_freq).real[::-1]/npts_max
             tr = Trace(data=ur)
             tr.stats.delta = self.dt
             self.rstream.append(tr)