アルゴリズム。ライブラリも当然あるとは思うが中身を知っておいたほうがよい。

緯度1度あたりの距離

地球の半径 = r
地球の円周 = 2 x PI x r
緯度1度あたりの距離 = 2 x PI x r / 360

これに具体値を入れると

2 x 3.14 x 6378 / 360 = 111.260666667

となり、経度1度あたりの距離となる。

経度は緯度によって1度あたりの距離が変わる(極に向かうとすぼまる)。赤道上では緯度と一致するので 111km になる。

地球を北極側から見ると緯線がなす(緯度違いで同心円に見える)円の半径は地球の半径を斜辺とする直角三角形の底辺になる。

地球の半径 = r
緯線がなす円の半径 = r x cos(緯度)
緯線がなす円の円周 = 2 x PI x r x cos(緯度)
経度1度あたりの距離 = 2 x PI x r x cos(緯度) / 360

(lat1, lon1), (lat2, lon2) 間の距離はこうなる

√(((lat1 - lat2) x 2 x PI x r / 360)^2 + ((lon1 - lon2) x 2 x PI x r x cos(lat1) / 360)^2)

具体値を入れるとこんな感じ

√(((lat1 - lat2) x 111)^2 + ((lon1 - lon2) x 111 x cos(lat1))^2)

緯度を日本付近に限定すると cos(35d) で 0.8 で

√((111(lat1 - lat2))^2 + (89(lon1 - lon2))^2)

こんな感じになるだろう

日本の最北端は北緯 45.5572222222222222 最南端は 20.4252777777777778 ということで cos は 0.70 と 0.93 ということになってかなり差があるのでこのへんの0.8で限定するのは誤差的には問題ありだな。

数学としては球面上に描画された三角形を扱う「球面三角法」というモノがあるのでこれを使う。 この三角法の特徴は辺の長さを球の中心角として扱うところである。

北極を頂点とした三角形を描くと北極の角は経度の差として表され、その角を挟む2辺の長さは直角から緯度を引いた角度として表される。 ここで2点間の距離はこの三角形の底辺として表現されることになる。

通常の三角形で二辺と挟む角から底辺を求める場合は「余弦定理」を使うが、球面三角法でも余弦定理が存在し底辺の長さ、球面三角法では中心角を求める事ができる。 中心角さえ分かってしまえば球の半径から長さが求まる

地球は赤道方向に長く潰れた饅頭のような形をしていると言われるが、その潰れは非常に微々たるモノで、用途にもよるだろうが、ほぼ球として扱って近似できる。