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

気象庁の数値予報モデルの計算結果は、気象庁サイト(数値予報天気図)様々な種類の天気図が公表されているものの、モノクロだったり欲しい時刻の天気図がなかったりします。
ということで、数値予報の生データ(GPV)を京都大学生存圏研究所のサイトにある数値予報GPVをダウンロードしてきて、GMTで自分好みの図を作図してみようと思います。

気象庁の数値予報データ(GPV:Grid Point Value)はGRIB2という形式のバイナリデータで、そのままテキストで読み取ることは出来ないため、何らかの方法でデコードする必要があります。NOAAがwgrib2というツールを公開してくれているのですが、色々と面倒くさそうだったので、私個人が慣れているGDAL¹⁾GDALはGRIB2にも対応しているをC#から触ってデコードしてみようと思います。

早速、GRIB2形式のGPVをダウンロードしてきて、下記のコードを使ってGDALをC#から使ってどんな情報が入っているのかを覗いてみます。

//冒頭にはこれが書かれているものとします
using System.IO;
using GDAL = OSGeo.GDAL;

/*中略*/
const string sGdalDir = @"C:\Program Files\GDAL"; //GDALのインストールディレクトリ
const string sInGRIB2 = @"ダウンロードしたGPVファイルのフルパス";

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);
        Console.WriteLine("ds.GetDescription:{0}", ds.GetDescription());

        //バンドの内容
        for (int i = 1; i <= ds.RasterCount; i++)
        {
            Console.WriteLine("Band{0}, GetDescription:{1}", i, ds.GetRasterBand(i).GetDescription());
        }

        //止める
        Console.ReadLine();
    }
}

//GDAL有効化【AnyCPUだとバグるので注意】
public static void InitGdal()
{
    const string sGdalDir = @"C:\Program Files\GDAL";
    string sVal;

    //PATH
    sVal = Environment.GetEnvironmentVariable("PATH") + ";" +
    sGdalDir + ";" +
    Path.Combine(sGdalDir, "csharp") + ";" +
    Path.Combine(sGdalDir, "gdal-data") + ";" +
    Path.Combine(sGdalDir, "gdalplugins") + ";" +
    Path.Combine(sGdalDir, "projlib");
    Environment.SetEnvironmentVariable("PATH", sVal);

    //GDAL_DATA
    Environment.SetEnvironmentVariable("GDAL_DATA", Path.Combine(sGdalDir, "gdal-data"));
    GDAL.Gdal.SetConfigOption("GDAL_DATA", Path.Combine(sGdalDir, "gdal-data"));

    //GDAL_DRIVER(Plugins)
    Environment.SetEnvironmentVariable("GDAL_DRIVER_PATH", Path.Combine(sGdalDir, "gdalplugins"));
    GDAL.Gdal.SetConfigOption("GDAL_DRIVER_PATH", Path.Combine(sGdalDir, "gdalplugins"));

    //PROJ_LIB
    Environment.SetEnvironmentVariable("PROJ_LIB", Path.Combine(sGdalDir, "projlib"));
    GDAL.Gdal.SetConfigOption("PROJ_LIB", Path.Combine(sGdalDir, "projlib"));

    GDAL.Gdal.SetConfigOption("GDAL_CACHEMAX", "100000");
    GDAL.Gdal.SetConfigOption("CPL_TMPDIR", @".\");

    GDAL.Gdal.AllRegister();
}

上記のコードを実行した結果がこちら。

GRIB2に含まれる様々な要素(気圧・風向風速・気温etc)が多数のバンドに分かれて保存されているようです。
Dataset.GetDescriptionメソッドの結果を見ると、気圧面・高度の種類が記載されているようですが、予報時刻や何の要素が記録されたバンドなのかはわかりません。
色々と試行錯誤してみたところ、予報時刻や気象要素はGDAL.Band.GetMetadataメソッドで読み取ることが出来るようです。
そこで、前述のコードにメタデータを読み取る部分を追加してみます。追加した部分のコードと、追加したコードによる出力結果は下記の通りでした。

//バンドの内容
for (int i = 1; i <= ds.RasterCount; i++)
{
    Console.WriteLine("Band{0}, GetDescription:{1}", i, ds.GetRasterBand(i).GetDescription());
    string[] ss = ds.GetRasterBand(i).GetMetadata("");
    for (int j = 0; j < ss.Length; j++) //追加した部分
    {
        Console.WriteLine("Band{0}, Metadata[{1}]:{2}", i, j, ss[j]);
    }
}

この内容を見ると、気象要素は「GRIB_COMMENT」に、初期値時刻と予報時刻は「GRIB_REF_TIME」と「GRIB_VALID_TIME」にそれぞれUNIXタイムで記録されていました。また、予報時間は「GRIB_FORECAST_SECONDS」に秒単位で記録されていました。

ひとまず、これでGRIB2からGPVの値を読み取る手順はなんとなくわかってきました。

次の記事で、GRIB2の情報をもう少しスッキリと読み取るクラスを使って、バンドの情報を読み取ってみたいと思います。
次の記事→ [C# GDAL 気象]数値予報モデルGPVを読み取る(その２)