GeophyAI

Python与地球物理数据处理

0%

OBSpy天然地震I

发表于 分类于 Valine：
地震数据处理包obspy简易教程(主要针对天然地震数据的获取、下载、台站数据显示等)。

天然地震

地震数据中Channel描述

Channel Description
EHZ/EHN/EHE Short Period 100 sps
BHZ/BHN/BHE Broad Band 20 sps
LHZ/LHN/LHE Long Period 1 sps
VHZ/VHN/VHE Very Long Period 0.1 sps
BCI Broad Band Calibration Signal
ECI Short Period Cal
LOG Console Log
ACE Administrative Clock Error
LCQ 1hz Clock Quality
OCF Opaque Configuration File

天然地震数据处理包obspy教程

obspy教程，以Tohoku-Oki:2011-03-11日本大地震为例

台站数据访问方法

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
import obspy
from obspy.clients.fdsn import Client

c_event = Client("IRIS")

# Event time.
event_time = obspy.UTCDateTime("2011-03-11T05:46:23.2")

# Get the event information. The temporal and magnitude constraints make it unique
cat = c_event.get_events(starttime=event_time - 10, endtime=event_time + 10,
                         minmagnitude=9)
cat.plot()
c = Client("IRIS")
# Download station information at the response level!
inv = c.get_stations(network="II", station="TLY", location="*", channel="BH?",
                     starttime=event_time - 60, endtime=event_time + 3600,
                     level="response")
inv.plot()
# Download 3 component waveforms.
# The unit of startime and endtime is seconds.
st = c.get_waveforms(network="II", station="TLY", location="*",
                     channel="BH?", starttime=event_time - 60,
                     endtime=event_time + 3600)
# Build-in plot method
st[0].plot()
# Get the data
traces = st[0].data

以上代码中，cat、inv和st的plot方法分别绘制了地震事件、台站位置信息以及波形信息。
cat
inv
trace

利用MassDownloader下载数据

台站数据下载教程(obspy.clients.fdsn.mass_downloader)

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
import obspy
from obspy.clients.fdsn.mass_downloader import CircularDomain, \
    Restrictions, MassDownloader

origin_time = obspy.UTCDateTime(2011, 3, 11, 5, 47, 32)

# Circular domain around the epicenter. This will download all data between
# 70 and 90 degrees distance from the epicenter. This module also offers
# rectangular and global domains. More complex domains can be defined by
# inheriting from the Domain class.
domain = CircularDomain(latitude=37.52, longitude=143.04,
                        minradius=70.0, maxradius=90.0)

restrictions = Restrictions(
    # Get data from 5 minutes before the event to one hour after the
    # event. This defines the temporal bounds of the waveform data.
    starttime=origin_time - 5 * 60,
    endtime=origin_time + 3600,
    # You might not want to deal with gaps in the data. If this setting is
    # True, any trace with a gap/overlap will be discarded.
    reject_channels_with_gaps=True,
    # And you might only want waveforms that have data for at least 95 % of
    # the requested time span. Any trace that is shorter than 95 % of the
    # desired total duration will be discarded.
    minimum_length=0.95,
    # No two stations should be closer than 10 km to each other. This is
    # useful to for example filter out stations that are part of different
    # networks but at the same physical station. Settings this option to
    # zero or None will disable that filtering.
    minimum_interstation_distance_in_m=10E3,
    # Only HH or BH channels. If a station has HH channels, those will be
    # downloaded, otherwise the BH. Nothing will be downloaded if it has
    # neither. You can add more/less patterns if you like.
    channel_priorities=["HH[ZNE]", "BH[ZNE]"],
    # Location codes are arbitrary and there is no rule as to which
    # location is best. Same logic as for the previous setting.
    location_priorities=["", "00", "10"])

# No specified providers will result in all known ones being queried.
mdl = MassDownloader()
# The data will be downloaded to the ``./waveforms/`` and ``./stations/``
# folders with automatically chosen file names.
mdl.download(domain, restrictions, mseed_storage="waveforms",
             stationxml_storage="stations")

以上代码会将时间origin_time之前5分钟、之后1小时以内，区域范围为圆形的domain（以经纬度以及环形区域为限定）的数据下载到当前目录下的./waveform文件夹中，台站数据则存储在./stations中。

本地文件读取(以.mseed为例)及显示

演示三通道数据读取与显示

1
2
3
4
5
6
7
8
9
import obspy
from obspy.core import read
path = r'./waveforms/'
tc = read(path+'AG.HHAR.00.HHE__20110311T054232Z__20110311T064732Z.mseed')
tc += read(path+'AG.HHAR.00.HHN__20110311T054232Z__20110311T064732Z.mseed')
tc += read(path+'AG.HHAR.00.HHZ__20110311T054232Z__20110311T064732Z.mseed')
print(tc)
# Save plot to file
tc.plot(outfile='test.png')

输出：
1
2
3
4
5
6
7
3 Trace(s) in Stream:
AG.HHAR.00.HHE | 2011-03-11T05:42:32.000000Z - 2011-03-11T06:47:32.000000Z | 100.0 Hz, 390001 samples
AG.HHAR.00.HHN | 2011-03-11T05:42:32.000000Z - 2011-03-11T06:47:32.000000Z | 100.0 Hz, 390001 samples
AG.HHAR.00.HHZ | 2011-03-11T05:42:32.000000Z - 2011-03-11T06:47:32.000000Z | 100.0 Hz, 390001 samples
"""
当流中的Traces过多时，print(stream)只会打印部分Trace，此时可通过print(stream.__str__(extended=True))来打印流中的所有Trace。
"""

基础绘图:

1
2
# Basic plotting
tc.plot()

tc

自定义绘图:

1
2
3
4
5
# Customized Plots
dt = tc[0].stats.starttime
tc.plot(color='red', number_of_ticks=7,
       tick_rotation=5, tick_format='%I:%M %p',
       starttime=dt + 60*60, endtime=dt + 60*60 + 120)

tc

使用Matplotlib进行自定义绘图：

1
2
3
4
5
6
7
8
# Custom Plotting using Matplotlib
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(tc[0].times("matplotlib"), tc[0].data, "b-")
ax.xaxis_date()
fig.autofmt_xdate()
plt.show()

tc

按小时绘制一整天道记录(文件中必须含有整天的数据)：

1
2
tc = read(path+'BW.ALTM..EHE__20120101T000000Z__20120102T000000Z.mseed')
tc.plot(type='dayplot')

tc

plot中显示地震事件

1
2
3
4
5
6
7
from obspy import read
st = read("./BW.ALTM..EHE__20121207T000000Z__20121208T000000Z.mseed")
#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})

tc

地震台站Google地图及信息查询

Network Codes

地震网络编码查询：https://www.fdsn.org/networks/?initial=1

Global Seismograph Network - IRIS/IDA

FDSN Network Information

FDSN code II Network name Global Seismograph Network - IRIS/IDA (GSN)
Start year 1986 Operated by Scripps Institution of Oceanography (SIO)
End year - Deployment region Global
Network Website http://ida.ucsd.edu/

Stations in this Network

Station Code Station Name Latitude Longitude Data Center(s)
AAK Ala Archa, Kyrgyzstan 42.6375 74.4942 IRISDMC
ABKT Alibek, Turkmenistan 37.9304 58.1189 IRISDMC
ABPO Ambohimpanompo, Madagascar -19.018 47.229 IRISDMC
ALE Alert, NU, Canada 82.5033 -62.35 IRISDMC
ARTI Arti, Russia 56.3879 58.3849 IRISDMC
ARU Arti, Russia 56.4302 58.5625 IRISDMC
ASCN Butt Crater, Ascension Island -7.9327 -14.3601 IRISDMC
BFO Black Forest Observatory, Schiltach, Germany 48.3301 8.3296 IRISDMC
BORG Borgarfjordur, Asbjarnarstadir, Iceland 64.7474 -21.3268 IRISDMC
BORK Burabay, Kazakhstan 53.0461 70.3184 IRISDMC
BRVK Borovoye, Kazakhstan 53.0581 70.2828 IRISDMC
CMLA Cha de Macela, Sao Miguel Island, Azores 37.7637 -25.5243 IRISDMC
COCO West Island, Cocos (Keeling) Islands -12.1901 96.8349 IRISDMC
DGAR Diego Garcia, Chagos Islands, Indian Ocean -7.4121 72.4525 IRISDMC
EFI Mount Kent, East Falkland Island -51.6753 -58.0637 IRISDMC
ERM Erimo, Hokkaido Island, Japan 42.015 143.1572 IRISDMC
ESK Eskdalemuir, Scotland, UK 55.3167 -3.205 IRISDMC
FFC Flin Flon, Canada 54.725 -101.9783 IRISDMC
GAR Garm, Tajikistan 39.0052 70.3328 IRISDMC
HOPE Hope Point, South Georgia Island -54.2836 -36.4879 IRISDMC
IASL Albuquerque Seismological Laboratory, N.Mex. 34.946201 -106.456703 IRISDMC
IBFO Black Forest Observatory, Schiltach, Germany 48.331902 8.3311 IRISDMC
JTS Las Juntas de Abangares, Costa Rica 10.2908 -84.9525 IRISDMC
KAPI Kappang, Sulawesi, Indonesia -5.0142 119.7517 IRISDMC
KDAK Kodiak Island, Alaska, USA 57.7828 -152.5835 IRISDMC
KIV Kislovodsk, Russia 43.9562 42.6888 IRISDMC
KURK Kurchatov, Kazakhstan 50.7154 78.6202 IRISDMC
KWAJ Kwajalein Atoll, Pacific Ocean 8.8019 167.613 IRISDMC
KWJN Gagan, Kwajalein Atoll, Marshall Islands 9.2873 167.5369 IRISDMC
LVZ Lovozero, Russia 67.8979 34.6514 IRISDMC
MBAR Mbarara, Uganda -0.6019 30.7382 IRISDMC
MSEY Mahe, Seychelles -4.6737 55.4792 IRISDMC
MSVF Monasavu, Fiji -17.7448 178.0528 IRISDMC
NIL Nilore, Pakistan 33.6506 73.2686 IRISDMC
NNA Nana, Peru -11.9875 -76.8422 IRISDMC
NRIL Norilsk, Russia 69.5049 88.4414 IRISDMC
NVS Novosibirsk, Russia 54.8404 83.2346 IRISDMC
OBN Obninsk, Russia 55.1146 36.5674 IRISDMC
PALK Pallekele, Sri Lanka 7.2728 80.7022 IRISDMC
PFO Pinon Flat, California, USA 33.6092 -116.4553 IRISDMC
RAYN Ar Rayn, Saudi Arabia 23.5225 45.5032 IRISDMC
RPN Rapanui, Easter Island, Chile -27.1267 -109.3344 IRISDMC
SACV Santiago Island, Cape Verde 14.9702 -23.6085 IRISDMC
SHEL Horse Pasture, St. Helena Island -15.9594 -5.7455 IRISDMC
SIMI Simiganj, Tajikistan 38.6585 69.0083 IRISDMC
SUR Sutherland, South Africa -32.3797 20.8117 IRISDMC
TAU Hobart, Tasmania, Australia -42.9099 147.3204 IRISDMC
TLY Talaya, Russia 51.6807 103.6438 IRISDMC
UOSS Univ. of Sharjah, Sharjah, United Arab Emirates 24.9453 56.2042 IRISDMC
WRAB Tennant Creek, NT, Australia -19.9336 134.36 IRISDMC
XBFO Black Forest Observatory, Schiltach, Germany 48.3301 8.3296 IRISDMC
XPF Pinon Flat, California, USA 33.6092 -116.4533 IRISDMC
XPFO Pinon Flat, California, USA 33.6107 -116.4555 IRISDMC