近年のスマートフォン普及の勢いは目覚ましく、日本においてはほとんどすべての人がこれを持ち歩く社会となりました。 そのため、私が開発・改修を手掛けるシステムの中にも、端末が備えるGPSモジュールによって取得された利用者の位置情報を利用するものが多々出てくるようになっています。 位置情報の利用方法は様々ですが、その中で最も大きな割合を占めるのが、システムに登録された施設 (店舗やATM等) と利用者の間の距離を計算して近いものを提案するというもの。 単純な直線距離ではなく、道路や線路などに沿った実際の移動距離を求める場合もありますが、その場合も連続する線分 (折れ線) で近似された経路の長さを計算することになるので、結局は二点間の直線距離の計算に行き着きます。

2地点の座標 (λ_A, φ_A), (λ_B, φ_B) (λ は緯度、φ は経度) からこれらの間の距離を求める計算はそれほど難しいものではありません。 まず、地球を、原点 o を中心とする半径1の球体 と仮定します。 北極を (0, 0, 1)、ヌル島を (1, 0, 0) とすれば、A, B の座標はそれぞれ、

となります。 これらの位置ベクトル OA (x_A, y_A, z_A) と OB (x_B, y_B, z_B) が成す角の余弦値 cos χ は内積として

と計算されるので、χ そのものの値は逆余弦 cos^-1 を用いて、

と求めることができます。

2点ABを結ぶ単位球面上の最短距離すなわち測地線の長さは、全周 2π に対する角度 χ の比 χ/2π でなので、これを地球の半径 R 倍だけ拡大してやれば、実際の二点間の距離 δ が得られます。

以上のプロセスを JavaScript で書き下したならば、以下のようになるでしょうか:

//地球半径 (単位: m)
const EARTH_RADIUS = 6378100;

//2点 (srcLat, srcLong), (dstLat, dstLong) 間の距離を計算
function calcDistance(srcLat, srcLong, dstLat, dstLong){

    var srcX = Math.cos(srcLat)*Math.cos(srcLong);
    var srcY = Math.cos(srcLat)*Math.sin(srcLong);
    var srcZ = Math.sin(srcLat);

    var dstX = Math.cos(dstLat)*Math.cos(dstLong);
    var dstY = Math.cos(dstLat)*Math.sin(dstLong);
    var dstZ = Math.sin(dstLat);

    return EARTH_RADIUS*Math.acos(srcX*dstX + srcY*dstY + srcZ*dstZ);
}

ただし、関数 calcDistance の引数 srcLat, srcLong, dstLat, dstLong の単位は一般的に用いられる度ではなく、ラジアンであることに注意してください。

これはプログラミングにおいてはごく初歩的なトピックであるため、学生ならばともかく、私のような働き始めて十ン年というソフトウェア技術者が解説記事を書くような内容だとは思っていませんでした。 仮に職業プログラマを名乗っていながらこの程度の計算ができないならば、あまりにも程度が低く、はっきり言ってプロとして金を取れるレベルに達していません。

......と言いたいところなのですが、実はこれまでに改修を手掛けてきたシステムのうち、この「距離の計算」を正しく行えているものは半数に満たないというのが実際のところ。 これまでに見た中で特にひどかったのは以下のような実装です:

const DISTANCE_PER_DEG = 111318;

function calcDistance(srcLat, srcLong, dstLat, dstLong){

    var diffLat  = DISTANCE_PER_DEG *(dstLat  - srcLat);
    var diffLong = DISTANCE_PER_DEG*(dstLong - srcLong);

    return Math.hypot(diffLat, diffLong);
}

定数名 DISTANEC_PER_DEG (= "1度あたりの距離") を見た瞬間に「これは駄目だ」と分かるシロモノ。 このプログラムは地球表面を平面と見なしてユークリッド距離の計算を適用していますが、これは正距円筒図法で描かれた地図の上に定規を当てて距離を測っているのと同じ。

ご存知のように、経度1度あたりの距離は赤道上で最も長く、緯度が上がる (/下がる) につれて短くなっていき、極点においては 0 になります。 DISTANEC_PER_DEG の定義値 111,318 (単位はm) は赤道上のものであり、これを用いた計算はインドネシア, ケニア, エクアドルといった赤道直下の国でしか正しい答を返しません。 北緯36度前後に位置する日本 (東京) ではこの長さは赤道上のそれの約80.9%となり、その差は1度あたり 21,259m にもなります。 このズレは実用的にも許容しがたい大きさであり、こんな計算をするプログラムは「欠陥品」と言われても仕方がないでしょう。 DISTANEC_PER_DEG の値を 90,059 に変更してやれば誤差はかなり小さくなりますが、それも北緯36度 (あるいは南緯36度) 付近に限った話であり、他の地域にこれを適用すればやはり致命的なズレが生じます。

これほどまでに杜撰な計算による問題が表面化しなかったのは、このプログラムが「利用者と複数の登録地点間の "距離" を計算し、その中で最も近くにあるものを提示する」という用途に使われていたため。 意味があるのは具体的な距離ではなく、その大小関係であるため、複数回行われる距離計算が一様にズレるならば、とりあえずは "正しい" 動作になる、という具合です。

私が手掛けてきたシステムの殆どがこれと同様にそれほどの精度を要求されない性質のプログラムであったことは事実で、おそらくは、カーナビや今話題の接触確認アプリといった位置情報の精度がクリティカルなプログラムではきちんとした計算がなされているのでしょう。 (そうだと思いたい。) しかし、それは正しくない計算をしていい理由にはなりません。

これほど滅茶苦茶な計算手法があちこちで使われている理由はいったい何なのかと考えたとき、ひとつ思い当たることがあります。 「緯度経度 距離 計算」というキーワードで Google 検索を掛けてみてください。

この記事を書いている2021年2月25日現在、検索結果の1ページ目に出てくる範囲では、冒頭で解説した計算式に相当するものは2件 (そのうちの1件はPDF) しか出てきません。

d = r cos^-1( siny₁ sin y₂ + cosy₁ cosy₂ cosΔx )

2地点間の距離と方位角 - 高精度計算サイト

L = 6370 arccos( sinφ₁ sinφ₂ + cosφ₁ cosφ₂ cos( λ₁ − λ₂) )

緯度経度を用いた3つの距離計算方法 [PDF]

他のページでは、さらに複雑な数式ないしアルゴリズムが紹介されています。 もちろん、複雑だから (単純でないから) 間違いであるということはありません。 地球を真球として近似せずに、より実際に近い形状として扱うための補正を加えていけば、数式は必然的に複雑なものになっていくでしょう。 しかし問題なのは、これらのページの多くが、数式の図形的 (幾何学的) な意味を吟味・理解することなく、

と言わんばかりの、数式やプログラムをコピペで載せてお終い、という姿勢を取っていること。 個人が趣味で開設しているブログに載っているだけならまだマシですが、職業プログラマが Google 検索でこれらの情報にたどり着き、自分が業務として作成しているプログラムにこれをコピペして使い始めたとすれば (実際問題そのようなことはあちこちで起きている)、これは「カーゴ・カルト・プログラミング」に他なりません。

カーゴ・カルト・プログラミング (cargo cult programming) とは、コンピュータープログラミングにおいて、実際の目的には必要のないコードやプログラム構造を儀式的に含むことを特徴とする悪習である。 カーゴ・カルト・プログラミングは、プログラマが、自身が解決しようとしているバグや明らかな解決策を理解していないことを示す兆候である。

カーゴ・カルト・プログラミングは、目の前の問題について経験の浅いプログラマが、他の場所にあるプログラムコードを、その仕組みや、それが本当に必要かどうかを理解することなしに、別の場所にコピーするときに生じうる。

カーゴ・カルト・プログラミング - Wikipedia

前述した Google 検索結果に含まれるページには「ヒュベニ (Hubeny) の式」と呼ばれる計算手法がたびたび現れます。 (出典として適当なページが見当たらなかったので、気になる方は各自検索してください。)

• 2地点の緯度の平均が使われている
• ユークリッド距離の計算と同じ二乗和の平方根 sqrt( x² + y² ) が使われている
• 日本語以外の出典がほとんど見当たらない

などの点に多々思うところは多々ありますが、それらはひとまず措いて、この式の特徴は、地球の扁平率、すなわち、「地球が完全な球ではなく南北方向に少しつぶれている」という事実を考慮に入れていることでしょう。 確かに、この記事の冒頭で解説した計算方法はこの扁平率を考慮していないため不完全であり、こちらの「ヒュベニの式」の方が "正しい" ように見えるかもしれません。 ところで、「地球が南北方向につぶれた球である」ことをご存知の方のうち、「そのつぶれ具合はどの程度か」を把握している方はどれくらいの割合になるでしょうか。

しかし地球が完全な球形でないことは、誰でも知っている通りである。 地球の表面にはヒマラヤもあり、日本海溝もある。 そういういちじるしい小凹凸の他に、地球は南北に縮んだ楕円体になっていることは、中学校で習うとおりである。 さらに大学程度になると、それもまたちがっていて、ほんとうは疑似楕円体であるということを教わる。 ところが専門の地球物理学者にいわせると、地球の形は、それらのいずれでもなく「狐の色が狐色であるごとく、地球の形は、地球形である」という。

ところでこれらのいろいろな説明の中で、一番真に近い形をかけといわれると、結局第3図のようにコンパスで円を描くより仕方がないのである。 というわけは、次のような簡単な計算をしてみればすぐ分る。 この円は直径6センチあって、線の幅は0.2ミリである。 それでこの円を地球とみると、地球の直径1万3000キロを6センチに縮尺して描いたことになる。 この縮尺率で計算すると、線の幅0.2ミリは44キロに相当する。

ところでエヴェレストの高さは海抜8.9キロであり、海の一番深いところといわれるエムデン海溝が10.8キロの深さである。 それで地球上の凹凸の差の極限は19.7キロにすぎない。 線の幅の半分以下である。 従って地球表面の普通の山や海の凹凸を忠実に描いてみても、だいたいこの線の幅の10分の1程度ていどの凹凸になってしまうので、どうにも描きようがない。 それから地球は楕円体になっているといっても、球からのへだたりは案外少ない量であって、赤道面内の半径よりも、南北の半径が約22キロ短いだけである。それで楕円体といっても、この線の幅の半分ていどの長短があるにすぎない。(p43-44)

『科学の方法』/ 中谷宇吉郎

このように、機器や天候等によるGPS信号の受信状況への影響を気にした方が良さそうな程度の差。 「地球は完全な球体ではない」という知識を持つが故に、取るに足らない程度の差を解消するための複雑なモデルを構築してしまうというのは、「球形の牛」ジョークの真逆の状況であり、なんとも言えない可笑しさを感じます。

球形の牛 (spherical cow) は、現実に起こる複雑な現象を高度に単純化した科学的なモデルでとらえることを喩えた表現。 物理学者がしばしば、考えることのできる最も単純なかたちまで問題を縮小して計算を容易にするために、モデルを現実に適用することが難しくなってしまうことさえある、というおかしさをついている。

球形の牛 - Wikipedia

まず大地を球体と仮定してください。 更なる補正が必要かどうかは、それから考えましょう。