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

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

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

ヒュベニの公式

一般的に使われているヒュベニの公式から求めるものです。
国土地理院とかにある距離計算もこれで求めてるっぽいです。
2024年現在、国土地理院の計算式はBowring(1996)になっています。
参考：経緯度を用いた2地点間の測地線長、方位角を求める計算

【追記】
[2022/05/12]　バグ修正　＞@rwanda_go_tanさんバグ報告ありがとうございます。
[2023/11/01] composer化してみました。
[2024/01/04] 国土地理院の距離計算が異なっていたので本文修正 ＞ @kkddさんありがとうございます。

/**
 * ２地点間の距離(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が使ってます。
実際は地球は正確な球ではないので誤差が生じますが、計算式が簡素です。

/**
 * ２地点間の距離を求める
 *   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)));
}

測地線航海算法

/**
 * ２地点間の距離(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