[GDAL&C#]gdal_csharp.dllを使って二次元配列の情報をGeoTiffファイルに書き出す

2020/4/2 C#, gdal, GIS, OGR, プログラミング

/* 追記:2007 Oct.18*/

こちらの記事もご参照ください
[QGIS & Gdal]Geotiffファイルのピクセルサイズが一定でない? ~その1~

/* 追記ここまで */

一つ前のエントリ「[DEM]数値標高モデルを使った地形表現」のプログラムの中の話です。

C#を使ってGeoTiffファイルを作ろうと最初に試したのが、LibTiff.Netを使う手法。
Tiffファイルを作るところまでは上手くいったのですが、Tiffファイルに位置情報を書き込むところが上手くいきません。
具体的にはTiffクラスの「SetField」というメソッドで
結局、ソースを見たら、位置情報・ピクセルサイズのタグは(別メソッドで取得はできても)SetFieldメソッドでの編集は出来ないようなので諦めました。

次に試したのが、GDALのC#用ラッパーを使う方法。
GDALのインストールディレクトリの中に「csharp」というディレクトリがあって、中には「gdal_csharp.dll」なんて如何にもなDLLファイルがあります。
さらに調べてみたら、VB.Netでの事例1)VB.NETでgdal/ogrを使ってみる。 [Chapter 1]-OpenなGISのこともあるようなので、頑張ってみました。

まず、手始めにVisualStudioの参照設定で「gdal_csharp.dll」を参照します。
次に、GDAL関連の環境変数を設定して、GdalクラスのAllRegister()でGDAL実行の準備をします。
設定が必要な環境変数は、PATH(システム環境変数)、GDAL_DATA、GDAL_DRIVER_PATH、PROJ_LIB、GDAL_CACHEMAX、CPL_TMPDIRです。下記のコードで環境変数設定＆AllRegistrメソッドを実行しました。

//冒頭にこれが書いてある前提とします。
//System.Data名前空間のDataSetクラスとの衝突を回避するためこの表記です。
using GDAL = OSGeo.GDAL;

(中略)

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();

ここで、AnyCPUのままビルドした場合(or GDALと設定したプラットフォームが一致しない)に「System.TypeInitializationException」が出て、「’OSGeo.GDAL.GdalPINVOKE’ のタイプ初期化子が例外をスローしました。」と言われてしまうことがあるかと思います。この場合、GDALの64bit/32bitを確認して、構成マネージャで適切なプラットフォーム(GDALが32bitならx86、64bitならx64)を設定してあげれば動く(はずです)。

次に、本番の処理に入ります。GeoTiffにしたいデータが、左下原点の二次元配列(Double型、下図参照)に入っていた場合に、GDALをつかってのGeoTiffへの変換のコードを。

データが格納された配列のイメージ

早速、サンプルコードです。実際の処理のポイントは…
14行目：GeoTiffを作ると宣言
20行目：GeoTiffファイルを作成(外枠)
23・24行目:位置情報設定
26行目：データの中身を作成。20行目でバンド数を設定したうえで、引数を増やしながらこの行(dsOut.GetRasterBand)を繰り返すと、複数バンドの入ったGeoTiffが作れます。
28～38行目：左下原点の二次元配列を、左上原点の1次元配列に。最初から1次元配列でデータを作れば、この工程は不要になりますね。
41・42行目：バンドデータを書き込み
46行目：測地系を設定。Proj4形式JGD2000(EPSG:4612)を書こうとしたら失敗しました。
49行目：GeoTiffファイルを保存。

//System.Data名前空間のDataSetクラスとの衝突を回避するためこの表記です。
//20行目でOSGeo.GDAL.DataSetクラスを使っています。
using GDAL = OSGeo.GDAL;

static void WriteGeoTiff(
    double[,] dZ, //標高データ
    double Xmin, double Ymax, //画像の左上の座標
    double dX, double dY, //1ピクセルあたりのサイズ
    string Projection, //測地系(OGC WKTフォーマットで、Proj4でもOK？)
    string OutFile) //出力するGeoTiffファイルの名称
{

    //出力するファイルの形式
    GDAL.Driver drv = GDAL.Gdal.GetDriverByName("GTiff");
    
    int iW = dZ.GetLength(0);
    int iH = dZ.GetLength(1);
    
    //出力ファイル：ファイル名・幅・高さ・バンド数・ピクセルのデータ型、オプション(null)
    using (GDAL.Dataset dsOut = drv.Create(OutFile, iW, iH, 1, GDAL.DataType.GDT_Float64, null))
    {
        //位置情報(Xmin, ΔX/px, 0, Ymax, 0, ΔY/px)
        double[] dGeo = new double[] { Xmin, dX, 0d, Ymax, 0, dY };
        dsOut.SetGeoTransform(dGeo);
    
        using (GDAL.Band bd = dsOut.GetRasterBand(1))
        {
            double[] dVal = new double[dZ.GetLength(0) * dZ.GetLength(1)];
            for (int n = 0; n < dVal.Length; n++) { dVal[n] = double.NaN; }
            int j = 0;
            for (int y = iH - 1; 0 <= y; y = y - 1)
            {
                for (int x = 0; x < iW; x++)
                {
                    dVal[j] = dZ[x, y];
                    j++;
                }
            }
    
            //データ書き込み
            bd.WriteRaster(0, 0, iW, iH, dVal, iW, iH, 0, 0);
            bd.FlushCache();
        }
        //Projection
        if (Projection != "") { dsOut.SetProjection(Projection); }

        //データ書き込み
        dsOut.FlushCache();
    }
}

脚注