[GDAL/OGR]地物同士の距離＆ポリゴンに含まれるか(C#)

2019/10/4 C#, gdal, GIS, OGR, プログラミング, 地理・地質・ジオ

GDAL/OGRを使って色々な処理をしていたのですが、地物(ジオメトリ)同士の距離の計測や、ポリゴンに内包されるかを判別するメソッドがGDAL/OGRに有りました。
今まで、ゴリゴリとコードを書いていた苦労は一体＿|￣|○

地物同士の距離はOSGeo.OGR.Geometry.Distance(OSGeo.OGR.Geometry)に、
内包確認はOSGeo.OGR.Geometry.Contains(OSGeo.OGR.Geometry)に実装されていました。

サンプルとして、下図のようにポイントとポリゴンが有ったとします。ファイル名(レイヤ名)で匂わせているとおり、両者の測地系はUTM54(平面直角座標、m単位)で揃えています。
この際、ポイントの「id」には１～６の連番が振られている(図中に表示された値)ものとします。
点①はポリゴンの外、②・③はポリゴンを構成する頂点の上(QGISのスナップ機能を使用)、④・⑤はポリゴンの中、⑥はポリゴンの外周線の上(QGISのスナップ機能を使用)に有ります。

ポイントとポリゴン

ポイントとポリゴン

下記に地物同士の距離計測とポリゴンが点を内包しているかを取得するサンプルコードを書きました。
長々と書いていますが、49行目でポリゴンが点を含んでいるかをbool型で取得していることと、
50行目でポリゴンと点の距離を取得しているトコロが味噌です。

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

//中略
const string sPt = @"略\Point_epsg3100.shp";
const string sPg = @"略\Polygon_epsg3100.shp";
const string sOut = @"略\01_ContainDistance.csv";
const string sGdalDir = @"C:\Program Files\GDAL"; //GDALのインストールディレクトリ

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(sOut, false, Encoding.ASCII))
    using (OGR.DataSource dsPt = drv.Open(sPt, 0))    //点DataSource
    using (OGR.Layer lyPt = dsPt.GetLayerByIndex(0))  //点Layer
    using (OGR.DataSource dsPg = drv.Open(sPg, 0))    //ポリゴンDataSource
    using (OGR.Layer lyPg = dsPg.GetLayerByIndex(0))  //ポリゴンLayer
    {
        //出力するテキストファイルのヘッダ
        sw.WriteLine("id,Contain?,Distance");

        //先頭の地物に
        lyPt.ResetReading(); lyPg.ResetReading();

        //ポリゴン(最初の地物)
        OGR.Feature fPolygon = lyPg.GetFeature(0);
        OGR.Geometry gmPolygon = fPolygon.GetGeometryRef();

        //点の数
        long lCnt = lyPt.GetFeatureCount(0);

        for (long l = 0; l < lCnt; l++)
        {
            OGR.Feature fPoint = lyPt.GetFeature(l);
            OGR.Geometry gmPoint = fPoint.GetGeometryRef();
            bool bCont = gmPolygon.Contains(gmPoint);   //ポリゴンが点を含んでいるか？
            double dDist = gmPolygon.Distance(gmPoint); //ポリゴンと点の距離
            int iID = fPoint.GetFieldAsInteger("id");   //ポイントの属性「id」を取得

            //距離＆内包確認出力
            sw.WriteLine("{0},{1},{2:e3}", iID, bCont.ToString(), dDist); //結果を出力
        }
    }
}

上記のコードを実行した結果のテキストファイルの内容が下図の通りです。

プログラム実行結果(ポイントのid、内包有無、距離

Containsメソッドの結果は、ポリゴンの内側に作図した点④・⑤のみがTrueで、線上の②・③・⑥はFalseに鳴りました。
また、Distanceの結果は、頂点上に作図した②・③はゼロになりましたが、線上に作図した⑥はわずかですが正の値となりました。

さて、点の測地系を緯度経度(EPSG:4612)、ポリゴンは平面直角座標(EPSG:3100、UTM54)とした場合の結果は、全てContainsの結果は全てFalseで、Distanceの値は非常に大きな値が返ってきました。
このため、二つのジオメトリの測地系をそろえる必要があります。

OGRでは上記サンプルコードの45～55行目を下記のように書き換える(TransformToメソッド)ことで測地系を変更する事が出来ます。
下記のコードの味噌は４行目でポイントのジオメトリをClone()で複製しているトコロと、５行目のTransformToで測地系を変換しているトコロです。
TransformToメソッドの引数(OSGeo.OSR.SpatialReferenceクラス)はOSGeo.OGR.Geometry.GetSpatialReference()で取得しています。

for (long l = 0; l < lCnt; l++)
{
    OGR.Feature fPoint = lyPt.GetFeature(l);
    OGR.Geometry gmPoint2 = fPoint.GetGeometryRef().Clone(); //★ジオメトリを複製する★
    gmPoint2.TransformTo(gmPolygon.GetSpatialReference());   //★測地系をポリゴンに合わせる★
    bool bCont = gmPolygon.Contains(gmPoint2);   //ポリゴンが点を含んでいるか？
    double dDist = gmPolygon.Distance(gmPoint2); //
    int iID = fPoint.GetFieldAsInteger("id");

    //距離＆内包確認出力
    sw.WriteLine("{0},{1},{2:e3}", iID, bCont.ToString(), dDist);
}

上記で測地系を揃えた場合の実行結果は下図の通りで、Containsメソッドの結果に変動はありませんでしたが、Distanceの結果がごくわずかですが変動しています。
座標変換の誤差の出方によっては、Containsメソッドの結果にも影響があるかもしれませんので境界線上の地物の内包確認は注意が必要かもしれません。

測地系を揃えた場合の実行結果