GeoDistanceとその他の測地線距離算出式の精度
Mathematicaには2点の緯度と経度を与えて、その間の測地線距離を返す関数としてGeoDistanceがある。 しかしながら、ここで書かれているように、その精度には疑問が呈されているようだ。
他の有名な測地線距離の計算方法として、
- ヒュベニの式(カシミール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点間の距離の精度
地球表面上で一様分布となるよう、ランダムに2点を選び、その間の距離の精度を評価した。
上図は横軸がGeographicLibで計算した正確な距離、縦軸が各計算方法で計算した距離である。
y = xの直線上にあれば正確ということだが、GeoDistanceの"Robbin61"では1万kmを超えると実際よりも極端に短く計算してしまうようだ。
またヒュベニの式は実際よりも長く計算することが多く、そのような場合について調べてみると高緯度になるほど誤差が大きくなることが分かった。
同じ計算結果を誤差についてプロットしたものが下図である。また表は各計算方法の最大誤差である。
| 計算方法 | 誤差の絶対値 |
|---|---|
| GeoDistance(“Vincenty75”) | 7.119 * 10^-12 |
| GeoDistance(“Robbin61”) | 9.536 * 10^-1 |
| GeoDistance(“ExtendedNewtonRaphson”) | 1.247 * 10^-6 |
| ヒュベニの式 | 2.022 * 10^1 |
| 測量計算サイト計算式 | 6.682 * 10^-3 |
| 完全球体モデル | 6.410 * 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 * 10 ^ -12 | 500 | - |
| GeoDistance(“Robbin61”) | 2,538,452.677 040 755 | 1.261 * 10 ^ -3 | 672 | - |
| GeoDistance(“ExtendedNewtonRaphson”) | 2,541,656.533 178 208 | 1.991 * 10 ^ -7 | 516 | - |
| ヒュベニの式 | 2,541,657.194 841 664 | 6.121 * 10 ^ -8 | 152 | 0.36 |
| 測量計算サイト計算式 | 2,541,657.039 253 852 | 2.202 * 10 ^ -12 | 124 | 0.56 |
| 完全球体モデル | 2,544,345.085 547 101 | 1.058 * 10 ^ -3 | 12 | 0.32 |
この表から、ヒュベニの式は7桁程度の精度が出ているが、 Google Maps APIのDirections Serviceで得られる緯度と経度の精度は7〜8桁(小数点以下5桁)なので、 日本付近の緯度においてパスの長さを計算する用途では十分であることが分かった。
他の計算方法ではランダムな2点間の距離を計算する場合とほぼ同じ傾向となった。
計算時間の面では、GeoDistanceよりもヒュベニの式や測量計算サイトの計算式のほうが4〜5倍程度高速であり、完全球体モデルはさらに5倍程度高速であった。 ただし、ヒュベニの式、測量計算サイト計算式、完全球体モデルはCompileを使った場合、もとよりも40倍〜400倍程度高速となった。
結論
- MathematicaのGeoDistanceではデフォルトのVencenty75が11桁程度の精度が出て、最も良い。
- ヒュベニの式は日本付近の緯度であれば7桁程度の精度が出るが、極付近を含む場合は使えない。
- 完全球体として計算する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)