Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationEventAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
161
Help us understand the problem. What are the problem?

More than 1 year has passed since last update.

@chiyoyo

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

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
161
Help us understand the problem. What are the problem?