[C# GDAL 気象]数値予報モデルGPVを読み取る(その２)

前回の記事でGRIB2形式で提供される数値予報モデル(GPV)を、GDALを使って開いてみました。GRIB2ファイルには複数の時刻・高度(気圧面)の気象要素(風向風速や気温など)が記録されています。GDALではそれぞれを「バンド」として扱い、それぞれのバンドの時刻・高度・要素はバンドのDescriptionやMetaDataを見ることによって特定可能というところまでを前回の記事に記載しました。
しかし、MetaDataは文字列の配列になっているなどそのままでは扱いづらいため、バンドの情報をとりまとめるクラスを作ってみました。単にDescriptionとMetaDataを読み取って、バンドを特定するのに特に重要と思われる時刻(初期値・予報)、予報時間、高度(地上・気圧面)をプロパティに整数型や日付時刻型で保存するとともに、その他の情報もMetaDataからキーと値の形でDicitionaryに保存してみました。

using GDAL = OSGeo.GDAL;

public class cGrib2Info
{
    //フィールド
    Dictionary<string, string> _dcMeta;
    int _iFT, _iPres;
    DateTime _dtInit, _dtValid;

    //プロパティ
    /// <summary>
    /// 予報時間:初期値時刻からの経過時間(hrs)
    /// </summary>
    public int FT { get { return (_iFT); } }
    /// <summary>
    /// 気圧面:地上の場合は1013を返す
    /// </summary>
    public int Pres { get { return (_iPres); } }
    /// <summary>
    /// メタデータ
    /// </summary>
    public Dictionary<string, string> MetaData { get { return (_dcMeta); } }
    /// <summary>
    /// バンドの内容
    /// </summary>
    public string Content { get { return (_dcMeta["GRIB_COMMENT"]); } }
    /// <summary>
    /// 初期値時刻
    /// </summary>
    public DateTime Initial { get { return (_dtInit); } }
    /// <summary>
    /// 予報時刻
    /// </summary>
    public DateTime Valid { get { return (_dtValid); } }
    
    /// <summary>
    /// バンド情報のコンストラクタ
    /// </summary>
    /// <param name="bd">バンド</param>
    public cGrib2Info(GDAL.Band Band)
    {
        //GRIB以外の時には例外を出す
        if (Band.GetDataset().GetDriver().ShortName != "GRIB") { throw new Exception("Dataset isn't GRIB"); }

        //気圧面(地上の時は1013とする)
        string sDec = Band.GetDescription();
        if (sDec.Length < 154) { _iPres = 1013; }
        else
        {
            if (!int.TryParse(sDec.Substring(0, sDec.Length - 154), out _iPres))
            {
                _iPres = 1013;
            }
        }

        //メタデータ
        _dcMeta = new Dictionary<string, string>();
        string[] ss = Band.GetMetadata("");
        foreach (string s in ss)
        {
            string[] s2 = s.Split('=');
            _dcMeta.Add(s2[0].Trim(), s2[1].Trim());
        }

        //予報時間
        string sf = _dcMeta["GRIB_FORECAST_SECONDS"];
        ss = _dcMeta["GRIB_FORECAST_SECONDS"].Split(' ');
        _iFT = int.Parse(ss[0]) / 3600;

        //初期＆予報時刻
        _dtInit = new DateTime(1970, 1, 1, 0, 0, 0, 0);
        ss = _dcMeta["GRIB_REF_TIME"].Split(' ');
        _dtInit = _dtInit.AddSeconds(double.Parse(ss[0]));

        _dtValid = new DateTime(1970, 1, 1, 0, 0, 0, 0);
        ss = _dcMeta["GRIB_VALID_TIME"].Split(' ');
        _dtValid = _dtValid.AddSeconds(double.Parse(ss[0]));
    }
}

このクラスを使って、下記のコードで気象庁のGRIB2形式のGPVを読み取ってみたのが下記のサンプルコードです。
14行目のGDAL有効化の部分は前回の記事のサンプルコードをご参照ください。
25行目でバンドのメタデータ(初期値時刻・高度・要素)を取得するクラス(上記参照)のインスタンスを発生させて、続く26・27行目で気圧面・初期値時刻・予報時間・要素を表示させています。

using System.IO;
using GDAL = OSGeo.GDAL;
using Grib2;

//略

const string sGdalDir = @"C:\Program Files\GDAL";
const string sInGRIB2 = @"D:\hogehoge\Z__C_RJTD_20180420090000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin";
//const string sInGRIB2 = @"D:\hogehoge\ReadGRIB2\Z__C_RJTD_20180420090000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin";

static void Main(string[] args)
{
    //GDAL有効化
    InitGDAL();

    //とりあえず開いてみる
    using (GDAL.Dataset ds = GDAL.Gdal.Open(sInGRIB2, GDAL.Access.GA_ReadOnly))
    {
        //バンド数
        Console.WriteLine("ds.RasterCount:{0}", ds.RasterCount);

        //バンドの内容
        for (int i = 1; i <= ds.RasterCount; i++)
        {
            cGrib2Info cGInfo = new cGrib2Info(ds.GetRasterBand(i));
            Console.WriteLine("Band:{0},P:{1} Init:{2}, FT:{3}, Element:{4}",
                i, cGInfo.Pres, cGInfo.Initial.ToString("yyyy-MM-dd_HH:mm"), cGInfo.FT, cGInfo.Content);
        }
    }

    Console.ReadLine();
}

実際の読み取り例は下記のようになりました。

ds.RasterCount:190
Band:1,P:1013 Init:2018-04-20_09:00, FT:0, Element:Pressure reduced to MSL [Pa]
Band:2,P:1013 Init:2018-04-20_09:00, FT:0, Element:Pressure [Pa]
・・・中略・・・
Band:189,P:1013 Init:2018-04-20_09:00, FT:14, Element:01 hr Total precipitation [kg/(m^2)]
Band:190,P:1013 Init:2018-04-20_09:00, FT:14, Element:Downward short-wave radiation flux [W/(m^2)]

また、サンプル９行目(コメント行)を有効にして、高層(気圧面)のデータを読み取った場合の読み取り結果は下記のようになりました。

ds.RasterCount:552
Band:1,P:1000 Init:2018-04-20_09:00, FT:0, Element:Geopotential height [gpm]
Band:2,P:1000 Init:2018-04-20_09:00, FT:0, Element:u-component of wind [m/s]
・・・中略・・・
Band:551,P:100 Init:2018-04-20_09:00, FT:15, Element:Temperature [C]
Band:552,P:100 Init:2018-04-20_09:00, FT:15, Element:Vertical velocity (pressure) [Pa/s]

これで、欲しい時刻・要素・高度のGPVのみを抽出して、Geotiffで出力するなり、欲しい地点の情報をテキストで出力するなど出来そうです。
なお、GDALの公式サイト(2018.4.22閲覧)を見ると「GMT Compatible netCDF」の生成は出来るように書かれていますが、私の環境ではC#のAPIを使ってGMTのnetCDFファイルを作成することは出来ませんでした(ただし、いったんGeoTIFFで吐き出した後でgdal_translateしたらファイルの作成＆GMTでの作図が出来た)¹⁾例外:System.ApplicationException「GDALDriver::Create() … no create method implemented for this format.」が発生しました。。

補足
高層のGPV(気圧面)のDescriptionが下記のようになっています。気圧の桁数が妙に多くて、とりあえず先頭部分だけを読み取ることで回避しています。

975000000000000240000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000[Pa] ISBL="Isobaric surface"

一方で、オフィシャルサイトでの例は下記のようになっています。

Description = 100000[Pa] ISBL="Isobaric surface"

なにか、この辺りで情報をお持ちの方がいらっしゃいましたら、教えて頂けると幸いです。

GDALを2.4.0にアップデートしたら上記の問題は解決していました。
単にGDALのバグだったようです。