Obspy学习笔记

Posted on | In

1、UTC标准时间、地震数据读取、数据中心

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

from obspy.core import UTCDateTime

# 设置UTC世界标准时间
    print(UTCDateTime("2012-09-07T12:15:00"))
    
    print(UTCDateTime(2012, 9, 7, 12, 15))
    
    print(UTCDateTime(1347020100.0))
    
    print(UTCDateTime("2012-09-07T12:15:00+02:00")) # 更改时区
    
    time = UTCDateTime("2012-09-07T12:15:00")
    print(time.year) # 年份属性
    print(time.julday) # 属性
    print(time.timestamp) #
    print(time.weekday) #存疑 返回一周的第几天，周一为0开始
    
    print(time + 3600)   # 增加3600秒
    time2 = UTCDateTime(2012, 9, 7)
    print(time - time2)  # 返回秒数差

## 地震数据读取
from obspy import read

# 读取obspy实例网站http://examples.obspy.org/内的数据。
st = read("http://examples.obspy.org/RJOB_061005_072159.ehz.new")
print(st)
tr = st[0]  # st中只有一个示例

print("next\n")
print(tr)
## 访问元信息
print("next\n")
print(tr.stats)
print(tr.stats.station)  # 返回stats中的station信息
print(tr.stats.gse2.datatype) # 返回datatype信息

# 访问波形数据
print(tr.data)  # 对tr中的波形数据检索，返回值为数组
print(tr.data[0:3])
print(len(tr)) # 返回波形数据长度

# 数据浏览
st.plot()  # 快速浏览波形
"""

'''
# 波形绘制
from obspy.core import read
singlechannel = read('https://examples.obspy.org/COP.BHZ.DK.2009.050')
print(singlechannel)
threechannels = read("https://examples.obspy.org/COP.BHE.DK.2009.050")
threechannels += read("https://examples.obspy.org/COP.BHN.DK.2009.050")
threechannels += read("https://examples.obspy.org/COP.BHZ.DK.2009.050")
print(threechannels)

singlechannel.plot()  # 绘制singlechannel图像
# 自定义选项
dt = singlechannel[0].stats.starttime
singlechannel.plot(color='red', number_of_ticks=7,
     tick_rotation = 5, tick_format='%l:%M %p',
     starttime = dt + 60*60, endtime = dt+60*60+120,
                   outfile = 'singlechannel.png')
## outfile 将图像保存到系统中。

threechannels.plot(size=(800, 600))
# 自行设置图像分辨率大小，多个trace会同时显示

singlechannel.plot(type='dayplot')
# 一天的trace图像，设置dayplot参数可以打印一天的图像
'''

'''
# 绘制具有事件信息的记录
from obspy import read
st = read("http://examples.obspy.org/GR.BFO..LHZ.2012.108")
st.filter("lowpass", freq=0.1, corners=2)
st.plot(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})


# stream.plot(type='section')
# 绘制stream中的某一段记录
'''

'''
# 使用matplotlib自定义绘图
import matplotlib.pyplot as plt
from obspy import read

st = read("http://examples.obspy.org/RJOB_061005_072159.ehz.new")
tr = st[0]
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(tr.times("matplotlib"), tr.data, "b-")
ax.xaxis_date()
fig.autofmt_xdate()
plt.show()
'''

重要

数据中心检索数据

1 FDSN网络服务

import obspy.clients.fdsn

2 Arclink网络服务

3 IRIS网络服务

4 Earthworm wave服务

import obspy.clients.earthworm

5 NEIC网络服务

import obspy.clients.neic

6 Syngine服务

import obspy.clients.syngine

计算震中距

参考：https://blog.csdn.net/weixin_45693832/article/details/118416213
使用geodetics中的gps2dist_azimuth计算。
示例：

1
2
3
from obspy import geodetics
distance=geodetics.gps2dist_azimuth(lat1,lon1,lat2,lon2)
angle=geodetics.kilometers2degrees(distance*0.001)