166
Help us understand the problem. What are the problem?

posted at

updated at

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

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

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

ヒュベニの公式

一般的に使われているヒュベニの公式から求めるものです。
国土地理院とかにある距離計算もこれで求めてるっぽいです。
【追記】
2022/05/12 バグ修正 >@rwanda_go_tanさんバグ報告ありがとうございます。

/**
 * 2地点間の距離(m)を求める
 * ヒュベニの公式から求めるバージョン
 *
 * @param float $lat1 緯度(始点)
 * @param float $lon1 経度(始点)
 * @param float $lat2 緯度(終点)
 * @param float $lon2 経度(終点)
 * @param boolean $mode 測地系 true:世界(default) false:日本
 * @return float 距離(m)
 */
function distance($lat1, $lon1, $lat2, $lon2, $mode=true)
{
    // 緯度経度をラジアンに変換
    $radLat1 = deg2rad($lat1); // 緯度(始点)
    $radLon1 = deg2rad($lon1); // 経度(始点)
    $radLat2 = deg2rad($lat2); // 緯度(終点)
    $radLon2 = deg2rad($lon2); // 経度(終点)

    // 緯度差
    $radLatDiff = abs($radLat1 - $radLat2);

    // 経度差算(距離の短い方を採用する)
    $radLonDiff = abs($radLon1 - $radLon2);
    $radLonDiff = min($radLonDiff, 2 * pi() - $radLonDiff);

    // 平均緯度
    $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 緯度(始点)
 * @param float $lon1 経度(始点)
 * @param float $lat2 緯度(終点)
 * @param float $lon2 経度(終点)
 * @return float 距離(m)
 */
public static function google_distance($lat1, $lon1, $lat2, $lon2) {
    // 緯度経度をラジアンに変換
    $radLat1 = deg2rad($lat1); // 緯度(始点)
    $radLon1 = deg2rad($lon1); // 経度(始点)
    $radLat2 = deg2rad($lat2); // 緯度(終点)
    $radLon2 = deg2rad($lon2); // 経度(終点)

    $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 緯度(始点)
 * @param float $lon1 経度(始点)
 * @param float $lat2 緯度(終点)
 * @param float $lon2 経度(終点)
 * @return float 距離(m)
 */
function distance($lat1, $lon1, $lat2, $lon2)
{
    // 緯度経度をラジアンに変換
    $radLat1 = deg2rad($lat1); // 緯度(始点)
    $radLon1 = deg2rad($lon1); // 経度(始点)
    $radLat2 = deg2rad($lat2); // 緯度(終点)
    $radLon2 = deg2rad($lon2); // 経度(終点)

    $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

Register as a new user and use Qiita more conveniently

  1. You can follow users and tags
  2. you can stock useful information
  3. You can make editorial suggestions for articles
What you can do with signing up
166
Help us understand the problem. What are the problem?