前線など立体的に把握するための鉛直断面図

黒良龍太

2021年10月10日 10:41

今回は、鉛直断面図の作成についてです。

鉛直断面図を分析することで、ジェット気流から地上前線に連なる前線帯の構造や、大気安定度や鉛直シアに関わる空間的な変化などを読み取ることができ、平面図の天気図と合わせて活用することで、気象現象を立体的にイメージできます。

どのような現象を分析したいかで、この断面図に表示する物理量は変わります。例えば、雪なら温度が必須でしょう。前線帯の解析では温位や相当温位があればよく、気温は必須ではないでしょう。

コードの紹介は、これまでの説明してきた内容は省いて、下の目次の内容に絞りました。また最後に、Jupyter Notebook用のコードをダウンロードできるようにしています。表示する物理量を書き換えるなどして、試してみてください。

前線解析のための鉛直断面図

今回のサンプルは前線解析のために、温位(黒線)や相当温位(赤線)・東西風速(青線)の等値線、シェードでは収束発散、風を矢羽で表示することにします。

等温位線の傾きが大きい（垂直に近い）領域は、前線帯にあたる可能性があります。風の鉛直変化から暖気移流か寒気移流か判断ができ、前線が強化している状況かを定性的に診断することもできます。さらに、相当温位の鉛直変化から不安定領域を把握できます。これらを総合し、また平面図も活用してすることで、前線帯の位置や前線の詳細な状況を理解することができます。

今回作成する鉛直断面図のサンプルは、次の通りです。

2021年9月14日3時の南北鉛直断面図(43-25N,130E)
等温位線(黒実線)、等相当温位線(赤実線)、東西方向の等風速線(青線)、
発散（シェード）、風(矢羽)

複数の気圧面をまとめて読み込む

ここから作図するためのコードの説明です。これまでは1つの特定の気圧面データを読み込んできました。次に示す通り、grbs()関数では、l >= 300 とすることで、300hPa以上の気圧面を一括して読み込んでくれます。

## データ読み込み部
# 読み込むGRIB2形式GSMのファイル名
gsm_fn_t="Z__C_RJTD_{0:4d}{1:02d}{2:02d}{3:02d}0000_GSM_GPV_Rgl_FD{4:04d}_grib2.bin"
gr_fn= gsm_fn_t.format(i_year,i_month,i_day,i_hourZ,i_ft)
#
# データOpen
grbs = pygrib.open(data_fld + gr_fn)
#
# データ読み込み。 tagHpの等圧面から下部のデータを全て読み込みます。
grbHt = grbs(shortName="gh",typeOfLevel='isobaricInhPa',level=lambda l:l >= tagHp)
grbWu = grbs(shortName="u",typeOfLevel='isobaricInhPa',level=lambda l:l >= tagHp)
grbWv = grbs(shortName="v",typeOfLevel='isobaricInhPa',level=lambda l:l >= tagHp)
grbTm = grbs(shortName="t",typeOfLevel='isobaricInhPa',level=lambda l:l >= tagHp)
grbRh = grbs(shortName="r",typeOfLevel='isobaricInhPa',level=lambda l:l >= topRh)

三次元のデータセット作成

前回までは2次元のデータセットを作成しています。次のコードのように、levelを指定すれば、3次元のデータセットも作成できます。

この３次元のデータセットを利用して、これまで2次元と同じコードで物理量を新たに計算すると、全ての等圧面に対して計算されます。ここでは、露点温度と収束を計算しています。

## 露点温度や収束などを計算するために、データセットを作る
ds = xr.Dataset(
   {
       "Geopotential_height": (["level","lat", "lon"], cubeHt * units.meter),
       "temperature": (["level","lat", "lon"], cubeTm * units('K')),
       "relative_humidity": (["level","lat", "lon"], cubeRh * units('%')),
       "u_wind": (["level","lat", "lon"], cubeWu * units('m/s')),
       "v_wind": (["level","lat", "lon"], cubeWv * units('m/s')),
   },
   coords={
       "level": levels,
       "lat": lats,
       "lon": lons,
       "time": [grbHt[0].validDate],
   },
)
ds['Geopotential_height'].attrs['units'] = 'm'
ds['temperature'].attrs['units']='K'
ds['relative_humidity'].attrs['units']='%'
ds['u_wind'].attrs['units']='m/s'
ds['v_wind'].attrs['units']='m/s'
ds['level'].attrs['units'] = 'hPa'
ds['lat'].attrs['units'] = 'degrees_north'
ds['lon'].attrs['units'] = 'degrees_east'
#
## 必要な物理量を計算する　単位の指定を忘れないように
# 露点温度算出
ds['dewpoint_temperature'] = mpcalc.dewpoint_from_relative_humidity(
       ds['temperature'],ds['relative_humidity'])
ds['dewpoint_temperature'].attrs['units']='K'
#
# 収束・発散算出
ds['divergence'] = mpcalc.divergence(ds['u_wind'],ds['v_wind'])
ds['dewpoint_temperature'].attrs['units']='1/s'

鉛直断面図用データの作成

MetPyのcross_section()関数に断面図の両端２地点指定すると、三次元のデータセットから断面図用データセットを簡単に作成できます。

このデータセットからも新たな物理量を計算し、追加することもできます。下のコードでは、温位や相当温位、断面に平行・垂直成分の風速を計算しています。

## 断面図データセットの作成   
start = (45.0, 130.0)
end   = (25.0, 130.0)
cross = cross_section(ds, start, end)
#
## 描画するデータを算出し、データセットに加える
# 温位、相当温位、相対湿度、風を断面に平行・垂直成分を描画する
cross['Potential_temperature'] = mpcalc.potential_temperature(
   cross['level'] * units.hPa,
   cross['temperature'])
cross['Equivalent_Potential_temperature'] = mpcalc.equivalent_potential_temperature(
   cross['level'], cross['temperature'], cross['dewpoint_temperature'])
cross['u_wind'] = (cross['u_wind']).metpy.convert_units('knots')
cross['v_wind'] = (cross['v_wind']).metpy.convert_units('knots')
cross['t_wind'], cross['n_wind'] = mpcalc.cross_section_components(
   cross['u_wind'],cross['v_wind'])

鉛直断面図の描画

前回までの平面図のコードと同様に、シェードはcountourf()を、等値線はcontour()を、矢羽はberbs()を利用して描画します。

縦軸は、set_yscale('symlog')と指定して、ログスケールとします。

## 断面図描画
#
# 断面図左上に描画する500hPa面天気図の描画位置
map_x_y_width_height = [0.1255, 0.572, 0.18, 0.33]
#
# 図の大きさを指定                                                                         
fig = plt.figure(1, figsize=(16., 9.))
ax = plt.axes()
#
# シェード:収束・発散
clevs_div = [-4,-2,-1,1,2,4]
div_contour = ax.contourf(cross[ax_cross], cross['level'], cross['divergence']*1e5,
                        clevs_div, cmap=plt.cm.bwr, extend='both', alpha = 0.7)
div_colorbar = fig.colorbar(div_contour)               
#                                                                                          
# 温位(黒実線)とラベル 
theta_contour = ax.contour(cross[ax_cross], cross['level'],
                          cross['Potential_temperature'],
                          levels=np.arange(252, 450, 3), colors='k', linewidths=2)
theta_contour.clabel(theta_contour.levels[1::2], fontsize=8, colors='k', inline=1,
                    inline_spacing=8, fmt='%i', rightside_up=True, use_clabeltext=True)
#
# 相当温位(赤実線)とラベル                                 
etheta_contour = ax.contour(cross[ax_cross], cross['level'],
                           cross['Equivalent_Potential_temperature'],
                           levels=np.arange(252, 450, 3), colors='r', linewidths=1)
etheta_contour.clabel(etheta_contour.levels[1::2], fontsize=8, colors='r', inline=1,
                    inline_spacing=8, fmt='%i', rightside_up=True, use_clabeltext=True)
#
# 断面の平行成分の風速表示(青実線：マイナスは点線表示となる)
v_contour = ax.contour(cross[ax_cross], cross['level'],
                      cross['t_wind'],
                      levels=np.arange(-100, 100, 5), colors='b', linewidths=1)
v_contour.clabel(v_contour.levels[1::2], fontsize=8, colors='b', inline=1,
               inline_spacing=8, fmt='%i', rightside_up=True, use_clabeltext=True)
#                                                                                          
# 風の表示                                                                                 
#  矢羽を間引き指定                                                                        
wind_slc_vert = list(range(0, len(levels),1))
wind_slc_horz = slice(5, 100, 2)
if flag_wind == 0:
 ## Tangential / Normal Winds                                                             
 ax.barbs(cross[ax_cross][wind_slc_horz], cross['level'][wind_slc_vert],
        cross['t_wind'][wind_slc_vert, wind_slc_horz],
        cross['n_wind'][wind_slc_vert, wind_slc_horz], color='k')
else:
 ## U / V Winds                                                                           
 ax.barbs(cross[ax_cross][wind_slc_horz], cross['level'][wind_slc_vert],
        cross['u_wind'][wind_slc_vert, wind_slc_horz],
        cross['v_wind'][wind_slc_vert, wind_slc_horz], color='k')
#
# Y軸をlog指定、ラベル指定、目盛り指定                                                     
ax.set_yscale('symlog')                                                             
ax.set_yticklabels(np.arange(1000, 200, -100))
ax.set_ylim(cross['level'].max(), 100.0)
ax.set_yticks(np.arange(1000, 290, -100))
#                                                                                          
## 左上の500hPa高度の天気図表示用
#  定義:CRS and inset axes                                                             
data_crs = ds['Geopotential_height'].metpy.cartopy_crs
# 地図の表示位置指定                                                                                                                                        
ax_inset = fig.add_axes(map_x_y_width_height, projection=data_crs)
# 地図の表示範囲指定                                                                                                                 
ax_inset.set_extent(map_latlon)
# 500hPa高度Plot(高度線の指定)                                                             
ax_inset.contour(ds['lon'], ds['lat'], ds['Geopotential_height'].sel(level=500.),
                levels=np.arange(5100, 6000, 30), cmap='inferno')
# 断面の位置をPlot                                                                         
endpoints = data_crs.transform_points(ccrs.Geodetic(),
                *np.vstack([start, end]).transpose()[::-1])
ax_inset.scatter(endpoints[:, 0], endpoints[:, 1], c='k', zorder=2)
ax_inset.plot(cross['lon'], cross['lat'], c='k', zorder=2)
#
# 海岸線を描画                                                 
ax_inset.coastlines()
#
## タイトルや軸ラベル                                                                     
ax_inset.set_title('')
if flag_wind == 0:
 ax.set_title('GSM Cross-Section \u2013 {} to {} \u2013 Valid: {}\n'
            'Potential Temperature (K),EPT (K), Tangent Wind Speed (knots), '
            'Tangential/Normal Winds (knots), Divergence\n'
            'Inset: Cross-Section Path and 500 hPa Geopotential Height'.format(
                start, end, dt.strftime('%Y-%m-%d %H:%MZ')))
else:
 ax.set_title('GSM Cross-Section \u2013 {} to {} \u2013 Valid: {}\n'
            'Potential Temperature (K),EPT (K), Tangent Wind Speed (knots), '
            'Winds (knots), Divergence\n'
            'Inset: Cross-Section Path and 500 hPa Geopotential Height'.format(
                start, end, dt.strftime('%Y-%m-%d %H:%MZ')))
ax.set_ylabel('Pressure (hPa)')
ax.set_xlabel('Longitude (degrees east)')
div_colorbar.set_label('Divergence (1e-5 1/s)')
#
plt.show()

サンプルコード

これまで説明してきたコードは肝の部分を抜き出しています。鉛直断面図作成コード（Jupyter Notebook用）は次からダウンロードできます。

次回は、NetCDFデータの読み込みついて、ご紹介する予定です。

この記事が気に入ったらサポートをしてみませんか？