測地線距離計算式・計算ライブラリの精度評価

雑記帳

2地点の緯度経度を与えてその間の距離を求める計算式はいくつかあり、 GeoDistanceとその他の測地線距離算出式の精度 ではランダムな2点間の距離や、日本での運転経路データを用いて代表的な計算式の精度を評価したが、 短い距離の計算精度の評価と、対蹠点付近の計算精度の評価が不十分であった。

そこで、比較する計算式を追加した上で、GeographicLibのテストデータを用いて計算精度の再評価を行った。

比較する計算式

今回比較する計算式と計算ライブラリ、および出典を以下に示す。

計算に使用する回転楕円体のパラメータとしては、Haversine以外はいずれもWGS84の値($a = 6378137 \mathrm{m}$, $f = 1 / 298.257223563$)を使用し、Haversineでは平均半径($6371000.790009154 \mathrm{m}$)を使用した。

国土地理院計算式については再実装を行い、前回ゾーン2の計算で間違えていたのを修正した。

またGeoDistance(ExtendedNewtonRaphson)は計算が収束しないことがあるため、TimeConstrainedを使って0.1秒以内に計算が終了しない場合は強制終了して0.0とした。

GeographicLibは一括計算するプログラムをJavaで書き、それをMathematicaからJ/Link経由で呼び出した。 その他の計算式はMathematicaでCompileし、オプションとしてCompilationTarget -> "C", RuntimeAttributes -> {Listable}, Parallelization -> True]を指定したもので計算した。

その他、測地線航海算法(Geodesic Sailing)、 には「小野の公式」というものが見つかるが、2点が極めて近い or 対蹠点付近となる場合に0での除算が発生し、回避方法がなかったため今回の比較から除外した。

使用するテストデータ

精度検証に使用するデータはGeographicLibによって提供されているテストデータを使用した。

テストデータの仕様についてはGeographicLib: Geodesics on an ellipsoid of revolutionに記載されており、 ランダムな2点、対蹠点近辺、短距離、片方が極付近、両方が極付近、子午線沿い、卯酉線沿いなどの計算テストデータが含まれている。

テストデータにはGeodTest.dat($n=500000$)と、GeodTest-short.dat($n=10000$)があるが、

  • 500,000点をグラフに描画すると点が重なりすぎてプロットが見えなくなってしまう
  • GeoDistance(特にExtendedNewtonRaphson)の計算が非常に遅い

の2点の理由から、GeodTest-short.datを使ったデータを使用した。

計算結果

各々の計算式で、テストデータの正確な値との絶対誤差(m)と相対誤差を比較した。

両対数グラフであり、横軸が真の距離(m)、縦軸が絶対誤差(m)/相対誤差である。

絶対誤差 相対誤差

また、各計算式の絶対誤差と相対誤差の最大値、およびMathematicaとJavaScriptのそれぞれで計算にかかった時間を下表にまとめた。

計算式 絶対誤差(m) 相対誤差 計算時間(Mathematica) (ms) 計算時間(JavaScript) (ms)
GeographicLib $1.118 \times 10^{-8}$ $1.130 \times 10^{-10}$ $137.4$ $29.5$
Haversine $3.763 \times 10^4$ $5.600 \times 10^{-3}$ $4.0$ $1.5$
GSI $2.013 \times 10^{-4}$ $1.112 \times 10^{-10}$ $5.0$ $7.7$
Hubeny(Simple) $8.778 \times 10^6$ $8.433 \times 10^{-1}$ $3.7$ $0.9$
Hubeny(Full) $2.515 \times 10^6$ $2.071 \times 10^{-1}$ $3.7$ $1.6$
Vincenty $1.339 \times 10^5$ $6.705 \times 10^{-3}$ $18.7$ $36.6$
Lambert $2.464 \times 10^4$ $1.232 \times 10^{-3}$ $4.2$ $3.2$
Andoyer-Lambert $2.464 \times 10^4$ $1.232 \times 10^{-3}$ $3.9$ $2.9$
Andoyer-Lambert-Thomas $2.261 \times 10^5$ $1.130 \times 10^{-2}$ $3.8$ $1.5$
GeoDistance (Vincenty75) $1.336 \times 10^5$ $6.688 \times 10^{-3}$ $10,713.6$ -
GeoDistance (Robbin61) $2.000 \times 10^7$ $1.000 \times 10^0$ $14,187.7$ -
GeoDistance (ExtendNewtonRaphson) $2.000 \times 10^7$ $1.000 \times 10^0$ $22,268.9$ -

Mathematicaは呼び出しのオーバーヘッドが大きく、各計算式の計算量の比較にはJavaScriptの計算時間を使うのが良いと考えられる。

各計算式・ライブラリへの評価

  • GeographicLib
    • さすがGeographicLib。
    • 誤差15nmという謳い文句は本当で、どんな距離でも安定して誤差15nm以内の精度で計算される。
  • Haversine
    • 安定して0.6%以内に収まる。
    • 計算式も単純で実装しやすい。
  • GSI (国土地理院計算式)
    • 対蹠点を含め、安定して9桁〜10桁の精度。
    • 前の記事では「サイトでは異なる計算式が使われている可能性がある」としたが、サイトでもきちんと説明通りの計算が行われていると訂正する。
  • Hubeny(Simple)
    • 短距離なら精度が高く、100mまでなら7桁程度、10kmまでなら5桁程度の精度を出せる一方、100km以上になると精度が悪化して1桁程度になってしまう。
    • ルートデータやGPSログなど短い線分の集合に対しては、計算速度と精度のバランスが良く、カシミール3Dが採用したのも納得。
  • Hubeny(Full)
    • 簡易版ヒュベニの式と同様、短距離に強く、長距離になると精度が悪化する傾向は変わらないが、10km程度まで12桁の精度が出て、長距離になっても精度悪化が遅い。
    • MathematicaでCompileした場合、簡易式と計算速度がほとんど変わらなかったが、JavaScriptでは簡易式と比べて計算時間は2倍弱であった。これも速度と精度のバランスは良いと言える。
  • Vincenty
    • 対蹠点以外では11桁程度の精度が出る。対蹠点付近では2桁まで精度が落ちるので使えない。
    • 反復計算があるため計算速度は遅く、同じ反復計算を含むGSI計算式と比較すると、対蹠点以外での精度は1桁良い一方で速度は1/10〜1/3程度。
  • Lambert
    • 1m程度から10000km程度まで安定して5桁の精度だが、10000km以上では精度が落ち、特に対蹠点付近では2桁程度になってしまう。
  • Andoyer-Lambert
    • Lambertとほぼ同じだが、10m以下の精度が1桁程度悪くなっている。
  • Andoyer-Lambert-Thomas
    • Lambertに比べて、10000km以下では2桁程度精度が向上していて7桁程度あるが、やはり対蹠点付近ではLambertと同じく2桁弱になってしまう。
    • 速度も他のLambert系と比較して良く、ヒュベニの式のオリジナル版と同等の計算速度。
    • ヒュベニの式より広い範囲(10000km以下)で精度と速度のバランスが良い。
  • GeoDistance(Vincenty75)
    • 対蹠点以外では常に11桁の精度が出ている。対蹠点付近では2桁程度まで落ちるが、他のMethodより良い。
  • GeoDistance(Robbin61)
    • 10000km以内では2桁精度、10000km以上では使い物にならない。
  • GeoDistance(ExtendNewtonRaphson)
    • 対蹠点以外では常に5桁程度の精度。
    • 対蹠点付近ではそもそも計算が収束せず結果が返ってこないので使えない。

結論

  • 精度を求めるならGeographicLibを使うのが一番良い(各言語に移植されている)
  • ルートやGPSログなど、短距離の計算に限られる場合は、ヒュベニの式が精度と速度のバランスが良い(簡易版は100m以下、オリジナル版は10km以下)
  • 10km以上になる場合はAndoyer-Lambert-Thomasが速度と精度のバランスが良くなる
  • 速度よりも精度を優先し、かつ外部ライブラリを使わない(使えない)場合は、国土地理院の計算式が良い(ただ文書を見ても非常に実装しにくいのが難点)
  • MathematicaならMethodオプションはデフォルトのままが一番で、対蹠点以外では問題ない。もし対蹠点付近を含む計算をするならGeographicLibをJ/Linkで呼び出して使う(mathematica-geodesic)。

Appendix