GPSデータから距離を算出する(ヒュベニの距離計算式)Perl版

以下、昔の自分のブログからの転写(2008-02-15 21:40)からの転写(NEW!)

GoogleMapsなどを使ってGPSデータを500地点ぐらい取った(仕事でバイト君にやらせた)。

そのデータを携帯電話を使って現在地から一番近いところを探せたら便利だなと。それはデータを引き算すれば良いのだけど、じゃぁ、どのぐらい離れているのか?GPSデータからの計算方法が分からなかった。ネットを調べているとヒュベニの計算式を使えばよいらしい。

これは地球の丸み(楕円形)を計算にいれているらしい。そんな細かくする必要はないのだが、地球の丸みを計算に入れるなんて。。。。バカらしいぐらいに素敵な計算式だ。

#!/usr/bin/perl

# 東京銀座 数寄屋橋交差点
my $lat1="35.672558912630656";
my $lon1="139.7632384300232";

# 東京銀座 銀座四丁目交差点
my $lat2="35.67123411394145";
my $lon2="139.76503014564514";

print int(&distance($lat1, $lon1, $lat2, $lon2)) ."m\n";

# ヒュベニの距離計算式
sub distance {

my ($lat1, $lon1, $lat2, $lon2) = @_;
my ($pi, $rd, $dP, $dR, $P, $M, $N);

# 定数
$pi = 4 * atan2(1,1); # 円周率
$rd = $pi / 180; # [ラジアン/度]

# ラジアンに変換
$lat1 *= $rd;
$lon1 *= $rd;
$lat2 *= $rd;
$lon2 *= $rd;

# dP:2点の緯度差
$dP = abs($lat1 - $lat2);
# dR:2点の経度差
$dR = abs($lon1 - $lon2);

# P:2点の平均緯度
$P = $lat1 + (($lat2 - $lat1) / 2);

# M:子午線曲率半径 M=6334834/sqrt((1-0.006674*sin(P)*sin(P))^3)
$M = 6335439 / sqrt((1 - 0.006694 * (sin($P) ** 2))**3);

# N:卯酉線曲率半径 N=6377397/sqrt(1-0.006674*sin(P)*sin(P))
$N = 6378137 / sqrt(1 - 0.006694 * (sin($P) ** 2));

# D:2点間の距離(m)
# D=sqrt((M*dP)*(M*dP)+(N*cos(P)*dR)*(N*cos(P)*dR))
sqrt((($M * $dP) ** 2) + (($N * cos($P) * $dR) ** 2));

}

__END__

■ atan2
-πからπの範囲でY/Xのアークタンジェント(逆正接)を返す
返却値の角度は、-PI~PI (-PIを除く) の範囲のラジアン
Usage: atan2(1.5,1);

■ abs
引数EXPRで指定した数値の絶対値を返す
Usage: abs($a - $b);

■ sqrt
引数の平方根を返す
Usage: sqrt(144);

▼動作結果
% 218 m

また、わざわざ計算しなくても緯度の0.01度(36秒)はおよそ1.11km、経度の0.01度(36秒)はおよそ0.91kmになるので、それを使っておよそ何キロ内を調べる方が速いんだけどね。