WMO FM92 GRIB Edition 2 についてはこちらをご覧ください。
https://qiita.com/e_toyoda/items/ce7497e1a633b16f1ff1
そんでもってそれを解読というか一部だけ抽出するプログラムを作って、せっかくだからライブラリとしてブラッシュアップしているという話です。誰の役にも立たないかもしれないが、まあ未来の自分には。
前提
大気というのは時空間4次元を占めるので、格子データは普通はこの四次元中、典型的には経度・緯度・高度・時間を軸として張る空間に置かれます。GRIB Edition 2 はこの4次元（アンサンブルなどを含めればさらに高次元）のデータを扱うことができるのですが、二次元配列（通常は水平二次元）の集まりに分解して保存します。
GRIB Edition 2 は節の集まりです。節には IS, IDS, LS, GDS, PDS, DRS, BMS, DS, ES という種類があります。日本公定訳は通報式の訳文である指示節、識別節、地域使用節、格子系定義節、プロダクト定義節、資料表現節、ビットマップ節、資料節、終端節ですが、実際上は第0節〜第8節というほうが多いでしょう。これらはこのような繰り返し構造に従い配列されます。
IS IDS LS? (GDS (PDS DRS BMS? DS)+)+ (LS (GDS (PDS DRS BMS? DS)+)+)* ES
厳密さを犠牲にしてわかりやすく書くならば（あるいは FM92 の文書ふうに書くと）こうですね。
IS IDS (LS? (GDS (PDS DRS BMS? DS)+)+)+ ES
ここで ? は0回または1回の出現、+は1回以上の出現です。
要するには、IS, IDS, GDS, PDS, DRS, BMS は DS を修飾して、DS からみて直近の同種の節を見ればよいわけです。
ただし、ビットマップ節についてだけは「直前のビットマップ節と同じ」という特殊な短縮形式のBMSがあります。
GRIB はシーケンシャルファイルで、何個目のDSが何オクテット目（オクテットとはバイトと同じと思ってください）にあるかは読んでみなければわかりません。豊田の知るユースケースは「ｎ番目のDSを取得したい」ではなく「12時間予報の500hPa面気温を抽出したい」といったものです。ファイルの先頭から各節を読んでいって、ほしいものが見つかったら後段処理に投げつけるということになります。
「12時間予報の500hPa面気温」というような絞り込みに使われる情報は基本的にPDSに書かれていますが、ひとつだけ「気温」といったパラメタ（物理量など）を識別する3つのオクテットのうち1つはIS第7オクテットにあります。DSのオクテット列の形式はDRSやBMSに書かれています。なんだかんだでIS〜DSの全セットを揃えて後段処理に渡す必要があります。
そんなことを考えて、似た処理は何度もしそうなので、いちおう汎用GRIB2解読ライブラリとして書いたのが https://github.com/etoyoda/grib2png/blob/main/gribscan.c です。これを使って grib2png という可視化プログラムと、データ抽出・ファイル連結プログラム gribslim が作られています。
動作説明
grib2scan_by_filename()
通常の使い方であれば grib2scan_by_filename(const char *fnam) 関数が入口です。これは名前 fnam でファイルを開いてオクテット列 "GRIB" から始まる部分を探して見つかるたびにデコーダ grib2decode() に渡します。GRIB メッセージをファイルに格納するにあたり、電文としてのヘッダがついていることがあったり、複数の GRIB メッセージが連結されることがあるので、こういう動きが求められます。典型的な PNG とか JPEG みたいなわけにはいかないんですわ。
オクテット列 "GRIB" を発見するための有限状態オートマトンは説明しなくてもいいですよね。
grib2decode()
関数 grib2decode(FILE *fp, const char *locator) は stdio ストリーム fp に格納された GRIB メッセージをデコードします。第二引数 locator はエラーメッセージ用文字列です（grib2scan_by_filename はファイル名とオクテット位置を渡します）。ストリームを引数に取るようにしているので、ちょっと改造すれば標準入力を読むようにもできます。つまりこの中の処理で fseek() は使わず一方向に読み進めるのです。シーケンシャルファイルだからこんなことができるんですね。
この中では IS（16オクテット固定長）のうち、既に読んでしまった"GRIB"を除いた残り12オクテットを読み込んで、解読用の構造体 struct grib2secs (下記) を構築して、節構造の解読ルーチン grib2loopsecs() に処理を渡しています。GRIB2 の8種類の節のうち、最初の IS と最後の ES だけは固定長ですが、その間の IDS〜DS はすべて、最初の4オクテットに節の長さが書いてある可変長構造なので、それをナメて集めていくわけです。
typedef struct grib2secs {
// comes from IS (indicator section) of the GRIB message
unsigned discipline;
size_t msglen;
unsigned char *ids, *gds, *pds, *drs, *bms, *ds;
size_t idslen, gdslen, pdslen, drslen, bmslen, dslen;
void *omake;
} grib2secs_t;
唯一のエラーチェックがエディション番号で、grib2scan_by_filename() がGRIBメッセージではないオクテット列 "GRIB" が混在しているファイルを読んでしまった場合、ここでエディション番号 2 (バイナリの2です）が然るべきところにないので弾かれます。エディション0やエディション1は全く構造が違うので相手にしません。
grib2loopsecs()
GRIB2 の8種類の節のうち、IDS〜DS はすべて、最初の4オクテットに節の長さが書いてある可変長構造なので、それを順次読み取って処理します。なので、4オクテット読んで、そのバイナリ値（ビッグエンディアン）が指示する長さだけ読んで節にします。
IDS～DSの各節の構造の分析がこの関数の本題です。ループ回してざっくりこんな感じになっているわけですが（エラー処理は全略）、
grib2loopsecs(grib2secs_t *gsp, FILE *fp, const char *locator)
{
for (;;) {
secbuf = mymalloc(64*1024);
fread(secbuf, 1, 4, fp);
recl = ui4(secbuf); // 4オクテットの符号なし整数として変換
if (recl == 0x37373737uL) break; // "7777" だったら終了
if (recl > 64*1024) { // 足りんかったら拡張
secbuf = myrealloc(secbuf, recl);
}
fread(secbuf+4, 1, recl-4, fp);
switch (secbuf[4]) {
case 1:
if (gsp->ids) myfree(gsp->ids);
gsp->ids = secbuf;
gsp->idslen = recl;
break;
...
case 6:
if (secbuf[5] == 254) break; // 無視すべきBMS
if (gsp->bms) myfree(gsp->bms);
gsp->ids = secbuf;
gsp->idslen = recl;
break;
case 7:
// gsp->ds 既に破棄済み
gsp->ds = secbuf;
gsp->dslen = recl;
checksec7(gsp);
break;
}
}
}
各節の第5オクテット（節の長さの次、コードでは secbuf[4] ）がオクテット種類番号（IDS:1, LS:2, ... DS:7）なので、分類に応じて gsp->ids, gsp->gds, ... gsp->ds に保存していきます。同種のものが複数回見つかれば上書きをするので、ある種類の節について常に最後に出現したものが構造体 gsp に置かれるわけです。そして DS が見つかったところでコールバック関数 checksec7(gsp) が呼び出されるわけです。
ここらでメモリ管理について
ここらでメモリ管理の話をしましょう。動的に確保する必要がある可変長のメモリは、 *gsp とそこから指す各節になります。
*gsp のライフサイクルは簡単ですね。 grib2decode() 内の new_grib2secs() で確保して del_grib2secs() で破棄します。
各節のメモリは節を読む前に確保せねばなりませんが、その長さは先頭4オクテットを読んでみなければわかりません。現状実装ではとりあえず 64 Kbyte 確保して、だめなら拡張することにしています。DS や BMS のように巨大になりうるもの以外は 64 Kbyte におさまるので myrealloc() はしなくて済みます。まあここは後年チューニングする余地はあるかもです。
各節のメモリを破棄するタイミングがちょっとトリッキーです。
基本的に *gsp 構造体に節のデータを持たせておいて、要らなくなったら捨てるわけですから、上記要約コードの case 1: のように、同種の節が現れたら古いものを myfree() で破棄します。 case 6: の場合にはちょっとした特例があって第6オクテット（secbuf[5]）が値254 の場合に無視するようになっていますが、これは「直前のBMSと同じ」という意味で BMS 本体が入っていないので、以前のものを破棄してしまったら後続処理が困るからですね。
case 7: ではメモリを開放していません。ちょっとライブラリAPIとしては変態的ですが、gsp->ds のメモリを破棄するのはコールバック checksec7() の責任にしてあります。これはですね、DS が巨大で、複数の DS を見ながら処理するアプリケーションが後で現れる可能性を考慮してのことです。そんなにメモリケチケチしなくてもいいハズなんですけどね。まあ、後で変えるかもしれません。
サービスサブルーチン的なやつら
オクテット列から数値型への変換
extern long si4(const unsigned char *buf);
extern unsigned long ui4(const unsigned char *buf);
extern int si2(const unsigned char *buf);
extern unsigned ui2(const unsigned char *buf);
extern double float40(const unsigned char *buf);
extern double float32(const unsigned char *buf);
extern unsigned
unpackbits(const unsigned char *buf, size_t nbits, size_t pos);
GRIB はオクテット列のバイナリでデータを表現しますから、数値型との相互変換が必要です。
ただビッグエンディアンというだけでない注意事項がいろいろあるのでライブラリ化が必要です。
si4() が signed integer in 4 octets です。 ui4(), si2(), ui2() もその要領です。
float32() は IEEE 754 形式なのでまあ話はわかりやすいですが、 float40() は５オクテットの浮動小数点形式なので他ではまず見ないものですね。
いちばんめんどくさいunpackbits() はたとえば10ビット整数列とかをオクテット列に詰めてあるのを解読します。コードはやっていることの割に長いですが、頻用される10ビット、12ビットなどがバリバリハードコードされています。こういうののデバグは人生でそう何度もやりたくはないですね。
NuSDaSっぽいデータ分類整理の話
ここをちゃんと書いておかなきゃいけない気がするんだがちょっと夜も遅くなったので今日はここまで。