从零开始的Obspy

Posted on | In

参考

Obspy中文文档

读取地震数据

访问元信息

from obspy import read
waves = read('/home/zhiyuanshi/Documents/Pythonworks/IU_ADK_10.mseed')
print(waves)  # 输出波形名称

输出：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
IU.ADK.10.BH1 | 2014-10-28T00:00:00.019500Z - 2014-10-28T23:59:59.994500Z | 40.0 Hz, 3456000 samples
IU.ADK.10.BH2 | 2014-10-28T00:00:00.019500Z - 2014-10-28T23:59:59.994500Z | 40.0 Hz, 3456000 samples
IU.ADK.10.BHZ | 2014-10-28T00:00:00.019500Z - 2014-10-28T23:59:59.994500Z | 40.0 Hz, 3456000 samples
IU.ADK.10.HH1 | 2014-10-28T01:29:59.388393Z - 2014-10-28T01:34:00.158393Z | 100.0 Hz, 24078 samples
IU.ADK.10.HH2 | 2014-10-28T01:30:02.058394Z - 2014-10-28T01:33:59.158394Z | 100.0 Hz, 23711 samples
IU.ADK.10.HHZ | 2014-10-28T01:29:41.818393Z - 2014-10-28T01:34:00.778393Z | 100.0 Hz, 25897 samples
IU.ADK.10.LH1 | 2014-10-28T00:00:00.069500Z - 2014-10-28T23:59:59.069500Z | 1.0 Hz, 86400 samples
IU.ADK.10.LH2 | 2014-10-28T00:00:00.069500Z - 2014-10-28T23:59:59.069500Z | 1.0 Hz, 86400 samples
IU.ADK.10.LHZ | 2014-10-28T00:00:00.069500Z - 2014-10-28T23:59:59.069500Z | 1.0 Hz, 86400 samples
IU.ADK.10.VH1 | 2014-10-28T00:00:00.069500Z - 2014-10-28T23:59:50.069500Z | 0.1 Hz, 8640 samples
IU.ADK.10.VH2 | 2014-10-28T00:00:00.069500Z - 2014-10-28T23:59:50.069500Z | 0.1 Hz, 8640 samples
IU.ADK.10.VHZ | 2014-10-28T00:00:00.069500Z - 2014-10-28T23:59:50.069500Z | 0.1 Hz, 8640 samples
IU.ADK.10.VMU | 2014-10-28T00:00:00.000000Z - 2014-10-29T00:00:00.000000Z | 0.1 Hz, 8641 samples
IU.ADK.10.VMV | 2014-10-28T00:00:00.000000Z - 2014-10-29T00:00:00.000000Z | 0.1 Hz, 8641 samples
IU.ADK.10.VMW | 2014-10-28T00:00:00.000000Z - 2014-10-29T00:00:00.000000Z | 0.1 Hz, 8641 samples
print(len(waves))  # waves中含有15个波形

1
15

输出最前面的波形名称

tra1 = waves[0]
print(tra1)  # 输出最前面的波形名称

1
IU.ADK.10.BH1 | 2014-10-28T00:00:00.019500Z - 2014-10-28T23:59:59.994500Z | 40.0 Hz, 3456000 samples

输出元信息，即地震数据的台站信息等:

print(tra1.stats)  # 输出元信息，即地震数据的台站信息等

1
2
3
4
5
6
7
8
9
10
11
12
13
      network: IU
      station: ADK
     location: 10
      channel: BH1
    starttime: 2014-10-28T00:00:00.019500Z
      endtime: 2014-10-28T23:59:59.994500Z
sampling_rate: 40.0
        delta: 0.025
         npts: 3456000
        calib: 1.0
      _format: MSEED
        mseed: AttribDict({'dataquality': 'M', 'number_of_records': 14448, 'encoding': 'STEIM2', 'byteorder': '>', 'record_length': 512, 'filesize': 22335488})

print(tra1.stats.station)

1
ADK

访问波形数据

使用data关键词对实际波形数据进行检索

print(tra1.data)

1
[-4017 -3853 -3594 ...  1167   807  1231]
print(len(tra1.data))

1
3456000   # 与上面的波形中npts相对应
tra1.plot()  # 浏览数据

波形绘制

读取数据

from obspy.core import read

# 读取一个数据
signalchannel = read('/home/zhiyuanshi/Documents/Pythonworks/IU_ADK_10_B03.sac')
# 读取三个数据
threechannels = read('/home/zhiyuanshi/Documents/Pythonworks/IU_ADK_10_B01.sac')
threechannels += read('/home/zhiyuanshi/Documents/Pythonworks/IU_ADK_10_B02.sac')
threechannels += read('/home/zhiyuanshi/Documents/Pythonworks/IU_ADK_10_B03.sac')
print(signalchannel)

1
2
1 Trace(s) in Stream:
IU.ADK.10.BHZ | 2014-10-28T00:00:00.019500Z - 2014-10-28T23:59:59.994500Z | 40.0 Hz, 3456000 samples
print("\n")
print(threechannels)

1
2
3
4
3 Trace(s) in Stream:
IU.ADK.10.BH1 | 2014-10-28T00:00:00.019500Z - 2014-10-28T23:59:59.994500Z | 40.0 Hz, 3456000 samples
IU.ADK.10.BH2 | 2014-10-28T00:00:00.019500Z - 2014-10-28T23:59:59.994500Z | 40.0 Hz, 3456000 samples
IU.ADK.10.BHZ | 2014-10-28T00:00:00.019500Z - 2014-10-28T23:59:59.994500Z | 40.0 Hz, 3456000 samples

基本绘制与保存

signalchannel.plot() # 基本绘制，默认大小800*250，可以通过size属性调整.

size属性使用方法：size=(1920,1080)

signalsstarttime = signalchannel[0].stats.starttime
signalchannel.plot(outfile='signalchannelred.png', color='red', number_of_ticks=7, tick_rotation=5, tick_format='%I:%M %p',starttime = signalsstarttime+60*60, endtime=signalsstarttime+60*60+120)
signalchannel.plot(outfile='signalchannel.png')  # 保存波形文件

signalchannelred
signalchannel

绘制多通道波形

threechannels.plot(outfile='threechannels1920_1080.png', size=(1920,1080))

threechannels1920_1080

绘制一天的trace图像，one-day图像

signalchannel.plot(outfile='signalchannel_oneday.png', size = (3400,1440), type='dayplot')

signalchannel_oneday

signalchannel.filter('lowpass', freq=0.1, corners=2)
signalchannel.plot(outfile='signalchannel_oneday_filter',
                   type='dayplot', interval=60, right_vertical_labels=False,
                   vertical_scaling_range=5e3, one_tick_per_line=True,
                   color=['k', 'r', 'b', 'g'],  show_y_UTC_label=False,
                   events={'min_magnitude': 6.5}, size=(3400, 1440))

signalchannel_oneday_filter

绘制一段记录 (未完成，需要定义单位或者地理位置)

threechannels.plot(outfile='section_threechannels', type='section')

使用matplotlib自定义绘图

  • 疑问：signalchannel 需要取read到的内容的[0]项，具体未知为何如此。
1
2
3
4
5
6
7
import matplotlib.pyplot as plt
fig = plt.figure()  # 创建一个画布
ax = fig.add_subplot(1, 1, 1) # 创建组合图，如subplot()
ax.plot(signalchannel[0].times("matplotlib"), signalchannel[0].data, 'b-')
ax.xaxis_date()
fig.autofmt_xdate()
plt.show()

matplotlib绘图

滤波

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
import numpy as np
import matplotlib.pyplot as plt
import obspy
st = obspy.read("https://examples.obspy.org/RJOB_061005_072159.ehz.new")
st.write("st.sac", format="sac")
tr = st[0]

print(st)
# .RJOB..Z | 2005-10-06T07:21:59.850000Z - 2005-10-06T07:24:59.845000Z | 200.0 Hz, 36000 samples
print(tr)
# .RJOB..Z | 2005-10-06T07:21:59.850000Z - 2005-10-06T07:24:59.845000Z | 200.0 Hz, 36000 samples

tr_filt = tr.copy() # 复制tr内容，保留原始数据
tr_filt.filter('lowpass', freq=1.0, corners=2, zerophase=True)
# freq指低通滤波的拐点频率，zerophase指零相位

# 横坐标t，即time，
t = np.arange(0, tr.stats.npts / tr.stats.sampling_rate, tr.stats.delta)
# arange(初值，终止值，步长，数据类型)，此处，初始时间为0，终止值为数据点数/(采样率)?存疑，步长为delta，即数据点采样周期

# matplot绘图
plt.subplot(211)
plt.plot(t, tr.data, 'k')
plt.ylabel('Raw Data.')
plt.subplot(212)
plt.plot(t, tr_filt.data, 'k')
plt.ylabel('Lowpassed Data')
plt.xlabel('TIME[S]')
plt.suptitle(tr.stats.starttime)
plt.show()

filter
采用低通滤波，波形变得稀疏了。滤除了高频成分，其中高频成分大多数可能是噪音。

下采样

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
tr_new = tr.copy()
tr_new.decimate(factor=4, strict_length=False)
# 下采样时自动滤波，滤波角频率为 0.4*new sampling rate
tr_filt = tr.copy()
tr_filt.filter('lowpass', freq=0.4 * tr.stats.sampling_rate / 4.0)

t = np.arange(0, tr.stats.npts / tr.stats.sampling_rate, tr.stats.delta)
t_new = np.arange(0, tr_new.stats.npts / tr_new.stats.sampling_rate, tr_new.stats.delta)

plt.plot(t, tr.data, 'k', label='Raw', alpha=0.3)  # alpha指点的不透明度
plt.plot(t, tr_filt.data, 'b', label='Lowpassed', alpha=0.7)
plt.plot(t_new, tr_new.data, 'r', label='Lowpassed/Downsampled', alpha=0.7)
plt.xlabel('Time[s]')
plt.xlim(82, 83.5)
plt.suptitle(tr.stats.starttime)
plt.legend()
plt.show()

downsampling
未了解下采样为何，从图像上感觉相位变化了？

合并

合并绘制三个重叠的地震图。即这三个数据时间是重合的，而且是同一组数据。将其合并后时间最长的是最终结果。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
st = obspy.read("https://examples.obspy.org/dis.G.SCZ.__.BHE")
st += obspy.read("https://examples.obspy.org/dis.G.SCZ.__.BHE.1")
st += obspy.read("https://examples.obspy.org/dis.G.SCZ.__.BHE.2")
print(st[0])
print(st[1])
print(st[2])
# 输出：
# G.SCZ..BHE | 2004-01-03T08:16:09.070990Z - 2004-01-03T09:22:36.120990Z | 20.0 Hz, 79742 samples
# G.SCZ..BHE | 2004-01-03T08:19:52.371008Z - 2004-01-03T09:26:18.721008Z | 20.0 Hz, 79728 samples
# G.SCZ..BHE | 2004-01-03T08:27:25.571021Z - 2004-01-03T09:37:27.521021Z | 20.0 Hz, 84040 samples

st.sort(['starttime'])  # 排序
dt = st[0].stats.starttime.timestamp

ax = plt.subplot(4, 1, 1)
for i in range(3):
    plt.subplot(4, 1, i+1, sharex=ax)
    t = np.linspace(st[i].stats.starttime.timestamp - dt, st[i].stats.endtime.timestamp - dt, st[i].stats.npts)
    plt.plot(t, st[i].data)

# 合并数据 merge the data
st.merge(method=1)
plt.subplot(4, 1, 4, sharex=ax)
t = np.linspace(st[0].stats.starttime.timestamp - dt, st[0].stats.endtime.timestamp - dt, st[0].stats.npts)
plt.plot(t, st[0].data, 'r')
plt.show()

merge
第二个图的波形有些奇怪。震级饱和现象有些相似，但是它的振幅最大是1，其他的是20000，总之有些奇怪。
其它图像看起来很正常，最后一个图像时间轴填满了，看得出来是1、3图片的结合。合并就是将不同记录时间的同一地震波形合并为一个总的时间的地震波形，简而言之是时间上的合并。

FK分析

直角坐标

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
import matplotlib.pyplot as plt
import matplotlib.dates as mdates

import obspy
from obspy.core.util import AttribDict
from obspy.imaging.cm import obspy_sequential
from obspy.signal.invsim import corn_freq_2_paz
from obspy.signal.array_analysis import array_processing

# 导入数据
st = obspy.read('https://examples.obspy.org/agfa.mseed')
print(st)
# 设置五个通道的PAZ和坐标
st[0].stats.paz = AttribDict({
'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)],
'zeros': [0j, 0j],
'sensitivity': 205479446.68601453,
'gain': 1.0})
st[0].stats.coordinates = AttribDict({
'latitude': 48.108589,
'elevation': 0.450000,
'longitude': 11.582967})

st[1].stats.paz = AttribDict({
'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)],
'zeros': [0j, 0j],
'sensitivity': 205479446.68601453,
'gain': 1.0})
st[1].stats.coordinates = AttribDict({
'latitude': 48.108192,
'elevation': 0.450000,
'longitude': 11.583120})

st[2].stats.paz = AttribDict({
'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)],
'zeros': [0j, 0j],
'sensitivity': 250000000.0,
'gain': 1.0})
st[2].stats.coordinates = AttribDict({
'latitude': 48.108692,
'elevation': 0.450000,
'longitude': 11.583414})

st[3].stats.paz = AttribDict({
'poles': [(-4.39823 + 4.48709j), (-4.39823 - 4.48709j)],
'zeros': [0j, 0j],
'sensitivity': 222222228.10910088,
'gain': 1.0})
st[3].stats.coordinates = AttribDict({
'latitude': 48.108456,
'elevation': 0.450000,
'longitude': 11.583049})

st[4].stats.paz = AttribDict({
'poles': [(-4.39823 + 4.48709j), (-4.39823 - 4.48709j), (-
2.105 + 0j)],
'zeros': [0j, 0j, 0j],
'sensitivity': 222222228.10910088,
'gain': 1.0})
st[4].stats.coordinates = AttribDict({
'latitude': 48.108730,
'elevation': 0.450000,
'longitude': 11.583157})

# 仪器校正至拐点评率1HZ
paz1hz = corn_freq_2_paz(1.0, damp=0.707)
st.simulate(paz_remove='self', paz_simulate=paz1hz)
print(st)

# 矩阵处理
stime = obspy.UTCDateTime("20080217110515")
etime = obspy.UTCDateTime("20080217110545")
kwargs = dict(
    # slowness网格X min, X max, Y min, Y max, Slow Step
    sll_x=-3.0, slm_x=3.0, sll_y=-3.0, slm_y=3.0, sl_s=0.03,
    # sliding window properties 滑动窗口属性
    win_len=1.0, win_frac=0.05,
    # frequency properties 频率特性
    frqlow=1.0, frqhigh=8.0, prewhiten=0,
    # restrict output 限制输出
    semb_thres=-1e9, vel_thres=-1e9, timestamp='mlabday',
    stime=stime, etime=etime
)
out = array_processing(st, **kwargs)   # 存储输出

# 绘图
labels = ['rel.power', 'abs.power', 'baz', 'slow']

xlocator = mdates.AutoDateLocator()
fig = plt.figure()
for i, lab in enumerate(labels):  # enumerate 将labels以(0,'rel.power'),(1,'abs.power)...形式返回
    ax = fig.add_subplot(4, 1, i + 1)
    ax.scatter(out[:, 0], out[:, i + 1], c=out[:, 1], alpha=0.6,
    edgecolors='none', cmap=obspy_sequential)
    ax.set_ylabel(lab)
    ax.set_xlim(out[0, 0], out[-1, 0])
    ax.set_ylim(out[:, i + 1].min(), out[:, i + 1].max())
    ax.xaxis.set_major_locator(xlocator)

ax.xaxis.set_major_formatter(mdates.AutoDateFormatter(xlocator))

fig.suptitle('AGFA skyscraper blasting in Munich %s' % (stime.strftime('%Y-%m-%d'), ))
fig.autofmt_xdate()
fig.subplots_adjust(left=0.15, top=0.95, right=0.95, bottom=0.2, hspace=0)
plt.show()

FK

极坐标

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colorbar import ColorbarBase
from matplotlib.colors import Normalize
import obspy
from obspy.core.util import AttribDict
from obspy.imaging.cm import obspy_sequential
from obspy.signal.invsim import corn_freq_2_paz
from obspy.signal.array_analysis import array_processing

# Load data
st = obspy.read("https://examples.obspy.org/agfa.mseed")

# Set PAZ and coordinates for all 5 channels
st[0].stats.paz = AttribDict({
'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)],
'zeros': [0j, 0j],
'sensitivity': 205479446.68601453,
'gain': 1.0})
st[0].stats.coordinates = AttribDict({
'latitude': 48.108589,
'elevation': 0.450000,
'longitude': 11.582967})

st[1].stats.paz = AttribDict({
'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)],
'zeros': [0j, 0j],
'sensitivity': 205479446.68601453,
'gain': 1.0})
st[1].stats.coordinates = AttribDict({
'latitude': 48.108192,
'elevation': 0.450000,
'longitude': 11.583120})

st[2].stats.paz = AttribDict({
'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)],
'zeros': [0j, 0j],
'sensitivity': 250000000.0,
'gain': 1.0})
st[2].stats.coordinates = AttribDict({
'latitude': 48.108692,
'elevation': 0.450000,
'longitude': 11.583414})

st[3].stats.paz = AttribDict({
'poles': [(-4.39823 + 4.48709j), (-4.39823 - 4.48709j)],
'zeros': [0j, 0j],
'sensitivity': 222222228.10910088,
'gain': 1.0})
st[3].stats.coordinates = AttribDict({
'latitude': 48.108456,
'elevation': 0.450000,
'longitude': 11.583049})

st[4].stats.paz = AttribDict({
'poles': [(-4.39823 + 4.48709j), (-4.39823 - 4.48709j), (-2.105 + 0j)],
'zeros': [0j, 0j, 0j],
'sensitivity': 222222228.10910088,
'gain': 1.0})
st[4].stats.coordinates = AttribDict({
'latitude': 48.108730,
'elevation': 0.450000,
'longitude': 11.583157})

# Instrument correction to 1Hz corner frequency
paz1hz = corn_freq_2_paz(1.0, damp=0.707)
st.simulate(paz_remove='self', paz_simulate=paz1hz)

# Execute array_processing
kwargs = dict(
# slowness grid: X min, X max, Y min, Y max, Slow Step
sll_x=-3.0, slm_x=3.0, sll_y=-3.0, slm_y=3.0, sl_s=0.03,
# sliding window properties
win_len=1.0, win_frac=0.05,
# frequency properties
frqlow=1.0, frqhigh=8.0, prewhiten=0,
# restrict output
semb_thres=-1e9, vel_thres=-1e9,
stime=obspy.UTCDateTime("20080217110515"), etime=obspy.UTCDateTime("20080217110545")
)
out = array_processing(st, **kwargs)
# Plot
cmap = obspy_sequential
# make output human readable, adjust backazimuth to values between 0 and 360
t, rel_power, abs_power, baz, slow = out.T
baz[baz < 0.0] += 360
# choose number of fractions in plot (desirably 360 degree/N is an integer!)
N = 36
N2 = 30
abins = np.arange(N + 1) * 360. / N
sbins = np.linspace(0, 3, N2 + 1)
# sum rel power in bins given by abins and sbins
hist, baz_edges, sl_edges = \
    np.histogram2d(baz, slow, bins=[abins, sbins], weights=rel_power)

# transform to radian 转换到极坐标
baz_edges = np.radians(baz_edges)

# add polar and colorbar axes 加入极坐标轴和色标
fig = plt.figure(figsize=(8, 8))
cax = fig.add_axes([0.85, 0.2, 0.05, 0.5])
ax = fig.add_axes([0.10, 0.1, 0.70, 0.7], polar=True)
ax.set_theta_direction(-1)
ax.set_theta_zero_location("N")

dh = abs(sl_edges[1] - sl_edges[0])
dw = abs(baz_edges[1] - baz_edges[0])
# circle through backazimuth
for i, row in enumerate(hist):
    bars = ax.bar(x=(i * dw) * np.ones(N2),
        height=dh * np.ones(N2),
        width=dw, bottom=dh * np.arange(N2),
        color=cmap(row / hist.max()))
ax.set_xticks(np.linspace(0, 2 * np.pi, 4, endpoint=False))
ax.set_xticklabels(['N', 'E', 'S', 'W'])
# set slowness limits
ax.set_ylim(0, 3)
[i.set_color('grey') for i in ax.get_yticklabels()]
ColorbarBase(cax, cmap=cmap, norm=Normalize(vmin=hist.min(), vmax=hist.max()))
plt.show()

信号包络 seismogram envelopes

滤波与包络信号绘制图像

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
import numpy as np
import matplotlib.pyplot as plt

import obspy
import obspy.signal

st = obspy.read('https://examples.obspy.org/RJOB_061005_072159.ehz.new')
data = st[0].data
npts = st[0].stats.npts
samprate = st[0].stats.sampling_rate

# 滤波
st_filt = st.copy()
st_filt.filter('bandpass', freqmin=1, freqmax=3, corners=2, zerophase=True)

# 滤波后数据的包络
data_envelope = obspy.signal.filter.envelope(st_filt[0].data)

# matplotlib绘制
t = np.arange(0, npts / samprate, 1 / samprate)
plt.plot(t, st_filt[0].data, 'k')
plt.plot(t, data_envelope, 'k:')
plt.plot(t, st[0].data, 'r:')
plt.title(st[0].stats.starttime)
plt.ylabel('Filtered Data w/ Envelope')
plt.xlabel('Time [s]')
plt.xlim(80, 90)
plt.show()

红色是原始波形的绘制图像，经过1-3HZ带通滤波后，波形包络和波形呈现如黑色曲线所示。

绘制频谱 plotting spectrograms

1
2
3
import obspy
st = obspy.read("https://examples.obspy.org/RJOB_061005_072159.ehz.new")
st.spectrogram(log=True, title='BW.RJOB ' + str(st[0].stats.starttime))

零极点和频率响应 poles and zeros Frequency response

计算LE-3D/1S地震仪的频率响应

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
import numpy as np
import matplotlib.pyplot as plt

from obspy.signal.invsim import paz_to_freq_resp

poles = [-4.440 + 4.440j, -4.440 - 4.440j, -1.083 + 0.0j]
zeros = [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j]
scale_fac = 0.4

h, f = paz_to_freq_resp(poles, zeros, scale_fac, 0.005, 16384, freq=True)

# paz_to_freq_resp(极点，零点，，采样间隔，傅里叶变换点，freq)

plt.figure()
plt.subplot(121)
plt.loglog(f, abs(h))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude')

plt.subplot(122)
phase = 2 * np.pi + np.unwrap(np.angle(h))
# unwrap() 实现相位卷绕，但注意不要对没有意义的角度列表进行，否则会出现错误。
# 参见：https://blog.csdn.net/lishuhuakai/article/details/78812540
plt.semilogx(f, phase)  # 半对数函数，坐标轴x为对数坐标。此即为绘图函数，绘制函数图像。
plt.xlabel('Frequency [Hz]')
plt.ylabel('Phase [radian]')
# ticks and tick labels at multiples of pi

plt.yticks(
    [0, np.pi / 2, np.pi, 3 * np.pi / 2, 2 * np.pi],
    ['$0$', r'$\frac{\pi}{2}$', r'$\pi$', r'$\frac{3\pi}{2}$', r'$2\pi$'])
# yticks([0], [1]) 把列表0的值用列表1的值表示。
plt.ylim(-0.2, 2 * np.pi + 0.2)
# title, centered above both subplots
plt.suptitle('Frequency Response of LE-3D/1s Seismometer')
# make more room in between subplots for the ylabel of right plot
plt.subplots_adjust(wspace=0.3)   # 为子图之间的空间保留的宽度，平均轴宽的一部分
plt.show()

绘制沙滩球 Beachball

1
2
3
4
5
6
7
8
from obspy.imaging.beachball import beachball

mt = [0.91, -0.89, -0.02, 1.78, -1.55, 0.47]
beachball(mt, size=200, linewidth=1, facecolor='b')
mt2 = [150, 87, 1]
beachball(mt2, size=200, linewidth=2, facecolor='r')
mt3 = [-2.39, 1.04, 1.35, 0.57, -2.94, -0.94]
beachball(mt3, size=200, linewidth=2, facecolor='g')

但是这个球的参数mt，目前不知道是什么含义。其中红色沙滩球为转换断层机制解。

mt可能是moment tensor？即地震矩张量。



basemap plots 地图绘制

自定义投影的basemap plot

通过Inventory.plot() Catalog.plot()内置方法执行

在使用过程中出现mpl_toolkits.basemap 不能调用的情况。在命令行中输入：

pip install pyproj -i https://pypi.tuna.tsinghua.edu.cn/simple

出现不能调用的情况，需要安装basemap库。安装包源码：https://github.com/matplotlib/basemap

1
2
3
python -m pip install basemap-data
python -m pip install basemap-data-hires
python -m pip install basemap
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
from mpl_toolkits.basemap import Basemap

import numpy as np
import matplotlib.pyplot as plt
from obspy import read_inventory, read_events
# Set up a custom basemap, example is taken from basemap users' manual
fig1, ax = plt.subplots()

# setup albers equal area conic basemap
# lat_1 is first standard parallel.
# lat_2 is second standard parallel.
# lon_0, lat_0 is central point.
m = Basemap(width=8000000, height=7000000,
resolution='c', projection='aea',
lat_1=40., lat_2=60, lon_0=35, lat_0=50, ax=ax)
m.drawcoastlines()  # 绘制海岸线
m.drawcountries()  # 绘制国界
m.fillcontinents(color='wheat', lake_color='skyblue')
# draw parallels and meridians. 绘制纬线和经线
m.drawparallels(np.arange(-80., 81., 20.))
m.drawmeridians(np.arange(-180., 181., 20.))
m.drawmapboundary(fill_color='skyblue')
ax.set_title("Albers Equal Area Projection")

# we need to attach the basemap object to the figure, so that obspy knows about
# it and reuses it
# 把basemap的属性赋给fig1。
fig1.bmap = m

# now let's plot some data on the custom basemap:
# 在原地图上添加一些台站和事件信息。但是目前画出来是三张分开的图，第一张是上面的地图，
# 后面两张是台站和事件在全球地图上的位置。
inv = read_inventory()
inv.plot(fig=fig1, show=False)
cat = read_events()
cat.plot(fig=fig1, show=False, title="", colorbar=False)

plt.show()

这些代码绘制出来的图像是三个独立的图像，并没有绘制在一起。可能是BUG？

确定地点的带有Beachball的basemap plot

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
import gzip
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

from obspy.imaging.beachball import beach
# read in topo data (on a regular lat/lon grid)
# (SRTM data from: http://srtm.csi.cgiar.org/)
with gzip.open("srtm_1240-1300E_4740-4750N.asc.gz") as fp:
    srtm = np.loadtxt(fp, skiprows=8)
# origin of data grid as stated in SRTM data file header
# create arrays with all lon/lat values from min to max and
lats = np.linspace(47.8333, 47.6666, srtm.shape[0])
lons = np.linspace(12.6666, 13.0000, srtm.shape[1])
# create Basemap instance with Mercator projection
# we want a slightly smaller region than covered by our SRTM data
m = Basemap(projection='merc', lon_0=13, lat_0=48,
resolution="h",
llcrnrlon=12.75, llcrnrlat=47.69, urcrnrlon=12.95,
urcrnrlat=47.81)

# create grids and compute map projection coordinates for lon/lat grid
x, y = m(*np.meshgrid(lons, lats))

# Make contour plot 绘制网格
cs = m.contour(x, y, srtm, 40, colors="k", lw=0.5, alpha=0.3)
m.drawcountries(color="red", linewidth=1) # 绘制国界线，用红色表示

# Draw a lon/lat grid (20 lines for an interval of one degree)
m.drawparallels(np.linspace(47, 48, 21), labels=[1, 1, 0, 0],
fmt="%.2f",
dashes=[2, 2])
m.drawmeridians(np.linspace(12, 13, 21), labels=[0, 0, 1, 1],
fmt="%.2f",
dashes=[2, 2])

# Plot station positions and names into the map
# again we have to compute the projection of our lon/lat values
lats = [47.761659, 47.7405, 47.755100, 47.737167]
lons = [12.864466, 12.8671, 12.849660, 12.795714]
names = [" RMOA", " RNON", " RTSH", " RJOB"]
x, y = m(lons, lats)
m.scatter(x, y, 200, color="r", marker="v", edgecolor="k", zorder=3)
for i in range(len(names)):
    plt.text(x[i], y[i], names[i], va="top", family="monospace", weight="bold")

# Add beachballs for two events
lats = [47.751602, 47.75577]
lons = [12.866492, 12.893850]
x, y = m(lons, lats)
# Two focal mechanisms for beachball routine, specified as [strike, dip, rake]
focmecs = [[80, 50, 80], [85, 30, 90]]
ax = plt.gca()
for i in range(len(focmecs)):
    b = beach(focmecs[i], xy=(x[i], y[i]), width=1000, linewidth=1)
    b.set_zorder(10)
    ax.add_collection(b)
plt.show()

带有beachball的全球Basemap

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from obspy.imaging.beachball import beach

m = Basemap(projection='cyl', lon_0=142.36929, lat_0=38.3215, resolution='c')

m.drawcoastlines()
m.fillcontinents()  # 填充大陆
m.drawparallels(np.arange(-90., 120., 30.))
m.drawmeridians(np.arange(0., 420., 60.))
m.drawmapboundary()

x, y = m(142.36929, 38.3215)
focmecs = [0.136, -0.591, 0.455, -0.396, 0.046, -0.615]

ax = plt.gca()
b = beach(focmecs, xy=(x, y), width=10, linewidth=1, alpha=0.85)
b.set_zorder(10)
ax.add_collection(b)
plt.show()

走时和射线路径绘制

走时

绘制给定距离和相位的使用iasp91速度模型计算出的走时

1
2
3
4
5
6
7
from obspy.taup import plot_travel_times
import matplotlib.pyplot as plt

fig, ax = plt.subplots()
print(ax)
print(fig)
ax = plot_travel_times(source_depth=10, ax=ax, fig=fig, phase_list=['P', 'PP', 'S'], npoints=200)

笛卡尔射线路径

笛卡尔射线路径指直角坐标系下的二维射线路径。

1
2
3
4
5
from obspy.taup import TauPyModel

model = TauPyModel(model='iasp91')
arrivals = model.get_ray_paths(500, 140, phase_list=['PP', 'SSS'])
arrivals.plot_rays(plot_type='cartesian', phase_list=['PP', 'SSS'], plot_all=False, legend=True)

球形射线路径

使用plot_rays()方法，绘制参数为spherical。

1
2
3
4
5
6
from obspy.taup.tau import Arrivals
from obspy.taup import TauPyModel

model = TauPyModel(model='iasp91')
arrivals = model.get_ray_paths(500, 140, phase_list=['PP', 'SSS'])
arrivals.plot_rays(plot_type='spherical', phase_list=['PP', 'SSS'], plot_all=False, legend=True)

多距离射线路径

能够绘制多个震中距台站的射线路径。

1
2
3
4
5
from obspy.taup.tau import plot_ray_paths
import matplotlib.pyplot as plt

fig, ax = plt.subplots(subplot_kw=dict(polar=True))
ax = plot_ray_paths(source_depth=1000, ax=ax, fig=fig, phase_list=['P', 'PKP'], npoints=6)
  • npoints: Number of receivers to plot, 但是在图中看不出来含义。

高级示例

绘制了更多震相和台站的射线路径。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
import numpy as np
import matplotlib.pyplot as plt
from obspy.taup import TauPyModel
PHASES = [
    # Phase, distance
    ('P', 26),
    ('PP', 60),
    ('PPP', 94),
    ('PPS', 155),
    ('p', 3),
    ('pPcP', 100),
    ('PKIKP', 170),
    ('PKJKP', 194),
    ('S', 65),
    ('SP', 85),
    ('SS', 134.5),
    ('SSS', 204),
    ('p', -10),
    ('pP', -37.5),('s', -3),
    ('sP', -49),
    ('ScS', -44),
    ('SKS', -82),
    ('SKKS', -120)]
model = TauPyModel(model='iasp91')

fig, ax = plt.subplots(subplot_kw=dict(polar=True))

# Plot all pre-determined phases
for phase, distance in PHASES:
    arrivals = model.get_ray_paths(700, distance, phase_list=[phase])
    ax = arrivals.plot_rays(plot_type='spherical',
    legend=False, label_arrivals=True,
    plot_all=True, show=False, ax=ax)

# Annotate regions
ax.text(0, 0, 'Solid\ninner\ncore',
horizontalalignment='center',
verticalalignment='center',
bbox=dict(facecolor='white', edgecolor='none',
alpha=0.7))
ocr = (model.model.radius_of_planet -
(model.model.s_mod.v_mod.iocb_depth +
model.model.s_mod.v_mod.cmb_depth) / 2)
ax.text(np.deg2rad(180), ocr, 'Fluid outer core',
horizontalalignment='center',
bbox=dict(facecolor='white', edgecolor='none',
alpha=0.7))
mr = model.model.radius_of_planet - model.model.s_mod.v_mod.cmb_depth / 2
ax.text(np.deg2rad(180), mr, 'Solid mantle', horizontalalignment='center',
bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))
plt.show()