緯度経度で与えられる2地点間の距離を計算する PHP ライブラリ｜加藤正人｜デベロッパーズブログ

DEVELOPERS BLOGデベロッパーズブログ

HOME >
加藤正人 >
タイトル

2014/08/23

緯度経度で与えられる2地点間の距離を計算する PHP ライブラリ

緯度・経度で与えられる2つの地点間の距離を計算する PHP ライブラリを考えてみた。

緯度経度による2地点間の距離は、Hubeny の公式が良く用いられているので、これを利用する。 (以下 http://yamadarake.jp/trdi/report000001.html より)

Hubeny の公式

地点1の経度 (ラジアン)を x_1、緯度を y_1 (ラジアン) とする。同様に地点2の経度および緯度をそれぞれ x_2、y_2 とすると、この2地点間の距離 d (m) は次の式で計算できる。

d = \sqrt[]{\mathstrut (d_y M)^2+(d_x N \cos{\mu_y})^2}

ここで各変数は以下のようになる。

d_x = x_1 - x_2	(2地点の経度の差)
d_y = y_1 - y_2	(2地点の緯度の差)
\mu_y = \frac{y_1 + y_2}{2}	2点間の緯度の平均値
M = \frac{a(1-e^2)}{W^3}	子午線曲率半径
N = \frac{a}{W}	卯酉線曲率半径
W = \sqrt{\mathstrut 1 - e^2 \sin^2 \mu_y}

なお、上記で使用される各定数は以下の通り。

e = \sqrt{\mathstrut \frac{a^2 - b^2}{a^2}}	第一離心率
a	長半径 (赤道半径)
b	短半径 (極半径)

ここで a (赤道半径)、b (極半径) は使用する楕円体モデルにより規定されているので、その値 (m) を使用する。Google Maps や GPS であれば WGS84 楕円体であり、a = 6,378,137.000、b = 6,356,752.314 245 となる。

PHP での実装

コアになる機能は以下のように実装できる。

class HubenyDistance
{
	// ヒュベニの公式に基づき緯度経度で与えられた２点間の距離を計算する。
	const a = 6378137.000000;	// WGS84 楕円体モデルの地球赤道半径 (長半径) (m)
	const b = 6356752.314245; // WGS84 楕円体モデルの地球極半径 (短半径) (m)

	// $point1: array('lon' => float, 'lat' => float)
	// $point2: array('lon' => float, 'lat' => float)
	public static function calculate($point1, $point2)
	{
		$lat1 = $point1['lat'];
		$lon1 = $point1['lon'];
		$lat2 = $point2['lat'];
		$lon2 = $point2['lon'];

		$e2 = (pow(self::a, 2) - pow(self::b, 2))/pow(self::a, 2);	// 第一離心率の二乗値

		$x1 = deg2rad($lon1);	// 経度 ⇒ ラジアン変換
		$y1 = deg2rad($lat1);	// 緯度 ⇒ ラジアン変換
		$x2 = deg2rad($lon2);	// 経度 ⇒ ラジアン変換
		$y2 = deg2rad($lat2);	// 緯度 ⇒ ラジアン変換

		$dy = $y1 - $y2;	// 緯度の差
		$dx = $x1 - $x2;	// 経度の差

		$mu_y = ($y1 + $y2) / 2.0; 	// 緯度の平均値

		$W = sqrt(1.0 - ($e2 * pow(sin($mu_y), 2)));

		$N = self::a / $W;	// 卯酉線曲率半径

		$M = (self::a * (1 - $e2)) / pow($W, 3);	// 子午線曲率半径

		return sqrt(pow($dy * $M, 2) + pow($dx*$N*cos($mu_y), 2));	// ２点間の距離 (m);
	}
}

これを元に CakePHP のヘルパーやコンポーネントを作れば、距離計算の必要なサイトの構築が楽になるかもしれない。