GRIBデータの扱い② - 気象系技術メモ

前回の続き。
前回まとめたような手続きでバイナリダンプしたGRIBのデータをfortranで読み込んで解析します。

バイナリファイルの構造

こちらやこちらに詳しくまとまっています。

wgribでGRIBをバイナリダンプしたファイルには、

時間--変数--高度--緯度--経度
の順で辞書式に値が入っています。複数変数のデータがある場合、その格納順序はもとのGRIBファイルにおけるそれと同じです。

個人的には、変数ごとにファイルが独立している方が何かとやりやすいのですが（複数の解析で同じデータファイルを再利用することを考えると、ファイルがまとまっていると場合によっては使わない変数まで読み込む必要が生じるため）、その場合ファイルは「経度方向グリッド数×緯度方向グリッド数×鉛直層数×時間ステップ数」の配列であることになります。

プログラムに読み込む

したがって、GRIBから作ったバイナリファイル（1変数）をFortranなどのプログラムに配列形式で読み込むとすると、

program sample
  implicit none
  integer, parameter :: nlon = 288       !経度方向グリッド数
  integer, parameter :: nlat = 145       !緯度方向グリッド数
  integer, parameter :: nlev = 37        !鉛直層数
  integer, parameter :: ntim = 1         !時間ステップ数
  real(4) :: data(nlon, nlat, nlev, ntim)

  open(11, file=ifile, access="direct",status="old", form="unformatted", convert="little_endian", recl=4*nlon*nlat*nlev*ntim)
  read(11, rec=1) data　　　!nlon*nlat*nlev*ntimの配列dataにファイルの値を入れる
  close(11)

みたいな感じになります。

探査の方式には順番探査（access="sequential"）と直接探査（access="direct"）があり、順番探査がデータの初めと終わりを示す記号（record marker）を見てその間にあるデータを読み込む方式であるのに対し、直接探査はあらかじめ定められたデータの単位長さ（record length、変数reclで指定）ごとにデータを読み込む方式です。
つまり、直接探査の場合たとえばレコード長を4バイトとするならば、4バイトのデータを一単位として、4バイトずつ読み込んでいくことになります。
ひとつひとつの数値のデータサイズがすでに決まっているバイナリファイルでは直接探査の方が便利なので、ここではaccess="direct"としています。

上記サンプルコードではrecl=4*nlon*nlat*nlev*ntimとなっていました。ここでは読み込むファイルの数値がひとつあたり4バイトで、それが「経度方向グリッド数×緯度方向グリッド数×鉛直層数×時間ステップ数」の数だけ存在するので、データのサイズは

4バイト×経度方向グリッド数×緯度方向グリッド数×鉛直層数×時間ステップ数
ということになります。したがって、recl=4*nlon*nlat*nlev*ntimというのは「レコード長＝ファイル全体のサイズ」ということで、ファイルをまるごと一気に読み込んでしまうということです。

そのあとのread文のrec変数は「レコード番号」で、レコード長ごとに頭からデータを数えていったときのデータの番号です。今回はファイルをまるごと一つのデータ（配列）として読み込んでいるのでレコード番号は1しかありませんが、例えば「各時間の500hPa面の値だけ取り出したい」ような場合には、レコード長を

recl=4バイト×経度方向グリッド数×緯度方向グリッド数
として、

real(4) :: data(nlon, nlat, 1, ntim)  !層数=1

  open(11, file=ifile, access="direct",status="old", form="unformatted", convert="little_endian", recl=4*nlon*nlat)
  do t=1,ntim
  read(11, rec=(t-1)*nlev+22) data　　
  close(11)

とすることもできます。22というのが、鉛直層の格納順序で見た時の500hPaの順番（22番目）です（使うデータによっても異なるので要確認）。
この場合、レコード長＝レコード単位を水平グリッド数としていますから、時刻tにおける500hPa高度（22番目）のグリッドデータのレコード番号は(t-1)*nlev+22ということになります。

あと地味に厄介なのがエンディアンの指定で、リトルエンディアンとビッグエンディアンを間違えると値が発散したりしてしまうため、計算実行時には注意が必要です。逆にいえば、計算結果が明らかにおかしいときには計算ミスかエンディアンのミスが疑われます。

なお、前回紹介した方法でwgribを使ってバイナリに変換した場合、変換時のエンディアンはCPUに依存するそうです。そのため、ほかの人からもらったデータではエンディアンが違う場合があります。
用いるファイルのエンディアンがすべて共通している場合には、コンパイル時に

gfortran -fconvert=little(または big)-endian -frecord-marker=4 hoge.f90

のようにオプションをつけてまとめて指定してしまってもよいですが、ファイルごとにエンディアンが異なる場合は読み込み時のopen文で

open(11, file=ifile, access="direct",status="old", form="unformatted", convert="little_endian", recl=4*nlon*nlat*nlev*ntim)

のようにconvert="little_endian"などと個別に設定してあげる必要があります。

このようにしてデータを読み込めたら、あとは配列同士で適当に計算して欲しい値を作ってあげればOKです。

出力も読み込みと同じく、

open(12, file=ofile, form="unformatted", status="replace", access="direct", convert="little_endian", recl=4*nlon*nlat*nlev*ntim)
write(12, rec=1) output_array
close(12)

のように行います。

上記のようにして作成したデータをGrADsで描画するためには、データを書き込んだファイルとは別に「コントロールファイル」というものが必要です。
詳しくは次の記事をご覧ください。