JRA-55から天気図を作成する

これまでは、気象庁GSMのGRIB2形式のデータから天気図を作成してきました。今回は、JRA-55のデータを表示させるために、NetCDF形式のデータ読み込みについて説明します。描画のコードは前回までと変わらないため省いています。

京都大学のサイトでは、気象庁JRA-55のNetCDF形式のデータをユーザー申請するとダウンロードできるようになります。これを利用すると、1955年以降の天気図を描画できるようになります。今回は、この等圧面データの読み込みについてご紹介します。

最後に、以前の記事「ジェット気流と上層発散を把握できる天気図」で作成した同じ天気図のNetCDF版のサンプルコードを添付します。コードの説明は、このサンプルコードを一部切り取っています。実行させるためには、サンプルコードをご利用ください。

NetCDF形式データの読み込み

xarrayを使ってnetCDFデータを読み込みます。xarrayのopen_dataset()を使うと、NetCDF形式のデータをxarray形式のデータセットに直接変換できます。

次のコードでは、まず読み込むデータの時刻や等圧面、データ切り出し領域を指定します。その後、要素(高度、発散など)ごとのNetCDFデータを読み込み、それぞれにデータセットを作成し、指定された時刻・等圧面・領域のみを切り出したデータセットに変換しています。

京都大学からダウンロードしたファイル名は、{TYPE}_{YYYY}{MM}.nc となっています。{TYPE}は、HGT(高度), TMP(気温), DEPR(湿数), UGRD(風速東西成分), VGRD(風速南北成分), STRM(流線関数), VPOT(速度ポテンシャル), VVEL(鉛直速度), RELV(相対渦度)です。{YYYY}は西暦年4桁、{MM}は月2桁です。 

## JRA55 netcdf形式データの読み込み処理
import datetime
from dateutil import tz
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import metpy.calc as mpcalc
import metpy.constants as m_const
from metpy.units import units
import numpy as np
import xarray as xr
import scipy.ndimage as ndimage
import scipy.constants as s_const
import sys
import math
#
## 読み込みデータの指定
# JRA55の読み込む年月日時をUTCで与えます。
i_year =1961
i_month = 7
i_day = 15
i_hourZ = 18
# 気圧面を指定
i_pre=300
# データ切り出し領域
lat_cut=slice(80.0,-20.0)
lon_cut=slice(70.0,190.0)
# データ格納先フォルダー名
DataFd="./data/Jra55/"
#
## データ読み込み、データセット作成
# 年月の6桁の整数(文字列)の取得                                      
yyyymm='{:04d}{:02d}'.format(i_year,i_month)
#
## 月別のNetCDF File名
# 高度
HgtFn ='{}HGT_{}.nc'.format(DataFd,yyyymm)
#
## 高度                                                                    
ds = xr.open_dataset(HgtFn)
dataHgt = ds.metpy.parse_cf('HGT').squeeze()

このデータセットの属性は、データセット名.attrs["units"] などと表記して、アクセスすることができます。下のコードでは、単位を指定しています。

dataHgt.attrs['units'] = 'meter'

データ切り出し

データの切り出しでは、データセット名.isel()  あるいは sel()を使います。前者は次元の何番目データを、後者は次元の座標値を指定して切り出せます。

# time番号                                                                 
time_targ=(i_day - 1) * 4 + i_hourZ // 6
# デー切り出し
dataHgt = dataHgt.isel(time=time_targ)

# 気圧面を指定
i_pre=300
# データ切り出し領域
lat_cut=slice(80.0,-20.0)
lon_cut=slice(70.0,190.0)
# デー切り出し
dataHgt = dataHgt.sel(level=i_pre,
                       lat=lat_cut,
                       lon=lon_cut)

参考のために、データセットの次元の文字列、次元のデータサイズ、次元の配列の値について、それぞれのアクセスの方法を下のコードに示します。データセット名.dims、データセット名.sizes、などを使います。

# データセットへのアクセス方法について
## 時刻や座標などの次元の取得
print(dataHgt.dims)
#
# ('time', 'level', 'lat', 'lon')
#
## 次元ごとに、xarray形式の配列の長さを取得し、配列の最後の値を取り出す。
for dim in dataHgt.dims:
   size = dataHgt.sizes[dim]
   print(dim,'   length:',size,'   last value:',dataHgt[dim][size - 1].data)
#
# time    length: 124    last value: 1961-07-31T18:00:00.000000000
# leve    length: 37     last value: 1000.0
# lat     length: 145    last value: -90.0
# lon     length: 288    last value: 358.75

JRA-55データの時刻を確認する

時刻の取得については、データセット名.time.astype(datetime.datetime)により、int型のUNIX timeが得られます。一方で、pythonのdatetime.datetimeは、Float型であるために、1e-9をかけて変換しています。

## 時刻文字列化                                                            
UTC = tz.gettz("UTC")
dt1 = datetime.datetime.fromtimestamp(dataHgt.time.astype(datetime.datetime) * 1e-9, tz=UTC)
dt_str = (dt1.strftime("%HZ%d%b%Y")).upper()
dt_str2 = dt1.strftime("%Y%m%d%H")
print(dt_str)
print(dt_str2)
## 出力
# 18Z15JUL1961
# 1961071518

地衡風などの物理量の計算例

下のコードの通り、これまでの記事で説明した通り、データセットを作成すると、MetPyを使うと簡単に計算することができます。

### 計算                                                                   
# 地衡風                                                                   
ug,  vg  = mpcalc.geostrophic_wind(dataHgt)                               
# 非地衡風                                                   
uag, vag = mpcalc.ageostrophic_wind(dataHgt, dataUgrd, dataVgrd)
# 風速 m/s => knot                                                         
dataWS = mpcalc.wind_speed(dataUgrd, dataVgrd)

サンプルコード

ジェット気流解析するための天気図を作図する、Jupyter Notebook用のコードを次からダウンロードできます。JRA-55のNetCDFデータが必要です。

note7.ipynb

511598 Bytes

ファイルダウンロードについて

ダウンロード

次回は、大気の安定度を表現する天気図について説明する予定です。