PHP - 2点間距離の計算
ある地点からある地点までの正確な距離を計算する方法をいくつかのアルゴリズムで計算してみる。
実行結果
東京タワーからエッフェル塔までの距離 Haversine 9719.47381 Hubeny 11436.102719 Lambert-Andoyer 9835.7223634185
Haversineは、球面三角法を利用した計算。球であることが前提で、楕円ではないので地球にあてはめると誤差が出る。
Hubenyは、赤道が膨らん考慮されているが、距離が離れると誤差がひどくなる。
Lambert-Andoyerは、測地線航海算法を利用した計算。計算時間は最もかかるが一番正確。
ソースコード
<?php
/**
* 2点間の直線距離を求める(Lambert-Andoyer)
*
* @param float $lat1 始点緯度(十進度)
* @param float $lon1 始点経度(十進度)
* @param float $lat2 終点緯度(十進度)
* @param float $lon2 終点経度(十進度)
* @return float 距離(Km)
*/
function distance_lambert($lat1, $lon1, $lat2, $lon2) {
// WGS84
$A = 6378137.0; // 赤道半径
$F = 1 / 298.257222101; // 扁平率
// 扁平率 F = (A - B) / A
$B = $A * (1.0 - $F); // 極半径
$lat1 = deg2rad($lat1);
$lon1 = deg2rad($lon1);
$lat2 = deg2rad($lat2);
$lon2 = deg2rad($lon2);
$P1 = atan($B/$A) * tan($lat1);
$P2 = atan($B/$A) * tan($lat2);
// Spherical Distance
$sd = acos(sin($P1)*sin($P2) + cos($P1)*cos($P2)*cos($lon1-$lon2));
// Lambert-Andoyer Correction
$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);
// Geodetic Distance
$distance = $A * ($sd + $delta);
return $distance / 1000.0;
}
/**
* 2点間の直線距離を求める(Hubenyの公式)
*
* @param float $lat1 始点緯度(十進度)
* @param float $lon1 始点経度(十進度)
* @param float $lat2 終点緯度(十進度)
* @param float $lon2 終点経度(十進度)
* @return float 距離(Km)
*/
function distance_hubeny($lat1, $lon1, $lat2, $lon2) {
$radlat = deg2rad(abs($lat1 - $lat2));
$radlon = deg2rad(abs($lon1 - $lon2));
$average = deg2rad($lat1 + (($lat2 - $lat1) / 2));
$meridian = 0;
$prime = 0;
// 世界測地系
$temp = 1.0 - 0.00669438 * pow(sin($average),2);
$meridian = 6335439.0 / sqrt(pow($temp,3)); // 子午線曲率半径
$prime = 6378137.0 / sqrt($temp); // 卯酉線曲率半径
// ミリメートルまで
return round(sqrt(pow($meridian*$radlat,2) + pow($prime*cos($average)*$radlon,2)),3) / 1000;
}
/**
* 2点間の直線距離を求める(Haversineの公式)
*
* @param float $lat1 始点緯度(十進度)
* @param float $lon1 始点経度(十進度)
* @param float $lat2 終点緯度(十進度)
* @param float $lon2 終点経度(十進度)
* @return float 距離(Km)
*/
function distance_haversine($lat1, $lon1, $lat2, $lon2) {
$radius = 6371; // radius of earth in meters
$dlat = deg2rad($lat2 - $lat1);
$dlon = deg2rad($lon2 - $lon1);
$rlat1 = deg2rad($lat1);
$rlat2 = deg2rad($lat2);
$a = sin($dlat/2) * sin($dlat/2)
+ sin($dlon/2) * sin($dlon/2) * cos($rlat1) * cos($rlat2);
$c = 2 * atan2(sqrt($a), sqrt(1-$a));
return round($radius * $c,6);
}
?>