GRIB2形式データを処理するプログラムwgrib2（入門編）

はじめに

数値予報モデルGPVなど、多くの格子点形式の気象データでは、世界気象機関WMOが定めるGRIB2というフォーマットが多く使われています。GRIB2の構造の解説は他のサイトや公的な資料にお任せするとして、ここでは、実用的に、GRIB2形式データを処理するプログラムwgrib2について、普段私が良く使うオプションやテクニック等を紹介したいと思います。

wgrib2は、米国海洋大気庁（NOAA）が提供しているGRIB2処理プログラムです。気象業界ではとても広く利用されています。今回は、入門編として、wgrib2のインストール、GRIB2形式データの基本的な利用方法、デコード（プログラムで利用しやすい形式に変換したり、CSV形式等に出力したりする）方法について解説します。

検証に利用したwgrib2のバージョンは2.0.8、OS環境はLinux x86_64 4.18.11-200 (gcc 8.1.1)です。

インストール

コンパイルはIntel Compiler等有償コンパイラでも可能ですが、ここではgccを利用します。
wgrib2の提供サイトからパッケージを取得して、適当なディレクトリで展開します。

$ tar -zxvf wgrib2.tgz.v2.0.8
$ cd grib2


bash環境ではexport、csh環境ではsetenvで環境変数を設定します。

# sh/bash
$ export CC=gcc
$ export FC=gfortran


# csh/tcsh
$ setenv CC gcc
$ setenv FC gfortran


あとはmakeするだけです。

$ make

wgrib2ディレクトリにwgrib2コマンドが作成されていれば成功です。パスを通すなどしておきましょう。

# mv grib2 /usr/local
$ vi ~/.bashrc
---
export WGRIB2="/usr/local/grib2"
export PATH="${WGRIB2}/wgrib2:${PATH}"
---


GRIB2形式データの中身を確認する

まずは、wgrib2の引数にGRIB2形式データを指定して実行してみましょう。データ要素を見ることができます。

$ wgrib2 grib2.bin

各要素の説明は気象庁技術資料を参照して確認します。
例えば、MSM-GPV(地上:Z__C_RJTD_yyyymmddhh0000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin)を引数にして実行すると以下の要素が出力されます。

記号	物理量
PRMSL	海面更生気圧
PRES	地上気圧
UGRD	東西風速
VGRD	南北風速
TMP	気温
RH	相対湿度
LCDC	下層雲量
MCDC	中層雲量
HCDC	上層雲量
TCDC	全雲量
APCP	降水量
DSWRF	日射量

面表現	意味
mean sea level	平均海面
surface	地表面
10 m above ground	地上から高度10m
1.5 m above ground	地上から高度1.5m

時間表現	意味
anl	初期時刻
N hour fcst	N時間予報
N1-N2 hour acc fcst	N1時間からN2時間までの積算値
N1-N2 hour ave fcst	N1時間からN2時間までの平均値

-grid: 格子配列を確認する

-gridオプションで格子配列の情報を出力できます。
例：MSM-GPV

$ wgrib2 -grid Z__C_RJTD_20191226000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin
1.1:0:grid_template=0:winds(N/S):
lat-lon grid:(481 x 505) units 1e-06 input WE:NS output WE:SN res 48
lat 47.600000 to 22.400000 by 0.050000
lon 120.000000 to 150.000000 by 0.062500 #points=242905
1.2:0:grid_template=0:winds(N/S):
lat-lon grid:(481 x 505) units 1e-06 input WE:NS output WE:SN res 48
lat 47.600000 to 22.400000 by 0.050000
lon 120.000000 to 150.000000 by 0.062500 #points=242905
～（略）～


詳細は省きますが、格子点が481×505、緯度方向に0.05度間隔、経度方向に0.0625度で、北西端から東、北から南に向かって格納されています。
格納順は北から南ですが、出力は南から北になるので注意が必要です。

単精度バイナリで出力し、Fortranのdirect accessで読む

気象業界では科学計算の得意なFortranをよく使います。Fortranのdirect access（配列単位で読み出しを指定する）で読んでみましょう。
まずは単精度バイナリファイルへの出力です。

$ wgrib2 Z__C_RJTD_20191226000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin | grep "TMP" | wgrib2 -i -no_header Z__C_RJTD_20191226000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin -no_append -ieee output.bin

オプションなしのwgrib2で出力したGRIB2形式データの情報をgrepコマンドで気温（TMP）だけ抽出してパイプで渡すと、そのデータだけ処理することができます。
その他のオプションは以下の通りです。
-i：入力ファイルの指定
-no_header：direct access用ファイル指定
-ieee：big_endianのバイナリファイルを出力
出力ファイルの容量を確認すると15545920バイト（東西481格子*南北505格子*16時間ステップ*4バイト）となっています。

$ ls -la output.bin
-rw-rw-r-- 1 user user 15545920 12月 27 13:54 output.bin


Fortranで読むには、以下のように変数宣言とopen文、read文を使います。

REAL array(481,505,16)
open(10,file='output.bin',form='unformatted',access='direct',recl=481*505*16*4)
read(10,rec=1) array


このバイナリファイルoutput.binは、GrADSという描画ソフトウェアでも扱うことができます。このような気温の分布図も簡単に描けます。

GrADSのコントロールファイルの例：

$ vi grads.ctl
dset output.bin
options big_endian
undef 9.999E+20
title MSM-GPV temperature
xdef 481 linear 120.0 0.0625
ydef 505 linear 22.4 0.05
zdef 1 levels 1000.0
tdef 16 linear 00Z26DEC2019 1hr
vars 1
tmp 1 0 temperature
endvars


上図のGrADS描画コマンド：

> open grads.ctl
> set mpdset hires
> set t 1
> set gxout grfill
> d tmp


変換したデータが正しいかどうか確認するために、値を標準出力したり、図を描いたりすると良いと思います。

-csv: CSV形式で出力する

$ wgrib2 Z__C_RJTD_20191226000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin | grep "TMP" | grep ":1 hour fcst:" | wgrib2 -i Z__C_RJTD_20191226000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin -csv output.csv


全部出力するとファイルサイズが大きくなるので、例では気温の1時間予報値のみとしています。

$ cat output.csv
"2019-12-26 00:00:00","2019-12-26 01:00:00","TMP","1.5 m above ground",120,22.4,295.64
"2019-12-26 00:00:00","2019-12-26 01:00:00","TMP","1.5 m above ground",120.062,22.4,295.64
"2019-12-26 00:00:00","2019-12-26 01:00:00","TMP","1.5 m above ground",120.125,22.4,295.609
"2019-12-26 00:00:00","2019-12-26 01:00:00","TMP","1.5 m above ground",120.188,22.4,295.609
"2019-12-26 00:00:00","2019-12-26 01:00:00","TMP","1.5 m above ground",120.25,22.4,295.64
"2019-12-26 00:00:00","2019-12-26 01:00:00","TMP","1.5 m above ground",120.312,22.4,295.64
"2019-12-26 00:00:00","2019-12-26 01:00:00","TMP","1.5 m above ground",120.375,22.4,295.625
～（略）～


1列目が初期時刻、2列目が予報時刻、5列目が経度、6列目が緯度、7列目に気温の値が入っています。
経度はよく見ると、小数点第3位までになっており、正確に出力されていません。気温も丸めて出力されているようで、気温の小数点以下の値は有意でないので気にしなくてよいですが、他の要素を出力するときは注意した方がよさそうです。
CSV形式であればpythonはもちろん、Excelや様々なプログラムで処理することができますね。

まとめ

GRIB2形式データを処理するプログラムwgrib2を用いて、GRIB2形式データの中身を確認する方法、単精度バイナリファイルに出力してFortranで読む方法、CSV形式ファイルに出力する方法について解説しました。
ほとんどのGRIB2形式データをwgrib2で処理できますが、wgrib2のバージョンやGRIB2形式データの圧縮方式によって対応していない場合があります。（wgrib2対応バージョンをまとめてみようと考えています。）
pythonでは、GRIB2形式データを直接読むライブラリがあるようですね。多くの方は、そちらの方が興味があるのではないでしょうか。私も勉強して、また紹介できればと思います。