obspy

ObsPy - 地震学データ処理

クイックリファレンス

from obspy import read, UTCDateTime
from obspy.clients.fdsn import Client

# Read local file (MiniSEED, SAC, etc.)
st = read("data.mseed")
tr = st[0]                          # First trace
print(tr.stats)                     # Metadata

# Fetch from FDSN
client = Client("IRIS")
t = UTCDateTime("2023-02-06T01:17:00")
st = client.get_waveforms("IU", "ANMO", "00", "LHZ", t, t + 3600)
st.plot()

主なクラス

クラス	用途
Stream	複数の Trace オブジェクトを格納するコンテナ
Trace	データとメタデータを含む単一の波形
UTCDateTime	精密な時刻処理
Inventory	観測点／チャンネルのメタデータ
Catalog	地震イベント情報

基本操作

読み込みと確認

st = read("data.mseed")              # Auto-detect format
tr = st[0]
print(tr.stats.station, tr.stats.channel, tr.stats.sampling_rate)

フィルタリングと処理

st.detrend("demean")                 # Remove mean
st.detrend("linear")                 # Remove trend
st.taper(max_percentage=0.05)        # Taper edges
st.filter("bandpass", freqmin=0.1, freqmax=10.0)

FDSN から波形を取得

client = Client("IRIS")
t1 = UTCDateTime("2023-02-06T01:17:00")
st = client.get_waveforms("IU", "ANMO", "00", "LHZ", t1, t1 + 3600)
st.write("output.mseed", format="MSEED")

地震を検索

client = Client("USGS")
cat = client.get_events(
    starttime=UTCDateTime("2023-01-01"),
    endtime=UTCDateTime("2023-12-31"),
    minmagnitude=7.0
)
event = cat[0]
print(f"M{event.magnitudes[0].mag} at {event.origins[0].latitude}")

観測点のメタデータを取得

inv = client.get_stations(
    network="IU", station="ANMO",
    level="response"                 # Required for response removal
)

計器応答を除去

st = client.get_waveforms("IU", "ANMO", "00", "LHZ", t, t + 3600)
inv = client.get_stations(network="IU", station="ANMO", level="response")
st.remove_response(inventory=inv, output="VEL")  # VEL, DISP, or ACC

トリミングと選択

st.trim(UTCDateTime("2023-01-01"), UTCDateTime("2023-01-01T01:00:00"))
st_z = st.select(channel="*Z")       # Vertical only
st_bh = st.select(channel="BH*")     # BH channels

マージとギャップの処理

st.print_gaps()                      # Check for gaps
st.merge(method=1, fill_value="interpolate")

データの書き込み

st.write("output.mseed", format="MSEED")
st.write("output.sac", format="SAC")

エラー処理

from obspy.clients.fdsn.header import FDSNNoDataException

try:
    st = client.get_waveforms("IU", "ANMO", "00", "LHZ", t1, t2)
except FDSNNoDataException:
    print("No data available")

よくあるヒント

1. アーティファクトを避けるため、フィルタリングの前に必ずトレンド除去とテーパー処理を行います
2. 計器補正用の観測点を取得するときは level="response" を使用します
3. 処理の前に st.print_gaps() で ギャップを確認します
4. クエリでは ワイルドカードを使用できます：station="A*"、channel="BH?"
5. UTCDateTime は ISO 文字列、タイムスタンプ、datetime オブジェクトを受け付けます

代替ツールとの使い分け

ツール	最適な用途
obspy	地震学の包括的なワークフロー、FDSN データへのアクセス、波形処理
segyio	反射法地震探査データ用の高速かつ低レベルな SEG-Y ファイル I/O
scipy.signal	地震学固有の機能を必要としない汎用的な信号処理

地震学的なコンテキストが必要な場合は、obspy を使用します：FDSN クライアント、計器応答の除去、地震カタログ、標準的な地震データ形式などに適しています。

反射法地震探査で SEG-Y ファイルのみを扱う場合は、segyio を使用します。segyio は、インライン／クロスラインにアクセスする大規模な 3D ボリュームでより高速です。

汎用的なフィルタリングやスペクトル解析だけが必要で、地震学固有のメタデータやデータアクセスが不要な場合は、scipy.signal を使用します。

よくあるワークフロー

地震の波形を取得して処理する

- [ ] Initialize FDSN client: `Client("IRIS")`
- [ ] Search events with `client.get_events()` for target magnitude/region
- [ ] Fetch waveforms with `client.get_waveforms()` for desired stations
- [ ] Get station metadata with `level="response"` for instrument correction
- [ ] Preprocess: detrend, taper, remove instrument response
- [ ] Filter to frequency band of interest
- [ ] Trim to analysis window and export or plot

参考資料

• FDSN データセンター - 利用可能なデータセンターと機能
• チャンネルコード - チャンネル命名規則の詳細
• トラブルシューティング - よくある問題と解決策

スクリプト

• scripts/fetch_earthquake.py - 地震の波形を取得
• scripts/batch_process.py - 地震データファイルを一括処理