二地点の緯度経度から距離を得る必要が出てきたので、調べてみるとヒュベニの公式というのが簡単でよいみたい。
ということで、「やまだらけ」さんによるJavaの実装[1]を Perl に移植しました。
使っている測地系は GRS80(世界測地系)のみです。
■コード (latlong-distance.pl):
■実行例: サンプルデータは[1]より。
出力の単位は m(メートル)です。
ということで、「やまだらけ」さんによるJavaの実装[1]を Perl に移植しました。
使っている測地系は GRS80(世界測地系)のみです。
■コード (latlong-distance.pl):
#!/usr/bin/env perl
use strict;
use warnings;
my ($lat1, $lng1, $lat2, $lng2) = @ARGV;
my $d = calcDistHubeny($lat1, $lng1, $lat2, $lng2);
print "$d\n";
sub calcDistHubeny {
my ($lat1, $lng1, $lat2, $lng2) = @_;
my $PI = 3.14159265358979; # atan2(1,1)*4;
my $a = 6378137.000; # $GRS80
my $e2 = 0.00669438002301188; # $GRS80
my $mnum = 6335439.32708317; # $GRS80
my $my = (($lat1 + $lat2) / 2.0) * $PI / 180.0;
my $dy = ($lat1 - $lat2) * $PI / 180.0;
my $dx = ($lng1 - $lng2) * $PI / 180.0;
my $sin = sin($my);
my $w = sqrt(1.0 - $e2 * $sin * $sin);
my $m = $mnum / ($w * $w * $w);
my $n = $a / $w;
my $dym = $dy * $m;
my $dxncos = $dx * $n * cos($my);
return sqrt($dym * $dym + $dxncos * $dxncos);
}
■実行例: サンプルデータは[1]より。
./latlong-distance.pl 36.10056 140.09111 35.65500 139.74472 58502.4589312411 ./latlong-distance.pl 35.65500 139.74472 33.59532 130.36208 890233.06430987 ./latlong-distance.pl 35.802739 140.380034 35.785796 140.392265 2180.94846976372
出力の単位は m(メートル)です。
参考文献
- [1] 二地点の緯度・経度からその距離を計算する(日本は山だらけ~)
http://yamadarake.jp/trdi/report000001.html - [2] 測量計算(距離と方位角の計算)
http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/bl2stf.html
