[OGR/GDAL]Geopackageを試してみた(C#)

QGIS3でデフォルトのデータ形式がシェープファイル(ESRI ShapeFile)からGeopackageに変わりました。シェープファイルには、OGRの(C#というより.Netの)APIからだとデータに全角文字が含まれていると(cpgファイルを書いていても)文字化けしてしまうこと、シェープファイルは１つのデータセットに複数のファイル(最低、shp・shx・dbfの３ファイル)が必要で管理が面倒といった問題点があります。一方でGeopackageの場合には、１つのファイルで完結し文字化けも(今までは)発生しない¹⁾文字化け、多発してましたorzするので、これを機に移行してみたいと思った次第です。
ところが・・・、軽い処理でもなんか重たい。

ということで、正方形を1万個描画する処理に要する時間を、末尾に記載したプログラムで比較してみました。
その結果が、シェープファイル:0.18秒、Geopackage:29.30秒。

圧倒的にGeopackageが遅いという結果になりました(複数回実行して比較するまでもない大差！！)。
Geopackageでの作図処理中に「gpkg-journal」ファイルが出来ていた²⁾GeopackageはSQLiteのファイルで、SQLiteにて大量にInsertなどの書き込み処理をする場合にトランザクションの開始終了を明示しないと非常に遅くなる。ので、トランザクション処理が抜けているという訳ではなさそうですけど、倍半分どころか100倍以上の差がついてしまいました。
コレではちょっと実用になりません(´・ω・｀)
GDAL/OGRのアップデートで改善は見込めるのでしょうか・・・

OGR.DataSourceにトランザクションの開始と終了を明示するメソッドがあるのを思いだして実装しました。具体的にはサンプルコードの51行目でトランザクションを開始して、84行目で終了しています。
すると、Geopackageでの処理時間が大幅に向上して、トランザクションの明示なしだと30秒近かった処理時間が0.22秒へと大幅に短縮しました。
複数回実行してみたところ、シェープファイルよりは少し時間がかかっているようですが、差はわずかですし、他の利点を考えると実用上全く問題のない範囲となりそうです。

ちなみに、当方の環境はこんな感じ。
Windows10(Ver1803、Build17134.286)
GDAL2.3.0(gdal-203-1911-x64-core.msi、2018.10.11時点で2.3.2が最新)
CPU:Core i7-8700、メモリ:16GB

using System.IO;
using OGR = OSGeo.OGR;
using OSR = OSGeo.OSR;

//中略

const string sOutD = @"D:\テストディレクトリ"; //描画結果を出力するディレクトリ
const string sGdalDir = @"C:\Program Files\GDAL"; //GDALのインストールディレクトリ
const string s3100 = "PROJCS[\"JGD2000 / UTM zone 54N\",GEOGCS[\"JGD2000\",DATUM[\"Japanese_Geodetic_Datum_2000\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6612\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4612\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",141],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],AUTHORITY[\"EPSG\",\"3100\"],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]";
const double d = 100d; //描画する図形の大きさ

static void Main(string[] args)
{
    #region "OGRの準備(環境変数Pathの設定→Ogr.RegisterAll())"
    string 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);
    OGR.Ogr.RegisterAll();
    #endregion

    OSR.SpatialReference sref = new OSR.SpatialReference(s3100); //UTM54帯
    OGR.wkbGeometryType gtp = OGR.wkbGeometryType.wkbPolygon;  //出力するジオメトリタイプ(ポリゴン)

    OGR.Driver[] drv = new OGR.Driver[2];
    drv[0] = OGR.Ogr.GetDriverByName("ESRI Shapefile"); //シェープファイル
    drv[1] = OGR.Ogr.GetDriverByName("GPKG"); //Geopackage
    string[] sExt = new string[] { "shp", "gpkg" }; //ファイルの拡張子

    for (int i = 0; i < drv.Length; i++) //シェープファイル→Geopackageの順
    {
        //時刻計測用
        System.Diagnostics.Stopwatch sw = new System.Diagnostics.Stopwatch();
        sw.Start();

        //出力するファイル名(DataSource)
        string sOutF = Path.Combine(sOutD, string.Format("Test{0}.{1}", i, sExt[i]));
        using (OGR.DataSource dsOut = drv[i].CreateDataSource(sOutF, null))
        using (OGR.Layer lyOut = dsOut.CreateLayer("Mesh", sref, gtp, null))
        {
            OGR.FeatureDefn fdOut = lyOut.GetLayerDefn();
            OGR.FieldDefn[] fld = new OGR.FieldDefn[2];
            fld[0] = new OGR.FieldDefn("X", OGR.FieldType.OFTInteger);
            fld[1] = new OGR.FieldDefn("Y", OGR.FieldType.OFTInteger);
            for (int j = 0; j < fld.Length; j++) { lyOut.CreateField(fld[j], 1); }

            long lFid = 0;
            dsOut.StartTransaction(0); //トランザクションの開始を明示
            for (int iX = 0; iX < 100; iX++)
            {
                for (int iY = 0; iY < 100; iY++)
                {
                    double dX0 = iX * d;
                    double dY0 = iY * d;

                    using (OGR.Feature fOut = new OGR.Feature(fdOut))
                    using (OGR.Geometry gmOutP = new OGR.Geometry(gtp))
                    using (OGR.Geometry gmOutR = new OGR.Geometry(OGR.wkbGeometryType.wkbLinearRing))
                    {
                        //FID
                        fOut.SetFID(lFid); lFid++;

                        //ジオメトリ作成：リング(外側なので時計回り)
                        gmOutR.AddPoint(dX0, dY0, 0d);
                        gmOutR.AddPoint(dX0+d, dY0, 0d);
                        gmOutR.AddPoint(dX0+d, dY0+d, 0d);
                        gmOutR.AddPoint(dX0, dY0+d, 0d);
                        gmOutR.CloseRings();
                        gmOutP.AddGeometry(gmOutR);
                        fOut.SetGeometry(gmOutP);

                        //フィールド作成
                        fOut.SetField("X", iX);
                        fOut.SetField("Y", iY);

                        //レイヤへ
                        lyOut.CreateFeature(fOut);
                    }
                }
            }
            dsOut.CommitTransaction(); //トランザクションの終了を明示
        }

        sw.Stop();
        Console.WriteLine("処理時間({0}):{1:0.00}秒", sExt[i], sw.ElapsedMilliseconds / 1000d);
        sw.Reset();
    }
    Console.WriteLine("END");
    Console.ReadLine();
}