Pythonゼミ　第7回　FALMA三次元標定結果、散布図、等高線図

Last updated: 2021/08/28

1．FALMAの三次元標定結果

三次元標定結果の一例は以下のリンクでダウンロードしてください。
http://tingwu.info/pylab/data/python07_3DLoc1.dat

このファイルは一つ雷の標定結果です。ConTEXTで開いたら、最初の部分は下図になります。

一行目はこの雷の開始時間及び観測サイトの情報です。
二行目からは標定結果です。このファイルに全部2881行があるので、標定結果は2880個があります。一行（一つの標定結果）に左からは時間(ms), X(m), Y(m), Z(m)の結果です。すなわち、一行は一つ放射源の発生時刻と三次元位置情報です。雷は進展しながら、電磁波を放射するので、たくさんの放射源の時刻と位置を特定できれば、雷全体の三次元構造が分かります。

三次元の標定結果を示すには、基本的に以下の図を作ります。この図に五つの部分があります。 aはHeight (Z) - Time図、bはHeight - X図（南北方向に見てこの雷の断面図）、cは標定点数の高度分布、 dはX - Y図（平面図）、eはHeight (Z) - Y図（東西方向に見てこの雷の断面図）です。各図に一つの点は一つの放射源（標定点）を示す。そして、点の色は時間を示す（青は最初発生した点、赤は最後発生した点）。このような図を作れば、この雷の三次元構造及び進展の特徴は全部分かります。

2．標定結果散布図の作り方

例として、上図のb, d, eの部分だけの図を作ります。下図になります。全体の図は宿題として皆さんに作ってもらいます。

まず図の各部分の位置と大きさを決めるパラメータを設定する。

Figureを作って、座標を入れる

import matplotlib.pyplot as plt

#前回と同じように、できた図の形を見ながら、以下の各値を最適な値に少しずつ調整する
te = 0.03
be = 0.08
h2 = 0.55
vb = 0.06
h1 = 1-te-be-h2-vb

le = 0.11
re = 0.05
hb = 0.06
w1 = 0.52
w2 = 1-le-re-w1-hb

f = plt.figure(figsize=(7.5,6.5))

ax1 = f.add_axes([le, 1-te-h1, w1, h1])
ax2 = f.add_axes([le, be, w1, h2])
ax3 = f.add_axes([le+w1+hb, be, w2, h2])

標定結果を読み込む

xMat = []
yMat = []
zMat = []
tMat = []

fid = open('F:/Course/python/python07_3DLoc1.dat')
strLine = fid.readline()
while True:
    strLine = fid.readline()
    if len(strLine)<2:
        break
    tVal = float(strLine[0:10])
    xVal = float(strLine[11:18])/1000
    yVal = float(strLine[19:26])/1000
    zVal = float(strLine[27:32])/1000
    xMat.append(xVal)
    yMat.append(yVal)
    zMat.append(zVal)
    tMat.append(tVal)
fid.close()

三つの散布図を作る

ax1.scatter(xMat, zMat, marker='o', s=10, c=tMat, edgecolor='none', 
            alpha=0.5, cmap='jet')
ax2.scatter(xMat, yMat, marker='o', s=10, c=tMat, edgecolor='none', 
            alpha=0.5, cmap='jet')
ax3.scatter(zMat, yMat, marker='o', s=10, c=tMat, edgecolor='none', 
            alpha=0.5, cmap='jet')

sは点のサイズ
c=tMatはtMatの値で点の色を決める（c=zMatに変えてみて）
alphaは点の透明度（点の数が多いので、透明度を設定したほうがよりきれいになります）
cmapはカーラーマップ、例えば、’jet’の場合、色は青から赤まで変化する。 matplotlibにはたくさんのカーラーマップがあります。例えば、cmap=’spring’, ‘summer’, ‘autumn’, ‘winter’に設定して、結果を見てください。 FALMAの標定結果を表示する場合は’jet’のほうがいいです。

座標のlabelを設定する

ax1.set_ylabel('Height (km)')
ax2.set_xlabel('West-East (km)')
ax2.set_ylabel('South-North (km)')
ax3.set_xlabel('Height (km)')

三つの部分は同じデータを表示するので、同じ範囲に設定したほうがいい。
範囲の設定もできた図を見ながら、最適な値に調整する。
ただし、xlimMax-xlimMin=ylimMax-ylimMinに設定したほうがいい。

xlimMin = -27
xlimMax = -8
ylimMin = -14
ylimMax = 5
zlimMin = -0.5
zlimMax = 13

ax1.set_xlim([xlimMin, xlimMax])
ax1.set_ylim([zlimMin, zlimMax])
ax2.set_xlim([xlimMin, xlimMax])
ax2.set_ylim([ylimMin, ylimMax])
ax3.set_xlim([zlimMin, zlimMax])
ax3.set_ylim([ylimMin, ylimMax])

a b c を付けて、図を保存する。

ax1.text(0.03, 0.88, '(a)', transform = ax1.transAxes)
ax2.text(0.03, 0.93, '(b)', transform = ax2.transAxes)
ax3.text(0.04, 0.93, '(c)', transform = ax3.transAxes)

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

後ろのtransform = ax1.transAxesはどいう意味か考えてください。

3．等高線図

[データ]　岐阜市一年間の温度データ: python05-data

以下のような等高線図（正確に言うと塗りつぶし等高線図）を作ります。
横軸は時間（hour）、縦軸は月（month）、一時間×一か月ごとの温度の平均値は色で表示します。この図を見ると、温度一番高いのは7～8月の14～15時、一番低いのは1月の5～6時ということがすぐ分かります。
このような等高線図は普通の線グラフや棒グラフより情報量が多いので、よく使われています。

まずfigureを作って、二つaxを入れます。

import os
import matplotlib.pyplot as plt
import numpy as np

le = 0.08
re = 0.08
w2 = 0.05
hb = 0.03
w1 = 1-le-re-w2-hb

te = 0.1
be = 0.1
h = 1-te-be

f = plt.figure(figsize=(10, 6))
ax1 = f.add_axes([le, be, w1, h])
ax2 = f.add_axes([le+w1+hb, be, w2, h])

座標を作らなくてもカラーバーを付けることができるが、時々特別な場所にカラーバーを置きたい、カラーバーのサイズを細かく調整したい、などの場合は座標を作ったほうがいいです。

二つの2次元の行列を作ります。
横24（24時間）×縦12（12か月）
sumArrayは温度の和を保存する、countArrayはデータの個数を保存する。後ほど平均温度を計算する。

sumArray = np.zeros((12, 24))
countArray = np.zeros((12, 24))

データを読み込む

strRoot = 'F:/Course/python/python05-data/'

monthMat = os.listdir(strRoot)
monthMat.sort()
for n in range(12):
    strMonth = monthMat[n]
    monthVal = int(strMonth[4:6])
    
    strDir = strRoot+strMonth
    fileMat = os.listdir(strDir)
    fileMat.sort()
    for fileName in fileMat:
        fid = open(strDir+'/'+fileName)
        while True:
            strLine = fid.readline()
            if len(strLine)==0:
                break
            hourVal = int(strLine[0:2])
            val = float(strLine[6:12])
            
            sumArray[monthVal-1, hourVal] += val
            countArray[monthVal-1, hourVal] += 1

        fid.close()

平均温度を計算する（2次元の行列の割り算はnumpyのdivide関数で）
全体の最高最低温度も計算する（2次元行列の最大値と最小値の計算はnumpyのamaxとaminで）

aveArray = np.divide(sumArray, countArray)
maxTemp = np.amax(aveArray)
minTemp = np.amin(aveArray)

横軸（時間）と縦軸（月）の行列を作る

xMat = np.arange(24)
yMat = np.arange(1, 13)

contourf関数でプロット
contourfというのはcontour fill、すなわち、塗りつぶし等高線図です。
contourという関数もあるので、使ってみてください。

c = ax1.contourf(xMat, yMat, aveArray, 
                 levels = np.arange(int(minTemp), int(maxTemp)+2, 1), 
                 cmap='jet')

levels = np.arange(…) を levels = [5, 10, 15, 20, 25] に入れ替えてみて、levelsの役割を理解してください。

カラーバーを付ける

cb = plt.colorbar(c, cax = ax2)

括弧の中のcは行55のcです（このカラーバーは行55に作った等高線図の付属のカラーバーということを指定する）
cbはこのカラーバーの名前、このカラーバーを細かく設定するとき使う。
cax = ax2はこのカラーバーをax2に置くことを指定する。

labelなどを付けて図を保存する。
カラーバーのlabelを付けるときはax2じゃなくてcbのset_label関数を使う。

ax1.set_xlabel('Time (hour)')
ax1.set_ylabel('Month')
ax1.set_title('My first contour plot')
cb.set_label('Temperature ($^\circ$C)')

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

等高線図は地形図、放電密度図、レーダエコー図等によく使われます。以下の図は実際に論文に乗せた図です。データの処理はちょっと複雑なので、ここでは紹介しません。今後必要の場合、相談してください。

問題

1．資料最初の全体の三次元の図を作る。

2．以下のリンクの三次元標定結果を使います。
https://tingwu.info/pylab/data/python07_3DLoc2.dat

このファイルのすべての標定結果の平面図をプロットすると、以下の図になります。標定点は二ヶ所に分けられています。これは二つ雷が別の場所にほぼ同時に発生していたということです。こういう場合は、三次元標定結果の図を作るとき、二つ雷を分けて二つ図にプロットする必要があります。

二つ雷の標定点を分けて、以下のように二つの三次元図を作ってください。（一つプログラムで二つ図を作る必要はない）

雷1

雷2

3．5～10月と11～4月を分けて、一つのグラフに二つ等高線図を作る（それぞれ別のカラーバーを付ける）

補足

1．[3. 等高線図]に作った温度の等高線図について、図を保存するコードの前に以下のコードを入れれば、図のX軸の表示が変わる。

ax1.set_xticks(range(0, 25, 5))
strHourMat = [str(i).zfill(2)+':00' for i in range(0, 25, 5)]
ax1.set_xticklabels(strHourMat)

実は、X,Y軸の記載はどんなものに変更してもできます。例えば、さっきの三行を以下の三行に変えて実行してみてください。

axText = 'meaningless'
ax1.set_xticks(np.linspace(0, 24, len(axText), endpoint=False))
ax1.set_xticklabels(axText)

2．2次元行列の処理

hstackとvstackの違いおよび2次元行列の中の一部を取り出す方法を理解してください。

import numpy as np

x = np.arange(12)
x = x.reshape((3, 4))
print(x)

y = np.arange(6)
y = y.reshape((3, 2))
print(y)

z = np.hstack((x, y))
print(z)

a = z[:, 3:5]
print(a)
>>
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
[[0 1]
 [2 3]
 [4 5]]
[[ 0  1  2  3  0  1]
 [ 4  5  6  7  2  3]
 [ 8  9 10 11  4  5]]
[[ 3  0]
 [ 7  2]
 [11  4]]

import numpy as np

x = np.arange(12)
x = x.reshape((3, 4))

y = np.arange(8)
y = y.reshape((2, 4))

z = np.vstack((x, y))
print(z)

a = z[2:4, 1:3]
print(a)
>>
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]
 [ 0  1  2  3]
 [ 4  5  6  7]]
[[ 9 10]
 [ 1  2]]

Pythonゼミ 第7回 FALMA三次元標定結果、散布図、等高線図

1．FALMAの三次元標定結果

2．標定結果散布図の作り方

3．等高線図

問題

補足

Pythonゼミ　第7回　FALMA三次元標定結果、散布図、等高線図