日本は山だらけ〜
技術研究本部
報告書 第一号
平成二十一年八月十五日 公開
二地点の緯度・経度からその距離を計算する
更新：コードのライセンスを MIT ライセンスに変更しました。（平成二十四年五月十日）
概要
各山のページに、
その山の近くにある別の山のリストを載せようと思い、
二地点の緯度と経度から距離を求める方法について調査しました。
ヒュベニ (Hubeny) の公式を用いると簡単な計算で精度よく距離が求まることがわかりました。
ヒュベニの公式ではいくつか定数が出てきますが、旧日本測地系、世界測地系、
GPS で利用されている測地系の 3 種類についての定数を調べ、算出しました。
Java および R のプログラムと合わせて研究成果を報告します。
緯度経度からの距離計算
国土地理院・測地部のウェブサイト
[1] に測量に関する様々な情報が載っています。
特に
[2] には様々な計算式が載っており、
緯度・経度から距離を計算する式
[3] もあります。
しかしここにあるのは緯度・経度を平面直角座標という座標系に変換した上で距離を計算で、
少し複雑です。
インターネットを調べるとヒュベニ (Hubeny) の公式と呼ばれている計算式が多く使われているようで、
3D 地図ソフトウェアとして有名なカシミール 3D
[4] でも、
距離計算に使用されているようです。
カシミール 3D のマニュアルページ
[5] にヒュベニの公式が紹介されています。
このヒュベニの公式は様々なところで紹介されていますが、
残念ながら１次ソース的な文献や論文は見つけることができませんでした。
しかし実際に国土地理院・測地部の距離計算プログラムと比較すると、
カシミール 3D の式とほぼ同じ値になることからこれをヒュベニの公式として引用します。
カシミール 3D ではベッセル楕円体に基づく旧日本測地系用のパラメータが用いられています。
そこで、国土地理院・測地部で紹介されている計算式を用い、
現在の正式な測地系（世界測地系または日本測地系２０００）である GRS80 に基づくパラメータおよび、
GPS で使われる WGS84 に基づくパラメータの算出を行いました。
ただし GRS80 と WGS84 は実用上はほとんど差がありません。
ヒュベニの公式
ヒュベニの公式を Fig. 1 に示します。
a,
b,
e や
M の分子 (
a (1 -
e 2))
などは定数ですのであらかじめ計算しておくことが可能です。
Fig.1 ヒュベニの公式
定数
実際に距離を求めるのに必要なのは長半径
a および短半径
b だけです。
それぞれ以下のように定められています。実際には長半径
a と扁平率の逆数で定義されていますが、
ここでは短半径を与えられる定数として扱います。
以下は Wikipedia に記載
[6] の長半径・短半径・扁平率の逆数、
および統計処理言語 R による計算結果です。
| 名称 |
長半径 a (m) |
短半径 b (m) |
扁平率の逆数 |
(第一離心率 e) 2 |
a (1 - e 2) |
| ベッセル楕円体（旧日本測地系） |
6,377,397.155 |
6,356,079.000 000 |
299.152 813 000 |
0.006 674 360 610 282 97 |
6,334,832.106 632 54 |
| GRS80（世界測地系） |
6,378,137.000 |
6,356,752.314 140 |
298.257 222 101 |
0.006 694 380 023 011 88 |
6,335,439.327 083 17 |
| WGS84 (GPS) |
6,378,137.000 |
6,356,752.314 245 |
298.257 223 563 |
0.006 694 379 990 197 58 |
6,335,439.327 292 46 |
精度評価
上記公式の精度の評価をしたいと思います。
本当の距離が分かっている２地点との比較を行うのがもっとも正確に精度の評価ができそうですが、
なかなかそういうところがありません。
国土地理院のウェブ上で最初に触れた複雑な式での距離計算プログラムを公開しているので、
それとの比較は容易です。そこで例として使用されている東京・筑波間の距離を比較してみます。
| 東京 |
筑波 |
距離 (m) |
| 緯度 | 経度 |
緯度 | 経度 |
国土地理院 |
ヒュベニの公式 |
差 |
| 35.65500 | 139.74472 |
36.10056 | 140.09111 |
58,501.873 |
58,502.459 |
0.586 |
結果は 1m 以下の差となりました。十分な精度が出ていると言える思います。
もう少し距離の長い地点の比較をしてみます。
同じく東京から
福岡ドーム
（北緯 33.59532 度、東経 130.36208 度）への距離を比較してみます。
| 東京 |
福岡ドーム |
距離 (m) |
| 緯度 | 経度 |
緯度 | 経度 |
国土地理院 |
ヒュベニの公式 |
差 |
| 35.65500 | 139.74472 |
33.59532 | 130.36208 |
889,826.431 |
890,233.064 |
406.633 |
今度は 400 m 以上の差が出ました。しかし距離が 890 km ありますので、
そう考えるとわずか 400m とも言えます。これだけ離れて差が 1km 以下ですので、
両者は実用上十分同じと言えそうです。
距離が分かってしかも地図上にもはっきりと載っているものとして滑走路が使えるのではないかと思いつきました。成田空港の北側暫定滑走路は距離が 2180m ですので、
この距離を計算してみます。国土交通省のサイトの地図 [7] で 2180m
がどこからどこまでかおおよそ分かります。両端の部分はこの距離には入らないようです。
地形図では始点・終点の位置が分かりづらいですので Google Maps の航空写真上で調べると、
滑走路の北端はおよそ北緯 35.802739 度、東経 140.380034 度、
南端は北緯 35.785796 度、東経 140.392265 度です。
| 北側滑走路北端 |
北側滑走路南端 |
距離 (m) |
| 緯度 | 経度 |
緯度 | 経度 |
ヒュベニの公式 |
実際の距離 |
差 |
| 35.802739 | 140.380034 |
35.785796 | 140.392265 |
2180.948470 |
2180.000000 |
0.948470 |
ギリギリ差は 1m 以内になりました。2 km の滑走路、
しかも Google Maps 航空写真からの緯度経度算出でここまで近くなりましたので、
非常に正確に求まっているといえそうです。
実装
R および Java による実装を公開します。プログラムは「MIT ライセンス」に従うことを条件に自由に利用できます。
統計処理言語 R 版
R による実装です。定数を関数内で計算しているため効率がよくありません。
もともと下記の Java プログラムのデバッグ用の実装なのでこれでよしとします。
deg2rad <- function(x){ x * pi / 180 }
dist <- function(x1, y1, x2, y2, a = 6378137.0, b = 6356752.314140){
dy <- deg2rad(y1 - y2)
dx <- deg2rad(x1 - x2)
my <- deg2rad((y1 + y2) / 2)
e2 <- (a^2 - b^2) / a^2
Mnum <- a * (1 - e2)
W <- sqrt(1 - e2 * sin(my)^2)
M <- Mnum / W^3
N <- a / W
d <- sqrt((dy * M)^2 + (dx * N * cos(my))^2)
d
}
|
実行例。通常は緯度、経度の順で書きますが、この関数は経度 (x)、
緯度 (y) の順番ですのでご注意ください。
> dist(140.09111,36.10056,139.74472,35.65500)
[1] 58502.4589312406
|
Java 版
public class Coords {
public static final double BESSEL_A = 6377397.155;
public static final double BESSEL_E2 = 0.00667436061028297;
public static final double BESSEL_MNUM = 6334832.10663254;
public static final double GRS80_A = 6378137.000;
public static final double GRS80_E2 = 0.00669438002301188;
public static final double GRS80_MNUM = 6335439.32708317;
public static final double WGS84_A = 6378137.000;
public static final double WGS84_E2 = 0.00669437999019758;
public static final double WGS84_MNUM = 6335439.32729246;
public static final int BESSEL = 0;
public static final int GRS80 = 1;
public static final int WGS84 = 2;
public static double deg2rad(double deg){
return deg * Math.PI / 180.0;
}
public static double calcDistHubeny(double lat1, double lng1,
double lat2, double lng2,
double a, double e2, double mnum){
double my = deg2rad((lat1 + lat2) / 2.0);
double dy = deg2rad(lat1 - lat2);
double dx = deg2rad(lng1 - lng2);
double sin = Math.sin(my);
double w = Math.sqrt(1.0 - e2 * sin * sin);
double m = mnum / (w * w * w);
double n = a / w;
double dym = dy * m;
double dxncos = dx * n * Math.cos(my);
return Math.sqrt(dym * dym + dxncos * dxncos);
}
public static double calcDistHubeny(double lat1, double lng1,
double lat2, double lng2){
return calcDistHubeny(lat1, lng1, lat2, lng2,
GRS80_A, GRS80_E2, GRS80_MNUM);
}
public static double calcDistHubery(double lat1, double lng1,
double lat2, double lng2, int type){
switch(type){
case BESSEL:
return calcDistHubeny(lat1, lng1, lat2, lng2,
BESSEL_A, BESSEL_E2, BESSEL_MNUM);
case WGS84:
return calcDistHubeny(lat1, lng1, lat2, lng2,
WGS84_A, WGS84_E2, WGS84_MNUM);
default:
return calcDistHubeny(lat1, lng1, lat2, lng2,
GRS80_A, GRS80_E2, GRS80_MNUM);
}
}
public static void main(String[] args){
System.out.println("Coords Test Program");
double lat1, lng1, lat2, lng2;
lat1 = Double.parseDouble(args[0]);
lng1 = Double.parseDouble(args[1]);
lat2 = Double.parseDouble(args[2]);
lng2 = Double.parseDouble(args[3]);
double d = calcDistHubeny(lat1, lng1, lat2, lng2);
System.out.println("Distance = " + d + " m");
}
}
|
実行例。Java 版は緯度、経度の順になっています。
$ java Coords 36.10056 140.09111 35.65500 139.74472
Coords Test Program
Distance = 58502.45893124115 m
|
参考文献
[1] 国土地理院 測地部
http://vldb.gsi.go.jp/sokuchi/
[2] 国土地理院 測地部の計算式のページ
http://vldb.gsi.go.jp/sokuchi/surveycalc/algorithm/
[3] 国土地理院 測地部の計算式のページ その２
http://vldb.gsi.go.jp/sokuchi/surveycalc/algorithm/xy2st/xy2st.htm
[4] カシミール 3D
http://www.kashmir3d.com/
[5] カシミール 3D マニュアルページ ヒュベニの距離計算式
http://www.kashmir3d.com/kash/manual/std_siki.htm
[6] Wikipedia 測地学
http://ja.wikipedia.org/wiki/測地学
[7] 国土交通省の成田国際空港の施設計画
http://www.mlit.go.jp/hakusyo/mlit/h16/hakusho/h17/html/g2052120.html
Copyright © 2007-2012
やまだらけ 無断転載禁止 Since 2007/7/26
このウェブサイトについて
当ウェブサイトで使用しているテキスト、情報、画像等の権利は、明記してある場合を除き全て
やまだらけ に帰属します。無断転載は一切禁止します。
当データベースの情報に限り、複製・二次利用等は出典を明記すれば自由に行うことができます。
