JRA-55のデータを使ったエマグラム

2021年11月9日 21:19

再解析データ(JRA-55)からエマグラムを作図するコードのご紹介です。

JRA-55のエマグラム作成の動機

昭和39年7月山陰北陸豪雨の事例解析を行なっています。JRA-55で作成した断面図を解析していると、ポイントのエマグラムを確認したくなりました。前線の遷移層の高度、暖域内の大気の安定度など断面図からも読み取れますが、エマグラムの方がより正確に読み取れ、これを参考にすることで断面図解析をより正確に解析できます。

例えば、図１の鉛直断面図右下付近にある温暖前線の遷移層を、等温位線の間隔や風向の鉛直変化（暖気移流を示す時計回りの変化）から解析しました。この遷移層の高度解析をより確からしくするために、黄色の線のポイントのエマグラムを確認しましょう。

そこで、JRA-55のエマグラムのGPVのポイントのエマグラムを描画させるコードを作成し、次のエマグラムを作図しました。この図の仕様は、高層ゾンデ観測のエマグラムの作図を説明した記事「エマグラム」と同じです。

図２にあるJRA-55の前述のポイントのエマグラムを確認すると、前線の遷移層は700から900hPaでした。断面図のみから解析した遷移層は厚みが足りませんでした。

エマグラムを作図するためのデータ変換

紹介するコードは、高層ゾンデ観測データによる「エマグラム」の記事にあるコードを流用しています。ここでの説明は、JRA-55の１格子のデータを前述のコードと同じデータ形式への変換のみとします。これ以外の説明は、前述の記事をご確認ください。

下のコードは高度のNetCDF形式のデータを読み込み、指定された位置に格子が存在するか確認し、指定した格子の１次元配列を得ています。

## データ保存先
DataFd="./data/Jra55/"
#
## NetCDF ファイル名                                                                    
yyyymm='{:04d}{:02d}'.format(year,month)
HgtFn ='{}HGT_{}.nc'.format(DataFd,yyyymm)
#
## 読み込む高度を指定するための変数  
cut_hgt=10                 #    10:100hPa 17:300hPa以下の高度データ  cut_hgt >= 10
cut_hgt_Dep= cut_hgt - 10  # DEP 0:100hPa以下の高度でデータがあるため
#
## xarrayデータセットとして、NetCDF形式のJRA-55を読み込む
ds = xr.open_dataset(HgtFn)
dataHgt = ds.metpy.parse_cf('HGT').squeeze()
#
## 緯度経度の配列を取得
ary_lon=np.array(dataHgt["lon"])
ary_lat=np.array(dataHgt["lat"])
#
## 指定した緯度経度が、JRA-55データの格子と一致するか確認
i_lat = -1
i_lon = -1
if np.where(ary_lon == f_lon)[0].size == 0:
   print("error!")
   print(ary_lon)
else:
   i_lon = np.where(ary_lon == f_lon)[0][0]
if np.where(ary_lon == f_lat)[0].size == 0:
   print("error!")
   print(ary_lat)
else:
   i_lat = np.where(ary_lat == f_lat)[0][0]
if i_lat == -1 or i_lon == -1:
   print("exit!"）
＃
dataHgt = dataHgt.isel(time=time_targ, lat=i_lat, lon=i_lon)
dataHgt = dataHgt[cut_hgt:]      # 100hPa以下の高度データとする

エマグラムの入力は、単位つきの配列データとなっていて、
units.Quantity(numpy配列、単位)
とすることで、変換できます。ただし、JRA-55のデータ配列は上層から下層に向かう順番となっていて、求めるデータとは逆順になっているため、np.flip()を利用して、順番を反転させています。

# 気温のGPVポイントデータの読み込み
TmpFn ='{}TMP_{}.nc'.format(DataFd,yyyymm)
ds = xr.open_dataset(TmpFn)
dataTmp = ds.metpy.parse_cf('TMP').squeeze()
dataTmp = dataTmp.isel(time=time_targ, lat=i_lat, lon=i_lon)
dataTmp = dataTmp[cut_hgt:]    # 100hPa以下の高度データとする
#
# 露点差や風の読み込みコードは省略
#
# united arraysへ変換
p = units.Quantity(np.flip(np.array(dataTmp["level"])), units.hPa)
T = units.Quantity(np.flip(np.array(dataTmp)), units.K)
Td = T - units.Quantity(np.flip(np.array(dataDep)), units.K)
u = units.Quantity(np.flip(np.array(dataUgrd)) / 0.514444, units('knots'))
v = units.Quantity(np.flip(np.array(dataVgrd)) / 0.514444, units('knots'))
h = units.Quantity(np.flip(np.array(dataHgt)), units.meter)

サンプルコード

このサンプルコードは、./data/Jra55/にデータを配置すると実行できます。なお、MetPy1.1を利用しています。

必要になれば、GRIB2形式のGSMのエマグラム作成コードも作って、いつかここにアップしたいと思います。

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