Ruby - 地球上の2点間の距離をほぼ正確に計算！

Updated: 2011-10-28 00:20

最近、 Ruby での UNIX MBOX メールヘッダの検証をやってきましたが、今回は気分転換に違った話題です。

それは「Rubyを使って地球上の2点間の距離を出来るだけ正確に計算してみよう」というものです。 ※きっかけは、地元の原子力発電所と自宅の距離を正確に知りたかったから。

ご存知のように地球は完全な球体ではありません。赤道方向の半径が極方向の半径より若干大きいです。

地球が完全な球体なら、高校数学程度の知識で計算できます。ちなみに、球面三角法による算出方法は以下のとおり。

DISTANCE_EARTH_1

正確に計算するには、まず正確な赤道半径・極半径が必要になります。測地学 - Wikipedia、測地系 - Wikipedia によると以下の3種類の測地系（測量の前提になる考え方）というものが存在し、それぞれに長半径（赤道半径）・短半径（極半径）が設定されています。

測地系	長半径（赤道半径）	短半径（極半径）
Bessel	6377397.155 m	6356079.000000 m
GRS 80	6378137.000 m	6356752.314140 m
WGS 84	6378137.000 m	6356752.314245 m

これらの測地系のうち、「GRS 80」と「WGS 84」が歴史も新しく精度も良いようで、中でも「GRS 80」が一番精度が良いようです。（と言っても、それぞれの誤差はわずかなものです）

次に、計算方法ですが一番正確な方法は国土地理院のサイトに掲載の方法のような気がします。

距離と方位角の計算計算式

しかし、非常に複雑で簡単に実現できそうな気がしません。

そこで、さらに調べてみると「ヒュベニの公式」というものが存在することが判明。

カシミール / 計算式

詳しくはここでは説明しません（できません）が、他のサイト等も参考にしてまとめると以下のような公式になります。

DISTANCE_EARTH_2

参考までに、以下にこの公式を利用した距離計算のRubyスクリプトを掲載しておきます。

Ruby スクリプト

ご利用の環境によっては微修正が必要かもしれません。

File: calc_dist.rb

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
# -*- coding: utf-8 -*-
# 
# ２地点間の距離を計算（２地点の緯度・経度から）
#
# date          name            version
# 2011.10.22    mk-mode.com     1.00 新規作成
#
# Copyright(C)2011 mk-mode.com All Rights Reserved.
#---------------------------------------------------------------------------------
# 引数 : [第１] 測地系 ( 0:BESSEL , 1:GRS80 , 2:WGS84 )
#        [第２] 緯度１ (  -90.00000000 ～  90.00000000 )
#        [第３] 経度１ ( -180.00000000 ～ 180.00000000 )
#        [第４] 緯度２ (  -90.00000000 ～  90.00000000 )
#        [第５] 経度２ ( -180.00000000 ～ 180.00000000 )
#---------------------------------------------------------------------------------
# 備考 : 1:GRS80 が最も正確
#        島根原発１号機 ( 35.5382 132.9998 )
#        自宅           ( 99.9999 179.9999 )
#      : 完全な球だと想定すると
#        D = R * acos( sin(y1) * sin(y2) + cos(y1) * cos(y2) * cos(x2-x1) )
#++

class CalcDist
  # 使用方法
  USAGE =  "USAGE : ruby calc_disk.rb Type Lat1 Lon1 Lat2 Lon2\n"
  USAGE << "Type       : 0:BESSEL, 1:GRS80, 2:WGS84\n"
  USAGE << "Lat1, Lat2 :  -90.00000000 -  [+]90.00000000\n"
  USAGE << "Lon1, Lon2 : -180.00000000 - [+]180.00000000"

  # 定数 ( ベッセル楕円体 ( 旧日本測地系 ) )
  BESSEL_R_X  = 6377397.155000 # 赤道半径
  BESSEL_R_Y  = 6356079.000000 # 極半径
  
  # 定数 ( GRS80 ( 世界測地系 ) )
  GRS80_R_X   = 6378137.000000 # 赤道半径
  GRS80_R_Y   = 6356752.314140 # 極半径
  
  # 定数 ( WGS84 ( GPS ) )
  WGS84_R_X   = 6378137.000000 # 赤道半径
  WGS84_R_Y   = 6356752.314245 # 極半径
  
  # 定数 ( 測地系 )
  MODE = [ "BESSEL", "GRS-80", "WGS-84" ]
  
  # [CLASS] 引数
  class Arg
    # 引数チェック
    def check_arg()
      begin
        # 存在チェック
        return false unless ARGV.length == 5
        
        # 測地系タイプチェック
        unless ARGV[0] =~ /^[012]$/
          return false
        end
        
        # 緯度チェック
        unless ARGV[1] =~ /^[+|-]?(\d|[1-8]\d|90)(\.\d{1,8})?$/ &&
               ARGV[3] =~ /^[+|-]?(\d|[1-8]\d|90)(\.\d{1,8})?$/
          return false
        end

        # 経度チェック
        unless ARGV[2] =~ /^[+|-]?(\d{1,2}|1[0-7]\d|180)(\.\d{1,8})?$/ &&
               ARGV[4] =~ /^[+|-]?(\d{1,2}|1[0-7]\d|180)(\.\d{1,8})?$/
          return false
        end

        return true
      rescue => e
        # エラーメッセージ
        str_msg = "[EXCEPTION][" + self.class.name + ".check_arg] " + e.to_s
        STDERR.puts( str_msg )
        exit! 1
      end
    end
  end

  # [CLASS] 計算
  class Calc
    def initialize(mode, lat_1, lon_1, lat_2, lon_2)
      @mode  = mode
      @lat_1 = lat_1
      @lon_1 = lon_1
      @lat_2 = lat_2
      @lon_2 = lon_2
    end

    # 距離計算
    def calc_dist
      begin
        # 指定測地系の赤道半径・極半径を設定
        case @mode
          when 0
            r_x = BESSEL_R_X
            r_y = BESSEL_R_Y
          when 1
            r_x = GRS80_R_X
            r_y = GRS80_R_Y
          when 2
            r_x = WGS84_R_X
            r_y = WGS84_R_Y
        end

        # 2点の経度の差を計算 ( ラジアン )
        a_x = @lon_1 * Math::PI / 180.0 - @lon_2 * Math::PI / 180.0

        # 2点の緯度の差を計算 ( ラジアン )
        a_y = @lat_1 * Math::PI / 180.0 - @lat_2 * Math::PI / 180.0

        # 2点の緯度の平均を計算
        p = ( @lat_1 * Math::PI / 180.0 + @lat_2 * Math::PI / 180.0 ) / 2.0

        # 離心率を計算
        e = Math::sqrt( ( r_x ** 2 - r_y ** 2 ) / ( r_x ** 2 ).to_f )

        # 子午線・卯酉線曲率半径の分母Wを計算
        w = Math::sqrt( 1 - ( e ** 2 ) * ( ( Math::sin( p ) ) ** 2 ) )

        # 子午線曲率半径を計算
        m = r_x * ( 1 - e ** 2 ) / ( w ** 3 ).to_f

        # 卯酉線曲率半径を計算
        n = r_x / w.to_f

        # 距離を計算
        d  = ( a_y * m ) ** 2
        d += ( a_x * n * Math.cos( p ) ) ** 2
        d  = Math::sqrt( d )

        # 地球を完全な球とみなした場合
        # ( 球面三角法 )
        # D = R * acos( sin(y1) * sin(y2) + cos(y1) * cos(y2) * cos(x2-x1) )
        d_1  = Math::sin( @lat_1 * Math::PI / 180.0 ) 
        d_1 *= Math::sin( @lat_2 * Math::PI / 180.0 )
        d_2  = Math::cos( @lat_1 * Math::PI / 180.0 )
        d_2 *= Math::cos( @lat_2 * Math::PI / 180.0 )
        d_2 *= Math::cos( @lon_2 * Math::PI / 180.0 - @lon_1 * Math::PI / 180.0 )
        d_0  = r_x * Math::acos( d_1 + d_2 ).to_f

        return [ d, d_0 ]
      rescue => e
        # エラーメッセージ
        str_msg = "[EXCEPTION][" + self.class.name + ".calc_dist] " + e.to_s
        STDERR.puts( str_msg )
        exit 1
      end
    end
  end

  #### MAIN ####
  if __FILE__ == $0
    # 引数チェック( エラーなら終了 )
    obj_arg = Arg.new
    unless obj_arg.check_arg
      # エラーの場合、終了
      puts USAGE
      exit!
    end

    # 引数取得
    mode  = ARGV[0].to_i
    lat_1 = ARGV[1].to_f
    lon_1 = ARGV[2].to_f
    lat_2 = ARGV[3].to_f
    lon_2 = ARGV[4].to_f

    # 距離計算
    obj_calc = Calc.new( mode, lat_1, lon_1, lat_2, lon_2 )
    dist = obj_calc.calc_dist

    # 結果出力
    puts "Mode        : #{MODE[mode]}"
    puts "Latitude (1): #{sprintf("%13.8f",lat_1)} degrees"
    puts "Longitude(1): #{sprintf("%13.8f",lon_1)} degrees"
    puts "Latitude (2): #{sprintf("%13.8f",lat_2)} degrees"
    puts "Longitude(2): #{sprintf("%13.8f",lon_2)} degrees"
    puts "Distance = #{dist[0]} m ( #{dist[1]} m )"
  end
end

使用例

上記のスクリプトに測地系(0:Bessel, 1:GRS80, 2:WGS84 )と地点１の緯度・経度、地点２の緯度・経度を引数に指定して実行します。（1:GRS80 が一番精度が高いと思います。）

緯度・経度は度で指定。Ａ度Ｂ分Ｃ秒は度に変換して指定。( A + B / 60 + C / 3600 )
南緯・西経はマイナスで指定。
小数点以下は8桁まで指定可能。

> ruby calc_dist.rb 1 35.5382 132.9998 35.47222 133.050556
Mode        : GRS-80
Latitude (1):   35.53820000 degrees
Longitude(1):  132.99980000 degrees
Latitude (2):   35.47222000 degrees
Longitude(2):  133.05055600 degrees
Distance = 8648.30331225593 m ( 8666.19473153661 m )

※上記は、島根原子力発電所１号機の位置（Googleマップより取得）と島根県庁の位置の距離です。
※計算結果の ( )内の数値は、地球が測地系の長半径を半径とする完全な球体とみなした場合の球面三角法による計算結果。

このように、地球が完全な球体であると仮定した場合との差は、今回の場合 18 メートル位となります。

ちなみに、島根原子力発電所１号機と自宅の距離は 8082.850806768164 mでした。 ※緯度・経度は自宅がばれるので公開しません。距離から自宅がある円周上に存在することは判明しますが・・・

以上。

Twitter Facebook Tumblr LinkedIn

mk-mode.com

Ruby - 地球上の2点間の距離をほぼ正確に計算！

Ruby スクリプト

使用例

Sponsored Link

Comments

You May Also Enjoy

Debian 12 (bookworm) - Ruby 3.3 インストール（ソースビルド）！ (2024-05-18)

2024年4月 - OS・ブラウザ別アクセス状況！ (2024-05-01)

2024年3月 - OS・ブラウザ別アクセス状況！ (2024-04-30)

2024年2月 - OS・ブラウザ別アクセス状況！ (2024-03-01)