Location::GeoToolの代わりに使えるProj.4

タイトルはオコガマシスで、本当はLocation::GeoToolの方がProj.4の代替でしかないのですが、Location::GeoToolを使っている人への呼びかけということなのでこのタイトルで。

先日、Web-APIを連発している某社のオフィスを訪ねたのですが、中の人の話によりますと、Perlで開発しているAPI内での経緯度測地系変換には拙作のLocation::GeoToolを使っていただいている、という話を伺いました。
いや、作った者として、使っていただけるのはありがたいのですが、昔昔に作ったモジュール構成を踏襲しているのでパフォーマンス的にも問題あるだろうし、また超エンタープライズな用途で変換結果に全幅の信頼を置かれてしまうのも中々辛いものがあるなと（いや、別に信頼を置かれてはマズイ問題は見つかっているわけではないのですが）、いうことで、間違いの許されない超エンタープライズな用途であるのならば、Proj.4のようなより専門的な投影変換ツールを使った方がよいのでは、という示唆をさせてもらいました。
ですが先方はProj.4をご存じなかったようなので、この記事で簡単に紹介させてもらいます。

Proj.4はオープンソースの地図投影変換ツール/ライブラリで、様々な地理的座標系間の相互変換機能を提供してくれます。
オープンソースのGISツールとして有名な、PostGISやMapServer等の投影変換部分のバックエンドとしても機能している、実績のある枯れたツール/ライブラリです。
利用法としてはコマンドラインで使うほか、Cライブラリとしても使うことができるのですが、とりあえずコマンドラインツールとしての使い方をご紹介。

インストール方法としては、Proj.4のサイトからソースをダウンロードして、

./configure
make
make install

するだけで超簡単。特に必要なライブラリ等もなく、すぐインストールできます。

以下、Location::GeoToolとの比較で、いろいろできることを列挙します。

１．世界測地系から日本測地系への変換

Location::GeoToolではこんな感じ（コマンドラインツールとの比較なので、コマンドラインでの実行で示します）。

    bash-3.1# perl -MLocation::GeoTool -e '
    > my $loc = Location::GeoTool->create_coord("35.00000","135.00000","wgs84","degree");
    > print join("\t",$loc->datum_tokyo->array)."\n";'
    34.9968036251296        135.002780769243

これを行おうと思うと、Proj.4では投影法間の相互変換を行うコマンドcs2csを使って、こんな感じで行います。

    bash-3.1# cs2cs +proj=latlong +ellps=WGS84 +to +proj=latlong +ellps=bessel
        +towgs84=-148,507,681 +no_defs -r -s -f "%.8f" <<EOF
    > 35.00000 135.00000
    > EOF
    34.99680329     135.00278102 -51.96734185

大体誤差範囲内の変換結果に収まって安心しました(汗)。
上記で、青で示した部分の設定が、変換元の世界測地系を表す投影法の設定であり、「+to」以降の赤字で表した部分が変換先の日本測地系を表す投影法の設定です。
日本測地系側の「+towgs84」で設定している楕円体重心のシフトパラメータは、以前私がWeb2.0ワークショップで提示した値（-146.336,506.832,680.254）と異なっていますが、この辺はどの値が最も正しいというのははっきり言って判らないところがあるので、今回はLocation::GeoToolで採用しているシフトパラメータを用いました。
その他のオプションは、-r、-sがそれぞれ入力・出力での座標値の並び順をY,Xの順（つまり緯度、経度の順）にする設定、-fは出力のフォーマット指定です。
cs2csの出力の第3値は高度ですが、世界測地系と日本測地系は楕円体が異なる上に原点も異なるので、高度ははっきり言って議論しても仕方ありません（日本測地系では高度0だけど世界測地系だと高度-51m、とか言われても困りますよね）。

２．日本測地系から世界測地系への変換

Location::GeoToolではこんな感じ。

    bash-3.1# perl -MLocation::GeoTool -e '
    > my $loc = Location::GeoTool->create_coord("35.00000","135.00000","tokyo","degree");
    > print join("\t",$loc->datum_wgs84->array)."\n";'
    35.0031965891345        134.997218898985

Proj.4ではやはりcs2csを使って、

    bash-3.1# cs2cs +proj=latlong +ellps=bessel +towgs84=-148,507,681 +no_defs
        +to +proj=latlong +ellps=WGS84 -r -s -f "%.8f" <<EOF
    > 35.00000 135.00000
    > EOF
    35.00319625     134.99721915 51.99721753

まあ、要するに１．の場合の逆、を行えばよいわけですね。

３．2点間の距離と方位角を求める

世界測地系で北緯35度、東経135度と北緯40度、東経140度間の距離と方位角を求める場合、Location::GeoToolだとこうなります。

    bash-3.1# perl -MLocation::GeoTool -e '
    > my $loc = Location::GeoTool->create_coord("35.00000","135.00000","wgs84","degree");
    > my $dir = $loc->direction_point("40.00000","140.0000","wgs84","degree");
    > print $dir->direction."\t".$dir->distance."\n";'
    37.047375828035 709253.729428108

で、結果は方位角37.047度（北北東）、距離709253mですね。

Proj.4では2点間の距離等の変換を行うgeodコマンドを用いて、

    bash-3.1# geod +ellps=WGS84 -I -p -f "%.8f" <<EOF
    > 35.00000 135.00000 40.00000 140.00000
    > EOF
    37.04737580     220.09531204    709253.728

のようにします。
入力の「+ellps」が計測するにあたって準拠する楕円体、-Iオプションは、このコマンドが本来1点と距離・方位角を与えて別の点を得るコマンドなので、その逆変換を指定するオプションです。
-pは、方位角の出力に負値を使わず、0-360度の値で出力させるオプション。
出力の第1値と第3値が、求める始点⇒終点方向の方位角、及び距離（m）ですね。
これも誤差範囲くらいの結果に終わってよかった（汗汗）。

ちなみに、出力の第2値は、終点から始点を見た場合の方位角です。
Location::GeoToolでこれを得たい場合には、次のようにすればOKです。

    bash-3.1# perl -MLocation::GeoTool -e '

    > my $loc = Location::GeoTool->create_coord("35.00000","135.00000","wgs84","degree");
    > my $dir = $loc->direction_point("40.00000","140.0000","wgs84","degree");
    > print $dir->reverse->direction."\n";'
    220.095312067574

４．1点からの距離と方位角から、その先の1点を求める。

世界測地系で北緯35度、東経135度から方位角45度（北東）、距離1000000mの地点の経緯度を求める場合、Location::GeoToolだとこうなります。

    bash-3.1# perl -MLocation::GeoTool -e '
    > my $loc = Location::GeoTool->create_coord("35.00000","135.00000","wgs84","degree");
    > my $dir = $loc->direction_vector(45.00000,1000000);
    > print join("\t",$dir->to_point->array)."\n";'
    41.0881639279743        143.411403618375

で、結果は北緯41.088度、東経143.411度です。

Proj.4だと、やはりgeodコマンドを-Iオプションなしで使って、

    bash-3.1# geod +ellps=WGS84 -p -f "%.8f" <<EOF
    > 35.00000 135.00000 45.00000 1000000
    > EOF
    41.08816395     143.41140363    230.19682191

出力の第3値は、やはり終点から始点を見た際の方位角です。

という感じで、Proj.4を使うと、Location::GeoToolと同等の変換を、より高い信頼性で行う事ができます。

とは言え、Location::GeoToolにはこれ以外にも、

• 経緯度を表す各種フォーマット間の相互変換が出来る
• 方位角をその方向を表す表現（例：北北東、NW、等等）に変換できる
• 2点間の中点等も扱える
• GridLocator、Locapoint、DoCoMoのｉエリア等との連携もできる

といった蓄積や利点もあるので、本当はLocation::GeoTool内部の変換エンジンを、Proj.4にアップデートできれば一番いいんですけどね。
でも、Proj.4ライブラリを叩くXSコードを書くのは私には敷居タカス...コマンドラインを叩いて実現する形なら、すぐにもできそうだけど...。