本稿の方針
本稿では、地震波形データをPythonで読み込む方法についてご紹介します。
Pythonを使って地震波形データを処理したい方に参考になればと思います。特に、地震の勉強・研究・業務をしてきた方で、Pythonを今まで使ってこなかった方や、反対に、Pythonは使いこなせる方で地震波形データに新たに手を出すような方々を想定しています。ちょっとしたフィルタリングや描画ができるようになったり、データ解析（機械学習を含む）のための基礎となったりしてお役に立てればと思います。
ここでは、自分で試した範囲のことを書いておくつもりです。網羅的な説明はしません。また、Pythonの説明などはしません。
（2022年1月25日追記）限定公開で同僚に内容を確認してもらいましたが、それほど問題なさそうなので、全体公開に切り替えました。
環境
Python 3系を使うことにします。今回はPython 3.7.2とObsPy 1.2.1を使っています。
ObsPyとは
(https://docs.obspy.org/_static/obspy-logo.png)
ObsPyは、
マニュアル等は下記のサイトにあります：
https://docs.obspy.org/
https://github.com/obspy/obspy/wiki
ObsPyをインストールしてみよう
pipが入っていれば、pip install obspyとするだけで、インストールできます。場合によっては、pip3 install obspyと、Python 3用のpipを明示的に使う必要があるかもしれません。
ObsPyのマニュアルでは、Anacondaを使ってインストールすることが推奨されているようです。https://github.com/obspy/obspy/wiki#installation をご参照ください。
地震波形データの読み込み
import obspy
st = obspy.read("data.sac")
これだけです。データ形式も自動的に認識して、対応してくれます。ここで例にしたSACファイルの場合は、ファイル1つに1本しか地震波形データが入っていませんが、WIN形式（Windowsのことではありません！）やSEED形式など、複数の地震波形データが入っている場合も全部読み込んでくれます。
対応しているデータ形式
https://docs.obspy.org/packages/autogen/obspy.core.stream.read.html#supported-formats に列挙されています。
MSEED, SAC, SEGYといった国際的に広く使われている形式を初めとして、多数の形式に対応しています。日本でよく使われるWIN形式やK-NET形式にも対応しています。（ただ、WINファイルの読み込みは非常に遅いです。それについてはまたの機会に。）変わったところでは、音声用のWAV形式にも対応しています。
読み込んだ地震波形データの中身
StreamクラスとTraceクラス
obspy.read()の返り値はStreamというクラスです。これは、地震波形が入ったTraceクラスが要素となったリストになっています。Traceクラスには地震波形データそのもののほか、観測網・観測点名、観測点位置、チャネル名、開始・終了時刻、サンプリングレート（その逆数のサンプリング間隔も）、サンプル数といったメタデータが入っています。
SACファイルのように1本しか地震波形データが入っていない場合は、st[0]しかTraceがありませんが、複数の地震波形データを一度に読み込んだ場合は、たくさんのTrace st[0], st[1], st[2], ...が入っています。
メタデータの表示
メタデータはTraceのstats（Statsクラス）に格納されています。
print(st[0].stats)
などとすると、それらのメタデータが表示されます。st.statsではないので、ご注意ください。
network:
station: GS.KOC1
location:
channel: U
starttime: 2017-06-14T10:24:00.000000Z
endtime: 2017-06-14T10:24:59.990000Z
sampling_rate: 100.0
delta: 0.01
npts: 6000
calib: 1.0
_format: SAC
sac: AttribDict({'delta': 0.0099999998, 'depmin': -0.0005446882, 'depmax': 0.00024386689, 'scale': 1.0, 'b': 0.0, 'e': 59.989998, 'depmen': 9.4378308e-08, 'nzyear': 2017, 'nzjday': 165, 'nzhour': 10, 'nzmin': 24, 'nzsec': 0, 'nzmsec': 0, 'nvhdr': 6, 'npts': 6000, 'iftype': 1, 'iztype': 9, 'leven': 1, 'lpspol': 1, 'lovrok': 1, 'lcalda': 0, 'kstnm': 'GS.KOC1', 'kcmpnm': 'U'})
ちなみに、単に
st[0].stats
としてもメタデータを見られますが、整形されておらず、見づらい場合があります。
Stats({'sampling_rate': 100.0, 'delta': 0.01, 'starttime': UTCDateTime(2017, 6, 14, 10, 24), 'endtime': UTCDateTime(2017, 6, 14, 10, 24, 59, 990000), 'npts': 6000, 'calib': 1.0, 'network': '', 'station': 'GS.KOC1', 'location': '', 'channel': 'U', 'sac': AttribDict({'delta': 0.0099999998, 'depmin': -0.0005446882, 'depmax': 0.00024386689, 'scale': 1.0, 'b': 0.0, 'e': 59.989998, 'depmen': 9.4378308e-08, 'nzyear': 2017, 'nzjday': 165, 'nzhour': 10, 'nzmin': 24, 'nzsec': 0, 'nzmsec': 0, 'nvhdr': 6, 'npts': 6000, 'iftype': 1, 'iztype': 9, 'leven': 1, 'lpspol': 1, 'lovrok': 1, 'lcalda': 0, 'kstnm': 'GS.KOC1', 'kcmpnm': 'U'}), '_format': 'SAC'})
各メタデータには、アトリビュートの形(例えば、st[0].stats.npts）でも辞書の形(st[0].stats["npts"])でもアクセスできます。
主なメタデータの意味はObsPyのマニュアルに載っていますが、以下の通りです。
| 意味 | |
|---|---|
| network | 観測網の名前 |
| station | 観測点名 |
| channel | チャネル名 |
| starttime | データの開始時刻 |
| endtime | データの終了時刻 |
| sampling_rate | サンプリングレート（1秒間当たりのサンプル数） |
| delta | サンプリング間隔[秒] |
| npts | サンプル数 |
波形の描画
読み込んだ波形を描画するのは、単に
st[0].plot()
とするだけでOKです。st.plot()とすると、stに入っているすべての地震波形データが描画されます。観測点名やチャネル名なども表示してくれて、気が利いています。
波形データは、st[0].dataに1次元のnumpy.ndarrayの形で入っていますので、matplotlibを使って描画することもできます。それについては、また改めてご紹介することにします。