距離については，以前 GPS 経路データから求めた距離の誤差について というのを書いた。 その時は六甲全山縦走路の距離に関して，GNSS (GPS) で得られた経路の距離の誤差について考えたものだった。 今回は，距離の求め方について書こうと思う。ま，ほとんどどっかのサイトの受け売りやけどね。
地図上の距離を求めるには定規で測って縮尺を考えて距離を出す，というのも考えられるが，GNSS (GPS) で経路を記録していると，緯度＋経度，という情報が得られる。 そこで，2点間の緯度と経度の情報を使って距離を求める，というのが今回のお話し。
緯度の経度について
まず，距離を求める前に緯度と経度の話を書こう。緯度は英語で「Latitude」，経度は「Longitude」と言う。 なので，計算スクリプトで lat やら Lat やらとあれば「緯度」を示し，lng や Lng, 場合によっては lon や Lon とあれば「経度」を表していると考えて欲しい。 ちなみに「高度」は「Altitude」が使われることが多い（と思う）。
緯度と言えば「北緯 35度25分50秒」，経度と言えば「東経 135度15分44秒」みたいなのを思い浮かべるかもしれないが，分，秒を用いると計算がややこしくなる。 そこで，ここではこれ以降，緯度，経度は「度」で表し，細かいところは小数で表すことにしよう。 つまり「北緯 35.430556，東経 135.262222」みたいに表すことにする。
また，測地系としては「WGS84」を用いる。 測地系は，地球上の場所の座標であり，いわゆる緯度，経度である。しかし，その基準をどうとるか（地球をどんな楕円体で近似するか，原点をどこにとるか）などで異なるみたい（それ以上はわかっていないかも…）。今使われているのは GRS80 と呼ばれるもの（1980年に決められたもの。多分，微妙な補正が入ってる思うけど）である。しかし，GPS 業界では WGS84 が用いられている。GRS80 と WGS84 の違いは，扁平率の微妙な違いみたい。どちらも，楕円体の長半径は 6,378,137 m なのだが，扁平率が GRS80 は 1/298.257 222 101 なのに対して，WGS84の扁平率は 1/298.257 223 563 らしい。分母の小数点以下4桁目からが違っている。微妙やなぁ…。具体的に短半径を書くと，GRS80 が 6,356,752.314 140 m なのに対し，WGS84 が 6,356,752.314 245 m となる。こんな微妙な違いなので，実用上は同じとして問題ない。 なお，昔の（2002年4月1日以前の）日本の測地系は，日本独自の日本測地系（旧日本測地系、Tokyo Datum）というものが使われていたので，古い資料を見る際には注意が必要である。
座標を表す場合，緯度と経度の順が違う場合がある。 地理院地図では「緯度」が先みたいだが，私はあまり細かいことは気にしていない。 基本的に日本国内の位置情報を考えているので，100 を越えていたら経度（東経）で，20〜40 ぐらいなら緯度（北緯）と思っている。 計算のスクリプトを書く際には，どちらを先にするかを決めておかないといけないけどね。 当然，逆の順に入力された際の処理を加えておくのがよさげである。
ヒュベニの公式
まずはヒュベニの公式を書こう。 原理のようなものに関しては，ネット上で検索してみてほしい。 一応，Perl のスクリプトを書いておこう。 my $long_r = 6378137.000; # [m] 長半径 GRS80 と WGS84 で共通
# my $short_r = 6356752.314140; # [m] 短半径 GRS80
my $short_r = 6356752.314245; # [m] 短半径 WGS84
my $rishin = sqrt(($long_r**2 - $short_r**2)/$long_r**2); # 第一離心率
my $a_e_2 = $long_r*(1-$rishin**2); # a(1-e^2)
my $pi = 3.14159265358979; # Pi
$lat0 = $lat0*$pi/180; # in radian
$lng0 = $lng0*$pi/180; # in radian
$lat1 = $lat1*$pi/180; # in radian
$lng1 = $lng1*$pi/180; # in radian
my $d_lat = $lat1 - $lat0;
my $d_lng = $lng1 - $lng0;
my $ave_lat = ($lat1+$lat0)/2;
my $W = sqrt(1-$rishin**2*sin($ave_lat)*sin($ave_lat));
my $N = $long_r/$W;
my $M = $a_e_2/$W/$W/$W;
my $dist = ($d_lat*$M)**2 + ($d_lng*$N*cos($ave_lat))**2; # square of distance
$dist = sqrt($dist);
ここで，$lat0, $lng0 が始点の位置座標であり，$lat1, $lng1 が終点の座標である。
それをラジアンに変換して計算に用いている。
最後の $dist が求める2点間の距離となっている。測地線航海算法
次に測地線航海算法について書こう。 これも比較的簡単に距離が求められるものである。 ただし，計算に tan, arccos (acos) と arctan (atan) を用いるので，Perl の場合は事前に Math モジュールを読み込んでおかないといけない。use Math::Trig qw/tan atan acos/;距離を求める部分を以下に書いておこう。
my $long_r = 6378137.000; # [m] 長半径 GRS80 と WGS84 で共通
# my $short_r = 6356752.314140; # [m] 短半径 GRS80
my $short_r = 6356752.314245; # [m] 短半径 WGS84
my $pi = 3.14159265358979; # Pi
# 扁平率
my $FF = ($long_r - $short_r)/$long_r;
# 度(deg) → radian
$lat0 = $lat0*$pi/180; # in radian
$lng0 = $lng0*$pi/180; # in radian
$LAT1 = $LAT1*$pi/180; # in radian
$LNG1 = $LNG1*$pi/180; # in radian
# 化成精度（phi）
my $phi0 = atan($short_r/$long_r * tan($lat0));
my $Phi1 = atan($short_r/$long_r * tan($LAT1));
# 球面上の距離 (len)
my $XX = acos( sin($phi0)*sin($Phi1) + cos($phi0)*cos($Phi1) * cos($lng0 - $LNG1) );
# 距離補正量 (delta_rho)
my $delta_rho = $FF/8 * ( (sin($XX) - $XX) * (sin($phi0)+sin($Phi1))**2 / cos($XX/2)**2
- (sin($XX) + $XX) * (sin($phi0) - sin($Phi1))**2 / sin($XX/2)**2);
# 測地線距離 (m)
my $dist = $long_r * ($XX + $delta_rho);
最後の行の $dist が求める2点間の距離である。
ここで $LAT1, $LNG1 のように大文字を使っているが，これは単にパッと見て $lat0, $lng1 と区別しやすいようにしただけで，それ以上の意味はない。ヒュベニの公式と航海算法の比較
ここで，上記の2つの方法を地理院の提供している距離計算のページの結果と比べてみよう。 地理院のページを使う際には，測地系は「GRS80」とし，「十進法度単位」を選んでおこう。(1) 長距離の場合（1000 km 程度）
まずは南北方向の場合を見るために，千葉県銚子市のポートタワー付近と，北海道稚内の国道の終点との距離を見てみよう。 ここでは，念のために測地系として GRS80 と WGS84 の両方の結果を書いている。銚子 (35.740961,140.862912) -- 稚内 (45.417198,141.676276)
地理院：1076.697201 km
WGS84
ヒュベニ：1076.70987537832 km 相対誤差：0.001177 %
航海算法：1076.69573094928 km 相対誤差：-0.000136 %
---------------------------------------------------------
GRS80
ヒュベニ：1076.70987536532 km 相対誤差：0.001177 %
航海算法：1076.69573093629 km 相対誤差：-0.000136 %
上記の結果を見ると，まず，GRS80 と WGS84 は実用上は全く差がないことがわかる。 ここには書いていないが，他にもいろいろ比較したが，GRS80 と WGS84 の間にほぼ差はなかった。 そこで，これ以降は WGS84 の結果のみを記載しよう。 また，ヒュベニの公式と航海算法にも大差ないが，どちらかというと航海算法の方が誤差が少ないことがわかる（地理院のページの計算に比べて，やけどね）。 どうも距離が長い（100 km より十分長い？）と，航海算法の方が精度がいいみたい。
次に東西方向のケースとして，銚子と対馬の間の距離を見てみよう。
銚子 (35.740961,140.862912) -- 対馬 (34.202710,129.287673)
地理院：1070.069226 km
ヒュベニ：1070.73123142042 km 相対誤差：0.061865 %
航海算法：1070.06913218936 km 相対誤差：-8e-06 %
ここでは航海算法の方はほぼ誤差なしだが，ヒュベニの公式の方は誤差が少し増えている。 南北方向（銚子--稚内）だとあまり差がなかったのに，東西方向（銚子--対馬）で差が出るのは， ヒュベニの公式は距離が長くなると，東西方向の方がより誤差が大きいことを示している。
ここでは2例しか示していないが，原理的に距離が長いと，南北，東西に関わらず，ヒュベニの公式よりも航海算法の方が精度がいいみたい。
(2) 中距離の場合（10 km 〜 100 km 程度）
次に条件を変えて，10 km 〜 100 km 程度の距離の場合を見てみよう。 まずは南北方向。ここでは上記の東京都庁付近から，宇都宮の東武宇都宮駅前までの距離を計算させている。東京都庁--東武宇都宮駅
地理院：98.126294 km
ヒュベニ：98.126361841262 km 相対誤差：6.9e-05 %
航海算法：98.1261650663422 km 相対誤差：-0.000131 %
これを見るとヒュベニの公式と航海算法の間にあまり差はない（細かい違いはあるが，実用上は差がないレベル，という意味）。
次に東西方向として，東京都庁から甲州市塩山駅前までの距離を見てみる。
東京都庁--甲州市塩山駅
地理院：86.658714 km
ヒュベニ：86.659057410539 km 相対誤差：0.000396 %
航海算法：86.6587135069539 km 相対誤差：0 %
ここでもヒュベニの公式と航海算法にはあまり差がない。多少航海算法の方がよさげではあるが…。
念のために少し距離を短くしたものも示そう。
東京都庁--さいたま県庁
地理院：18.983221 km
ヒュベニ：18.9832215547376 km 相対誤差：2e-06 %
航海算法：18.983195965596 km 相対誤差：-0.000131 %
(3) 短距離の場合（100 m 以下程度）
次に極端に条件を変えて，100 m 以下程度の距離の場合を見てみよう。 まずは南北方向。ここでは上記の東京都庁付近の data（35.68944,139.691670）を用い，緯度の数字の小数点以下4桁目をいじって2点を作っている。東京都庁付近 南北方向（35.68944 <-> 35.68964）
地理院：0.022191 km (22 m)
ヒュベニ：0.0221906511971943 km 相対誤差：-0.001571 %
航海算法：0.0221908520376506 km 相対誤差：-0.000666 %
次に東西方向。今度は経度の数値の小数点以下4桁目をいじっている。
東京都庁付近 東西方向（139.691670 <-> 139.691770）
地理院：0.0090516 km (9 m)
ヒュベニ：0.00905158702795335 km 相対誤差：-0.000143 %
航海算法：0.00905143839182423 km 相対誤差：-0.001785 %
この場合，ヒュベニの公式と航海算法の間に大差はないように見える。
もう一桁距離を近づけてみると，以下のようになった。
東京都庁付近 南北方向（35.68944 <-> 35.68946）
地理院：0.0022191 km (2.2 m)
ヒュベニ：0.0022190650867564 km 相対誤差：-0.001573 %
航海算法：0.00221993984116535 km 相対誤差：0.037846 %
# 東京都庁付近 東西方向（139.691670 <-> 139.691690）
地理院：0.0018103 km (1.8 m)
ヒュベニ：0.00181031740559067 km 相対誤差：0.000961 %
航海算法：0.0018082906698376 km 相対誤差：-0.110994 %
これを見ると，極端に短距離の場合は，ヒュベニの公式の方が航海算法よりも精度がよさそうな感じに見える。 まぁ，2 m 程度となると，使っている GNSS (GPS) ロガーの精度よりよくなるので，そこまで必要なのか？という気もするが…。
ということで，今回のまとめとしては，たかだか2例ずつしか示せていないし，あくまで地理院のサイトの計算が正しいとしての結果だが，
(1) 長距離（100km 以上など）では，航海算法の方が精度がよさそう。
(2) 極端な短距離（100m 以下など）では，ヒュベニの公式の方が精度が良さそう。
ということになった。
考えてみれば，航海算法は遠くの港に行くための計算なので，あまり短距離は使わないから，長距離で精度がいいのは正しいことかもしれない。 しかし，私は登山で歩いた際の GNSS (GPS) のログから距離を出したいので，その場合はヒュベニの公式を使おうかな？と思ってる。
0 件のコメント:
コメントを投稿