忙しい人のために
distance_km = 6371 * acos(
cos(latA/180*pi) * cos((lonB - lonA)/180*pi) * cos(latB/180*pi) +
sin(latA/180*pi) * sin(latB/180*pi)
);
さて
緯度経度で表される位置情報を二つもらったとき、2地点間の距離を求めるにはどうすればいいでしょう。
地球は真球ではなく回転楕円体であることを真面目に反映しようとすると、非常に難しいアルゴリズムを使う必要があります。 Vincenty法といい、収束するまで計算を繰り返して精度を高めていくタイプのアルゴリズムで、地点によっては収束しないことがあるという厄介なもの。
さてこのVincenty法のことは知識としてだけ頭に入れておいて、地球を真球、半径6371kmの球と仮定しますと、距離を求めるのは非常に容易になります。使うのは、高校の代数幾何で習ったベクトルの知識。そう、これは「図形情報を扱うのに大事なことは、全部高校で教わった」の続編となります。
2地点間の距離をベクトル計算で求める
2地点間を最短距離で結ぶ球面上の線、大圏航路と呼んだり測地線と呼んだりするものは、真球上においては必ず、地球半径と同じ半径の円、大円の一部です。
大円の中だけで考えれば、これは立体図形の問題ではありません。扇形の弧の長さを求めるだけという平面の問題です。
扇形の半径は？ そう6371km。扇形のなす角度は？ それをθと呼ぶことにしてこれから求めていきます。θがラジアンで求まれば、弧の長さは6371・θkmと決定できますから。では、答えの方から順を追っていきましょう。
- cosθかsinθが求まれば、θは算出できます。コンピュータには
sin,cosの逆関数asin,acosが使えますから。 - cosθは簡単に求まります。内積です。
地球中心から見たA地点、B地点の座標をベクトルa, bとすると、その内積abは 【aの長さ】・【bの長さ】・cosθ ですから。そして、両地点の座標ベクトルを最初から長さ1に正規化しておけば、内積は cosθ そのものです。
内積を求める計算は簡単でした。(x,y,z)と(u,v,w)の内積だったらxu+yv+zwです。 - 地点は緯度経度で表されています。これを（地球中心が原点になり、ベクトルの長さが1であるような）三次元座標に変換しないといけません。
XYZ軸の置き方は地球固定座標系というのを使いましょう。X軸が経度0緯度0のギニア湾方向、Y軸が東経90度緯度0度のスマトラ沖方向。Z軸が北極方向。するとこうです。(cos東経・cos北緯, sin東経・cos北緯, sin北緯)
ただし、経度・緯度は通常ラジアンではないので、cos, sin関数に渡すにはラジアンに変換する必要があります。そこまでちゃんと書くとこう。(cos(東経/180*π)cos(北緯/180*π), sin(東経/180*π)cos(北緯/180*π), sin(北緯/180*π))
一本の式にまとめると
以上の考え方を一本の式にまとめると、こうなります。
地点A（北緯latA, 東経lonA）, 地点B(北緯latB, 東経lonB)間の距離は、
distance_km = 6371 * acos(
cos(lonA/180*pi) * cos(latA/180*pi) * cos(lonB/180*pi) * cos(latB/180*pi) +
sin(lonA/180*pi) * cos(latA/180*pi) * sin(lonB/180*pi) * cos(latB/180*pi) +
sin(latA/180*pi) * sin(latB/180*pi)
);
と記述できます。そしてこの式は簡単にできます。始点終点を並行移動して、片方が経度0にくるように持ってきてあげましょう。lonAとlonBからそれぞれlonAを引いて0と(lonB - lonA)にしてあげるのです。すると、cos 0 == 1とsin 0 == 0が出てくるので
distance_km = 6371 * acos(
cos(latA/180*pi) * cos((lonB - lonA)/180*pi) * cos(latB/180*pi) +
sin(latA/180*pi) * sin(latB/180*pi)
);
冒頭のこの式になりました。