本稿の方針

本稿では、地震波形データを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を使って描画することもできます。それについては、また改めてご紹介することにします。