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

# 方位角とは

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

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

```php
/**
 * 指定の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)

```php
/**
 * 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`,<br>経度`139.767125` | 緯度`40.712775`<br>経度`-74.005973` |

```php
$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";
```

```sh
方法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回でこの程度の実行時間ですので、ほとんど気にしなくても良さそうです。

## 正確さ

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

![map.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/416295/a8a7c629-4f2d-a93d-e0c6-0c60fca141eb.png)

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