Edited at

2つの座標間の距離を求める

2地点の緯度、経度を指定して、その直線距離を求める必要があったのでメモ書き。

<注意>

稀にソースをコピペしてブログ等に載せてる方がいらっしゃいますが、間違いを修正する事もありますので、せめてリンクくらいは貼って頂けると困る人が減るかと思います。


ヒュベニの公式

一般的に使われているヒュベニの公式から求めるものです。

国土地理院とかにある距離計算もこれで求めてるっぽいです。

/**

* 2地点間の距離(m)を求める
* ヒュベニの公式から求めるバージョン
*
* @param float $lat1 緯度1
* @param float $lon1 経度1
* @param float $lat2 緯度2
* @param float $lon2 経度2
* @param boolean $mode 測地系 true:世界(default) false:日本
* @return float 距離(m)
*/

function distance($lat1, $lon1, $lat2, $lon2, $mode=true)
{
// 緯度経度をラジアンに変換
$radLat1 = deg2rad($lat1); // 緯度1
$radLon1 = deg2rad($lon1); // 経度1
$radLat2 = deg2rad($lat2); // 緯度2
$radLon2 = deg2rad($lon2); // 経度2

// 緯度差
$radLatDiff = $radLat1 - $radLat2;

// 経度差算
$radLonDiff = $radLon1 - $radLon2;

// 平均緯度
$radLatAve = ($radLat1 + $radLat2) / 2.0;

// 測地系による値の違い
$a = $mode ? 6378137.0 : 6377397.155; // 赤道半径
$b = $mode ? 6356752.314140356 : 6356078.963; // 極半径
//$e2 = ($a*$a - $b*$b) / ($a*$a);
$e2 = $mode ? 0.00669438002301188 : 0.00667436061028297; // 第一離心率^2
//$a1e2 = $a * (1 - $e2);
$a1e2 = $mode ? 6335439.32708317 : 6334832.10663254; // 赤道上の子午線曲率半径

$sinLat = sin($radLatAve);
$W2 = 1.0 - $e2 * ($sinLat*$sinLat);
$M = $a1e2 / (sqrt($W2)*$W2); // 子午線曲率半径M
$N = $a / sqrt($W2); // 卯酉線曲率半径

$t1 = $M * $radLatDiff;
$t2 = $N * cos($radLatAve) * $radLonDiff;
$dist = sqrt(($t1*$t1) + ($t2*$t2));

return $dist;
}


球面三角法を利用したもの

GoogleMapAPIのgeometory.computeDistanceBetweenが使ってます。

実際は地球は正確な球ではないので誤差が生じますが、計算式が簡素です。

/**

* 2地点間の距離を求める
* GoogleMapAPIのgeometory.computeDistanceBetweenのロジック
* 浮動小数点の精度が足りないためGoogleより桁数が少ないかもしれません
*
* @param float $lat1 緯度1
* @param float $lon1 経度1
* @param float $lat2 緯度2
* @param float $lon2 経度2
* @return float 距離(m)
*/

public static function google_distance($lat1, $lon1, $lat2, $lon2) {
// 緯度経度をラジアンに変換
$radLat1 = deg2rad($lat1); // 緯度1
$radLon1 = deg2rad($lon1); // 経度1
$radLat2 = deg2rad($lat2); // 緯度2
$radLon2 = deg2rad($lon2); // 経度2

$r = 6378137.0; // 赤道半径

$averageLat = ($radLat1 - $radLat2) / 2;
$averageLon = ($radLon1 - $radLon2) / 2;
return $r * 2 * asin(sqrt(pow(sin($averageLat), 2) + cos($radLat1) * cos($radLat2) * pow(sin($averageLon), 2)));
}


測地線航海算法

/**

* 2地点間の距離(m)を求める
* 測地線航海算法バージョン
*
* @param float $lat1 緯度1
* @param float $lon1 経度1
* @param float $lat2 緯度2
* @param float $lon2 経度2
* @return float 距離(m)
*/

function distance($lat1, $lon1, $lat2, $lon2)
{
// 緯度経度をラジアンに変換
$radLat1 = deg2rad($lat1); // 緯度1
$radLon1 = deg2rad($lon1); // 経度1
$radLat2 = deg2rad($lat2); // 緯度2
$radLon2 = deg2rad($lon2); // 経度2

$A = 6378137.0; // 赤道半径
$B = 6356752.314140356; // 極半径
// $F = ($A - $B) / $A;
$F = 0.003352858356825; // 扁平率

$BdivA = 0.99664714164317; // $B/$A
$P1 = atan($BdivA * tan($radLat1));
$P2 = atan($BdivA * tan($radLat2));

$sd = acos(sin($P1)*sin($P2) + cos($P1) * cos($P2) * cos($radLon1 - $radLon2));

$cos_sd = cos($sd/2);
$sin_sd = sin($sd/2);
$c = (sin($sd) - $sd) * pow(sin($P1)+sin($P2),2) / $cos_sd / $cos_sd;
$s = (sin($sd) + $sd) * pow(sin($P1)-sin($P2),2) / $sin_sd / $sin_sd;
$delta = $F / 8.0 * ($c - $s);

return $A * ($sd + $delta);
}

参考

http://www2.nc-toyama.ac.jp/~mkawai/lecture/sailing/geodetic/geosail.html