GRIB2 data の 3D 表示 - 2 - 風を重ねてみた

　先日の続きで、風を 3D 表示してみました。matplotlib mplot3d には "Quiver" という矢筒を描く機能があったので楽勝かと思っていたのですが…。

１．mplot3d "Quiver" について

２．気象庁 MSM データ

３．Pygrib で緯度経度を指定する

４．データの型の確認

５．結果はこんな感じ

６．Code はこんな感じです

７．まとめ

１．mplot3d "Quiver" について

ax.quiver(X, Y, Z, U, V, W, **kwargs)

　X、Y、Z が座標、U、V、W が風の座標軸成分を代入するだけでグラフ中に風の矢羽根を描いてくれます。何でも出来るんですね。

２．気象庁 MSM データ

　気象庁の MSM データに含まれるのは、配信資料に関する技術資料や GRIB2 ファイルを扱った色々な方のサイトでも分かるように、気圧面高度、風の東西、南北水平成分（u, v m/s）および鉛直 p 速度（dP/dt Pa/s）、そして温度、相対湿度です。

　ということで、mplot3d の Z 成分をどうしようかなと思ったのですが、今回はとりあえず無視し、ざっくり「0」にしました。

　というのも、鉛直 p 速度から静力学平衡の仮定の下、試しに求めてみたのですが、水平成分よりも圧倒的に小さいですし、そもそも台風で静力学平衡を仮定するのも微妙ですし、今回は描画のお試し優先です。

３．Pygrib で緯度経度を指定すると、データの型が変わる

　これまでは地図に気象データを重ねていたので、描画する緯度経度はすべて地図を扱うライブラリで指定していました。今回は matplotlib のみを使っているので、Pygrib で読み込む時に指定しています。前回はやらなかったので気づかなかったのですが、今回風を描画するのに読み込むデータの step を指定した際、Error が返ってきます。というのもどうやら、

　「Pygrib で緯度経度を指定すると、データの型が変わるらしい…」

　なんでこうなるかな、と思いながら公式ページに助けを求めると、どうやら以下の文が説明に該当するようです。

　"If the grid type is unprojected lat/lon and a geographic subset is requested (by using the lat1,lat2,lon1,lon2 keywords), then 2-d arrays are returned, otherwise 1-d arrays are returned."

　残念ですが…、文を読んでも分かりません（泣）。でも、返されるデータの型が２通りあり得ることは分かったので、さっそくデータの中身を確認することとしました。

４．データの型の確認

　読み込んだ GRIB2 データを print したところ、以下のようでした。

①緯度経度を指定した場合

(array([[-0.5148735, -0.6711235, -0.7336235, ..., -2.2648735, -1.2336235, -0.5148735], ....,
 [ 3.2663765, 2.9226265, 2.9538765, ..., 3.7663765, 3.8601265,
 3.8913765]]), 
 array([[33. , 33. , 33. , ..., 33. , 33. , 33. ], ....,
 [24.1, 24.1, 24.1, ..., 24.1, 24.1, 24.1]]), 
 array([[121. , 121.125, 121.25 , ..., 131.75 , 131.875, 132. ], ....,
 [121. , 121.125, 121.25 , ..., 131.75 , 131.875, 132. ]]))

　緯度経度を指定すると、各気象データ、緯度、経度の順にタプルになっています。前回等圧面を描画した時には、このデータのみを見て、タプルの要素をそのまま使っていたので、特に問題がありませんでした。逆に言えば、このデータを numpy.[::10] 等で抽出しようとしたから Error なのでした。

② GRIB2 データをそのまま読み込んだ場合

[[ 3.7976265 3.8913765 3.7976265 ... -0.0773735 0.0788765 0.2038765]
 ...
 [ 8.3601265 8.4538765 8.4538765 ... -3.8586235 -3.9211235 -3.9211235]]

　緯度×経度の配列の形で、各気象データが格納されています。こちらの配列なら、step 指定の抽出が使えるわけです。

　要するに、風の矢羽根のためには、タプルの各要素にしてから step にすればいいということです。ライブラリにお任せではなく、自分が何を扱っているのかきちんと把握しましょう、ということですね。

５．結果はこんな感じ

f:id:ohigehige:20201005124502p:plain

　950 hPa 等圧面の風と、参考に概ね台風の中心気圧の、975 hPa 等高線を図示してみました。確かに台風の鉛直構造が窺えます。

６．Code はこんな感じです　

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
 
import pygrib
import numpy as np

# Set parameter 
dat = "0824"
tm ="06"
filenm = "Z__C_RJTD_2020"+dat+tm+"0000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin"
grbs = pygrib.open(filenm)
ms2kts = 0.5144

# Read data from GRIB2 file
req_data_1= grbs.select(forecastTime=0, level=950)[0]
wu_1 = grbs.select(forecastTime=0, level=950)[1]
wv_1 = grbs.select(forecastTime=0, level=950)[2]
height_2 = grbs.select(forecastTime=0, level=975)[0]

# Extract the data by LAT/LON
hgt_1 = req_data_1.data(22,34,121,132)
hgt_2 = height_2.data(22,34,121,132)
wu_1 = wu_1.data(22,34,121,132)
wv_1 = wv_1.data(22,34,121,132)

# Thin out wind data
x = hgt_1[2][::10, ::10]
y = hgt_1[1][::10, ::10]
z = hgt_1[0][::10, ::10]
u = wu_1[0][::10, ::10]
v = wv_1[0][::10, ::10]
w = u*0

fig = plt.figure(figsize=(14,7))
ax1 = fig.add_subplot(111, projection='3d')

ax1.quiver(x,y,z,u,v,w, length = 1, arrow_length_ratio=0.5, normalize=True)
ax1.contour(hgt_2[2], hgt_2[1], hgt_2[0], linewidth=0)

ax1.set_title("Typhoon #8 (Bavi) 950 hPa Wind\n (Geo-potential height reference line 975 hPa)", fontsize=14)
ax1.set_xlabel('\nLongitude\n (2 deg = about 222 km)', fontsize=12) 
ax1.set_ylabel('Latitude', fontsize=12) 
ax1.set_zlabel('\nGeo-potential Height\n (m)', fontsize=12) 
ax1.set_box_aspect((1,1,0.6))

fig.show()