[GDAL/OGR In CSharp] C#のOGRを使って図形の情報を読み取る(位置・属性)

前回までは、C#でGDAL/OGRを使ってシェープファイルを作図する記事を書いてまいりましたが、今回はシェープファイルの内容(ジオメトリ・属性)を読み取ってみたいと思います。

基本的には、Driver→DataSource→Layer→Feature→Geometryの順番に開いていき、ジオメトリや属性テーブルの内容を読み取っていきます。
早速ですが、サンプルコードを。

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

/*中略*/

const string sGdalDir = @"C:\Program Files\GDAL"; //GDALのインストールディレクトリ
const string sShp = @"D:\hoge\2_poly.shp"; //読み取るシェープファイル
const string sCsv = @"D:\hoge\2_poly_ct.csv"; //読み取った結果を書き込むCSV

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

    OGR.Driver drv = OGR.Ogr.GetDriverByName("ESRI Shapefile");

    //使用するファイルを開く
    using(StreamWriter sw = new StreamWriter(sCsv, false, Encoding.ASCII))
    using (OGR.DataSource ds = drv.Open(sShp, 0)) //SHP:DataSource
    using (OGR.Layer ly = ds.GetLayerByIndex(0))  //レイヤ
    {
        sw.WriteLine("Centroid_X,Centroid_Y,Attribute,GeomNo,PointNo");
        
        //地物をレイヤ先頭の地物にリセットする(?)
        //参考→ http://www.gdal.org/classOGRLayer.html#aad0f2cd7f0587584b8f382c6a913583c
        ly.ResetReading();
        
        long lCnt = ly.GetFeatureCount(0); //レイヤ内の地物の数
        for (long l = 0; l &lt; lCnt; l++)
        {
            OGR.Feature f = ly.GetFeature(l);
            OGR.Geometry g = f.GetGeometryRef(); //地物に登録されたジオメトリ
            
            string sAtt = f.GetFieldAsString("fStr");
            //OGR.Geometry gc = g.Centroid(); //コメントを外すと地物の重心を取得(穴の空いたポリゴンの重心をとっているみたい)
            //sw.WriteLine("{0},{1},{2}", gc.GetX(0), gc.GetY(0), sAtt);
            for (int i = 0; i &lt; g.GetGeometryCount(); i++) //ジオメトリの中のジオメトリを見ていく
            {
                OGR.Geometry gg = g.GetGeometryRef(i);
                for (int j = 0; j &lt; gg.GetPointCount(); j++) //ジオメトリの構成点を見ていく
                {
                    sw.WriteLine("{0},{1},{2},{3},{4}", gg.GetX(j), gg.GetY(j), sAtt, i, j); //座標出力
                }
            }
        }
    }
}

44・45行目がコメントアウトされていますが、ここを外すとポリゴンの重心を出力します。
ポリゴンは穴が空いていますが、穴が空いた状態の重心を出力しているようです。

このコードでは、ポリゴンを構成するリング(閉じた線)の座標をすべて読み取っています。
穴あきポリゴンのため、1つの地物にリングは二つあります。
この二つのリングを識別するためにCSVファイルではGeomNo(サンプルソースの変数名:i)を出力しています。

出力されたCSVファイル(冒頭)

もとのシェープファイルと、構成点位置・属性を出力したCSVファイルをQGISで図化したのが下図です。ここではGeomNo=0が赤、GeomNo=1が青になっています。

ポリゴンシェープファイルの読み取り結果

毎度のオマケ:OGR.Ogr.RegisterAll()で「System.TypeInitializationException」が出たときは、構成マネージャーでx68かx64に変更(64bit版GDALを入れていればx64、32bit版ならx86)すると幸せになれるでしょう。