Pythonゼミ　第8回　FALMA波形データ、複数の図を一括作る

Last updated: 2023/04/07

データダウンロード：
https://tingwu.info/pylab/data/python08-data.zip

解凍したら、下図のように12個のファイルがあります。

これはFALMAの12か所の観測サイト記録した同じ時刻の1秒間の波形データです。

ファイル名各部分の意味

例えば、AKKL1548422986.bz2
AKK: サイトネーム　（各ファイルのサイトネームは違います）
L: LF帯観測データ　（Hを付けるファイルも存在します。それはVHF帯あるいはHF帯の観測データです）
1548422986: Unix time （データの時刻。各ファイルの時刻は同じなので、同じ放電の波形を記録してるということが分かります。 Unix timeの意味は自分で調べてください。普通の時間との変換方法は後で紹介します。）
bz2：拡張子。1種の圧縮ファイルです。（pythonでは直接bz2ファイルを読み込むことができるので、解凍する必要はないです。）
一つのファイルに1秒間のデータが保存されています。サンプリングレートは１MHz、すなわち、1秒間に10⁶個の点が記録されます。

波形データと標定結果との関係

FALMAは複数の観測サイトに構成されます。各サイトに雷の電磁波を受信するアンテナがあります。そのアンテナと記録装置を使って、各サイトは雷の電磁波を受信できます。受信された電磁波は波形データになります。
各サイトの場所は違うので、雷との距離も違います。だから、雷が放射する電磁波は各サイトに到達する時刻にも少し差があります。その時間差を利用して、雷の発生時刻と場所を計算できます。計算の結果は標定結果になります。
FALMAの詳細は以下のページと論文を参考してください。
https://www1.gifu-u.ac.jp/~lrg/falma.html
https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2018GL077628

1. 1サイトのデータ（HIML1548422986.bz2）だけをプロットする

import bz2
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['agg.path.chunksize'] = 20000

1行目：bz2ファイルを処理するために、bz2モジュールをインポートする必要がある。
4行目：プロットする点は多いので、時々メモリエラーが出ます。この行を加えたら、メモリエラーを防ぐことができます。

f = plt.figure(figsize=(9,4))
ax = f.add_axes([0.12, 0.12, 0.82, 0.81])

座標を作る。

fid = bz2.BZ2File('F:/Course/python/python08-data/HIML1548422986.bz2')
dataMat = np.frombuffer(fid.read(), np.int16)
fid.close()

データを読み込む。FALMAの波形データは実はバイナリデータです。バイナリデータの読み込む方法はデータ保存のフォーマットによって違います。 FALMAの波形データは非常に簡単なフォーマットで保存されるので、この三行だけで読み込めます。 FALMAのすべての波形データはこのように読み込めます（具体的な意味は分からなくでも大丈夫です）。
読み込んだ後、1秒間の波形データはnumpy行列dataMatに保存される。前述のように、1秒間10⁶個の点がある。 dataMatの長さは10⁶ということを確認してください。

tMat = np.linspace(0, 1000, 1000000)

x軸（時間軸）の行列です。
y軸（dataMat）と同じ長さ（1000000）、最初の点の値は0、最後の点の値は1000（ミリ秒）です。すなわち、時間軸は0ミリ秒から1000ミリ秒までを示します。
もし時間軸の単位は秒の場合は、tMatはどう作りますか自分で考えてください。

ax.plot(tMat, dataMat)
ax.set_xlim([tMat[0], tMat[-1]])
ax.set_xlabel('Time (ms)')
ax.set_ylabel('Digital Unit')

plt.savefig('F:/Course/python/08/test1.png', dpi=300)

プロットして保存する。

1秒間の波形です。波形の変化は雷放電による電場の変化を示します。電場変化の単位はV/mですが、V/mを求めるためにはアンテナを校正する必要があります。ここはdigital unitで示しています。Digital unitはV/mと比例しますが、実際V/mはどのぐらいか校正しないと分からないです。

上図の中に、一番大きいなパルスは帰還雷撃の波形です。次は、帰還雷撃波形のピークの前後0.5 msの波形をプロットします。

import bz2
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['agg.path.chunksize'] = 20000


f = plt.figure(figsize=(9,4))
ax = f.add_axes([0.12, 0.12, 0.82, 0.81])


fid = bz2.BZ2File('F:/Course/python/python08-data/HIML1548422986.bz2')
dataMat = np.frombuffer(fid.read(), np.int16)
fid.close()

tMat = np.linspace(0, 1000, 1000000)

#以下の3行だけがさっきと違います
maxIndex = dataMat.argmax()
begIndex = maxIndex-500
endIndex = maxIndex+500

ax.plot(tMat[begIndex:endIndex], dataMat[begIndex:endIndex])
ax.set_xlim([tMat[begIndex], tMat[endIndex]])
ax.set_xlabel('Time (ms)')
ax.set_ylabel('Digital Unit')

plt.savefig('F:/Course/python/08/test2.png', dpi=300)

2. 複数の図を一括作る

データpython08-dataの中に、12サイトの波形データがあるので、12サイトの波形図を一括作る

以下はできた12個の図です。

import os
import bz2
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['agg.path.chunksize'] = 20000

f = plt.figure(figsize=(9,4))
ax = f.add_axes([0.12, 0.12, 0.82, 0.81])

#フォルダー中のデータファイル名の行列を作る
dataDir = 'F:/Course/python/python08-data/'
fileMat = os.listdir(dataDir)
fileMat.sort()

#図の保存先。存在してなければ、作る
picDir = 'F:/Course/python/08/testPic/'
if os.path.exists(picDir)==False:
    os.mkdir(picDir)

#tMatはすべての図は同じなので、ループの外に置く
tMat = np.linspace(0, 1000, 1000000)

for fileName in fileMat:
    fid = bz2.BZ2File(dataDir+fileName)
    dataMat = np.frombuffer(fid.read(), np.int16)
    fid.close()

    ax.clear() #plotする前にまず座標をclearする
    ax.plot(tMat, dataMat)
    ax.set_xlim([tMat[0], tMat[-1]])
    ax.set_xlabel('Time (ms)')
    ax.set_ylabel('Digital Unit')
    
    plt.savefig(picDir+fileName[0:14]+'.png', dpi=300)

問題

1．12サイトの波形を下図のように1つのグラフにプロットする。各波形の右上に当該波形の観測サイトの名前を表示する。（この図を見ると、各サイトの波形は似ているということが分かります）

2．各サイトの波形の84～85ミリ秒の部分（帰還雷撃）をプロットする。（この図を見ると、各サイトがこの期間雷撃を記録した時刻には少し差があることが分かります。）

3．データpython08-dataの中に、各サイトのデータにおいて、以下の図を作ってください。上は一秒間全体の波形、下は最大値を中心に1msの波形です。

一つプログラムでpython08-dataの中の12個のデータの図を作ってください。

できた12個の図

補足：Unix timeの変換

import datetime

unixTimeVal = 1548422986
strNatTime = datetime.datetime.fromtimestamp(unixTimeVal).strftime('%Y%m%d %H:%M:%S')
print(strNatTime)
>>
20190125 22:29:46

import datetime

unixTimeVal = 1548422986
unixTimeVal += 60
strNatTime = datetime.datetime.fromtimestamp(unixTimeVal).strftime('%Y%m%d%H%M%S')
print(strNatTime)
>>
20190125223046

import datetime
import time

strNatTime = '20200612103030'
unixSec = time.mktime(datetime.datetime.strptime(strNatTime, '%Y%m%d%H%M%S').timetuple())
print(unixSec)
>>
1591925430.0

Pythonゼミ 第8回 FALMA波形データ、複数の図を一括作る

1. 1サイトのデータ（HIML1548422986.bz2）だけをプロットする

2. 複数の図を一括作る

問題

補足：Unix timeの変換

Pythonゼミ　第8回　FALMA波形データ、複数の図を一括作る