pygribでオープンできないGRIB2形式の気象データをnetCDFに変換して可視化する
2024.5 一部修正
この記事で使用しているデータ(1kmメッシュ全国合成レーダーGPV)について、気象庁からの配信が終了しました。
詳細は気象庁資料をご覧ください。
～ 10分毎１kmメッシュ全国合成レーダーGPVの配信終了 ～
（配信資料に関する技術情報第 568 号関連）
これは技術の進歩によって、データ配信間隔が10分から5分に変更となったり、メッシュ間隔が1kmから250mになったりしたためです。
気象庁によるとしばらくの間、250m間隔メッシュを低解像度化した1kmメッシュが提供されるそうです。
これに伴い、元の投稿での以下が変更になっています。
・ファイル名
ファイル名途中の _Prr10lv_が、_Prr05lv_に変更
・パラメータ名
Netcdfファイルでの指定パラメータがvar0_1_201_surfaceからvar0_1_203_surfaceに変更
はじめに
面的立体的な広がりを持った気象データのフォーマットは「GRIB2」形式が業界標準です。
GRIB2形式については@e_toyodaさんによる詳しい紹介があります。
WMO（世界気象機関）の格子データ形式GRIB2について
GRIB2形式のファイルからデータを読み出すための様々なツールがあります。
pythonならばpygribというモジュールを利用できます。
ところが、全国合成レーダーや推計気象分布などは全国1km間隔でのデータであることからデータ量が大きく、GSMやMSMとは違った、日本気象庁が独自に指定している圧縮形式が用いられています。
これがpygribでサポートされていないようで、オープンできないデータが存在しています。
＊全国合成レーダは以下に解説があります。
１kmメッシュ全国合成レーダーGPV
このデータを用いたDeep LearningプログラムをKerasで作りたいので、pythonで扱えるようにしたいところです。
そこで気象庁フォーマットをサポートしているwgrib2というツールを使用して、pythonで扱えるnetCDF形式に変換しました。
pygribでオープンできない全国合成レーダ
全国合成レーダは京都大学生存圏研究所からダウンロードしました。
このデータファイルをpygribでオープンすると以下のエラーが表示されます。
>>> import pygrib
>>> grb=pygrib.open("Z__C_RJTD_20240331230000_RDR_JMAGPV_Ggis1km_Prr05lv_ANAL_grib2.bin")
ECCODES ERROR : Unable to find template productDefinition from grib2/template.4.50008.def
ECCODES ERROR : Unable to find template productDefinition from grib2/template.4.50008.def
これは気象庁が使用している定義がpygribの中に存在していないためと思われます。
wgrib2を用いた変換
wgrib2
wgrib2については、@ysomeiさんが下記に纏めておられます。
アメリカの「大気海洋庁（NOAA）」が提供しているツールです。
今更ですが、気象データを解析してみた
色々なオプションがあって様々なことができるのですが、例えば以下のようにGRIBファイルの内容をCSVに出力するとことができます。
できたCSVファイルをのぞくと、データを指定するラベルにvar0_1_201とsurfaceが入っています。
できたCSVファイルをのぞくと、データを指定するラベルにvar0_1_203とsurfaceが入っています。
最後の３要素が、経度、緯度、データ（「データ代表値」）です。
> wgrib2 -V GRIB_formatted_file -csv tmp.csv
> cat tmp.csv
略
"2024-03-31 23:00:00","2024-03-31 23:00:00","var0_1_203","surface",141.419,46.0125,0
"2024-03-31 23:00:00","2024-03-31 23:00:00","var0_1_203","surface",141.431,46.0125,0
"2024-03-31 23:00:00","2024-03-31 23:00:00","var0_1_203","surface",141.444,46.0125,0
"2024-03-31 23:00:00","2024-03-31 23:00:00","var0_1_203","surface",141.456,46.0125,0
"2024-03-31 23:00:00","2024-03-31 23:00:00","var0_1_203","surface",141.469,46.0125,0
後略
合成レーダのデータ形式（データ代表値とレベル値）
合成レーダの「データ代表値」は降雨強度であるmm/hが単位です。
元のGRIB2形式のファイルには、直接降雨強度の数値がセットされてはおらず、「レベル値」という0から251の値がセットされています。
例えば「レベル値」が「37」というのは、降雨強度が「7.0mm/h以上、7.5mm/h未満」ということを意味しています。
その「データ代表値」は「7.75mm/h」であるということです。
つまり、「レベル値」から「データ代表値」への変換テーブルが存在していることになりますが、wgrib2はそこもやってくれるようです（wgrib2を解読した訳ではありませんが）。
NetCDF形式への変換
wgrib2を用いて、GRIB2からnetCDF形式へ変換するには下記のオプションを使います。
wgrib2 GRIB_formatted_file -netcdf GRIB_formatted_file.nc
これでnetCDF形式ファイルが出来上がります。
ただし、圧縮が解除されるためか、猛烈にファイルサイズが大きくなりますので注意が必要です。
103424 Z__C_RJTD_20240331230000_RDR_JMAGPV_Ggis1km_Prr05lv_ANAL_grib2.bin
34454836 Z__C_RJTD_20240331230000_RDR_JMAGPV_Ggis1km_Prr05lv_ANAL_grib2.bin.nc
pythonでnetCDF4を用いたデータの読み込み
このファイルはpythonモジュールnetCDF4を用いると下記のように取り扱うことができます。
from netCDF4 import Dataset
(中略）
nc = Dataset( fn , 'r' )
xx = nc.variables['longitude'][:] #
yy = nc.variables['latitude' ][:] #
data = nc.variables['var0_1_203_surface'][:]
lons, lats = np.meshgrid( xx , yy )
_save_file = savefilenamelist[ptr_sfile]
drawmap( data[0,:,:] , levels_rai , precip_colormap , _save_file ,_flag_show , "RAIN" , norm)
ファイルに含まれているデータを特定するために、var0_1_203_surfaceを指定しています。
wgrib2でGRIBファイルの中身を参照した際には、var0_1_203とsurfaceに分かれて見えていましたが、NetCDFファイルに変換した際にvar0_1_203_surfaceというラベルに変換されるようです。
下記のように、netcdfファイルを情報出力させると、variablesの最後にこのラベルが出てきています。
>>> from netCDF4 import Dataset
>>> nc=Dataset("www.nc")
>>> nc
<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF3_CLASSIC data model, file format NETCDF3):
Conventions: COARDS
History: created by wgrib2
GRIB2_grid_template: 0
dimensions(sizes): latitude(3360), longitude(2560), time(1)
variables(dimensions): float64 latitude(latitude), float64 longitude(longitude), float64 time(time), float32 var0_1_203_surface(time,latitude,longitude)
groups:
可視化
可視化するところのpythonソースです。
matplotlibを用いています。
手前味噌ではありますが、可視化のあれこれは、気象データをもとに「天気図っぽい前線」を機械学習で描いてみる（2) 入力となる気象データ可視化の話に記載しております。
可視化するカラーマップに白色のレンジがあるので、下記の可視化では地図そのものに少し色をつけています。
色のついた地図とコンターマップの重ね合わせの順番をzorderで指定しています。
def drawmap( values , levels , cmap , _save_filename , flag_show , name ,norm) :
fig,ax = plt.subplots(figsize=(5,5))
plt.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.02)
m = Basemap(projection='stere', llcrnrlat=25, urcrnrlat=47, llcrnrlon=125, urcrnrlon=150, lat_0=60, lon_0=140, resolution='i' )
m.drawparallels(np.arange(-80.,81.,10.))
m.drawmeridians(np.arange(-180.,181.,10.))
m.drawcoastlines()
m.drawmapboundary(fill_color='#eeeeee' , zorder=-1)
m.fillcontinents(color='#ceceee' , zorder=-1)
x , y = m(lons, lats)
m.contourf( x , y , values , levels , cmap=cmap , zorder=0 , norm=norm )
fig.savefig(_save_filename)
if flag_show == 'yes' :
plt.show()
plt.close()
下図のような画像が作成できます(2019年9月8日15時UTCのデータ）。
##まとめ
GRIB2形式だけどpygribで扱えないファイルをnetCDF形式に変換してpythonから扱えるようにした方法を記載しました。