Zenn
🌏

wintools: PythonによるWINフォーマット地震波形の読み書きモジュール

2024/09/23に公開
Python
ライブラリ
地震学
tech

※あまりにも雑な記事なので，時間あるときに写真など増やします．


柔軟性・コード作成効率が高いPythonで，WINに特化したライブラリを作成した．
WINシステムをインストールしていない環境でも，Pythonユーザーであれば容易にデータを読み込み，プロット，処理できる．
応用すればwin2sacやsac2winとしての活用もできる．


例: WINのサンプルデータの読み込み結果のプロット

背景

知識

  • WINシステム:
    WINフォーマットのデータを扱うためのC, Fortranのプログラム群．
  • WINフォーマット:
    日本国内で実質的な共通規格として用いられている，地震波形データのフォーマット．

現状の問題点

  • Pythonの地震波形処理ライブラリObsPyによるWINの読み込みは不十分．
    • 年の処理(90年代のデータが「2090年」のように誤解釈される)
    • チャンネルテーブル情報を付加できない
      • 観測点名，成分名を扱えない
      • 振幅がnビット整数値のままで，物理量（速度や加速度）として扱えない．
    • WINフォーマットでの保存ができない

作成したライブラリ

wintoolsという名前にしているが，GitHubに同名のライブラリ（Windows関連）が山ほどあるので考え直した方が良いかも．

特徴

  • WINシステムと同じ年の処理(yy -> yyyy)方法を採用
  • チャンネルテーブルと同時にデータを読み込み
  • 複数の連続したファイルの同時読み込み・自動連結
  • WINフォーマットでのデータ書き込み
  • 読み込みはObsPyと同等
  • ObsPyとの相互変換（チャンネル名，成分名，観測点名の対応に悩み中）
  • hypomhに必要なstan, seisファイルなどの作成，読み込みにも対応．

公開予定

結構雑なコードでバグも多そうだし，一般公開予定なし．
それを承知したうえで興味のある方はX(twitter)などでご相談ください．
（必ずしも提供できるとは限りませんが...）

実例

winのサンプルデータを読み込んだ例を示す．

読み込み

チャンネルテーブルを合わせて与えることで，
観測点名などの情報が保持される．

import wintools
import time

fp = "./etc/991109.064607"
chtbl = "./etc/991109.064607.ch"

ts = time.perf_counter()
dat = wintools.read(
    fp,
    chtbl,
    )
te = time.perf_counter()
print("read time: ", te - ts)

出力:

> INFO|19:44:16|__readwin__| Given 1 files.
> INFO|19:44:16|__readwin__| Loading...
> INFO|19:44:17|read_chtbl| 19
> INFO|19:44:17|read_chtbl| 19
read time:  0.31686024599912344

読み込んだdatのデータ概要は次のようになっている．

    0: 0200 (ASO-U)	| fs: 100.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    1: 0201 (ASO-N)	| fs: 100.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    2: 0202 (ASO-E)	| fs: 100.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    3: 0206 (GNZ-U)	| fs: 100.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    4: 0207 (GNZ-N)	| fs: 100.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    5: 0208 (GNZ-E)	| fs: 100.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    6: 0209 (GNZ-wU)	| fs: 20.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    7: 020A (GNZ-wN)	| fs: 20.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    8: 020B (GNZ-wE)	| fs: 20.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    9: 020C (KRO-U)	| fs: 100.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   10: 020D (KRO-N)	| fs: 100.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   11: 020E (KRO-E)	| fs: 100.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   12: 0218 (KBH-U)	| fs: 100.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   13: 0219 (KBH-N)	| fs: 100.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   14: 021A (KBH-E)	| fs: 100.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   15: 0234 (NIK-U)	| fs: 100.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   16: 0235 (NIK-N)	| fs: 100.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   17: 0236 (NIK-E)	| fs: 100.0 Hz, unit: m/s	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)

一応読み込みの確認をした結果

ObsPy（v1.4.1）で読み込んだ結果と比較したところ，
波形の値は完全に一致していました．
恐らく正しく読み込めていると思います．

プロット

calibrateメソッドで波形の振幅を物理量に変換してからプロット．
（整数値のままでもプロットは可能）

dat.copy().calibrate().plot() #整数から物理量に変換


プログラムで読み込んだ波形のプロット．ASO-Uの絶対値振幅は最大で約1.8\times10^{-6} m/sで，後に示すWINシステム上での振幅と一致する．

ここで比較用にWINシステム上での波形も見てみる．
波形，振幅ともに正しく処理できているようだ．

WINシステム上での表示．ASO-Uの絶対値振幅は最大で約1.8\times10^{-6} m/s．

ObsPyとの相互変換

wintoolsからobspyへ

次のように1行でObsPyのStreamへ変換できます．

st = dat.to_obspy()

結果:

18 Trace(s) in Stream:
.ASO..U  | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.990000Z | 100.0 Hz, 6200 samples
.ASO..N  | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.990000Z | 100.0 Hz, 6200 samples
.ASO..E  | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.990000Z | 100.0 Hz, 6200 samples
.GNZ..U  | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.990000Z | 100.0 Hz, 6200 samples
.GNZ..N  | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.990000Z | 100.0 Hz, 6200 samples
.GNZ..E  | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.990000Z | 100.0 Hz, 6200 samples
.GNZ..wU | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.950000Z | 20.0 Hz, 1240 samples
.GNZ..wN | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.950000Z | 20.0 Hz, 1240 samples
.GNZ..wE | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.950000Z | 20.0 Hz, 1240 samples
.KRO..U  | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.990000Z | 100.0 Hz, 6200 samples
.KRO..N  | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.990000Z | 100.0 Hz, 6200 samples
.KRO..E  | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.990000Z | 100.0 Hz, 6200 samples
.KBH..U  | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.990000Z | 100.0 Hz, 6200 samples
.KBH..N  | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.990000Z | 100.0 Hz, 6200 samples
.KBH..E  | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.990000Z | 100.0 Hz, 6200 samples
.NIK..U  | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.990000Z | 100.0 Hz, 6200 samples
.NIK..N  | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.990000Z | 100.0 Hz, 6200 samples
.NIK..E  | 1999-11-09T06:45:47.000000Z - 1999-11-09T06:46:48.990000Z | 100.0 Hz, 6200 samples

obspyからwintoolsへ

wintoolsはwinフォーマットでの書き出しもできるので，
これを使えばobspyで読み込んだsacなどのファイルをwinへ変換できるようになると思います．
まだ実際には試していませんが，wintosacの逆，sactowinとして機能できるはずです．
ただし，チャンネル番号は適切に与えなければなりません．

win = wintools.from_obspy(st)

結果:

    0: None (ASO-U)	| fs: 100.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    1: None (ASO-N)	| fs: 100.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    2: None (ASO-E)	| fs: 100.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    3: None (GNZ-U)	| fs: 100.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    4: None (GNZ-N)	| fs: 100.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    5: None (GNZ-E)	| fs: 100.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    6: None (GNZ-wU)	| fs: 20.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    7: None (GNZ-wN)	| fs: 20.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    8: None (GNZ-wE)	| fs: 20.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
    9: None (KRO-U)	| fs: 100.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   10: None (KRO-N)	| fs: 100.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   11: None (KRO-E)	| fs: 100.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   12: None (KBH-U)	| fs: 100.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   13: None (KBH-N)	| fs: 100.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   14: None (KBH-E)	| fs: 100.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   15: None (NIK-U)	| fs: 100.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   16: None (NIK-N)	| fs: 100.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)
   17: None (NIK-E)	| fs: 100.0 Hz, unit: None	| 1999/11/09T06:45:47 - 1999/11/09T06:46:48 (62.0s)

Discussion