diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 41cdba6e..9d12d857 100644 --- a/.github/workflows/test.yml +++ b/.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 diff --git a/seispy/__init__.py b/seispy/__init__.py index 319196c1..416a26a6 100644 --- a/seispy/__init__.py +++ b/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 diff --git a/seispy/core/depmodel.py b/seispy/core/depmodel.py index 8a02805f..463bfda3 100644 --- a/seispy/core/depmodel.py +++ b/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): diff --git a/seispy/io.py b/seispy/io.py index 5545f31f..f4305bbd 100644 --- a/seispy/io.py +++ b/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() diff --git a/seispy/scripts.py b/seispy/scripts.py index cb8c905c..4dc6cdac 100644 --- a/seispy/scripts.py +++ b/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='///', 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='[/]', default=None) + parser.add_argument('-r', help='Box range with min and max latitude and logitude, omited when specify \'-d\' ', + metavar='///', default=None) + parser.add_argument('-p', help='Focal depth, optional for max depth', + metavar='[/]', 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='///', 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='///', 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)) \ No newline at end of file diff --git a/seispy/seisfwd.py b/seispy/seisfwd.py index a99d205c..7d2de7a4 100644 --- a/seispy/seisfwd.py +++ b/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) diff --git a/setup.py b/setup.py index 01348fc2..f2ab2368 100644 --- a/setup.py +++ b/setup.py @@ -25,6 +25,7 @@ 'obspy>=1.2.1', 'pyside6>=6.2.0', 'scikits.bootstrap>=1.0.0', + 'numba', 'pyproj'], entry_points={'console_scripts': ['gen_rayp_lib=seispy.psrayp:gen_rayp_lib', 'prf=seispy.scripts:prf', @@ -42,7 +43,9 @@ 'ccp3d=seispy.scripts:ccp3d', 'rfharmo=seispy.scripts:rfharmo', 'get_pierce_points=seispy.scripts:get_pierce_points', - 'veltxt2mod=seispy.modcreator:veltxt2mod']}, + 'veltxt2mod=seispy.modcreator:veltxt2mod', + 'get_stations=seispy.scripts:get_stations', + 'get_events=seispy.scripts:get_events',]}, zip_safe=False, classifiers=['Programming Language :: Python', 'Programming Language :: Python :: 3.9', diff --git a/test/test_case02.py b/test/test_case02.py index bf0738de..2d05bc02 100644 --- a/test/test_case02.py +++ b/test/test_case02.py @@ -16,6 +16,8 @@ def test_download(): def test_sub01(): rfs = RFStation('ex-rfani/SC.LTA') rfs.jointani(3, 8, tlen=3.5, stack_baz_val=10, weight=[0.4, 0.4, 0.2]) + rfs.ani.plot_stack_baz() + rfs.ani.plot_polar() def test_sub02(): @@ -24,6 +26,8 @@ def test_sub02(): rfs.psrf2depth() rfs.psrf_1D_raytracing() rfs.slantstack() + rfs.slant.syn_tps(['P410s', 'P660s']) + rfs.slant.plot() def test_sub03(): diff --git a/test/test_case05.py b/test/test_case05.py index aacbcbf9..a416fcef 100644 --- a/test/test_case05.py +++ b/test/test_case05.py @@ -1,6 +1,7 @@ from seispy.core.depmodel import DepModel from seispy.seisfwd import SynSeis from seispy.rfcorrect import RFStation +from seispy.geo import skm2srad import pytest import numpy as np @@ -16,5 +17,23 @@ def test_sub01(): rfsta.moveoutcorrect(ref_rayp=0.06, dep_range=np.arange(100)) +def test_sub02(): + h = [30, 0] + vp = [6.5, 8.04] + vs = [3.6, 4.48] + rho = [2.67, 3.2] + rayp = 0.06 + rrayp = skm2srad(rayp) + dt = 0.1 + depth = np.arange(100) + basemodel = DepModel.read_layer_model(depth, h, vp, vs, rho=rho) + basemodel.plot_model(show=False) + basemodel.tpds(rrayp, rrayp) + basemodel.tpppds(rrayp, rrayp) + basemodel.tpspds(rrayp) + ssfwd = SynSeis(basemodel, rayp, dt, npts=1200) + ssfwd.run_fwd() + + if __name__ == '__main__': test_sub01()