地震波形データをPythonで読み込んでみよう

# 本稿の方針
本稿では、地震波形データをPythonで読み込む方法についてご紹介します。

Pythonを使って地震波形データを処理したい方に参考になればと思います。特に、地震の勉強・研究・業務をしてきた方で、Pythonを今まで使ってこなかった方や、反対に、Pythonは使いこなせる方で地震波形データに新たに手を出すような方々を想定しています。ちょっとしたフィルタリングや描画ができるようになったり、データ解析（機械学習を含む）のための基礎となったりしてお役に立てればと思います。

ここでは、自分で試した範囲のことを書いておくつもりです。網羅的な説明はしません。また、Pythonの説明などはしません。

（2022年1月25日追記）限定公開で同僚に内容を確認してもらいましたが、それほど問題なさそうなので、全体公開に切り替えました。

## 環境
Python 3系を使うことにします。今回はPython 3.7.2とObsPy 1.2.1を使っています。

---

# ObsPyとは
[![alt](https://docs.obspy.org/_static/obspy-logo.png)](https://docs.obspy.org/)
(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 をご参照ください。

---

# 地震波形データの読み込み
```python
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`クラス）に格納されています。

```python
print(st[0].stats)
```
などとすると、それらのメタデータが表示されます。`st.stats`ではないので、ご注意ください。

```console
         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'})
```

ちなみに、単に

```python
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のマニュアル](https://docs.obspy.org/packages/autogen/obspy.core.trace.Stats.html#obspy.core.trace.Stats)に載っていますが、以下の通りです。

|       |意味    |
|:------|:------|
|network|観測網の名前|
|station|観測点名|
|channel|チャネル名|
|starttime|データの開始時刻|
|endtime|データの終了時刻|
|sampling_rate|サンプリングレート（1秒間当たりのサンプル数）|
|delta|サンプリング間隔[秒]|
|npts|サンプル数|

---


## 波形の描画
読み込んだ波形を描画するのは、単に

```python
st[0].plot()
```
とするだけでOKです。`st.plot()`とすると、`st`に入っているすべての地震波形データが描画されます。観測点名やチャネル名なども表示してくれて、気が利いています。
![plot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/242916/d883c1ba-1e5a-ef74-0139-e3cbeec0a7a2.png)

波形データは、`st[0].data`に1次元の`numpy.ndarray`の形で入っていますので、matplotlibを使って描画することもできます。それについては、また改めてご紹介することにします。