はじめに
「緯度経度と平面直角座標の相互変換をPythonで実装する」を参考に Grasshopper の GhPython で平面直角座標を緯度経度に変換してみた。@sw1227 さん大変参考になりました！ありがとうございます！
変えたところ
GhPython では numpy を使えないので・・・
・math で使える関数に置き換えた
・np.array をリストに置き換えた
・np.dot を for文で書き直した
といったあたりを変更
GhPythonで作成
Type hint はすべて float にする
コードはこちら
import math
""" 平面直角座標を緯度経度に変換する
- input:
(x, y): 変換したいx, y座標[m]
(phi0_deg, lambda0_deg): 平面直角座標系原点の緯度・経度[度]（分・秒でなく小数であることに注意）
- output:
lat: 緯度[度]
lon: 経度[度]
"""
# 平面直角座標系原点をラジアンに直す
phi0_rad = math.radians(phi0_deg)
lambda0_rad = math.radians(lambda0_deg)
# 補助関数
def A_array(n):
A0 = 1 + (n**2)/4. + (n**4)/64.
A1 = - (3./2)*( n - (n**3)/8. - (n**5)/64. )
A2 = (15./16)*( n**2 - (n**4)/4. )
A3 = - (35./48)*( n**3 - (5./16)*(n**5) )
A4 = (315./512)*( n**4 )
A5 = -(693./1280)*( n**5 )
return [A0, A1, A2, A3, A4, A5]
def beta_array(n):
b0 = None
b1 = (1./2)*n - (2./3)*(n**2) + (37./96)*(n**3) - (1./360)*(n**4) - (81./512)*(n**5)
b2 = (1./48)*(n**2) + (1./15)*(n**3) - (437./1440)*(n**4) + (46./105)*(n**5)
b3 = (17./480)*(n**3) - (37./840)*(n**4) - (209./4480)*(n**5)
b4 = (4397./161280)*(n**4) - (11./504)*(n**5)
b5 = (4583./161280)*(n**5)
return [b0, b1, b2, b3, b4, b5]
def delta_array(n):
d0 = None
d1 = 2.*n - (2./3)*(n**2) - 2.*(n**3) + (116./45)*(n**4) + (26./45)*(n**5) - (2854./675)*(n**6)
d2 = (7./3)*(n**2) - (8./5)*(n**3) - (227./45)*(n**4) + (2704./315)*(n**5) + (2323./945)*(n**6)
d3 = (56./15)*(n**3) - (136./35)*(n**4) - (1262./105)*(n**5) + (73814./2835)*(n**6)
d4 = (4279./630)*(n**4) - (332./35)*(n**5) - (399572./14175)*(n**6)
d5 = (4174./315)*(n**5) - (144838./6237)*(n**6)
d6 = (601676./22275)*(n**6)
return [d0, d1, d2, d3, d4, d5, d6]
def show_angle(deg):
""" 小数点の角度[deg]を度,分,秒で表記 """
d = int(math.floor(deg))
m = int(math.floor((deg%1) * 60))
s = ( ((deg%1)*60) % 1 ) * 60
return """ {0}°{1:02d}'{2}" """.format(d, m, s) # 分は10の位を0埋めする
# 定数 (a, F: 世界測地系-測地基準系1980（GRS80）楕円体)
m0 = 0.9999
a = 6378137.
F = 298.257222101
# (1) n, A_i, beta_i, delta_iの計算
n = 1. / (2*F - 1)
A_array = A_array(n)
beta_array = beta_array(n)
delta_array = delta_array(n)
# (2), S, Aの計算
A_ = ( (m0*a)/(1.+n) )*A_array[0]
S_ = ( (m0*a)/(1.+n) ) * (A_array[0]*phi0_rad + sum([A_array[j]*(math.sin(2*phi0_rad*j)) for j in range(1,6)]))
# (3) xi, etaの計算
xi = (x + S_) / A_
eta = y / A_
# (4) xi', eta'の計算
xi2 = xi - sum([beta_array[j]*math.sin(2*xi*j)*math.cosh(2*eta*j) for j in range(1,6)])
eta2 = eta - sum([beta_array[j]*math.cos(2*xi*j)*math.sinh(2*eta*j) for j in range(1,6)])
# (5) chiの計算
chi = math.asin( math.sin(xi2)/math.cosh(eta2) )
# (6) 緯度(latitude), 経度(longitude)の計算
latitude = chi + sum([delta_array[j] * math.sin(2*chi*j) for j in range(1,7)])
longitude = lambda0_rad + math.atan( math.sinh(eta2)/math.cos(xi2) )
lat = show_angle(math.degrees(latitude))
lon = show_angle(math.degrees(longitude))