これは何ですか
最近LeafletによるタイルベースのWeb地図で天気図っぽいものを色々作っています。たとえばこれですね。
地点データは一通りできたので、次は格子データなんですが、まずはデータ形式から考えなければいけません。
Web friendlyというかLeaflet/JavaScriptで楽に扱えるデータ形式として、地点データは細かいことを置いとくとJSONで決まりといえそうですが、格子データは自明ではありません。
模索から始まるので、答のコードから示しはじめても意味が伝わりません。いろいろの配慮をまずは共有するところからしましょう。という話です。
入力形式
GISC Tokyo ウェブサイトから1.25度格子GSM解析・予報値が取得できます。たとえばこんなかんじ…とはいっても常時取得の開発運用に手間がかかるのですがまあそれは別の機会に。かつては WMO DDB という名前でやっていたものの後継です。
これはGRIB Edition 1です。もういまどきEdition 1のデータサービスは再解析とここくらいしかありません。逆に言えば処理プログラムはGRIBの豊饒な機能すべてに備えずとも、決め打ちでガッと書いて恥じないでいいとも言えます。
すべてに備えなくてもいいんですが、いくつか鬱陶しいものに備える必要があって、Thinned Gridと呼ばれる変形経緯度格子かつ全球８分割電文なので始末が悪い。まあこれを1.25×1.25度格子の全球に張り合わせるくらいはあらかじめしておいた、中間加工データでないと後が不便でたまりません。ちなみに素直に全球を1.25×1.25度格子で覆うとその大きさは
(180°÷1.25° + 1)×(360°÷1.25° + 1) = 289×145 = 41905 [個格子]
となります。１ピクセル４バイトならば164 KBですが、圧縮がかかる形式ならば40 KBくらいになるでしょうか。手頃な大きさといえそうです。
その二次元面データがたくさんあるんです。次の4次元の掛け算になります。
- 参照時刻（初期時刻）：毎日4回（00Z, 06Z, 12Z, 18Z）
- 予報時間：三日半先まで6時間毎（0すなわち解析,6,12,18,24,30,36,42,48,54,60,66,72,78,84h）および12Z初期値に限り八日先まで12時間毎（96,108,120,132,144,156,168,180,192h）、つごう最大24種類
- 鉛直位置：最大18種類（地上ならびに1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10hPa）
- 要素：おおむね6種類、最大11種類（ジオポテンシャル（地上を欠く）、東西風、南北風、気温、相対湿度（300hPa以下に限る）、鉛直風（300hPa以下に限りまた地上を欠く）、渦度（500hPaのみ）、速度ポテンシャル及び流線関数（850・250hPaのみ）、海面更正気圧（地上のみ）、総降水量（地上のみ、解析以後の積算値のため解析を欠く））
なにしろこの掛け算なのでまじめにやると1000を軽く超えます。まあそうは言っても高層天気図 https://www.jma.go.jp/jp/metcht/suuchi.html が全部作られているかというとそうでもなく、注目ポイントは決まってくるのでそこを絞り込むと100くらいになるでしょうか。
出力形式
画像フォーマット（PNG, JPGとか）しかないでしょう。気象屋さんならnetCDFを挙げたくなるでしょうが、僕は違うと思うんだ。中間データとはいえ、ウェブに張り込めないと専用ソフトウェアなしに中身が検査できない。そういうのはだめだと思うんです。
画像なら Leaflet の ImageOverlay クラス で貼り付けられます。リファレンスを読むと、経緯度boundsを与えれば貼り付けられることがわかります。
もちろんメルカトル投影でなければなりません。まあざっくりWeb Mercator空間（緯度南北85度以下）内に限定すると、ピクセル数が2倍ですね。
南北に隣接したラインで相関が高いとなると、ひょっとするとJPEG圧縮が効くかもしれませんが、基本はPNGで考えましょう。
次の迷いどころは、格子間隔をタイル標準に合わせて全球２冪分割にすべきかですね。ちょっと迷ったのですが、赤道円288分割ならば補間によって「なまる」のが防げるので、289x289ピクセルにしておきましょうかね。たしかに289格子にしても南北方向には補間がされざるを得ず、また熱帯外では thinned grid によって288より少ない緯度円格子数からの補間がされてしまうのですが、結果を無駄に大きくしないことを重視して。
補間は、絶対「なまる」のを防ぐと思えば高尚な方法もあるんでしょうけれど、出力格子ごとに（入力格子からの距離による）重みと入力格子値×重みを積算して、最後に割るなら入力を1回読むだけでできます。ワークメモリが足りていてI/Oが遅いという前提ではこれが速いでしょう。
数値フォーマットはさしあたり国土地理院の標高タイルの方式にしてみましょう。24ビット符号付整数（欠損地 0x800000 を除外して -8,388,607～+8,388,607 という余裕のある空間なので、「現実的なスケール、オフセットなし」の固定的運用でもオーバーフローを気にしなくてすみそうです。
（風ベクトルを色またはバンドで表現といったことが考えられなくもないですが、まあちょっと発展課題ということで）
|パラメタ番号|要素|単位|想定値域|
|:---------|:---|:---|:---|:------|
|2|海面更正気圧|0.1 hPa|800 0～1100 0|
|7|ジオポテンシャル高度|gpm|-400～100 000|
|11|気温|0.1 K|150 0～400 0|
|14|相当温位|0.1 K|150 0～400 0|
|18|湿数|0.1 K|150 0～400 0|
|33|風の u 成分|0.1 m/s|-300 0～+300 0|
|34|風の v 成分|0.1 m/s|-300 0～+300 0|
|35|流線関数|検討中||
|36|速度ポテンシャル|検討中||
|39|鉛直速度|0.01 Pa/s|-100 00～100 00|
|43|相対渦度|1e-6 /s|-1000～1000|
|52|相対湿度|%|0～100|
|61|総降水量|kg/m2|0～1024|
UIと構成管理
さて、（参照時刻×予報時間×鉛直層×要素）の組み合わせから、少数の出力したいものを選び、その設定を変更しやすくしておきたいわけです。経験上、独自文法の設定ファイルを要求すると使うのがものすごく億劫になってしまうので、なるべくコマンドラインにしておきたいでしょう。
とすると、任意個の出力ファイルを列挙させて、そのファイル名から（あるいはファイル名に付随して）設定を読み取らせるのがよさそうです。
さしあたり全球289x289格子にしますが、いずれ z/x/y.png 型 256x256 格子のファイル名に拡張するので念頭に置きます。
任意個引数は最後にという黄金ルールからすると、入力ファイル名は１つに限定して先頭に置くのでしょう。
プログラムは最初に引数を全チェックして、収集すべき変数を列挙し、Thinned GRID分のメモリを確保してから入力ファイルをスキャンして格子補間することになります。コマンドライン引数の処理順が記載順でなく気味が悪いですが、黄金ルールを優先しましょう。
相当温位や風速（風ベクトル絶対値）を出力できるようにするには、必要な入力をThinned GRIDで確保してから変換することになるのでしょう。
メタデータ
ファイル内にメタデータを記載すべきか。
まずは記述場所。いくらでも複雑なPNGを作ることができますが、まあ穏健にテキストコメントを入れることにするなら、
せいぜいPNG specification のkeywordsくらいです。
記述内容ですが、だいたい次の３通り
- ジオリファレンス: z/x/y.png 型ならば不要。全球PNGもべつに要らない気がする。
- 単位: 標高タイルで cm 単位だと明示していない例をみれば、不要な気もするので、human readable description として。
- what's this
とすると例えばこんなかな。
- Title: c34p3/2020-01-28T00Z/2020-01-28T18Z/SURF/T
- Description: T [0.1K] on SURF at 18 hours since 2020-01-28T00Z (valid 2020-01-28T18Z) from process 3 centre 34
- Software: https://github.com/etoyoda/grib1tile