ObsPyを触ってみる その2

ObsPy Documentation (1.4.0) — ObsPy 1.4.0 documentation

もう少しObsPyで遊んでみよう。

WIN32形式をObsPyで読み込む
1. SACに変換する
2. WINに変換する
連続ウェーブレット変換をしてみる
1. 地震波形を取得する
2. 連続ウェーブレット変換を行う

WIN32形式をObsPyで読み込む

WIN32形式をObsPyで読み込むには、SACやWINに変換して、読み込ませます。

support reading WIN32 format by seisman · Pull Request #1692 · obspy/obspy

This PR adds basic support for reading WIN32 format. It's a little similar to obspy.io.datamark. It works under Python 3...

SACに変換する

これは説明済みです。

SACのインストール

IRISが提供するSACのインストール方法について説明しています。

WINに変換する

WIN32形式からWIN形式に変換するには、win32toolsのw32tow1を用います。

w32tow1 [入力ファイル] [出力ファイル]

たとえば、0101_202005030000_60.cntというWIN32形式のファイルを変換するには、

w32tow1 0101_202005030000_60.cnt 0101_202005030000

とすればよい。簡単。

連続ウェーブレット変換をしてみる

Continuous Wavelet Transform — ObsPy 1.4.0 documentation

ここでは焼岳周辺で発生している、深部低周波地震の波形を連続ウェーブレット変換してみよう。

ObsPyでは、なんと標準で連続ウェーブレット変換ができてしまう。

地震波形を取得する

というわけで早速Hi-netから、地震波形を取得します。
おすすめの観測点は、N.FRKHとN.KOKH。

今回は2020年5月3日0時0分から45分間の波形を取得します。

HinetPyなら、

from HinetPy import Client, win32
client = Client("UserID", "Password")
client.
data, ctable = client.get_continuous_waveform('0101', 202005030000, 45)
win32.extract_sac("0101_202005030000_45.cnt", "0101_20200503.ch")

で波形取得とSACに変換ができます。

連続ウェーブレット変換を行う

早速変換してみよう。今回はN.KOKH.N.SACを変換してみる。

import numpy as np
import matplotlib.pyplot as plt

import obspy
from obspy.imaging.cm import pqlx
from obspy.signal.tf_misfit import cwt

st = obspy.read("./N.KOKH.N.SAC")
tr = st[0]
npts = tr.stats.npts
dt = tr.stats.delta
t = np.linspace(0, dt * npts, npts)
f_min = 1
f_max = 50

scalogram = cwt(tr.data, dt, 8, f_min, f_max)

fig = plt.figure()
ax = fig.add_subplot(111)

x, y = np.meshgrid(
    t,
    np.logspace(np.log10(f_min), np.log10(f_max), scalogram.shape[0]))

ax.pcolormesh(x, y, np.abs(scalogram), cmap=pqlx)
ax.set_xlabel("Time after %s [s]" % tr.stats.starttime)
ax.set_ylabel("Frequency [Hz]")
ax.set_yscale('log')
ax.set_ylim(f_min, f_max)
fig.savefig("hogehoge.png")

この結果が、

となります。

ObsPyを触ってみるその2