In [11]:
st.plot(show=False)
Out[11]:
## EARTH INVERSION
## Utpal Kumar
## Plot record section for synthetic seismograms
## Import libraries
import obspy
import instaseis
import matplotlib.pyplot as plt
import numpy as np
from obspy.core import UTCDateTime
%matplotlib inline
plt.style.use('seaborn')
!20s_PREM_ANI_FORCES/
/usr/local/bin/bash: line 1: 20s_PREM_ANI_FORCES/: Is a directory
db = instaseis.open_db('./20s_PREM_ANI_FORCES/') #requires precomputed green functions
print(db)
ReciprocalInstaseisDB reciprocal Green's function Database (v7) generated with these parameters: components : vertical and horizontal velocity model : prem_ani attenuation : True dominant period : 20.000 s dump type : displ_only excitation type : dipole time step : 4.869 s sampling rate : 0.205 Hz number of samples : 370 seismogram length : 1796.8 s source time function : errorf source shift : 34.085 s spatial order : 4 min/max radius : 5700.0 - 6371.0 km Planet radius : 6371.0 km min/max distance : 0.0 - 180.0 deg time stepping scheme : newmark2 compiler/user : gfortran 4.9.1 by lion on dhcp-10-181-12-69.dynamic.eduroam.mwn.de directory/url : 20s_PREM_ANI_FORCES size of netCDF files : 598.2 MB generated by AxiSEM version 60945ec at 2014-10-23T15:35:34.000000Z
## Specify receivers
recvr = instaseis.Receiver(latitude=0., longitude=0., network='SN', station='XYZ')
print(recvr)
Instaseis Receiver: Longitude : 0.0 deg Latitude : 0.0 deg Network : SN Station : XYZ Location :
src = instaseis.Source(
latitude=10.,
longitude=15.,
depth_in_m = 1500.,
m_rr = 4.71e+17,
m_tt = 3.81e+15,
m_pp =-4.74e+17,
m_rt = 3.99e+16,
m_rp =-8.05e+16,
m_tp =-1.23e+17,
origin_time=UTCDateTime(2020, 1, 1, 0, 0)
)
print(src)
Instaseis Source: Origin Time : 2020-01-01T00:00:00.000000Z Longitude : 15.0 deg Latitude : 10.0 deg Depth : 1.5e+00 km km Moment Magnitude : 5.73 Scalar Moment : 4.96e+17 Nm Mrr : 4.71e+17 Nm Mtt : 3.81e+15 Nm Mpp : -4.74e+17 Nm Mrt : 3.99e+16 Nm Mrp : -8.05e+16 Nm Mtp : -1.23e+17 Nm
st = db.get_seismograms(
source = src,
receiver = recvr,
components = 'ZRT',
dt = 0.05,
)
st
3 Trace(s) in Stream: SN.XYZ..BXZ | 2020-01-01T00:00:00.000000Z - 2020-01-01T00:28:24.200000Z | 20.0 Hz, 34085 samples SN.XYZ..BXR | 2020-01-01T00:00:00.000000Z - 2020-01-01T00:28:24.200000Z | 20.0 Hz, 34085 samples SN.XYZ..BXT | 2020-01-01T00:00:00.000000Z - 2020-01-01T00:28:24.200000Z | 20.0 Hz, 34085 samples
st.plot(show=False)
fig, ax = plt.subplots(figsize=(10, 6))
for depth in np.arange(100, 600, 100):
print(depth)
src = instaseis.Source(latitude = 0.0,
longitude=0.,
depth_in_m = depth * 1e3,
m_rr = 4.71e+17,
m_tt = 3.81e+15,
m_pp =-4.74e+17,
m_rt = 3.99e+16,
m_rp =-8.05e+16,
m_tp =-1.23e+17, )
st = db.get_seismograms(
source = src,
receiver = recvr,
components = 'ZRT',
dt = 0.05,)
plt.plot(st[0].data,
label=f"Depth: {depth/1000.}")
plt.legend()
100 200 300 400 500
## Slight modified after the instaseis documentation
from obspy.taup import TauPyModel
from collections import defaultdict
m = TauPyModel(model="ak135")
# some paramters
depth_in_km = 150.
mindist = 0.
maxdist = 180.
numrec = 30
fmin = 0.008
fmax = 0.05
component = "Z"
phases = ["P", "PP", "Pdiff", "S", "SS", "PS"]
colors = ["r", "r", "r", "b", "b", "g"]
# define instaseis source
src = instaseis.Source.from_strike_dip_rake(
latitude=0., longitude=0.,
depth_in_m=depth_in_km * 1e3,
M0=1e+21, strike=32., dip=62., rake=90.)
# storage for traveltimes
distances = defaultdict(list)
ttimes = defaultdict(list)
fig, ax = plt.subplots(figsize=(10, 6))
# loop over distances
for dist in np.linspace(mindist, maxdist, numrec):
# define receiver
rec = instaseis.Receiver(latitude=0, longitude=dist)
# generate seismogram, filter and plot
tr = db.get_seismograms(source=src, receiver=rec, components=[component])[0]
tr.filter('highpass', freq=fmin)
tr.filter('lowpass', freq=fmax)
tr.normalize()
plt.plot(tr.times(), tr.data * 5 + dist, color="black")
# get traveltimes
arrivals = m.get_travel_times(distance_in_degree=dist,
source_depth_in_km=depth_in_km,
phase_list=phases)
for arr in arrivals:
distances[arr.name].append(dist)
ttimes[arr.name].append(arr.time)
# plot traveltimes
for color, phase in zip(colors, phases):
plt.scatter(ttimes[phase], distances[phase], s=40, color=color)
plt.xlim(tr.times()[0], tr.times()[-1])
plt.ylim(mindist - 3, maxdist + 3)
plt.xlabel('time / s')
plt.ylabel('epicentral distance / degree')
plt.show()
