Top | 日本の山岳データベース | GPS 入門 | 装備品 | 本棚 | 技術研究本部 | ブログ || 旧ウェブサイト
日本は山だらけ〜 技術研究本部 報告書 第一号 平成二十一年八月十五日 公開

二地点の緯度・経度からその距離を計算する

更新:コードのライセンスを MIT ライセンスに変更しました。(平成二十四年五月十日)

概要

各山のページに、 その山の近くにある別の山のリストを載せようと思い、 二地点の緯度と経度から距離を求める方法について調査しました。 ヒュベニ (Hubeny) の公式を用いると簡単な計算で精度よく距離が求まることがわかりました。 ヒュベニの公式ではいくつか定数が出てきますが、旧日本測地系、世界測地系、 GPS で利用されている測地系の 3 種類についての定数を調べ、算出しました。 Java および R のプログラムと合わせて研究成果を報告します。

緯度経度からの距離計算

国土地理院・測地部のウェブサイト [1] に測量に関する様々な情報が載っています。 特に [2] には様々な計算式が載っており、 緯度・経度から距離を計算する式 [3] もあります。 しかしここにあるのは緯度・経度を平面直角座標という座標系に変換した上で距離を計算で、 少し複雑です。 インターネットを調べるとヒュベニ (Hubeny) の公式と呼ばれている計算式が多く使われているようで、 3D 地図ソフトウェアとして有名なカシミール 3D [4] でも、 距離計算に使用されているようです。 カシミール 3D のマニュアルページ [5] にヒュベニの公式が紹介されています。 このヒュベニの公式は様々なところで紹介されていますが、 残念ながら1次ソース的な文献や論文は見つけることができませんでした。 しかし実際に国土地理院・測地部の距離計算プログラムと比較すると、 カシミール 3D の式とほぼ同じ値になることからこれをヒュベニの公式として引用します。 カシミール 3D ではベッセル楕円体に基づく旧日本測地系用のパラメータが用いられています。 そこで、国土地理院・測地部で紹介されている計算式を用い、 現在の正式な測地系(世界測地系または日本測地系2000)である GRS80 に基づくパラメータおよび、 GPS で使われる WGS84 に基づくパラメータの算出を行いました。 ただし GRS80 と WGS84 は実用上はほとんど差がありません。
▲ Top

ヒュベニの公式

ヒュベニの公式を Fig. 1 に示します。 a, b, eM の分子 (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

▲ Top

精度評価

上記公式の精度の評価をしたいと思います。 本当の距離が分かっている2地点との比較を行うのがもっとも正確に精度の評価ができそうですが、 なかなかそういうところがありません。 国土地理院のウェブ上で最初に触れた複雑な式での距離計算プログラムを公開しているので、 それとの比較は容易です。そこで例として使用されている東京・筑波間の距離を比較してみます。

東京 筑波 距離 (m)
緯度経度 緯度経度 国土地理院 ヒュベニの公式
35.65500139.74472 36.10056140.09111 58,501.873 58,502.459 0.586

結果は 1m 以下の差となりました。十分な精度が出ていると言える思います。 もう少し距離の長い地点の比較をしてみます。 同じく東京から福岡ドーム (北緯 33.59532 度、東経 130.36208 度)への距離を比較してみます。

東京 福岡ドーム 距離 (m)
緯度経度 緯度経度 国土地理院 ヒュベニの公式
35.65500139.74472 33.59532130.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.802739140.380034 35.785796140.392265 2180.948470 2180.000000 0.948470

ギリギリ差は 1m 以内になりました。2 km の滑走路、 しかも Google Maps 航空写真からの緯度経度算出でここまで近くなりましたので、 非常に正確に求まっているといえそうです。
▲ Top

実装

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
▲ Top

参考文献

[1] 国土地理院 測地部 http://vldb.gsi.go.jp/sokuchi/
[2] 国土地理院 測地部の計算式のページ http://vldb.gsi.go.jp/sokuchi/surveycalc/algorithm/
[3] 国土地理院 測地部の計算式のページ その2 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
▲ Top



Copyright © 2007-2012   やまだらけ   無断転載禁止   Since 2007/7/26   このウェブサイトについて
当ウェブサイトで使用しているテキスト、情報、画像等の権利は、明記してある場合を除き全て やまだらけ に帰属します。無断転載は一切禁止します。 当データベースの情報に限り、複製・二次利用等は出典を明記すれば自由に行うことができます。