GSMの500hPa面天気図作成 2
この第２回では、第１回で準備した気象庁全球モデルの格子点データセットを図にするコードを紹介します。
等値線やハッチ描画、地図の指定、グリッド描画などの参考になると思います。ただし、今回は、気象庁の数値予報天気図にある低気圧や高気圧のマークなどの表示について紹介できていません。第3回にご紹介します。
なお、第3回にJupyter Notebook用のコードを公開しています。適宜ご利用ください。
描画の準備
次のコードで、表示する図に必要なパラメーター、地図の図法や表示範囲、表示する等値線の高度や渦度の値などを指定します。
## 表示用の指定
# 図法指定 気象庁の数値予報資料と同じ図法に設定
proj = ccrs.Stereographic(central_latitude=60, central_longitude=140)
latlon_proj = ccrs.PlateCarree() # 緯度経度の処理用に正距円筒図法も使う
#
# 地図の描画範囲指定
areaAry = [108, 156, 17, 55] # 極東
#areaAry = [115, 151, 20, 50] # 日本付近
#
# 等値線の間隔を指定
levels_ht =np.arange(4800, 6000, 60) # 高度を60m間隔で実線
levels_ht2=np.arange(4800, 6000, 300) # 高度を300m間隔で太線
levels_vr =np.arange(-0.0002, 0.0002, 0.00004) # 渦度4e-5毎に等値線
#
# 渦度のハッチの指定
levels_h_vr = [0.0, 0.00008, 1.0] # 0.0以上で 灰色(0.9), 8e-5以上で赤
colors_h_vr = ['0.9','red']
alpha_h_vr = 0.3 # 透過率を指定
#
# 緯度・経度線の指定
dlon,dlat=10,10 # 10度ごとに
#
## タイトル文字列用
# 予想時間を得る
ft_hours=int(i_ft/100) * 24 + int(i_ft%100)
# 初期時刻の文字列
dt_i = grbHt.analDate
dt_str = (dt_i.strftime("%H00UTC%d%b%Y")).upper()
dt_str2 = dt_i.strftime("%Y%m%d%H")図のサイズなどの指定
図の大きさや余白などを指定するコードは次の通りです。
## 図のSIZE指定inch
fig = plt.figure(figsize=(10,8))
#
## 余白設定
plt.subplots_adjust(left=0, right=1, bottom=0.06, top=0.98)
#
## 作図開始
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.set_extent(areaAry, latlon_proj)等値線やハッチの描画
次のコードで、相対渦度のハッチ描画、高度の等値線(太線、実線、特定の高度5820m(茶色)と5400m(青色))描画を行なっています。等値線の値表示も指定しています。しかし、値が表示されない等値線が出ることがあります。バグなのか仕様なのかわかっていません。
# 500hPa 相対渦度のハッチ 0.0以上:着色 0.8*10**-4以上:赤
cn_relv_hatch2 = ax.contourf(ds['lon'], ds['lat'], ds['vorticity'],
levels_h_vr, colors=colors_h_vr,
alpha=alpha_h_vr, transform=latlon_proj)
# 5000hPa 相対渦度 実線 0.00004毎 負は破線
cn_relv = ax.contour(ds['lon'], ds['lat'], ds['vorticity'],
levels_vr, colors='black', linewidths=1.0, transform=latlon_proj)
# 500hPa 等高度線 実線 step1:60m毎
cn_hgt = ax.contour(ds['lon'], ds['lat'], ds['Geopotential_height'],
colors='black',
linewidths=1.2, levels=levels_ht, transform=latlon_proj )
ax.clabel(cn_hgt, levels_ht, fontsize=15, inline=True, inline_spacing=5,
fmt='%i', rightside_up=True)
# 500hPa 等高度線 太線 step1:300m毎
cn_hgt2= ax.contour(ds['lon'], ds['lat'], ds['Geopotential_height'],
colors='black',
linewidths=1.5, levels=levels_ht2, transform=latlon_proj)
ax.clabel(cn_hgt2, fontsize=15, inline=True, inline_spacing=0,
fmt='%i', rightside_up=True)
# 500hPa 5820gpm高度線 茶色
cn5820 = ax.contour(ds['lon'], ds['lat'], ds['Geopotential_height'],
colors='brown',linestyles='dashdot',
linewidths=1.2, levels=[5820], transform=latlon_proj )
# 500hPa 5400gpm高度線 青色
cn5400 = ax.contour(ds['lon'], ds['lat'], ds['Geopotential_height'],
colors='blue',linestyles='dashdot',
linewidths=1.2, levels=[5400], transform=latlon_proj )海岸線やグリッド、タイトル描画
海岸線や緯度線、経度線、タイトルの描画などのコードを、次に示します。
## 海岸線
ax.coastlines(resolution='50m',)
## グリッド
xticks=np.arange(0,360.1,dlon)
yticks=np.arange(-90,90.1,dlat)
gl = ax.gridlines(crs=ccrs.PlateCarree()
, draw_labels=False
, linewidth=1, alpha=0.8)
gl.xlocator = mticker.FixedLocator(xticks)
gl.ylocator = mticker.FixedLocator(yticks)
#
## Title
fig.text(0.5,0.01,
"GSM FT{0:d} IT:".format(ft_hours)+dt_str+" 500hPa Height(m),VORT" ,
ha='center',va='bottom', size=18)
# 出力
plt.show()このコードで出力される図は次の通りとなります。
画像をファイル出力するには、次のコードを利用します
## ファイル出力
out_fn="gsm_{0}UTC_FT{1:03d}_50.png".format(dt_str2,ft_hours)
plt.savefig(out_fn)第1回と第2回のコードについて、パラメータなど変更すれば、気圧面の気温などの等値線なども比較的簡単に作成できます。試してみてください。
次回は、極大値や極小値の位置にマークをつけるために必要なコードを紹介します。
この記事が気に入ったらサポートをしてみませんか？