[GDAL/OGR]座標系を変換する (CoordinateTransformation)

2019/7/10 C#, gdal, GIS, OGR, プログラミング

OGRを使って座標変換(緯度経度から平面直角座標など)する方法について、備忘録を兼ねて。

今まで、OGRを使って座標変換をするときに、今まではポイントのジオメトリを作った後に、TransformToメソッドを使って座標変換をしていました。

だいたい、こんな感じのコードです。

OSGeo.OGR.Geometry geom = new OSGeo.OGR.Geometry(OSGeo.OGR.wkbGeometryType.wkbPoint);
geom.AssignSpatialReference(rfSrc); //rfSrc: 変換前の測地系(OSGeo.OSR.SpatialReference)
geom.AddPoint(X, Y, 0d);
geom.TransformTo(rfDest); //rfDest: 変換先の測地系(OSGeo.OSR.SpatialReference)
X = geom.GetX(0);
Y = geom.GetY(0);

変換した結果をGISデータとして残したい場合にはこの方法なんですけど、座標だけを変換したい場合には処理が猛烈に遅くなります。例えばラスタデータの全てのピクセルの位置を変換するような場合には、無駄に時間を過ごすことになります。

何か良い方法は無いのかとOSGeo.OSR名前空間の中を眺めたら「CoordinateTransformation」といういかにもなクラスがあったので試してみました。

CoordinateTransformationクラスのコンストラクタは、「CoordinateTransformation(SpatialReference src, SpatialReference dst)」となっていて、第一引数が変換元(変換前)の座標系定義、第二引数が変換先(変換後)の座標系定義となっています。
ちなみに、座標定義(SpatialReferenceクラス)のコンストラクタの引数は、OGC_WKTフォーマットで、SpatialReferenceなどでEPSGから調べることが出来ます。

早速、多数地点の座標系を変換してみました。
コードはこんな感じです。

//冒頭にこれが書かれている前提です
using System.IO;
using OSGeo.OSR;

//座標変換結果を出力
const string sOut = "Pos.csv";

//変換元:WGS84(EPSG:4326)、変換先:JGD2000のUTM・54帯(EPSG:3100)
const string sSrc = Osr.SRS_WKT_WGS84; //メジャーなのは定義がある
const string sDest = "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]]";

//GDALのインストール先
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(); //ここではOGRは使ってない
    #endregion

    //変換前後の測地系定義
    SpatialReference rfSrc = new SpatialReference(sSrc);
    SpatialReference rfDest = new SpatialReference(sDest);

    //CoordinateTransformation→OSGeo.OSR名前空間
    CoordinateTransformation Trans =  new CoordinateTransformation(rfSrc, rfDest);

    //緯度経度から平面直角座標(UTM・54帯、JGD2000)
    using (StreamWriter sw = new StreamWriter(sOut, false, Encoding.ASCII))
    {
        sw.WriteLine("Lon,Lat,X,Y");
        for (double dLon = 138; dLon < 141.1d; dLon += 0.1d)
        {
            for (double dLat = 35; dLat < 38.1d; dLat += 0.1d)
            {
                double[] dPos = new double[] { dLon, dLat, 0d };
                Trans.TransformPoint(dPos); //TransformPoint(double[] inout);

                //dPos:の値が変換後(UTM・54帯)の座標に更新されている
                sw.WriteLine("{0:f1},{1:f1},{2:f1},{3:f1}", dLon, dLat, dPos[0], dPos[1]);
            }
        }
    }
}

いろいろと前振りなどが長くなっていますが、コードの本体は33行目と44行目。
33行目でCoordinateTransformationクラスのインスタンス「Trans」を作成して、44行目でTransを使って座標系を変換しています。

気になる変換の誤差は、ほぼゼロでした。
ただし、WGS84からJGD2000なので地球楕円体の形自体がほとんど変わっていません。そのため、地球楕円体の形そのものが変わる場合などには、誤差が出るかも知れませんのでお気をつけください。