Ruby - 地球楕円体上の2地点間中心角の計算!

更新日時:


地球楕円体上の2地点と地球中心がなす中心角を Ruby で計算してみました。

単純に2点の直交座標を計算後、2ベクトルのなす角を計算するだけ。(他にも算出方法はありますが)

0. 前提条件

  • LMDE 3 (Linux Mint Debian Edition 3; 64bit) での作業を想定。
  • Ruby 2.6.4 での作業を想定。

1. 計算手順

  1. 地点1の緯度・経度を , 地点2の緯度・経度を とする。
  2. 地点1の緯度 ,経度 , 地点2の緯度 ,経度 を直交座標(ECEF(Earth Centerd Earth Fixed)座標) に変換する。
  3. 地球中心から地点1, 地球中心から地点2という2つのベクトルとみなし、2ベクトルのなす角 を計算する。

ECEF 座標の計算については過去記事を参照。(今回の計算のキモとなる部分)

2. Ruby スクリプトの作成

  • 地球楕円体は GRS-80, WGS-84, BESSEL を考慮。
  • Shebang ストリング(1行目)では、フルパスでコマンド指定している。(当方の慣習

File: calc_central_angle.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
#! /usr/local/bin/ruby
#---------------------------------------------------------------------------------
# 2地点の緯度・経度から中心角を計算
# * 2地点の緯度(Beta)/経度(Lambda)/楕円体高(Height)から
#   地球中心から2地点への線分のなす中心角を求める。
# * CalcCentralAngle クラスが実行クラス
# * CentralAngle クラスが計算の本質部分
#
#   Date          Author          Version
#   2019.08.29    mk-mode.com     1.00 新規作成
#
# Copyright(C) 2019 mk-mode.com All Rights Reserved.
#---------------------------------------------------------------------------------
# コマンドライン引数:
#   第1: 緯度-1
#   第2: 経度-1
#   第3: 緯度-2
#   第4: 経度-2
#   第5: 地球楕円体("GRS80"|"WGS84"|"BESSEL") (optional)
#   * 緯度は、東経が +, 西経が - で指定。
#   * 経度は、北緯が +, 南緯が - で指定。
#   * 2地点の高度(標高)は 0m に限定。
#   * 地球楕円体省略時は "GRS80" とみなす。
#---------------------------------------------------------------------------------
#
class CalcCentralAngle
  USAGE = "[USAGE] ./calc_central_angle.rb LAT_1 LNG_1 LAT_2 LNG_2 [ELLIPSOID]"

  # Initialization
  # * コマンドライン引数の取得
  def initialize
    if ARGV.size < 4
      puts USAGE
      exit
    end
    args = ARGV[0, 4].map(&:to_f)
    @pt_1 = args[0, 2]
    @pt_2 = args[2, 2]
    @el   = "GRS80"
    @el   = ARGV[4] if ARGV[4]
  end

  # Execution
  def exec
    puts "      POINT-1: %8.4f° %9.4f°" % @pt_1
    puts "      POINT-2: %8.4f° %9.4f°" % @pt_2
    puts "    ELLIPSOID: #{@el}"
    obj = CentralAngle.new(*((@pt_1 + @pt_2) << @el))
    ca = obj.calc_central_angle
    puts "CENTRAL-ANGLE: %8.4f°" % ca
  rescue => e
    msg = "[#{e.class}] #{e.message}\n"
    e.backtrace.each { |tr| msg << "\t#{tr}\n" }
    $stderr.puts msg
    exit 1
  end
end

class CentralAngle
  # 地球楕円体(長半径, 扁平率の逆数)
  ELLIPSOID = {
    "GRS80"  => [6_378_137.0  , 298.257_222_101],
    "WGS84"  => [6_378_137.0  , 298.257_223_563],
    "BESSEL" => [6_377_397.155, 299.152_813_000]
  }
  PI_180 = Math::PI / 180.0

  def initialize(lat_1, lng_1, lat_2, lng_2, ellipsoid="GRS80")
    @lat_1, @lng_1 = lat_1, lng_1
    @lat_2, @lng_2 = lat_2, lng_2
    el  = ELLIPSOID[ellipsoid]
    @a  = el[0]
    f   = 1 / el[1]
    @e2 = f * (2 - f)
  end

  # 中心角の計算
  # * 2ベクトルのなす角とみなして計算
  #
  # @return: 中心角 (Unit: °)
  def calc_central_angle
    x_1, y_1, z_1 = blh2ecef(@lat_1, @lng_1)
    x_2, y_2, z_2 = blh2ecef(@lat_2, @lng_2)
    return Math.acos(
      (x_1 * x_2 + y_1 * y_2 + z_1 * z_2) /
      (Math.sqrt(x_1 * x_1 + y_1 * y_1 + z_1 * z_1) *
       Math.sqrt(x_2 * x_2 + y_2 * y_2 + z_2 * z_2))
    ) / PI_180
  rescue => e
    raise
  end

private
  # BLH座標 -> ECEF直交座標 変換
  #
  # @param        lat: BLH  座標 緯度 (Unit: °)
  # @param        lng:           経度 (Unit: °)
  # @param         ht:           高度 (Unit: m)
  # @return [x, y, z]: ECEF 座標      (Unit: m)
def blh2ecef(lat, lng, ht=0.0)
    s = Math.sin(lat * PI_180)
    n = @a / Math.sqrt(1.0 - @e2 * s * s)
    x = (n + ht) \
      * Math.cos(lat * PI_180) \
      * Math.cos(lng * PI_180)
    y = (n + ht) \
      * Math.cos(lat * PI_180) \
      * Math.sin(lng * PI_180)
    z = (n * (1.0 - @e2) + ht) \
      * Math.sin(lat * PI_180)
    return [x, y, z]
  rescue => e
    raise
  end
end

CalcCentralAngle.new.exec if __FILE__ == $0

3. Ruby スクリプトの実行

地点1の緯度・経度、地点2の緯度・経度、オプションで地球楕円体(GRS80|WGS84|BESSEL)を指定して実行。

$ ./calc_central_angle.rb 35.0 135.0 40.0 140.0
      POINT-1:  35.0000°  135.0000°
      POINT-2:  40.0000°  140.0000°
    ELLIPSOID: GRS80
CENTRAL-ANGLE:   6.3792°

$./calc_central_angle.rb 35.0 135.0 40.0 140.0 WGS84
      POINT-1:  35.0000°  135.0000°
      POINT-2:  40.0000°  140.0000°
    ELLIPSOID: WGS84
CENTRAL-ANGLE:   6.3792°

$./calc_central_angle.rb 35.0 135.0 40.0 140.0 BESSEL
      POINT-1:  35.0000°  135.0000°
      POINT-2:  40.0000°  140.0000°
    ELLIPSOID: BESSEL
CENTRAL-ANGLE:   6.3792°

$ ./calc_central_angle.rb 35.0 135.0 40.0 135.0
      POINT-1:  35.0000°  135.0000°
      POINT-2:  40.0000°  135.0000°
    ELLIPSOID: GRS80
CENTRAL-ANGLE:   4.9912°

$ ./calc_central_angle.rb 0.0 135.0 0.0 140.0
      POINT-1:   0.0000°  135.0000°
      POINT-2:   0.0000°  140.0000°
    ELLIPSOID: GRS80
CENTRAL-ANGLE:   5.0000°

$./calc_central_angle.rb 0.0 -45.0 0.0 135.0
      POINT-1:   0.0000°  -45.0000°
      POINT-2:   0.0000°  135.0000°
    ELLIPSOID: GRS80
CENTRAL-ANGLE: 180.0000°

$ ./calc_central_angle.rb 0.0 -135.0 0.0 135.0
      POINT-1:   0.0000° -135.0000°
      POINT-2:   0.0000°  135.0000°
    ELLIPSOID: GRS80
CENTRAL-ANGLE:  90.0000°

$ ./calc_central_angle.rb 90.0 135.0 90.0 135.0
      POINT-1:  90.0000°  135.0000°
      POINT-2:  90.0000°  135.0000°
    ELLIPSOID: GRS80
CENTRAL-ANGLE:   0.0000°

当初、地図上の正確な位置計算に利用しようと考えていたが、結局、別の方法を採用することにしたので、今回の方法は使い道がなくなった。

しかし、高校数学の復習にはなった。

以上。

 

Sponsored Link

 

コメントする