GeoDistanceとその他の測地線距離算出式の精度

Mathematicaには2点の緯度と経度を与えて、その間の測地線距離を返す関数としてGeoDistanceがある。 しかしながら、ここで書かれているように、その精度には疑問が呈されているようだ。

---

17 Mar. 2022 追記) 続編記事を投稿

---

他の有名な測地線距離の計算方法として、

• ヒュベニの式(カシミール3Dが採用しているが、英語圏では情報が見つからない)
• 国土地理院の測量計算サイトの計算式のついてのドキュメントを実装したもの
• 完全な球体とみなして計算(Google EarthやGoogle Maps API v3のcomputeDistanceBetweenが採用、haversine formula)

が見つかったので、GeoDistanceのMethodオプション3種類("Vincenty75", "Robbin61", "ExtendedNewtonRaphson")と上記の3種類の計算方法の合計6種類について精度を評価してみた。

正確な距離としては、GeographicLib (Wikipediaによると誤差15nmらしい)をJLink経由で使って計算したもの採用し、準拠楕円体はWGS84とした。

ただし、完全な球体とした場合の半径はGoogle EarthやGoogle Mapsにあわせて6378137 mとした。

特定の2点間の距離の精度

まずは北緯35度東経135度の地点から、北緯40度東経140度の地点までの距離を計算した。

計算方法	距離(m)	誤差
GeographicLib	$709,253.729\ 434\ 529\ 0$	-
GeoDistance(“Vincenty75”)	$709,253.729\ 435\ 369\ 6$	$5.465 \times 10^{-12}$
GeoDistance(“Robbin61”)	$708,476.302\ 141\ 391\ 7$	$1.096 \times 10^{-3}$
GeoDistance(“ExtendedNewtonRaphson”)	$709,253.691\ 659\ 883\ 5$	$5.327 \times 10^{-8}$
ヒュベニの式	$709,527.686\ 179\ 327\ 0$	$3.863 \times 10^{-4}$
測量計算サイト計算式	$709,253.729\ 433\ 037\ 8$	$2.103 \times 10^{-12}$
完全球体モデル	$710,213.967\ 093\ 551\ 9$	$1.354 \times 10^{-3}$

GeoDistanceのデフォルトVincenty75は他のMethodに比べて精度がよく、次いでExtendedNewtonRaphson、Robbin61の順だった。

ヒュベニの式は4桁程度であり、完全球体モデルから1桁程度しか精度が増えていなかった。

測量計算サイトの計算式は最も精度が良く、11桁の精度を持っていた。

ランダムな2点間の距離の精度

地球表面上で一様分布となるよう、ランダムに2点を選び、その間の距離の精度を評価した。

上図は横軸がGeographicLibで計算した正確な距離、縦軸が各計算方法で計算した距離である。

y = xの直線上にあれば正確ということだが、GeoDistanceのRobbin61では1万kmを超えると実際よりも極端に短く計算してしまうようだ。

またヒュベニの式は実際よりも長く計算することが多く、そのような場合について調べてみると測地線の緯度変化が一様でない(東京～ロンドンでは測地線が北極付近を通過するように、緯度が上がってから下がる場合など)に誤差が大きくなることが分かった。

同じ計算結果を誤差についてプロットしたものが下図である。また表は各計算方法の最大誤差である。

計算方法	誤差の絶対値
GeoDistance(“Vincenty75”)	$7.119 \times 10^{-12}$
GeoDistance(“Robbin61”)	$9.536 \times 10^{-1}$
GeoDistance(“ExtendedNewtonRaphson”)	$1.247 \times 10^{-6}$
ヒュベニの式	$2.022 \times 10^{1}$
測量計算サイト計算式	$6.682 \times 10^{-3}$
完全球体モデル	$6.410 \times 10^{-3}$

以上の結果からわかることをまとめると、

• MathematicaのGeoDistanceで最も精度が良いのはデフォルトのVincenty75で、常に11桁程度の精度を持っている。
• Robbin61は1万km以内でも2桁程度で、1万kmを超えると実用にならない(誤差95%)。
• ExtendedNewtonRaphsonは5桁程度の精度。
• 測地線の緯度変化が一様でない場合はヒュベニの式は精度が出ない(誤差2000%)。
• 測量計算サイトの計算式は1万kmまでは10桁程度の精度が出ていたが、1万kmを超えると2桁程度(誤差0.7%)まで一気に落ちる。
• 完全球体として計算する、Google EarthやGoogle Mapsの計算は2桁の精度(誤差0.6%)。

ただし、測量計算サイトの計算式については、測量計算サイトのフォームに入力して計算させると常にほぼ正しい結果が返ってくるので、ドキュメントと異なる計算が行われているのではないかと考えられる。 反復計算の収束閾値を変えてみても結果はほぼ変わらなかった。

日本付近におけるパスの長さを計算する場合の精度と計算時間

次に、もう少し実用的(ここでいう実用的とは専門家ではなく一般人にとって)な場合を考えて、 鹿児島県の枕崎から北海道の稚内までの運転経路の道のりの長さを計算させてみた。

Google Mapsで枕崎駅から稚内駅まで一般道路のみの経路を計算させ、その結果を400点に削減したものを用いた。 使用したデータはmakurazaki_wakkanai_reduced.gpxである。

計算方法	計算された距離(m)	誤差	計算時間 (ms)	計算時間(Compile使用時) (ms)
GeographicLib	$2,541,657.039\ 259\ 448$	-	$212$ (JLink使用)	-
GeoDistance(“Vincenty75”)	$2,541,657.039\ 264\ 922$	$2.154 \times 10^{-12}$	$500$	-
GeoDistance(“Robbin61”)	$2,538,452.677\ 040\ 755$	$1.261 \times 10^{-3 }$	$672$	-
GeoDistance(“ExtendedNewtonRaphson”)	$2,541,656.533\ 178\ 208$	$1.991 \times 10^{-7 }$	$516$	-
ヒュベニの式	$2,541,657.194\ 841\ 664$	$6.121 \times 10^{-8 }$	$152$	$0.36$
測量計算サイト計算式	$2,541,657.039\ 253\ 852$	$2.202 \times 10^{-12}$	$124$	$0.56$
完全球体モデル	$2,544,345.085\ 547\ 101$	$1.058 \times 10^{-3 }$	$12$	$0.32$

この表から、ヒュベニの式は7桁程度の精度が出ており、特定の2点間の距離を計算した場合に比べて精度が良かった。 これはヒュベニの式は比較的近い2地点においては精度が高くなることを示しており、 折れ線パスの長さを計算する場合は比較的精度が良くなることを示している。

Google Maps APIのDirections Serviceで得られる緯度と経度の精度は7〜8桁(小数点以下5桁)なので、 日本付近の緯度においてパスの長さを計算する用途では十分であることが分かった。

他の計算方法ではランダムな2点間の距離を計算する場合とほぼ同じ傾向となった。

計算時間の面では、GeoDistanceよりもヒュベニの式や測量計算サイトの計算式のほうが4〜5倍程度高速であり、完全球体モデルはさらに5倍程度高速であった。 ただし、ヒュベニの式、測量計算サイト計算式、完全球体モデルはCompileを使った場合、もとよりも40倍〜400倍程度高速となった。

結論

• MathematicaのGeoDistanceではデフォルトのVencenty75が11桁程度の精度が出て、最も良い。
• ヒュベニの式は日本付近の緯度で、パスの長さの計算など2地点間が近い場合であれば7桁程度の精度が出るが、極付近を含む場合や2地点が離れた場合は精度が悪化する。
• 完全球体として計算するGoogle EarthやGoogle Mapsは2桁程度の精度(誤差0.6%程度)。
• 測量計算サイトはドキュメントと異なる計算が行われていると考えられる。 -> 続編記事にて修正

追記

• GeographicLibを参考にMathematicaで測地線距離を計算するパッケージを作ってみた。→Google Code: mathematica-geodesic
• Google Codeが終了するのでGitHubに移転→mathematica-geodesic (Mar. 13 2015)
• いつのまにか本家GeographicLibにMathematicaでの実装としてリンクされている(!)
• 国土地理院のサイトのURLが変更されたので、リンク先を変更 (Jan. 17 2016)
• 続編記事を投稿。一部の結論が変更。 (Mar. 17 2022)