緯度・経度から２地点間の距離を求める方法（PHP・MySQL）

最近、緯度・経度から２地点間の距離を求める必要があったので、その際に調べた内容を記録しておきます。

緯度・経度の数値から距離を求めるには「ヒュベニの公式」という計算式を使うそうです。
（地球は球体なので、単純に緯度経度を平面に置き換えて距離を求めようとすると、湾曲している分の誤差が発生します。）

詳しい解説は他のサイトをご覧いただくとして、「ヒュベニの公式」は次のような計算式になります。

上記のうち、定数：ａ（長半径）・ｂ（短半径）は使用する「測地系」によって異なりますが、例えばGoogleマップなどでは「世界測地系（WGS84）」が使われているので、「ａ= 6,378,137.000（m）」「ｂ= 6,356,752.314 245（m）」になります。

これらをもとに、次のようなコード（PHP）で緯度・経度から２地点間の距離を求めることができます。

$a = 6378137.000000;  // 世界測地系（WGS84）の赤道半径 (長半径) (m)
$b = 6356752.314245;  // 世界測地系（WGS84）の極半径 (短半径) (m)

// 第一離心率の二乗値
$e2 = ( pow( $a, 2) - pow( $b, 2 )) / pow( $a, 2 );

// 緯度・経度をラジアンに変換
$x1 = deg2rad( $lng1 );
$y1 = deg2rad( $lat1 );
$x2 = deg2rad( $lon2 );
$y2 = deg2rad( $lat2 );

// 緯度の差
$dy = $y1 - $y2;
// 経度の差
$dx = $x1 - $x2;

// 緯度の平均値
$mu_y = ( $y1 + $y2 ) / 2.0;

$W = sqrt( 1.0 - ( $e2 * pow( sin( $mu_y ), 2 )));

// 卯酉線曲率半径
$N = $a / $W;
// 子午線曲率半径
$M = ( $a * ( 1 - $e2 )) / pow( $W, 3 );

// ２点間の距離 (m)
$d = sqrt( pow( $dy * $M, 2 ) + pow( $dx * $N * cos( $mu_y ), 2 ));

また、緯度・経度のデータが格納されたMySQLデータベースがあり、ある地点（緯度・経度）から近い順に並び替えて取得するには、次のようなSQLを使います。

select no, lng, lat, ( 6371 * acos( cos( radians( LAT1 )) * cos( radians( lat )) * cos( radians( lng ) - radians( LNG1 )) + sin( radians( LAT1 )) * sin( radians( lat )))) as distance from TABLE order by distance;