[DotSpatial]ポイントがポリゴンに含まれているか判別する③

DotSpatialでポイントとポリゴンの重なり判定について、今回はポイントとポリゴンのシェープファイルの測地系が違った場合にどうなるのかを検討してみます。使用するコードは前回と同じですが、９行目で指定するポイントシェープファイルと８行目のポリゴンシェープファイルでは座標系が違っています。

using DotSpatial.Data;
using DotSpatial.Topology;
using DotSpatial.Projections;

/* 中略 */

//読み取るシェープファイル
const string sPolygon = @"C:\hoge\Polygon.shp";
const string sPoint = @"C:\hoge\Point_4326.shp";

//判定結果はテキストファイルへ
const string sOut = @"C:\hoge\Result.txt";

//処理時間計測用
System.Diagnostics.Stopwatch swTime = new System.Diagnostics.Stopwatch();

//シェープファイルを開く・準備する
using(PolygonShapefile psPoly = new PolygonShapefile(sPolygon))
using(PointShapefile psPt = new PointShapefile(sPoint))
using(StreamWriter sw = new StreamWriter(sOut, false, Encoding.UTF8))
{
    //Prjファイル読み取り
    psPoly.ReadProjection();
    psPt.ReadProjection();
    
    //読み取った結果を表示
    Console.WriteLine("Poly :{0}", psPoly.Projection);
    Console.WriteLine("Point:{0}", psPt.Projection);

    //手法①
    //DotSpatial.Topology.Polygon.Contains(Point)
    swTime.Start();
    foreach (Feature fPoly in psPoly.Features)
    {
        Polygon pg = new Polygon(fPoly.Coordinates);
        foreach (Feature fPt in psPt.Features)
        {
            Point pt = new Point(fPt.Coordinates[0]);
            if (pg.Contains(pt)) //重なり判定
            {
                sw.WriteLine("Method1 Polygon:{0} Point:{1}",
                    fPoly.DataRow[0].ToString(),
                    fPt.DataRow[0].ToString());
            }
        }
    }
    Console.WriteLine("Method1:{0}", swTime.ElapsedMilliseconds);

    //手法②
    //DotSpatial.Topology.Polygon.IsWithinDistance(Point, double)
    swTime.Reset();
    swTime.Start();
    foreach (Feature fPoly in psPoly.Features)
    {
        Polygon pg = new Polygon(fPoly.Coordinates);
        foreach (Feature fPt in psPt.Features)
        {
            Point pt = new Point(fPt.Coordinates[0]);
            if (pg.IsWithinDistance(pt, 0)) //重なり判定
            {
                sw.WriteLine("Method2 Polygon:{0} Point:{1}",
                    fPoly.DataRow[0].ToString(),
                    fPt.DataRow[0].ToString());
            }
        }
    }
    Console.WriteLine("Method2:{0}", swTime.ElapsedMilliseconds);

    //手法③
    //DotSpatial.Data.Feature.IsWithinDistance(Feature, double)
    swTime.Reset();
    swTime.Start();
    foreach (Feature fPoly in psPoly.Features)
    {
        foreach (Feature fPt in psPt.Features)
        {
            if (fPoly.IsWithinDistance(fPt, 0)) //重なり判定
            {
                sw.WriteLine("Method3 Polygon:{0} Point:{1}",
                    fPoly.DataRow[0].ToString(),
                    fPt.DataRow[0].ToString());
            }
        }
    }
    Console.WriteLine("Method3:{0}", swTime.ElapsedMilliseconds);
}

そのため、コンソールの標準出力の１行目と２行目に表示されるProjectionの内容が異なっています。

ポリゴンはUTM座標系(JGD2000・53帯)なのに対してポイントは緯度経度(WGS84)と座標系も、座標の単位(度・ｍ)も異なっています。
この場合、「IsWithinDistance」メソッドの第二引数の距離ってどうするんだろうという問題が真っ先に思い浮かぶのですが、とりあえず実行してみました。

・・・実行してみました、がポリゴン内にあるポイントは無いという結果が帰ってきました。
(Result.txtが空)

ということで、DotSpatialでポリゴンとポイントの重なり関係の判定などを行う際には、測地系を揃えておかないといけないようです。

測地系はShapefile(PolygonShapefile・PointShapefile etc)のReprojectメソッドで変更することが出来ます(下記サンプルの11・12行目)。
引数のProjectionInfoクラスはDotSpatial.Projectionsにあります。

//シェープファイルを開く・準備する
using(PolygonShapefile psPoly = new PolygonShapefile(sPolygon))
using(PointShapefile psPt = new PointShapefile(sPoint))
using(StreamWriter sw = new StreamWriter(sOut, false, Encoding.UTF8))
{
    //Prjファイル読み取り
    psPoly.ReadProjection();
    psPt.ReadProjection();
　　
    //測地系を変換する
    psPoly.Reproject(ProjectionInfo.FromEpsgCode(3100));
    psPt.Reproject(ProjectionInfo.FromEpsgCode(3100));

    //後略

参考にさせて頂いたサイト

DotSpatial.Projections