【MetPy】気象庁でも利用している (らしい) 前線解析のやり方【気象予報】

【MetPy】気象庁でも利用している (らしい) 前線解析のやり方【気象予報】

はじめに


気象予報にとって,前線の位置を解析したり,前線付近でこの後対流活動が活発化するか否かを判定する作業は極めて重要な位置を占めています。実際に,気象予報士試験では,資料を用いてこの判定をさせる問題が出題されます。

気象庁が行う前線解析の基本的な考え方は,予報技術に関する資料集のページにある「アジア太平洋地上天気図における前線解析」という資料にとてもわかり易くまとめられています (私が受験するときにも欲しかった…)。

この資料にあるように,基本的には,わたしたちのような非気象庁職員でも,気象庁から配信される数値予報天気図等を用いて前線の解析することができます。

ところが,前掲の資料の24ページから25ページに,気になる内容が書かれています。

前線の解析においては、前節で述べた様々な観測データに加えて、数値予報結果から客観的に前線と考えられる場所を計算した「客観前線」プロダクトも参考資料として利用している。

残念ながら,ここにある資料はわれわれ素人には配信されていません。今回は,ちょっとズルいなと思ったので,このような客観解析をすることができる pythonのプログラムを作成してみたいと思います!

動作検証済み環境

macOS Monterey 12.6.2, Python 3.9.15, MetPy 1.3.1

気象庁が行っている (らしい) 前線の客観解析


目標は,資料にあるような図を作ることです。

出典:気象庁「アジア太平洋地上天気図における前線解析」25ページ

気象庁予報課では,GSMによる925hPaおよび950hPaの相当温位を利用したTFPと850hPaの気温,地上の等圧線を重ね合わせたものを利用しているようです。

TFPとは,Thermal Front Parameter (サーマル・フロント・パラメタ)の頭文字で,次のように定義されます:

$$ \mathrm{TFP} = -\nabla_\mathrm{h} |\nabla_\mathrm{h}\tau|\cdot\frac{\nabla_\mathrm{h}\tau}{|\nabla_\mathrm{h}\tau|} $$

ここに,$\tau$ は何らかの熱力学的変数です。具体的には,気温や相当温位が使われます。$\nabla_\mathrm{h}$ は水平方向の微分演算子を意味します。

「相当温位を利用したTFP」とは,「熱力学変数 $\tau$ として,相当温位 $\theta_\mathrm{e}$ を使う」という意味です。

また,ジェット軸と前線帯との対応関係を確認するために,GSMによる300hPaの風速場分布図と500hPaの高度場分布図を重ね合わせたものを用いているようです。

GSMの結果は利用できないではないですが,これまでのプログラムを利用してGFSを使って似たようなプロット,つまり,

  • 925hPaの相当温位を利用したTFPと850hPaの気温,地上の等圧線を重ね合わせた図
  • 300hPaの風速と500hPaの高度,地上の等圧線を重ね合わせた図

を作ってみたいと思います。

なお,GFSの読み込みには,pygribが必要です。pygribのインストールについては,以前の記事を参照してください。

TFPの計算


TFPの定義式を見るとわかるように,TFPを計算するには,相当温位を計算し,さらに,その水平傾度を計算しなければなりません。

MetPyにはこのような気象要素を計算するのに必要な関数が揃っています。それらを組み合わせて計算します。MetPyのインストールについては,以前の記事を参照してください。

相当温位の計算

相当温位は,metpy.calc.equivalent_potential_temperature(pressure, temperature, dewpoint) で計算することができます。pressure (気圧) と temperature (温度) は GFS を読み込むと得られます。

dewpoint は露点温度で,metpy.calc.dewpoint_from_relative_humidity(temperature, relative_humidity) を使うと,temperature (温度) と relative_humidity (相対湿度) は GFS から得られます。

水平傾度の計算

水平傾度を計算するには,metpy.calc.gradient(f, deltas=(dx, dy)) を利用します。f には水平傾度を計算したい量 (ここでは相当温位 (のアレイ)) を,deltas で指定するのは f がサンプルされている間隔 (のアレイ) です。

f として,GFSが出力する緯度経度のグリッドのデータを用いるので,deltas で指定する間隔は mpcalc.lat_lon_grid_deltas(longitude, lattitude) で得ることができます。

実装


以上述べたことをPythonで実装すると,次のようになります。以下のコードをコピーしてplot-front-analysis.pyという名前で,適当なディレクトリ(ここでは,~/Desktop/labcode/python/front-analysis とします)に保存します。

import datetime

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


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

# 最大,最小の緯度と経度
lon_min = 60
lon_max = 230
lat_min = -10
lat_max = 80


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

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

    # 925hPaの温度
    rhum_925, _, _ = gribin.select(name='Relative humidity', level=925)[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)

    # 地上気圧
    mslp, _, _ = gribin.select(name='Pressure reduced to MSL')[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)

    # 850hPaの温度
    tmpr_850, _, _ = gribin.select(name='Temperature', level=850)[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)

    # 300hPaの風
    uwnd_300, _, _ = gribin.select(name='U component of wind', level=300)[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)
    vwnd_300, _, _ = gribin.select(name='V component of wind', level=300)[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)

    # 500hPaの高度
    hght_500, _, _ = gribin.select(name='Geopotential Height', level=500)[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)

    
# ----- ガウシアンフィルタで平滑化する
tmpr_925 = ndimage.gaussian_filter(tmpr_925, sigma=1.0)
rhum_925 = ndimage.gaussian_filter(rhum_925, sigma=1.0)
mslp = ndimage.gaussian_filter(mslp, sigma=1.0)
tmpr_850 = ndimage.gaussian_filter(tmpr_850, sigma=1.0)

uwnd_300 = ndimage.gaussian_filter(uwnd_300, sigma=1.0)
vwnd_300 = ndimage.gaussian_filter(vwnd_300, sigma=1.0)
hght_500 = ndimage.gaussian_filter(hght_500, sigma=1.0)


# ----- 単位をつける
tmpr_925 = tmpr_925 * units('K')
lat = lat * units('degrees_north')
lon = lon * units('degrees_east')
rhum_925 = rhum_925 * 0.01 * units('')
mslp = mslp * units('Pa')
tmpr_850 = tmpr_850 * units('K')

uwnd_300 = uwnd_300 * units('m/s')
vwnd_300 = vwnd_300 * units('m/s')
hght_500 = hght_500 * units('m')


# ----- MetPyで相当温位に基づくTFPを計算する
# まず,相当温位を計算する
dewt_925 = mpcalc.dewpoint_from_relative_humidity(tmpr_925, rhum_925)
ept_925 = mpcalc.equivalent_potential_temperature(925*units('hPa'), tmpr_925, dewt_925)

# 次に,相当温位の水平傾度の大きさを計算する
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
grad_ept_925 = np.array(mpcalc.gradient(ept_925, deltas=(dy, dx)))
mgntd_grad_ept_925 = np.sqrt(grad_ept_925[0]**2 + grad_ept_925[1]**2)

# 相当温位の水平傾度の大きさの水平傾度を計算する
grad_mgntd_grad_ept_925 = np.array(mpcalc.gradient(mgntd_grad_ept_925, deltas=(dy, dx)))

# tfpを計算する
tfp_925 = np.zeros_like(ept_925)
for i in range(ept_925.shape[0]):
    for j in range(ept_925.shape[1]):
        tfp_925[i, j] = -np.dot(grad_mgntd_grad_ept_925[:, i, j], grad_ept_925[:, i, j]/mgntd_grad_ept_925[i, j])


#------------------------------------------------------------------------------
# プロット
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")}-front-analysis.png'

fig.text(0.10, 0.97, f"NOAA GFS (resolution 0.5 deg)", fontsize=12)
fig.text(0.10, 0.95, f"initial: {dtm_ini_gfs.strftime('%Y/%m/%d %H')}UTC", fontsize=12)
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=12)

# 地図用の変換
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.87, 0.50, 0.02, 0.40])
ax02 = fig.add_axes([0.10, 0.05, 0.80, 0.40], projection=crs_map)
cb02 = fig.add_axes([0.87, 0.05, 0.02, 0.40])


# ----- プロット その1 500hPa 気温
# 海岸線を描画する
ax01.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5, color="blue")

# 地図を描く緯度経度の範囲の設定
ax01.set_extent([107.0, 173.0, 5.0, 75.0], crs_data)

# 図の上にタイトルをつける
fig.text(0.15, 0.905, "T850 & EPT925TFP & PSEA", fontsize=12)

# TFP
cmap = mpl.colors.ListedColormap(['greenyellow', 'yellow', 'gold', 'orange', 'red'])
cmap.set_over('magenta')
cmap.set_under('white')

clevs_tfp = [0.5, 1.0, 1.5, 2.0, 4.0, 8.0]
norm = mpl.colors.BoundaryNorm(clevs_tfp, cmap.N)
cf01 = ax01.contourf(lon, lat, tfp_925*1.0e10, 
                   clevs_tfp, norm=norm, cmap=cmap, extend='both',
                   transform=crs_data)

mpl.colorbar.ColorbarBase(cb01, cf01, ticks=clevs_tfp,
                          orientation='vertical', 
                          label='EPT925TFP (K/100km$^2$)',
                          ).minorticks_off()

# 地上気圧
clevs_mslp = np.arange(800, 1200, 4)
ax01.contour(lon, lat, mslp.to(units('hPa')), clevs_mslp, 
             colors='black', linestyles='solid', linewidths=[1.25, 0.75, 0.75, 0.75, 0.75],
             transform=crs_data)

# 850hPa気温
clevs_tmpr_850 = np.arange(-273, 60, 3)
ax01.contour(lon, lat, tmpr_850.to(units('degC')), clevs_tmpr_850, 
             colors='red', linestyles='solid', linewidths=[1.25, 0.75, 0.75, 0.75, 0.75], 
             transform=crs_data)

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


# ----- プロット その2 300hPaのIsotac分布図,500hPaの高度場および等圧線の重ね合わせ図
# 海岸線を描画する
ax02.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5, color="blue")

# 地図を描く緯度経度の範囲の設定
ax02.set_extent([107.0, 173.0, 5.0, 75.0], crs_data)

# 図の上にタイトルをつける
fig.text(0.15, 0.455, "ISOTAC300 & Z500 & PSEA", fontsize=12)

# 地上気圧
clevs_mslp = np.arange(800, 1200, 4)
ax02.contour(lon, lat, mslp.to(units('hPa')), clevs_mslp, 
             colors='black', linewidths=[1.25, 0.75, 0.75, 0.75, 0.75],
             transform=crs_data)

# 300hPa風
wndsp_300 = mpcalc.wind_speed(uwnd_300, vwnd_300).to(units.knots) # wind speed

cmap = mpl.colors.ListedColormap(['darkgrey', 'dodgerblue', 'greenyellow', 'yellow', 'orange', 'magenta'])
cmap.set_over('white')
cmap.set_under('white')

clevs_wndsp_300 = [40.0, 60.0, 80.0, 100.0, 120.0, 140.0, 160.0]
norm = mpl.colors.BoundaryNorm(clevs_wndsp_300, cmap.N)

cf02 = ax02.contourf(lon, lat, wndsp_300, clevs_wndsp_300, norm=norm,
                     cmap=cmap, extend='both',
                     transform=crs_data)

mpl.colorbar.ColorbarBase(cb02, cf02, 
                          orientation='vertical', 
                          label=u'wind speed (knots)'
                          ).minorticks_off()

# 500hPa高度
clevs_hght_500 = np.arange(60.0, 8000.0, 60.0)
cs02 = ax02.contour(lon, lat, hght_500, clevs_hght_500, 
                    colors='purple', linewidths=[1.25, 0.75, 0.75, 0.75],
                    transform=crs_data)

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

#------------------------------------------------------------------------------
# 図を保存する
plt.savefig(out_fig_path, transparent=False)

実行する前に

NOAAが運用する気象予報モデルGSMの出力データをダウンロードする必要があります。ダウンロードの仕方は,以前の記事を参照してください。

今回は,例として2023年4月6日00UTC (日本時間9時) の予報時間0時間 (つまり,初期時刻で解析値) を使って解析・プロットしています。

ファイル名は,16行目にあるように,gfs_4_20230406_0000_000.grb2 で,プログラムと同じディレクトリ (~/Desktop/labcode/python/front-analysis) に置いてあります。

プログラムを実行する

ターミナルを開き,

$ cd ~/Desktop/labcode/python/front-analysis

と入力し,ディレクトリを移動します。あとは以下のコマンドを実行するだけです。( $マークは無視してください)

$ python plot-front-analysis.py

実行結果

~/Desktop/LabCode/python/front-analysis にplot-gfs-ini2023040600-fcst2023040600-front-analysis.pngというファイルができて,以下のようになっていれば成功です!

上の図には,850hPaの等温線 (赤線), 925hPa面の相当温位を利用したTFPの分布 (カラープロット),地上の等圧線 (黒線) がプロットしてあります。

下の図には,300hPa面の等風速分布 (カラープロット),500hPaの等高度線 (紫線),地上の等圧線 (黒線) がプロットしてあります。

日本海から九州を通って東シナ海,華南に伸びるTFPの大きい帯状の部分と,日本のはるか東のTFPの大きい帯状の領域に目が行くと思います。

ちなみに,気象庁から発表された2023年4月6日00UTC (日本時間9時) の天気図は以下のようになっていました。日本海に中心をもつ低気圧から前線が解析され,日本のはるか東にも前線が解析されています。TFPや強風帯と対応が良いことがわかります (詳しい人は,閉塞点と強風軸の対応に目が行くかもしれません)。

2023年4月6日00UTC (日本時間9時) の天気図 (気象庁ホームページ「過去の天気図」より)

コードの解説


上に書いたソースコードの解説をしていきます。

import datetime

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

まず,必要なモジュールをインポートします。個々でエラーが出る場合は,モジュールをインストールしてください。

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

# 最大,最小の緯度と経度
lon_min = 60
lon_max = 230
lat_min = -10
lat_max = 80


つぎに,読み込むデータのファイルパス (現在のディレクトリにあるときはそのファイル名) を指定します。データを読み込む緯度と経度の範囲も個々で指定しています。

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

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

    # 925hPaの温度
    rhum_925, _, _ = gribin.select(name='Relative humidity', level=925)[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)

    # 地上気圧
    mslp, _, _ = gribin.select(name='Pressure reduced to MSL')[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)

    # 850hPaの温度
    tmpr_850, _, _ = gribin.select(name='Temperature', level=850)[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)

    # 300hPaの風
    uwnd_300, _, _ = gribin.select(name='U component of wind', level=300)[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)
    vwnd_300, _, _ = gribin.select(name='V component of wind', level=300)[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)

    # 500hPaの高度
    hght_500, _, _ = gribin.select(name='Geopotential Height', level=500)[0].data(lat1=lat_min, lat2=lat_max, lon1=lon_min, lon2=lon_max)

    
# ----- ガウシアンフィルタで平滑化する
tmpr_925 = ndimage.gaussian_filter(tmpr_925, sigma=1.0)
rhum_925 = ndimage.gaussian_filter(rhum_925, sigma=1.0)
mslp = ndimage.gaussian_filter(mslp, sigma=1.0)
tmpr_850 = ndimage.gaussian_filter(tmpr_850, sigma=1.0)

uwnd_300 = ndimage.gaussian_filter(uwnd_300, sigma=1.0)
vwnd_300 = ndimage.gaussian_filter(vwnd_300, sigma=1.0)
hght_500 = ndimage.gaussian_filter(hght_500, sigma=1.0)

# ----- 単位をつける
tmpr_925 = tmpr_925 * units('K')
lat = lat * units('degrees_north')
lon = lon * units('degrees_east')
rhum_925 = rhum_925 * 0.01 * units('')
mslp = mslp * units('Pa')
tmpr_850 = tmpr_850 * units('K')

uwnd_300 = uwnd_300 * units('m/s')
vwnd_300 = vwnd_300 * units('m/s')
hght_500 = hght_500 * units('m')


この部分では,データを読み込んで,ガウシアンフィルタをかけて平滑化し,単位を付けている部分です。以前の記事と同様です。

# ----- MetPyで相当温位に基づくTFPを計算する
# まず,相当温位を計算する
dewt_925 = mpcalc.dewpoint_from_relative_humidity(tmpr_925, rhum_925)
ept_925 = mpcalc.equivalent_potential_temperature(925*units('hPa'), tmpr_925, dewt_925)

# 次に,相当温位の水平傾度の大きさを計算する
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
grad_ept_925 = np.array(mpcalc.gradient(ept_925, deltas=(dy, dx)))
mgntd_grad_ept_925 = np.sqrt(grad_ept_925[0]**2 + grad_ept_925[1]**2)

# 相当温位の水平傾度の大きさの水平傾度を計算する
grad_mgntd_grad_ept_925 = np.array(mpcalc.gradient(mgntd_grad_ept_925, deltas=(dy, dx)))

# tfpを計算する
tfp_925 = np.zeros_like(ept_925)
for i in range(ept_925.shape[0]):
    for j in range(ept_925.shape[1]):
        tfp_925[i, j] = -np.dot(grad_mgntd_grad_ept_925[:, i, j], grad_ept_925[:, i, j]/mgntd_grad_ept_925[i, j])


この部分は,MetPyの関数を使ってTFPを計算しているところです。

まず露点温度 dewt_925,それから相当温位 ($\theta_\mathrm{e}$) ept_925 を計算します。

次に,相当温位の水平傾度 ($\nabla_\mathrm{h}\theta_\mathrm{e}$) grad_ept_925 を計算し,その大きさ ($|\nabla_\mathrm{h}\theta_\mathrm{e}|$) mgntd_grad_ept_925 も計算します。

そして,その水平傾度 ($\nabla_\mathrm{h}|\nabla_\mathrm{h}\theta_\mathrm{e}|$) を計算し,最後に内積を np.dot(a, b) とってTFP tfp_925 を計算しています。

#------------------------------------------------------------------------------
# プロット
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")}-front-analysis.png'

fig.text(0.10, 0.97, f"NOAA GFS (resolution 0.5 deg)", fontsize=12)
fig.text(0.10, 0.95, f"initial: {dtm_ini_gfs.strftime('%Y/%m/%d %H')}UTC", fontsize=12)
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=12)

# 地図用の変換
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.87, 0.50, 0.02, 0.40])
ax02 = fig.add_axes([0.10, 0.05, 0.80, 0.40], projection=crs_map)
cb02 = fig.add_axes([0.87, 0.05, 0.02, 0.40])

# ----- プロット その1 500hPa 気温
# 海岸線を描画する
ax01.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5, color="blue")

# 地図を描く緯度経度の範囲の設定
ax01.set_extent([107.0, 173.0, 5.0, 75.0], crs_data)

# 図の上にタイトルをつける
fig.text(0.15, 0.905, "T850 & EPT925TFP & PSEA", fontsize=12)

# TFP
cmap = mpl.colors.ListedColormap(['greenyellow', 'yellow', 'gold', 'orange', 'red'])
cmap.set_over('magenta')
cmap.set_under('white')

clevs_tfp = [0.5, 1.0, 1.5, 2.0, 4.0, 8.0]
norm = mpl.colors.BoundaryNorm(clevs_tfp, cmap.N)
cf01 = ax01.contourf(lon, lat, tfp_925*1.0e10, 
                   clevs_tfp, norm=norm, cmap=cmap, extend='both',
                   transform=crs_data)

mpl.colorbar.ColorbarBase(cb01, cf01, ticks=clevs_tfp,
                          orientation='vertical', 
                          label='EPT925TFP (K/100km$^2$)',
                          ).minorticks_off()

# 地上気圧
clevs_mslp = np.arange(800, 1200, 4)
ax01.contour(lon, lat, mslp.to(units('hPa')), clevs_mslp, 
             colors='black', linestyles='solid', linewidths=[1.25, 0.75, 0.75, 0.75, 0.75],
             transform=crs_data)

# 850hPa気温
clevs_tmpr_850 = np.arange(-273, 60, 3)
ax01.contour(lon, lat, tmpr_850.to(units('degC')), clevs_tmpr_850, 
             colors='red', linestyles='solid', linewidths=[1.25, 0.75, 0.75, 0.75, 0.75], 
             transform=crs_data)

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

# ----- プロット その2 300hPaのIsotac分布図,500hPaの高度場および等圧線の重ね合わせ図
# 海岸線を描画する
ax02.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5, color="blue")

# 地図を描く緯度経度の範囲の設定
ax02.set_extent([107.0, 173.0, 5.0, 75.0], crs_data)

# 図の上にタイトルをつける
fig.text(0.15, 0.455, "ISOTAC300 & Z500 & PSEA", fontsize=12)

# 地上気圧
clevs_mslp = np.arange(800, 1200, 4)
ax02.contour(lon, lat, mslp.to(units('hPa')), clevs_mslp, 
             colors='black', linewidths=[1.25, 0.75, 0.75, 0.75, 0.75],
             transform=crs_data)

# 300hPa風
wndsp_300 = mpcalc.wind_speed(uwnd_300, vwnd_300).to(units.knots) # wind speed

cmap = mpl.colors.ListedColormap(['darkgrey', 'dodgerblue', 'greenyellow', 'yellow', 'orange', 'magenta'])
cmap.set_over('white')
cmap.set_under('white')

clevs_wndsp_300 = [40.0, 60.0, 80.0, 100.0, 120.0, 140.0, 160.0]
norm = mpl.colors.BoundaryNorm(clevs_wndsp_300, cmap.N)

cf02 = ax02.contourf(lon, lat, wndsp_300, clevs_wndsp_300, norm=norm,
                     cmap=cmap, extend='both',
                     transform=crs_data)

mpl.colorbar.ColorbarBase(cb02, cf02, 
                          orientation='vertical', 
                          label=u'wind speed (knots)'
                          ).minorticks_off()

# 500hPa高度
clevs_hght_500 = np.arange(60.0, 8000.0, 60.0)
cs02 = ax02.contour(lon, lat, hght_500, clevs_hght_500, 
                    colors='purple', linewidths=[1.25, 0.75, 0.75, 0.75],
                    transform=crs_data)

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

#------------------------------------------------------------------------------
# 図を保存する
plt.savefig(out_fig_path, transparent=False)


残りは,プロットをしているところです。以前の記事と同様です。コメントしておくべきところとしては,非等間隔の範囲に対して,指定した色をつけてカラープロットしている部分です。

例えば,上図は [0.5, 1.0, 1.5, 2.0, 4.0, 8.0] という等間隔でない範囲の境界となる数値が与えられて,その間 (5つに区切られる) を ['greenyellow', 'yellow', 'gold', 'orange', 'red'] で塗り,上限 (上の場合は8.0) を超える部分には ‘magenta’ を,下限を下回る場合は ‘white’ で塗るように指定しています。

これには,matplotlib.colors.ListedColormap()set_over()set_under()matplotlib.colors.BoundaryNorm() を組み合わせるとうまくできます。

最後に


今回は,気象庁が前線の客観解析に使っている (らしい) TFPやと風速分布図のプロットをしてみました。

MetPyには実装するには面倒くさい気象要素の計算や微分の計算のための関数が数多く用意されています。他にも,重要な気象学的なパラメタを簡単に計算しプロットすることができます。

今回実装したプログラムから出力される資料を基に自分で前線を解析してみて,気象庁のプロの予報官が解析した前線と比較してみると面白いかもしれません。

参考


この記事を執筆するに当たり,以下の資料を参考にしました。

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です