二点間測地線距離（地球上）：各計算方法の誤差評価 #JavaScript

こんにちは。
回転楕円体面上（地球上）の二点間測地線距離 (geodesic distance) を計算する各方法の最大相対誤差 (maximum relative accuracy) の距離依存性を調べ比較しました（評価結果は下記チャート）。

方法を選定したい方のためには、簡素な計算のうちでは精度が良い Gauss-mid-latitude もしくは Andoyer あたりが候補だろうと思いました（下記のコメント）。
詳細評価のための処理および各計算法は「測地線距離計算式・計算ライブラリの精度評価」（330K INFO）からそっくり拝借して使いました（ありがとうございます）。
なお実装の際の注意点は haversine 型の距離計算を精度上使うことが重要です¹。

入力データ

入力データは楕円体面上に稠密に二点を発生させ（15万対）、その際 GeographicLib（の GeodSolve）を利用し距離を与えています。距離 10,000 km 以上では、"GeodTest.dat" (GeographicLib) も使いました。

各計算法²

Spherical: Menelaus of Alexandria (about 100); Great-circle distance - Wikipedia; 完全な球面近似であり回転楕円体面を考慮しない形; 通称 haversine
Flat-surface: Newton, Cassini (1713); 起源探究は断念しました; Geographical distance - Wikipedia; 一次近似（平面近似）; 日本での通称は Hubeny
Gauss-mid-latitude: Gauss (1843, 1846); Lambert and Swick (1935); Rapp (1991) §6.4; Geographical distance - Wikipedia; 回転楕円体面を考慮した球面近似の形; 「回転楕円体面に沿う最短距離の式」 - wikipedia
Hubeny (1954, 1959)
Andoyer: Forsynth (1895), Andoyer (1932)
Andoyer-Lambert: Thomas (1965)³
Andoyer-Lambert-Thomas: Thomas (1970)³
(Pittman: Pittman (1986); Deakin and Hunter (2007); not tested)
Bowring: Bowring (1996); used by Geospatial Information Authority of Japan (GSI) based on this method, 距離と方位角の計算　計算式（国土地理院）
Vincenty: Vincentry (1975, 1976); 回転楕円体上の測地線及び航程線の算出について（海洋情報部研究報告） - 国立国会図書館デジタルコレクション
Karney: Karney (2011); GeographicLib

なお今回私自身で書き下して用いたソース（下記）はこれらのうち spherical, flatsurface, gaussmidlatitude, andoyer の4つだけです（どれも簡素な計算方法）。

今回の評価結果に対するコメントおよび注意点は（厳密ではありませんが）、

spherical は最大相対誤差 0.5%。距離依存性がほぼなく、対蹠点まで安定している。
なおどの式も誤差は緯度依存性も持つが今回の評価結果には表現されない（例えば sphericalを除き、赤道上の二点間に対しては式の上では誤差が無くなる）。
flatsurface と hubeny は経緯度変数による局地点展開形に由来し、特異点を持つことが要注意です。使用条件を誤らなければ近似精度は良いです。
- 誤差の緯度依存性が大きい。
- もしも測地線が極上を通る場合の最大相対誤差は、距離には依存せず定数 $\frac{\pi}{2}-1$ に達する（この最大誤差はチャート上に示していない）。
- 180度経線を跨ぐ際にも特異性を持つが、対処は容易です（下記コード）。
gaussmidlatitude の誤差は 1,000 km で 8 ppm、10,000 km で 0.3 %（spherical に逆転される）。10,000 km 以上は対蹠点まで安定している。
andoyer⁴ の誤差は、10,000 km 以下では常に 15 ppm。ただし 10,000 km 以上で急激に増大する（対蹠点で0.3 %）。
andoyerlambert の誤差は、10,000 km 以下では常に 1.5 ppm⁵。ただし 10,000 km 以上で急激に増大する（対蹠点で0.3 %）。
andoyerlambertthomas の誤差は、10,000 km 以下では常に 3 ppb。ただし 10,000 km 以上で急激に増大する（対蹠点で0.3 %）。
hubeny の誤差は、使用条件を誤らなければ、絶対誤差が 10 nm 以下（karney と同等）。
vencenty と bowring の誤差は安定して小さいが、対蹠点直近で誤差が急増する。計算速度はそれほど速くない。
bowring は実装が容易ではなく、加えて NaN となる計算結果の回避や誤差急増対策が不明。
- これらへ対処可能であれば有望と思われる。
karney の誤差は、今回の検討では絶対誤差が常に 10 nm 以下。対蹠点まで安定している。計算速度は速くない。

総合すると、Gauss-mid-latitude が、簡素な計算法を求めるが多少近似精度を気にする方の常用のための一候補に見えます（spherical や flatsurface に替えて）。ただし 100 km 以上での誤差は 1 ppm を超え増大していきます。これを気にする方は、Andoyer あたりが次の候補になるかと思います。

評価結果

ソース

spherical.js

function spherical(lat1, lon1, lat2, lon2){
  "use strict";
  const a_adjusted = 6367000;  // Earth radius
  const degree = Math.PI / 180;
  const sin = Math.sin;
  const cos = Math.cos;
  const asin = Math.asin;
  const sqrt = Math.sqrt;

  const latdiffhalf = (lat1 - lat2) / 2 * degree;
  const londiffhalf = (lon1 - lon2) / 2 * degree;
  const lat = (lat1 + lat2) / 2 * degree; // middle

  return 2 * a_adjusted * asin(hypot(  // arc_half = asin(chord_half)
    sin(londiffhalf) * cos(lat),
    cos(londiffhalf) * sin(latdiffhalf)
  ));
}

flatsurface.js

function flatsurface(lat1, lon1, lat2, lon2){
  "use strict";
  const a = 6378137;  // GRS80
  const f = 1 / 298.257223563;  // WGS84
  const e2 = f * (2 - f);
  const degree = Math.PI / 180;
  const abs = Math.abs;
  const sin = Math.sin;
  const cos = Math.cos;
  const sqrt = Math.sqrt;
  const hypot = Math.hypot;
  
  const londiff = abs(lat1) == 90 || abs(lat2) == 90 ? 0 :  // at the north and south poles
                 (((lon1 - lon2 + 180) % 360) - 180) * degree;  // across the 180th meridian
  const latdiff = (lat1 - lat2) * degree;
  const lat = (lat1 + lat2) / 2 * degree; // middle
  const sinlat = sin(lat);
  const coslat = cos(lat);
  const n2 = 1 / (1 - e2 * sinlat ** 2);
  const n = sqrt(n2);  // prime vertical radius of curvature
  const m_by_n = (1 - e2) * n2;  // meridian ratio

  return a * n * hypot(
    londiff * coslat, 
    latdiff * m_by_n
  );
}

gaussmidlatitude.js

function gaussmidlatitude(lat1, lon1, lat2, lon2){
  "use strict";
  const a = 6378137;  // GRS80
  const f = 1 / 298.257223563;  // WGS84
  const e2 = f * (2 - f);
  const degree = Math.PI / 180;
  const sin = Math.sin;
  const cos = Math.cos;
  const asin = Math.asin;
  const sqrt = Math.sqrt;
  const hypot = Math.hypot;

  const latdiffhalf = (lat1 - lat2) / 2 * degree;
  const londiffhalf = (lon1 - lon2) / 2 * degree;
  const lat = (lat1 + lat2) / 2 * degree; // middle
  const sinlat = sin(lat);
  const coslat = cos(lat);
  const n2 = 1 / (1 - e2 * sinlat ** 2);
  const n = sqrt(n2);  // prime vertical radius of curvature
  const m_by_n = (1 - e2) * n2;  // meridian ratio	

  return 2 * n * a * asin(hypot(  // arc_half = asin(chord_half)
    sin(londiffhalf) * coslat,
    cos(londiffhalf) * sin(latdiffhalf * m_by_n)
  ));
}

andoyer.js

function andoyer(lat1, lon1, lat2, lon2){
  "use strict";
  const a = 6378137;  // GRS80
  const f = 1 / 298.257223563;  // WGS84
  const degree = Math.PI / 180;
  const sin = Math.sin;
  const cos = Math.cos;
  const asin = Math.asin;
  const hypot = Math.hypot;
  
  const lat = (lat1 + lat2) / 2 * degree; // middle
  const latdiffhalf = (lat1 - lat2) / 2 * degree;
  const londiffhalf = (lon1 - lon2) / 2 * degree;

  const sin_lat = sin(lat);
  const cos_lat = cos(lat);
  const sin_latdiffhalf = sin(latdiffhalf);
  const cos_latdiffhalf = cos(latdiffhalf);
  const d = 2 * asin(hypot(  // 0-th order distance; arc_half = asin(chord_half)
    sin(londiffhalf) * cos_lat,
    cos(londiffhalf) * sin_latdiffhalf
  ));
  const c = (3 * sin(d) - d) * (sin_lat * cos_latdiffhalf / cos(d / 2)) ** 2;
  const s = (3 * sin(d) + d) * (cos_lat * sin_latdiffhalf / sin(d / 2)) ** 2;
  const delta_1st = f * (c - s) / 2;  // 1st order
  return a * (d + delta_1st);
}

もしも Spherical law of cosines（余弦の球面法則）を用いると短距離での計算精度が劣化します。世の中の実装をいくつか調べてみますと、少なからず見られました（例えば、Boost.Geometry、Forsynth, Andoyer, Lambert, Thomas 系実装）。 ↩
なお文献一覧は、"A geodesic bibliography" (GeographicLib) もよくまとまっています。 ↩
この Thomas (1965) は、Forsynth (1895), Andoyer (1932); (Marie Henri Andoyer - Wikipedia), Lambert: Lambert (1942); (Geographical distance - Wikipedia) の先行研究に基づき、一次近似を再整理したもの。Thomas (1970) は、二次近似を正しく求めたもの。 ↩ ↩²
Boost.Geometry がデフォルト計算法として採用。 ↩
parametric latitude 使用による精度向上の反面、計算速度が若干低下します。 ↩

二点間測地線距離（地球上）：各計算方法の誤差評価

入力データ

各計算法2

コメント

評価結果

ソース

各計算法²