OTENKI Data Science

TechBlog

K.Sakurai

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 fcstN時間予報
N1-N2 hour acc fcstN1時間からN2時間までの積算値
N1-N2 hour ave fcstN1時間から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形式データを直接読むライブラリがあるようですね。多くの方は、そちらの方が興味があるのではないでしょうか。私も勉強して、また紹介できればと思います。
 

この記事を書いた人

K.Sakurai

気象予報士/技術士(応用理学)/防災士 
総合気象数値計算システムSACRA、データ提供システムCOSMOS及びお天気データサイエンスの開発を一から携わる。
WRF-5kmモデル、虹予報、1㎞メッシュ雨雪判別予測データを開発。
趣味はバドミントンと登山。バドミントンのシャトルは気温と湿度で最適な号数を選びたいと常々思っている。登山は、関西地方の低山でトレーニングしながら、年に1回くらいは高山に挑戦。

前の記事へ 天気予報の変化
記事一覧へ
次の記事へ GPVから最も近い格子点を抽出する方法
  1. お天気データサイエンスHOME
  2. 技術情報
  3. 技術ブログ