Created
February 24, 2019 04:16
-
-
Save komasaru/f4693296ad38b562bd90067858ef1560 to your computer and use it in GitHub Desktop.
Ruby script to convert ECEF to WGS84(BLH) coordinate.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#! /usr/local/bin/ruby | |
#--------------------------------------------------------------------------------- | |
# ECEF -> BLH 変換 | |
# : ECEF(Earth Centered Earth Fixed; 地球中心・地球固定直交座標系)座標を | |
# WGS84 の緯度(Latitude)/経度(Longitude)/楕円体高(Height)に変換する。 | |
# | |
# Date Author Version | |
# 2019.02.19 mk-mode.com 1.00 新規作成 | |
# | |
# Copyright(C) 2019 mk-mode.com All Rights Reserved. | |
#--------------------------------------------------------------------------------- | |
# | |
class ConvEcef2Blh | |
USAGE = "[USAGE] ./ecef2blh.rb X Y Z" | |
PI_180 = Math::PI / 180.0 | |
# WGS84 座標パラメータ | |
A = 6378137.0 # a(地球楕円体長半径(赤道面平均半径)) | |
ONE_F = 298.257223563 # 1 / f(地球楕円体扁平率=(a - b) / a) | |
B = A * (1.0 - 1.0 / ONE_F) # b(地球楕円体短半径) | |
E2 = (1.0 / ONE_F) * (2 - (1.0 / ONE_F)) | |
# e^2 = 2 * f - f * f | |
# = (a^2 - b^2) / a^2 | |
ED2 = E2 * A * A / (B * B) # e'^2= (a^2 - b^2) / b^2 | |
# Initialization | |
# * コマンドライン引数の取得 | |
def initialize | |
if ARGV.size < 3 | |
puts USAGE | |
exit | |
end | |
@ecef = ARGV[0, 3].map(&:to_f) | |
end | |
# Execution | |
def exec | |
puts "ECEF: X = %12.3fm" % @ecef[0] | |
puts " Y = %12.3fm" % @ecef[1] | |
puts " Z = %12.3fm" % @ecef[2] | |
blh = ecef2blh(@ecef) | |
puts "--->" | |
puts "BLH: LATITUDE(BETA) = %12.8f°" % blh[0] | |
puts " LONGITUDE(LAMBDA) = %12.8f°" % blh[1] | |
puts " HEIGHT = %12.3fm" % blh[2] | |
rescue => e | |
msg = "[#{e.class}] #{e.message}\n" | |
e.backtrace.each { |tr| msg << "\t#{tr}\n" } | |
$stderr.puts msg | |
exit 1 | |
end | |
private | |
# ECEF -> BLH 変換 | |
# | |
# @param ecef: ECEF 座標 [x, y, z] | |
# @return blh: BLH 座標 [latitude, longitude, height] | |
def ecef2blh(ecef) | |
x, y, z = ecef | |
n = lambda do |x| | |
s = Math.sin(x * PI_180) | |
A / Math.sqrt(1.0 - E2 * s * s) | |
end | |
p = Math.sqrt(x * x + y * y) | |
theta = Math.atan2(z * A, p * B) / PI_180 | |
lat = Math.atan2( | |
z + ED2 * B * Math.sin(theta * PI_180)**3, | |
p - E2 * A * Math.cos(theta * PI_180)**3 | |
) / PI_180 | |
lon = Math.atan2(y, x) / PI_180 | |
ht = (p / Math.cos(lat * PI_180)) - n.call(lat) | |
return [lat, lon, ht] | |
rescue => e | |
raise | |
end | |
end | |
ConvEcef2Blh.new.exec if __FILE__ == $0 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment