PHPで二地点間の方位角を求める

方位角とは

ある地図上の一点から、目標のもう一点に対して最短経路で結んだときの方角を 方位角 と言います。
例えば、東京から見てニューヨークの方位角は およそ北北東 となります。
これをPHPで計算して求めてみます。
ここでは数行で簡単に求める方法と、 Vincenty法の逆解法 を使う方法の2種類で試してみます。
緯度経度は 世界測地系 の 十進法表記 を使いますので、 日本測地系 や 度分秒表記 を使っている場合は事前に変換してください。

簡単に求める方法 (方法A)

/**
 * 指定の2地点間の地点1からの方位角(ラジアン)を求める。
 *
 * @param  float $lat1 緯度1 (十進法度単位)
 * @param  float $lng1 経度1 (十進法度単位)
 * @param  float $lat2 緯度2 (十進法度単位)
 * @param  float $lng2 経度2 (十進法度単位)
 * @return float       方位角 (ラジアン)
 */
function azimuth($lat1, $lng1, $lat2, $lng2) {
    $earth_minor_radius = 6356752.0;
    $earth_major_radius = 6378137.0;

    $lat1   = deg2rad($lat1);
    $lng1   = deg2rad($lng1);
    $lat2   = deg2rad($lat2);
    $lng2   = deg2rad($lng2);

    // 測地緯度から化成緯度を求める
    $rlat1  = atan($earth_minor_radius / $earth_major_radius * tan($lat1));
    $rlat2  = atan($earth_minor_radius / $earth_major_radius * tan($lat2));

    $dlng   = $lng1 - $lng2;

    return -atan2(sin($dlng), cos($rlat1) * tan($rlat2) - sin($rlat1) * cos($dlng));
}

Vincenty法の逆解法を使う方法 (方法B)

/**
 * Vincenty法の逆解法により、指定の2地点間の地点1からの方位角(ラジアン)を求める。
 *
 * @param  float $lat1 緯度1 (十進法度単位)
 * @param  float $lng1 経度1 (十進法度単位)
 * @param  float $lat2 緯度2 (十進法度単位)
 * @param  float $lng2 経度2 (十進法度単位)
 * @return float       方位角 (ラジアン)
 */
function azimuth_vincenty($lat1, $lng1, $lat2, $lng2) {
    $earth_minor_radius = 6356752.0;
    $earth_major_radius = 6378137.0;
    $earth_flattening   = 0.0033528106811823;
    $acceptable_error   = 1E-12;

    $lat1   = deg2rad($lat1);
    $lng1   = deg2rad($lng1);
    $lat2   = deg2rad($lat2);
    $lng2   = deg2rad($lng2);

    // 測地緯度から化成緯度を求める
    $rlat1     = atan($earth_minor_radius / $earth_major_radius * tan($lat1));
    $rlat2     = atan($earth_minor_radius / $earth_major_radius * tan($lat2));
    $sin_rlat1 = sin($rlat1);
    $cos_rlat1 = cos($rlat1);
    $sin_rlat2 = sin($rlat2);
    $cos_rlat2 = cos($rlat2);

    // 経度差を求める
    $dlng      = $lng1 - $lng2;

    $lambda      = $lng1 - $lng2;
    $prev_lambda = $lambda + 1;
    $sin_sigma   = null;
    $cos_sigma   = null;
    $cos2_alpha  = null;
    $cos_2_sigma = null;
    // ラムダ値を計算する。 (値が収束するまで繰り返す)
    for ($i = 0; $i<100 && abs($lambda - $prev_lambda) > $acceptable_error; $i++) {
        $sin_lambda  = sin($lambda);
        $cos_lambda  = cos($lambda);
        $sin_sigma   = sqrt(pow($cos_rlat2 * $sin_lambda, 2) + pow($cos_rlat1 * $sin_rlat2 - $sin_rlat1 * $cos_rlat2 * $cos_lambda, 2));
        $cos_sigma   = $sin_rlat1 * $sin_rlat2 + $cos_rlat1 * $cos_rlat2 * $cos_lambda;
        $sigma       = atan2($sin_sigma, $cos_sigma);
        $sin_alpha   = $cos_rlat1 * $cos_rlat2 * $sin_lambda / ($sin_sigma === (float) 0 ? 1 : $sin_sigma);
        $cos2_alpha  = 1 - $sin_alpha * $sin_alpha;
        $cos_2sigma  = $cos_sigma - 2 * $sin_rlat1 * $sin_rlat2 / ($cos2_alpha === (float) 0 ? 1 : $cos2_alpha);
        $c           = $earth_flattening / 16.0 * $cos2_alpha * (4 + $earth_flattening * (4 - 3 * $cos2_alpha));
        $tmp         = $sigma + $c * $sin_sigma * ($cos_2sigma + $c * $cos_sigma * (-1 + 2 * pow($cos_2sigma, 2)));
        $prev_lambda = $lambda;
        $lambda      = $dlng + (1 - $c) * $earth_flattening * $sin_alpha * $tmp;
    }

    // 必要な値を事前に計算
    $u2    = $cos2_alpha * (pow($earth_major_radius, 2) - pow($earth_minor_radius, 2)) / pow($earth_minor_radius, 2);
    $a     = 1 + $u2 / 16384 * (4096 + $u2 * (-768 + $u2 * (320 - 175 * $u2)));
    $b     = $u2 / 1024 * (256 + $u2 * (-128 + $u2 * (74 - 47 * $u2)));

    // デルタ値を求める
    $tmp   = 1.0 / 6.0 * $b * $cos_2sigma * (-3 + 4 * pow($sin_sigma, 2)) * (-3 + 4 * pow($cos_2sigma, 2));
    $delta = $b * $sin_sigma * ($cos_2sigma + 0.25 * $b * ($cos_sigma * (-1 + 2 * pow($cos_2sigma, 2)) - $tmp));
    return -atan2($cos_rlat2 * sin($lambda), $cos_rlat1 * $sin_rlat2 - $sin_rlat1 * $cos_rlat2 * cos($lambda));
}

簡単に求める方法と比べて、かなり難解なコードになりました。

結果

東京 (東京駅)から、ニューヨーク (市庁舎) への方位角を計算した結果

東京	ニューヨーク
緯度35.681236, 経度139.767125	緯度40.712775 経度-74.005973

$tokyo   = [35.681236, 139.767125];
$newyork = [40.712775, -74.005973];
echo "方法A: " . azimuth($tokyo[0], $tokyo[1], $newyork[0], $newyork[1]) . "\n";
echo "方法B: " . azimuth_vincenty($tokyo[0], $tokyo[1], $newyork[0], $newyork[1]) . "\n";

方法A: 0.4398738241381
方法B: 0.43852844993238

それぞれの関数で返ってくる結果は、方位角のラジアン値となります。
これは、真北を0とし、時計回りの方角を正の数としたものです。3.14...が真南ということになります。
度数に変換する場合はrad2deg関数を使ってください。その結果が以下となります。

簡単に求める方法	Vincenty法を使った方法
25.202913641393度	25.12582937754度

2つの方法の違い

実行時間

Vincenty法 は、途中収束するまで繰り返す処理などがあり、簡単な方式よりは処理時間がかかりそうです。
10000回計算を繰り返したときの実行時間を比べてみました。

簡単に求める方法	Vincenty法を使った方法
0.02219秒	0.09495秒

確かに、Vincenty法を使った方法の方が4倍程度時間がかかっています。
とは言え、10000回でこの程度の実行時間ですので、ほとんど気にしなくても良さそうです。

正確さ

東京から、先程求めた方位角の方角に、事前に別の方法で計算しておいた二地点間の距離だけ移動した地点の座標をプロットしてみます。

右上の赤いマーカーが本来の座標 (ニューヨーク市庁舎)です。
簡単に求める方法で得られた方位角を使ってプロットすると、左下の青いマーカーの位置になってしまいました。
ニューヨークの範囲内には収まっていますが、日本とアメリカとなると、少々ずれが大きいです。
それに対して、Vincenty法を使った方法では、赤いマーカーとほぼ寸分違わず重なりました。