Pythonの練習6
数値計算のoutput等でよく使用されるnetcdf形式のファイルを図面にする練習について書く．
今回は過去に計算した気象モデルWRFによる2018年台風21号の中心気圧に関する図を作成する．
Python上では，時間・行・列という感じでデータが読み取られるので，LANDMASKやLON，LATの情報を読む際には，”0″を指定している（実際にはどの変数も1時刻なので0でよい）．
最終的に読んだ変数達を用いて2次元のcontourfとcontour図にしておしまい．
import numpy as np
from netCDF4 import Dataset
import matplotlib.pyplot as plt
# WRFOUTファイルのパス
wrfout_file = "./wrfout_d02_2018-09-04_04_00_00"
# WRFOUTファイルを開く
wrfout = Dataset(wrfout_file, "r")
# データの取得
p = wrfout.variables["PSFC"][:] # 気圧 [Pa]
H = wrfout.variables["HGT"][:] # ジオポテンシャル高度
T = wrfout.variables["T2"][:] # 2m高度気温
MASK = wrfout.variables["LANDMASK"][0,:,:] #Land mask
slp=p*(1-0.0065*H/(T+0.0065*H))**(-5.257)
slp=slp/100
#大きさが1の行列を削除する（あってもなくてもよいが，あった方がclearなコードになる）
slp = np.squeeze(slp)
# 緯度と経度の取得，最初から１列目を０とすると２列目・３列目だけを指定して取れる
lat = wrfout.variables["XLAT"][0, :, :] # 緯度
lon = wrfout.variables["XLONG"][0, :, :] # 経度
# プロット
plt.figure(figsize=(10, 8))
plt.contour(lon, lat, MASK, colors='black')
plt.contourf(lon, lat, slp, np.linspace(950,1000,100), cmap="jet", extend='both')
plt.colorbar(label="Sea level pressure",ticks=np.linspace(950,1000,5)) #範囲を指定して，その間を５個の数値で表現
plt.xlim(130,138)
plt.ylim(30,35.8)
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.title("2D Pressure Distribution")
plt.grid()
plt.savefig("test-wrfout.png")
plt.show()
# WRFOUTファイルを閉じる
wrfout.close()
