【気象データ】数値予報モデルが出力したGPVデータの取得・読み込み・描画【GRIB】

$ pip3 install pygrib
$ pip3 install metpy

import datetime

import pygrib
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.patheffects as path_effects
from metpy.units import units
import numpy as np
from scipy.ndimage import gaussian_filter


#------------------------------------------------------------------------------
# 読み込むデータのファイルパス
fname_gfs = "gfs_4_20221224_0000_012.grb2"

# 最大，最小の緯度と経度
lon_min = 90
lon_max = 180
lat_min = 10
lat_max = 70

#------------------------------------------------------------------------------
# NOAA GFS 予報データの読み込みと平滑化など
# データの読み込み
print("*** GFSデータを読み込みます：", fname_gfs)
with pygrib.open(fname_gfs) as gribin:
    # 500 hPa 温度
    tmpr_500, lat, lon = gribin.select(name='Temperature', level=500)[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)
    
    # 初期値の日付と予報時間
    dtm_ini_gfs = gribin.select(name='Temperature', level=500)[0].analDate
    dtm_fcst_gfs = gribin.select(name='Temperature', level=500)[0].validDate
    
    # 850 hPa 温度
    tmpr_850, _, _ = gribin.select(name='Temperature', level=850)[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)


# 予報時間を計算
fcst_hrs = dtm_fcst_gfs - dtm_ini_gfs

# ガウシアンフィルタで平滑化する
tmpr_500 = gaussian_filter(tmpr_500, sigma=1.0)
tmpr_850 = gaussian_filter(tmpr_850, sigma=1.0)

# 単位をつける
lat = lat * units('degrees_north')
lon = lon * units('degrees_east')
tmpr_500 = tmpr_500 * units('K')
tmpr_850 = tmpr_850 * units('K')


#------------------------------------------------------------------------------
# プロット
fig = plt.figure(figsize=(210/25.4, 294/25.4), dpi=100)
out_fig_path = f'plot-gfs-ini{dtm_ini_gfs.strftime("%Y%m%d%H")}-fcst{dtm_fcst_gfs.strftime("%Y%m%d%H")}-500-850-temperature.png'

fig.text(0.10, 0.97, f"NOAA GFS (resolution 0.5 deg)", fontsize=14)
fig.text(0.10, 0.95, f"initial: {dtm_ini_gfs.strftime('%Y/%m/%d %H')}UTC", fontsize=14)
fig.text(0.10, 0.93, f"forecast: {dtm_fcst_gfs.strftime('%Y/%m/%d %H')}UTC (initial + {fcst_hrs/datetime.timedelta(hours=1)}hrs)", fontsize=14)

# 地図用の変換
crs_map = ccrs.NorthPolarStereo(central_longitude=140.0, true_scale_latitude=60.0) # 気象庁「日々の天気図」に近いもの

# データ用の変換
crs_data = ccrs.PlateCarree()

# プロット枠の設定
ax01 = fig.add_axes([0.10, 0.50, 0.80, 0.40], projection=crs_map)
cb01 = fig.add_axes([0.80, 0.50, 0.01, 0.40])
ax02 = fig.add_axes([0.10, 0.05, 0.80, 0.40], projection=crs_map)
cb02 = fig.add_axes([0.80, 0.05, 0.01, 0.40])


# ----- プロット その1 500hPa 気温
# 海岸線や国境を描画する
ax01.add_feature(cfeature.COASTLINE, linewidth=0.5, color="green")
ax01.add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':', color="green")

# 地図を描く緯度経度の範囲の設定
ax01.set_extent([117.0, 156.0, 18.0, 61.0], crs_data)

# 500 hPa 気温
clevs_tmpr_500 = np.arange(-54, 15, 3)
cs01 = ax01.contour(lon, lat, tmpr_500.to(units.degC), clevs_tmpr_500, 
                    colors='black', transform=crs_data, linewidths=0.5, linestyles='dashed')
clbls01 = ax01.clabel(cs01, fmt='%d')

plt.setp(clbls01, path_effects=[path_effects.withStroke(linewidth=2, foreground="w")])

# 雪の目安の等温線を塗る
ax01.contour(lon, lat, tmpr_500.to(units.degC), [-36.0, -30.0],
             colors=['purple', 'blue'], transform=crs_data, linewidths=[3.0], linestyles='solid', alpha=0.5)


# cflevs_tmpr_500 = np.arange(-24, 0, 6) # 夏用
cflevs_tmpr_500 = np.arange(-54, 3, 3) # 冬用
norm = mpl.colors.BoundaryNorm(cflevs_tmpr_500, 256)

cf = ax01.contourf(lon, lat, tmpr_500.to(units.degC),
                   cflevs_tmpr_500, norm=norm,
                   cmap='PuOr_r', transform=crs_data, extend='both', alpha=0.5)
mpl.colorbar.ColorbarBase(cb01, cf, ticks=cflevs_tmpr_500,
                          orientation='vertical', 
                          label=u'temperature at 500hPa ($^\circ$C)',
                          ).minorticks_off()

# グリッドを描く
ax01.gridlines(xlocs=np.arange(100, 180, 10), 
               ylocs=np.arange(10, 90, 10),
               linestyle='dashed', color='green', linewidth=0.5)


# ----- プロット その2 850hPa 気温
# 海岸線や国境を描画する
ax02.add_feature(cfeature.COASTLINE, linewidth=0.5, color="green")
ax02.add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':', color="green")

# 地図を描く緯度経度の範囲の設定
ax02.set_extent([117.0, 156.0, 18.0, 61.0], crs_data)

# 850 hPa 気温
clevs_tmpr_850 = np.arange(-45, 30, 3)
cs02 = ax02.contour(lon, lat, tmpr_850.to(units.degC), clevs_tmpr_850, 
                    colors='black', transform=crs_data, linewidths=0.5, linestyles='dashed')
clbls02 = ax02.clabel(cs02, fmt='%d')

plt.setp(clbls02, path_effects=[path_effects.withStroke(linewidth=2, foreground="w")])

# 雪の目安の等温線を塗る
ax02.contour(lon, lat, tmpr_850.to(units.degC), [-12.0, -6.0],
             colors=['purple', 'blue'], transform=crs_data, linewidths=[3.0], linestyles='solid', alpha=0.5)

# cflevs_tmpr_850 = np.arange(-24, 24, 6) # 夏用
cflevs_tmpr_850 = np.arange(-33, 33, 3) # 冬用
norm = mpl.colors.BoundaryNorm(cflevs_tmpr_850, 256)

cf = ax02.contourf(lon, lat, tmpr_850.to(units.degC),
                   cflevs_tmpr_850, norm=norm,
                   cmap='PuOr_r', transform=crs_data, extend='both', alpha=0.5)
mpl.colorbar.ColorbarBase(cb02, cf, ticks=cflevs_tmpr_850,
                          orientation='vertical', 
                          label=u'temperature at 850hPa ($^\circ$C)',
                          ).minorticks_off()

# グリッドを描く
ax02.gridlines(xlocs=np.arange(100, 180, 10), 
               ylocs=np.arange(10, 90, 10),
               linestyle='dashed', color='green', linewidth=0.5)


#------------------------------------------------------------------------------
# 図の保存
plt.savefig(out_fig_path, transparent=False)

#------------------------------------------------------------------------------
# 読み込むデータのファイルパス
fname_gfs = "gfs_4_20221224_0000_048.grb2"

$ cd Desktop/LabCode/python/gfs-temperature

$ python3 plot-gfs-tempeature.py

import datetime

import pygrib
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.patheffects as path_effects
from metpy.units import units
import numpy as np
from scipy.ndimage import gaussian_filter

#------------------------------------------------------------------------------
# 読み込むデータのファイルパス
fname_gfs = "gfs_4_20221224_0000_048.grb2"

# 最大，最小の緯度と経度
lon_min = 90
lon_max = 180
lat_min = 10
lat_max = 70

#------------------------------------------------------------------------------
# NOAA GFS 予報データの読み込みと平滑化など
# データの読み込み
print("*** GFSデータを読み込みます：", fname_gfs)
with pygrib.open(fname_gfs) as gribin:
    # 500 hPa 温度
    tmpr_500, lat, lon = gribin.select(name='Temperature', level=500)[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)
    
    # 初期値の日付と予報時間
    dtm_ini_gfs = gribin.select(name='Temperature', level=500)[0].analDate
    dtm_fcst_gfs = gribin.select(name='Temperature', level=500)[0].validDate
    
    # 850 hPa 温度
    tmpr_850, _, _ = gribin.select(name='Temperature', level=850)[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)


# 予報時間を計算
fcst_hrs = dtm_fcst_gfs - dtm_ini_gfs

# ガウシアンフィルタで平滑化する
tmpr_500 = gaussian_filter(tmpr_500, sigma=1.0)
tmpr_850 = gaussian_filter(tmpr_850, sigma=1.0)

# 単位をつける
lat = lat * units('degrees_north')
lon = lon * units('degrees_east')
tmpr_500 = tmpr_500 * units('K')
tmpr_850 = tmpr_850 * units('K')

#------------------------------------------------------------------------------
# プロット
fig = plt.figure(figsize=(210/25.4, 294/25.4), dpi=100)
out_fig_path = f'plot-gfs-ini{dtm_ini_gfs.strftime("%Y%m%d%H")}-fcst{dtm_fcst_gfs.strftime("%Y%m%d%H")}-500-850-temperature.png'

fig.text(0.10, 0.97, f"NOAA GFS (0p50)", fontsize=14)
fig.text(0.10, 0.95, f"initial:{dtm_ini_gfs.strftime('%Y/%m/%d %H')}UTC", fontsize=14)
fig.text(0.10, 0.93, f"forecast:{dtm_fcst_gfs.strftime('%Y/%m/%d %H')}UTC (initial + {fff_gfs}hrs)", fontsize=14)

# 地図用の変換
crs_map = ccrs.NorthPolarStereo(central_longitude=140.0, true_scale_latitude=60.0) # 気象庁「日々の天気図」に近いもの

# データ用の変換
crs_data = ccrs.PlateCarree()

# プロット枠の設定
ax01 = fig.add_axes([0.10, 0.50, 0.80, 0.40], projection=crs_map)
cb01 = fig.add_axes([0.80, 0.50, 0.01, 0.40])
ax02 = fig.add_axes([0.10, 0.05, 0.80, 0.40], projection=crs_map)
cb02 = fig.add_axes([0.80, 0.05, 0.01, 0.40])

# ----- プロット その1 500hPa 気温
# 海岸線や国境を描画する
ax01.add_feature(cfeature.COASTLINE, linewidth=0.5, color="green")
ax01.add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':', color="green")

# 地図を描く緯度経度の範囲の設定
ax01.set_extent([117.0, 156.0, 18.0, 61.0], crs_data)

# 500 hPa 気温
clevs_tmpr_500 = np.arange(-54, 15, 3)
cs01 = ax01.contour(lon, lat, tmpr_500.to(units.degC), clevs_tmpr_500, 
                    colors='black', transform=crs_data, linewidths=0.5, linestyles='dashed')
clbls01 = ax01.clabel(cs01, fmt='%d')

plt.setp(clbls01, path_effects=[path_effects.withStroke(linewidth=2, foreground="w")])

# 雪の目安の等温線を塗る
ax01.contour(lon, lat, tmpr_500.to(units.degC), [-36.0, -30.0],
             colors=['purple', 'blue'], transform=crs_data, linewidths=[3.0], linestyles='solid', alpha=0.5)


# cflevs_tmpr_500 = np.arange(-24, 0, 6) # 夏用
cflevs_tmpr_500 = np.arange(-54, 3, 3) # 冬用
norm = mpl.colors.BoundaryNorm(cflevs_tmpr_500, 256)

cf = ax01.contourf(lon, lat, tmpr_500.to(units.degC),
                   cflevs_tmpr_500, norm=norm,
                   cmap='PuOr_r', transform=crs_data, extend='both', alpha=0.5)
mpl.colorbar.ColorbarBase(cb01, cf, ticks=cflevs_tmpr_500,
                          orientation='vertical', 
                          label=u'temperature at 500hPa ($^\circ$C)',
                          ).minorticks_off()

# グリッドを描く
ax01.gridlines(xlocs=np.arange(100, 180, 10), 
               ylocs=np.arange(10, 90, 10),
               linestyle='dashed', color='green', linewidth=0.5)

# ----- プロット その2 850hPa 気温
# 海岸線や国境を描画する
ax02.add_feature(cfeature.COASTLINE, linewidth=0.5, color="green")
ax02.add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':', color="green")

# 地図を描く緯度経度の範囲の設定
ax02.set_extent([117.0, 156.0, 18.0, 61.0], crs_data)

# 850 hPa 気温
clevs_tmpr_850 = np.arange(-45, 30, 3)
cs02 = ax02.contour(lon, lat, tmpr_850.to(units.degC), clevs_tmpr_850, 
                    colors='black', transform=crs_data, linewidths=0.5, linestyles='dashed')
clbls02 = ax02.clabel(cs02, fmt='%d')

plt.setp(clbls02, path_effects=[path_effects.withStroke(linewidth=2, foreground="w")])

# 雪の目安の等温線を塗る
ax02.contour(lon, lat, tmpr_850.to(units.degC), [-12.0, -6.0],
             colors=['purple', 'blue'], transform=crs_data, linewidths=[3.0], linestyles='solid', alpha=0.5)

# cflevs_tmpr_850 = np.arange(-24, 24, 6) # 夏用
cflevs_tmpr_850 = np.arange(-33, 33, 3) # 冬用
norm = mpl.colors.BoundaryNorm(cflevs_tmpr_850, 256)

cf = ax02.contourf(lon, lat, tmpr_850.to(units.degC),
                   cflevs_tmpr_850, norm=norm,
                   cmap='PuOr_r', transform=crs_data, extend='both', alpha=0.5)
mpl.colorbar.ColorbarBase(cb02, cf, ticks=cflevs_tmpr_850,
                          orientation='vertical', 
                          label=u'temperature at 850hPa ($^\circ$C)',
                          ).minorticks_off()

# グリッドを描く
ax02.gridlines(xlocs=np.arange(100, 180, 10), 
               ylocs=np.arange(10, 90, 10),
               linestyle='dashed', color='green', linewidth=0.5)

#------------------------------------------------------------------------------
# 図の保存
plt.savefig(out_fig_path, transparent=False)