从数据中心申请地震数据

通过 obspy.clients.fdsn 模块可以访问任何支持 FDSN Web Services 协议的网络服务器,功能强大,基本上可以满足科研需求,比如下载波形数据、地震目录和台站元数据等。

警告

由于数据中心和 Web 服务可能会有变化,所以本节的相关示例未来可能会失效,当出现这种情况的时候请自行参阅最新文档。

第一步,初始化 client 对象,总是需要执行这一步:

from obspy.clients.fdsn import Client
client = Client("IRIS")  # 这里既可以用数据中心的简写,也可以用 URL 地址
Copy to clipboard

显示当前可用的数据中心列表:

from obspy.clients.fdsn.header import URL_MAPPINGS
for key in sorted(URL_MAPPINGS.keys()):
    print("{0:<7} {1}".format(key,  URL_MAPPINGS[key]))
Copy to clipboard

Out:

BGR     http://eida.bgr.de
EMSC    http://www.seismicportal.eu
ETH     http://eida.ethz.ch
GEONET  http://service.geonet.org.nz
GFZ     http://geofon.gfz-potsdam.de
ICGC    http://ws.icgc.cat
INGV    http://webservices.ingv.it
IPGP    http://ws.ipgp.fr
IRIS    http://service.iris.edu
ISC     http://isc-mirror.iris.washington.edu
KNMI    http://rdsa.knmi.nl
KOERI   http://eida.koeri.boun.edu.tr
LMU     http://erde.geophysik.uni-muenchen.de
NCEDC   http://service.ncedc.org
NIEP    http://eida-sc3.infp.ro
NOA     http://eida.gein.noa.gr
ODC     http://www.orfeus-eu.org
ORFEUS  http://www.orfeus-eu.org
RASPISHAKE http://fdsnws.raspberryshakedata.com
RESIF   http://ws.resif.fr
SCEDC   http://service.scedc.caltech.edu
TEXNET  http://rtserve.beg.utexas.edu
UIB-NORSAR http://eida.geo.uib.no
USGS    http://earthquake.usgs.gov
USP     http://sismo.iag.usp.br
Copy to clipboard

下载波形数据

通过 get_waveforms() 方法从服务器下载波形数据时,可以添加关键字参数自定义申请数据。

例如,添加参数 filename="1.mseed" 后将会把申请的数据保存在本地目录中的 1.mseed ,返回对象的结果为 None

设置 attach_response=True 将为波形数据添加仪器响应信息。

以下示例申请了美国 IU 台网 ANMO 和 AFI 台站 LHZ 分量从 2010-02-27 06:45 (UTC) 开始的 60 分钟连续波形数据,结果作为 Stream 对象返回。如果想同时发送多个请求可以使用 get_waveforms_bulk() 方法。

from obspy import UTCDateTime
t = UTCDateTime("2010-02-27T06:45:00.000")
st = client.get_waveforms("IU", "ANMO,AFI", "00", "LHZ", t, t + 60 * 60)  # 多个参数使用逗号分隔
print(st)
#st.plot(size=(800,300))
Copy to clipboard

Out:

2 Trace(s) in Stream:
IU.AFI.00.LHZ  | 2010-02-27T06:45:00.069536Z - 2010-02-27T07:44:59.069536Z | 1.0 Hz, 3600 samples
IU.ANMO.00.LHZ | 2010-02-27T06:45:00.069538Z - 2010-02-27T07:44:59.069538Z | 1.0 Hz, 3600 samples
Copy to clipboard

当然也可以使用通配符作为参数:

st = client.get_waveforms("IU", "A*", "1?", "LHZ", t, t + 5)
print(st)
#st.plot()
Copy to clipboard

Out:

3 Trace(s) in Stream:
IU.ADK.10.LHZ  | 2010-02-27T06:45:00.069536Z - 2010-02-27T06:45:04.069536Z | 1.0 Hz, 5 samples
IU.AFI.10.LHZ  | 2010-02-27T06:45:00.069536Z - 2010-02-27T06:45:04.069536Z | 1.0 Hz, 5 samples
IU.ANMO.10.LHZ | 2010-02-27T06:45:00.069538Z - 2010-02-27T06:45:04.069538Z | 1.0 Hz, 5 samples
Copy to clipboard

申请波形数据时添加仪器响应信息:

t = UTCDateTime("2012-12-14T10:36:01.6Z")
st = client.get_waveforms("TA", "E42A", "*", "BH?", t+300, t+400,
                          attach_response=True)
st.remove_response(output="VEL")   # 去除仪器响应
print(st)
#st.plot(size=(800,400))
Copy to clipboard

Out:

3 Trace(s) in Stream:
TA.E42A..BHE | 2012-12-14T10:41:01.624998Z - 2012-12-14T10:42:41.599998Z | 40.0 Hz, 4000 samples
TA.E42A..BHN | 2012-12-14T10:41:01.624998Z - 2012-12-14T10:42:41.599998Z | 40.0 Hz, 4000 samples
TA.E42A..BHZ | 2012-12-14T10:41:01.624998Z - 2012-12-14T10:42:41.599998Z | 40.0 Hz, 4000 samples
Copy to clipboard

利用 get_waveforms_bulk() 方法可以同时提交多个申请, 符合要求的 bulk 有以下形式:

  • 多个列表项组成的列表,每一个列表项必须包含 network, station, location, channel, starttime and endtime

  • 包含有效 request 的字符串

  • 包含有效 request 的文件名

  • 包含有效 request 的已打开文件的句柄

多个列表项组成的列表:

client = Client("IRIS")
t1 = UTCDateTime("2010-02-27T06:30:00.000")
t2 = t1 + 1
t3 = t1 + 3
bulk = [("IU", "ANMO", "*", "BHZ", t1, t2),
          ("IU", "AFI", "1?", "BHE", t1, t3),
            ("GR", "GRA1", "*", "BH*", t2, t3)]
st = client.get_waveforms_bulk(bulk)
print(st)
#st.plot()
Copy to clipboard

Out:

5 Trace(s) in Stream:
GR.GRA1..BHE   | 2010-02-27T06:30:01.040000Z - 2010-02-27T06:30:02.990000Z | 20.0 Hz, 40 samples
GR.GRA1..BHN   | 2010-02-27T06:30:01.040000Z - 2010-02-27T06:30:02.990000Z | 20.0 Hz, 40 samples
GR.GRA1..BHZ   | 2010-02-27T06:30:01.039999Z - 2010-02-27T06:30:02.989999Z | 20.0 Hz, 40 samples
IU.ANMO.00.BHZ | 2010-02-27T06:30:00.019538Z - 2010-02-27T06:30:00.969538Z | 20.0 Hz, 20 samples
IU.ANMO.10.BHZ | 2010-02-27T06:30:00.019538Z - 2010-02-27T06:30:00.994538Z | 40.0 Hz, 40 samples
Copy to clipboard

包含有效 request 的字符串:

bulk = 'quality=B\n' + \
         'longestonly=false\n' + \
         'IU ANMO * BHZ 2010-02-27 2010-02-27T00:00:02\n' + \
         'IU AFI 1? BHE 2010-02-27 2010-02-27T00:00:04\n' + \
         'GR GRA1 * BH? 2010-02-27 2010-02-27T00:00:02\n'
st = client.get_waveforms_bulk(bulk)
print(st)
#st.plot()

# 包含 request 的文件:
#st = client.get_waveforms_bulk("request.txt", attach_response=True)
#st.remove_response(output="VEL")
Copy to clipboard

Out:

5 Trace(s) in Stream:
GR.GRA1..BHE   | 2010-02-27T00:00:00.040000Z - 2010-02-27T00:00:01.990000Z | 20.0 Hz, 40 samples
GR.GRA1..BHN   | 2010-02-27T00:00:00.039999Z - 2010-02-27T00:00:01.989999Z | 20.0 Hz, 40 samples
GR.GRA1..BHZ   | 2010-02-27T00:00:00.039999Z - 2010-02-27T00:00:01.989999Z | 20.0 Hz, 40 samples
IU.ANMO.00.BHZ | 2010-02-27T00:00:00.019500Z - 2010-02-27T00:00:01.969500Z | 20.0 Hz, 40 samples
IU.ANMO.10.BHZ | 2010-02-27T00:00:00.019500Z - 2010-02-27T00:00:01.994500Z | 40.0 Hz, 80 samples
Copy to clipboard

下载地震事件

通过 get_events 方法可以从服务器申请地震数据,申请结果返回 Catalog 对象。

eventid 申请地震:

client = Client("IRIS")
cat = client.get_events(eventid=609301)
print(cat)
Copy to clipboard

Out:

1 Event(s) in Catalog:
1997-10-14T09:53:11.070000Z | -22.145, -176.720 | 7.8 mw
Copy to clipboard

用起止时间和震级等条件筛选地震,然后利用 Catalogplot() 方法绘制震中分布图:

t1 = UTCDateTime("2001-01-07T00:00:00")
t2 = UTCDateTime("2001-01-07T03:00:00")
cat = client.get_events(starttime=t1, endtime=t2, minmagnitude=4,
                        catalog="ISC")
print(cat)
cat.plot(projection="ortho")
Copy to clipboard
3 events (2001-01-07 to 2001-01-07) - Color codes depth, size the magnitude

Out:

3 Event(s) in Catalog:
2001-01-07T02:55:59.290000Z |  +9.801,  +76.548 | 4.9 mb
2001-01-07T02:35:35.170000Z | -21.291,  -68.308 | 4.4 mb
2001-01-07T00:09:25.630000Z | +22.946, -107.011 | 4.0 mb
/usr/share/miniconda3/envs/seisnote/lib/python3.8/site-packages/obspy/imaging/maps.py:762: MatplotlibDeprecationWarning: The 'cmap' parameter to Colorbar has no effect because it is overridden by the mappable; it is deprecated since 3.3 and will be removed two minor releases later.
  cb = Colorbar(cm_ax, scatter, cmap=colormap,

<Figure size 640x480 with 2 Axes>
Copy to clipboard

下载台站数据

通过 get_stations 方法可以从服务器申请台站数据,申请结果返回 Inventory 对象:

inventory = client.get_stations(network="IU", station="A*",
                                 starttime=t1,
                                 endtime=t2, level="response")  # 申请所有台站的所有分量和仪器响应文件
print(inventory)
Copy to clipboard

Out:

Inventory created at 2021-02-03T15:46:58.000000Z
        Created by: IRIS WEB SERVICE: fdsnws-station | version: 1.1.47
                    http://service.iris.edu/fdsnws/station/1/query?starttime=2001-01-07...
        Sending institution: IRIS-DMC (IRIS-DMC)
        Contains:
                Networks (1):
                        IU
                Stations (3):
                        IU.ADK (Adak, Aleutian Islands, Alaska)
                        IU.AFI (Afiamalu, Samoa)
                        IU.ANMO (Albuquerque, New Mexico, USA)
                Channels (88):
                        IU.ADK..BC0, IU.ADK.00.BC0, IU.ADK.00.BHZ, IU.ADK.00.BHN,
                        IU.ADK.00.BHE, IU.ADK.00.LHZ, IU.ADK.00.LHN, IU.ADK.00.LHE,
                        IU.ADK.00.UHZ, IU.ADK.00.UHN, IU.ADK.00.UHE, IU.ADK.00.VE1,
                        IU.ADK.00.VHZ, IU.ADK.00.VHN, IU.ADK.00.VHE, IU.ADK.00.VK1,
                        IU.ADK.00.VMZ, IU.ADK.00.VMN, IU.ADK.00.VME, IU.ADK.10.HLZ,
                        IU.ADK.10.HLN, IU.ADK.10.HLE, IU.ADK.10.LLZ, IU.ADK.10.LLN,
                        IU.ADK.10.LLE, IU.AFI..BC0, IU.AFI..BC1, IU.AFI.00.BHZ,
                        IU.AFI.00.BHN, IU.AFI.00.BHE, IU.AFI.00.LHZ, IU.AFI.00.LHN,
                        IU.AFI.00.LHE, IU.AFI.00.UHZ, IU.AFI.00.UHN, IU.AFI.00.UHE,
                        IU.AFI.00.VE1, IU.AFI.00.VHZ, IU.AFI.00.VHN, IU.AFI.00.VHE,
                        IU.AFI.00.VK1, IU.AFI.00.VMZ, IU.AFI.00.VMN, IU.AFI.00.VME,
                        IU.AFI.10.BHZ, IU.AFI.10.BHN, IU.AFI.10.BHE, IU.AFI.10.HHZ,
                        IU.AFI.10.HHN, IU.AFI.10.HHE, IU.AFI.10.LHZ, IU.AFI.10.LHN,
                        IU.AFI.10.LHE, IU.AFI.20.HLZ, IU.AFI.20.HLN, IU.AFI.20.HLE,
                        IU.AFI.20.LLZ, IU.AFI.20.LLN, IU.AFI.20.LLE, IU.ANMO..BC0,
                        IU.ANMO..BC1, IU.ANMO.00.BHZ, IU.ANMO.00.BH1, IU.ANMO.00.BH2,
                        IU.ANMO.00.LHZ, IU.ANMO.00.LH1, IU.ANMO.00.LH2, IU.ANMO.00.UE1,
                        IU.ANMO.00.UHZ, IU.ANMO.00.UH1, IU.ANMO.00.UH2, IU.ANMO.00.VE1,
                        IU.ANMO.00.VHZ, IU.ANMO.00.VH1, IU.ANMO.00.VH2, IU.ANMO.00.VK1,
                        IU.ANMO.10.BHZ, IU.ANMO.10.BH1, IU.ANMO.10.BH2, IU.ANMO.10.HHZ,
                        IU.ANMO.10.HH1, IU.ANMO.10.HH2, IU.ANMO.10.LHZ, IU.ANMO.10.LH1,
                        IU.ANMO.10.LH2, IU.ANMO.10.VHZ, IU.ANMO.10.VH1, IU.ANMO.10.VH2
Copy to clipboard

利用 Inventoryplot() 方法绘制台站分布图:

inventory.plot(projection="local")  # 绘制台站分布
Copy to clipboard
plot 5 retrieving data

Out:

<Figure size 640x480 with 1 Axes>
Copy to clipboard

利用 Inventoryplot_response() 方法绘制仪器响应文件:

inventory.plot_response(min_freq=0.01, output="DISP", station="ADK", channel="BH?")
Copy to clipboard
plot 5 retrieving data

Out:

<Figure size 640x480 with 2 Axes>
Copy to clipboard

Total running time of the script: ( 0 minutes 10.383 seconds)

Gallery generated by Sphinx-Gallery