schedule2019-06-06

三角球面法を使った2点間距離の算出 | php

PHPで三角球面法を使って2点間距離を算出する関数を作成しました。

精度は正確には測っていません。

三角球面法

球面三角法 - Wiki

基本的な三角関数を組み合わせて、球面上での辺や角の関係を扱う分野です。 平面にもある余弦定理や正弦定理もまとまっています。 地球は丸いため、球面上の計算が必要です。

球面上での2点間距離の式

緯度経度を用いた3つの距離計算方法 pdfの式(1)を基にしています。

赤道半径をRR、2地点の緯度経度をそれぞれ(緯度ϕ1\phi_1, 経度λ1\lambda_1)、(緯度ϕ2\phi_2, 経度λ2\lambda_2)として、2点間距離LLとすると以下の式になります。 距離の単位は赤道半径に依ります。

L=Rarccos(sinϕ1sinϕ2+cosϕ1cosϕ2cos(λ1λ2))L = R \arccos{( \sin \phi_1 \sin \phi_2 + \cos \phi_1 \cos \phi_2 \cos{(\lambda_1 - \lambda_2)})}

phpでのソースコード

クラーク胸像と札幌駅間の直線距離を計算してみました。

phpの三角関数はラジアンに変換する必要があります。

/**
 * 2点間距離[m]を算出する。
 * 引数の緯度経度の単位は度。
 *
 * @param float $lat1_deg 地点1の緯度
 * @param float $lng1_deg 地点1の経度
 * @param float $lat2_deg 地点2の緯度
 * @param float $lng2_deg 地点2の経度
 * @return float
 */
function calcDistance($lat1_deg, $lng1_deg, $lat2_deg, $lng2_deg){
    // 赤道半径[m]
    $pole_radius = 6378137.0;

    // ラジアンに変換
    $lat1 = deg2rad($lat1_deg);
    $lng1 = deg2rad($lng1_deg);
    $lat2 = deg2rad($lat2_deg);
    $lng2 = deg2rad($lng2_deg);

    // 経度差
    $lngDiff = $lng1 - $lng2;
    $dist = sin($lat1)*sin($lat2) + cos($lat1)*cos($lat2)*cos($lngDiff);
    
    return $pole_radius * acos(min(max($dist,-1.0),1.0));
}


// クラーク胸像
$lat1_deg = 43.070694;
$lng1_deg = 141.343490;

// 札幌駅
$lat2_deg = 43.068806;
$lng2_deg = 141.350699;

$distance = calcDistance($lat1_deg, $lng1_deg, $lat2_deg  ,$lng2_deg);

echo $distance; // 622.78135287911

GoogleMapでの計測が623.40mでしたので差は1m以内でした。(クリックした位置がずれた可能性もあります。)

map

おわりに

より高い精度はヒュベニの公式やライブラリで計算されると良いかと。

参考

JavaScriptでのヒュベニと球面三角法の関数、GeoLocationAPIでの計算が載っています。