![見出し画像](https://assets.st-note.com/production/uploads/images/65304540/rectangle_large_type_2_2f6be4fc4573ce27b57ef179c836f430.png?width=1200)
JRA-55のデータを使ったエマグラム
再解析データ(JRA-55)からエマグラムを作図するコードのご紹介です。
JRA-55のエマグラム作成の動機
昭和39年7月山陰北陸豪雨の事例解析を行なっています。JRA-55で作成した断面図を解析していると、ポイントのエマグラムを確認したくなりました。前線の遷移層の高度、暖域内の大気の安定度など断面図からも読み取れますが、エマグラムの方がより正確に読み取れ、これを参考にすることで断面図解析をより正確に解析できます。
例えば、図1の鉛直断面図右下付近にある温暖前線の遷移層を、等温位線の間隔や風向の鉛直変化(暖気移流を示す時計回りの変化)から解析しました。この遷移層の高度解析をより確からしくするために、黄色の線のポイントのエマグラムを確認しましょう。
![画像1](https://assets.st-note.com/production/uploads/images/65300695/picture_pc_7609811d3c5f3d4e4b69bbe6a9ecc703.png?width=1200)
そこで、JRA-55のエマグラムのGPVのポイントのエマグラムを描画させるコードを作成し、次のエマグラムを作図しました。この図の仕様は、高層ゾンデ観測のエマグラムの作図を説明した記事「エマグラム」と同じです。
図2にあるJRA-55の前述のポイントのエマグラムを確認すると、前線の遷移層は700から900hPaでした。断面図のみから解析した遷移層は厚みが足りませんでした。
![画像2](https://assets.st-note.com/production/uploads/images/65301105/picture_pc_9e5f1fa3949e2332088ad0b7fde29960.png)
エマグラムを作図するためのデータ変換
紹介するコードは、高層ゾンデ観測データによる「エマグラム」の記事にあるコードを流用しています。ここでの説明は、JRA-55の1格子のデータを前述のコードと同じデータ形式への変換のみとします。これ以外の説明は、前述の記事をご確認ください。
下のコードは高度のNetCDF形式のデータを読み込み、指定された位置に格子が存在するか確認し、指定した格子の1次元配列を得ています。
## データ保存先
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のエマグラム作成コードも作って、いつかここにアップしたいと思います。
この記事が気に入ったらサポートをしてみませんか?