Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Ruby script to calculate a central angle betweetn 2 points on earth.
#! /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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.