東京タワーとエッフェル塔の距離
PostGISをインストールした（詳細は“PostGISのインストール”）。目的は地球上の2点間の距離と方向を計算できるプログラムが欲しかったから。で、早速遊んでみた。
“東京タワーとエッフェル塔はどれくらい離れているか？”で試してみた。
まず、東京タワーとエッフェル塔の経度・緯度が必要だ。（余談だが、日本では“緯度・経度”の方が一般的だが、PostGISとか英語のページとか見ていると“経度・緯度”が一般的のようだ。）便利なページもあるもので Geocoding というページで“東京タワー”とか“エッフェル塔”と入れると出てくる。検索結果は以下の通りとなっている。
東京タワー 座標(WGS84): 緯度 35度39分30.996秒(35.65861), 経度 139度44分43.609秒(139.745447) エッフェル塔 座標(WGS84): 緯度 48度51分29.538秒(48.858205), 経度 2度17分39.692秒(2.294359)
早速、この値を使ってPostGISに計算させてみる。
adsaria@postgis:~$ psql -U adsaria -d gistestdb -h localhost
Password for user adsaria:
psql (8.4.2)
SSL connection (cipher: DHE-RSA-AES256-SHA, bits: 256)
Type "help" for help.
gistestdb=> SELECT
ST_Distance_Spheroid(
ST_GeomFromText('POINT(139.745447 35.65861)'),
ST_GeomFromText('POINT(2.294359 48.858205)'),
'SPHEROID["GRS 80", 6378137, 298.257222101]'
) ;
st_distance_spheroid
----------------------
9743082.98115865
(1 row)
gistestdb=>約9,743km。
では、国土地理院のホームページではどうなっているか。“距離と方位角の計算”のページで計算してみる。
約9,743kmで同じ。PostGISと国土地理院のページでは“1mm”の違いである。
では、Google Earthでは？ Google Earthの“定規”を使って測る。画面の下にマウスのある緯度・経度が表示されるので、上記の緯度・経度に合わせて測定すると、
おや、11kmも違う。約１万キロに対して10キロだから1/1000程度の話なのだが。例えて言えば、1m先のコップの位置が1mm遠いという話と同じレベルなのだけど。
Google Earthが使っている地球楕円体の定義と国土地理院が使っている地球楕円体の定義は異なっている。Google EarthはWGS 84という定義を使っているぽい。国土地理院のページはGRS 80とBesselという楕円体を使って計算できるが、上の例はGRS 80を使っている。WGS 84とGRS 80は扁平率が僅かに異なる。とはいっても有効数字で8桁は同じで9桁目で違う程度で1億分の一程度の違いしかない。
PostGISでは任意の楕円体を使って計算できる。先ほどの例は国土地理院と同じGRS 80を使って計算したが念のためWGS 84を使って計算すると以下のようになる。
gistestdb=> SELECT
ST_Distance_Spheroid(
ST_GeomFromText('POINT(139.745447 35.65861)'),
ST_GeomFromText('POINT(2.294359 48.858205)'),
'SPHEROID["WGS 84", 6378137, 298.257223563]'
) ;
st_distance_spheroid
----------------------
9743082.98109636
(1 row)
gistestdb=>WGS84とGRS 80では1万キロで約0.06mmの違いしかない（ただし、測地線長は2点間の相対的な位置関係で変わってくるので同じ1万キロでも常に誤差が同じではない）。12桁目で違いが出てくる程度であり、演算誤差による違いと言っても良いくらいの話だ（“浮動小数点の精度の確認”参照）。ということはGoogle Earthの“定規”は何か別の理由で誤差が出ていることになる。ただ、11kmと随分他のソフトウェアと異なる結果を出しているので単なるアルゴリズムの違いとは思えない。東京・パリ間を高度10000m（＝約10km）で飛行した場合の“実際の”航跡なのだろうか？ PostGIS等で使っている2点間の距離の計算方法（アルゴリズム）の原型は1960年頃には考案されているようでFORTRANで実装されてる（“SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY”、“PC Software Download - INVERSE and FORWARD”）。従ってGoogleだけが大きく誤差があるというのはちょっと奇妙な気がする。
ところで、実は私も自分でプログラムを組んで計算してみたのだが、私のプログラムでは“10153958.9492”mであり、410kmも違う。“公共測量教程 測量計算 改訂版”という本を参考にプログラムを組んでみた。測地線を計算する単純な式をそのままアルゴリズム化しただけなので、2点間の距離が近い時は良い精度が出るが、距離が離れると精度が悪くなる。（私のプログラムでも数十キロの距離であれば国土地理院のページと同じ結果になる。）距離が離れても精度を保つためには、微分値の最小を反復計算をして求める等と色々な方があるらしい。（余談だが、“公共測量教程 測量計算 改訂版”は余りお勧めできない。三角関数の最初の定義のところで、cosA = c/b というコサインの定義の2行下に sinA=cos(90-A)=c/b と書いてある。cosA = sinA ということになってしまう。こうなってくると例題で答えが合わない場合に“本に書かれた数式が本当に正しいのか？”と懐疑的になってしまう。改訂版であるにも関わらず、だ。それに説明無しで専門用語が出てきたりするので理解するのに時間がかかってしまった。）
さて、PostGISの計算で2点間の距離を十分正確に計算できることが分かった。さて、では2点間の方位は？と思ってPostGISのマニュアルを見てみてたが、無い！ 2次元平面上の方位角を求める azimuth（アジマス）という関数はあるのだが楕円体上の方位角を求める関数が見当たらない。上の例にもあるように、2点間の距離と同時に方位角も重要なファクターであり、国土地理院のホームページでもGoogle Earthでも出てくる。PostGISのソースコードも眺めてみたが、lwgeom_spheroid.cの中には距離を求める関数distance_ellipse_calculationはあるが、方位角を求める関数がない。GISのアプリケーションを使う人達って、楕円体（もしくは球体）上の方位角って使わないのかな。
仕方ないので、他のソフトウェアが無いか調べてみた。
次にたどり着いたのは、“Geo-Ellipsoid-1.12”というPerlで書かれたクラスライブラリ。簡単なテストプログラムを作って走らせると、次のような感じ。
adsaria@postgis:~/Geo-Ellipsoid-1.12/test$ ./mytest.pl
Enter test_ellipsoid
Using Geo::Ellipsoid version 1.12
WGS84 ellipsoid values:
Equatorial radius = 6378137.0000000000
Polar radius = 6356752.3142451793
Flattening = 0.0033528106647475
Eccentricity = 0.0818191908426215
Reciprocal flattening = 298.25722356
Here = [35.658610000000,139.745446944444]
There = [48.858205000000,2.294358888889]
Range = 9743082.9823 m., bearing = 333.5616 deg.
East = -4337958.2060 m., north = 8724092.1936 m.
There2 = [1.303168565570,227.689275034013]
adsaria@postgis:~/Geo-Ellipsoid-1.12/test$しかし、私はPerlは使ったことなし、どうしよう。Cのライブラリあたりがあれば良いのだが。PostGISのソースか、Geo-EllipsoidのPerlスクリプトを参考に自前のライブラリを作るしかないか、と考え始めていたところ、“ひょっとしてPostGISの動作環境としてインストールしたGEOSかProj4にそういった機能はないのか”と思い調べてみた。
どうもGEOSには方位角まで求めるクラスは見当たらなかったが、Proj4にはありました。Linuxのコマンドとして。invgeodというコマンドがまさにそれ。
adsaria@postgis:~$ invgeod +ellps=WGS84 <<EOF 35.65861 139.745447 48.858205 2.294359 EOF -26d26'18.071" 33d19'38.911" 9743082.958 adsaria@postgis:~$
灯台下暗し。
ちなみに、始点と終点を与えてその2点間の距離と方位角を求めるのが“逆問題（inverse problem）”と呼ばれ、始点に方位角と距離（測地線長）を与えて終点を求めるのが直接問題（forward problem）”という言うらしい。なので“invgeod”。その逆（つまり直接問題を解く方）は“geod”で、実行すると、
adsaria@postgis:~$ geod +ellps=WGS84 <<EOF 35.65861 139.745447 -26d26'18.071" 9743082.958 EOF 48d51'29.539"N 2d17'39.694"E 33d19'38.912" adsaria@postgis:~$
となる。（東京タワーの緯度・経度から方位角と測地線長を与えて、エッフェル塔の緯度・経度を求めている。）
Proj4にはlibproj-devという開発用のパッケージがあるので、そこにgeod/invgeodの関数があるかなと期待したのだけど、PROJ.4の“PROJ.4 - geod program”にあるように、geodのアルゴリズムはライブラリには含まれていないということで、少々がっかりした。もともとPROJ.4は地球楕円体と2次元の地図平面との座標変換システムとして作られているので、測地計算の機能はあくまで補助的な位置づけなのだろう。また、“Gerald writes that geod is based upon a poorer algorithm”とあるように、余り優秀なアルゴリズムではない、というのも理由のようだ。しかし、国土地理院の計算結果と1万kmで22cmしか違わないのであれば、十分だと思う。
当面はCからfork/execでgeodを呼び出して使えばいいかな。
こうして幾つかのソフトウェアで測地計算をすると違いが出てくる。Google Earthでは11kmもの誤差がある。しかし、国土地理院の方がより正確でGoogleの方が誤差が大きいとは（私には）言いきれない。こうなってくると、どれが正しいのか分からない。あとは実用の範囲で問題なければ良いということになってしまう。私のプログラムは4倍精度の浮動小数で書いてみたが、所詮、そこまで計算精度を高めても元のアルゴリズムが悪ければ意味がない。それに（無重力中の“水の塊”のように）地球自身も絶えず変形している筈だし、大陸プレートだって動いている。それを考えれば倍精度の10桁精度でも十分過ぎる位だし、アプリケーションによっては単精度の有効桁数7桁でも十分かもしれない。（1万キロで10メートル程度の誤差。ただし実際には反復計算で誤差の積み上がりがあるので精度はかなり落ちるが。今度暇な時に測地計算で単精度と倍精度の比較をしてみよう。）
なんかPostGISまでインストールする必要はなかったが、まぁ、PostGISはPostGISで色々と使えそうだから良いかな。
