はじめに
4年前には「GhPythonで平面直角座標を緯度経度に変換する」をやったけど今回は逆に緯度経度を平面直角座標に変換したくなった
前回と同じく「緯度経度と平面直角座標の相互変換をPythonで実装する」を参考にしてます。@sw1227 さんありがとうございます！
変えたところ
GhPython では numpy を使えないので・・・
・math で使える関数に置き換えた
・np.array をリストに置き換えた
・np.dot を for文で書き直した
といったあたりを変更
GhPythonで作成
Type hint はすべて float にする
コードはこちら
import math
""" 緯度経度を平面直角座標に変換する
- input:
(lat_deg, lon_deg): 変換したい緯度・経度[度]（分・秒でなく小数であることに注意）
(lat0_deg, lon0_deg): 平面直角座標系原点の緯度・経度[度]（分・秒でなく小数であることに注意）
- output:
x: 変換後の平面直角座標[m]
y: 変換後の平面直角座標[m]
"""
# 緯度経度・平面直角座標系原点をラジアンに直す
phi_rad = math.radians(lat_deg)
lambda_rad = math.radians(lon_deg)
phi0_rad = math.radians(lat0_deg)
lambda0_rad = math.radians(lon0_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 alpha_array(n):
a0 = None # dummy
a1 = (1./2)*n - (2./3)*(n**2) + (5./16)*(n**3) + (41./180)*(n**4) - (127./288)*(n**5)
a2 = (13./48)*(n**2) - (3./5)*(n**3) + (557./1440)*(n**4) + (281./630)*(n**5)
a3 = (61./240)*(n**3) - (103./140)*(n**4) + (15061./26880)*(n**5)
a4 = (49561./161280)*(n**4) - (179./168)*(n**5)
a5 = (34729./80640)*(n**5)
return [a0, a1, a2, a3, a4, a5]
# 定数 (a, F: 世界測地系-測地基準系1980（GRS80）楕円体)
m0 = 0.9999
a = 6378137.
F = 298.257222101
# (1) n, A_i, alpha_iの計算
n = 1. / (2*F - 1)
A_array = A_array(n)
alpha_array = alpha_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) lambda_c, lambda_sの計算
lambda_c = math.cos(lambda_rad - lambda0_rad)
lambda_s = math.sin(lambda_rad - lambda0_rad)
# (4) t, t_の計算
t = math.sinh( math.atanh(math.sin(phi_rad)) - ((2*math.sqrt(n)) / (1+n))*math.atanh(((2*math.sqrt(n)) / (1+n)) * math.sin(phi_rad)) )
t_ = math.sqrt(1 + t*t)
# (5) xi', eta'の計算
xi2 = math.atan(t / lambda_c) # [rad]
eta2 = math.atanh(lambda_s / t_)
# (6) x, yの計算
x = A_ * (xi2 + sum(
alpha_array[n] * math.sin(2 * xi2 * n) * math.cosh(2 * eta2 * n)
for n in range(1, 6))) - S_
y = A_ * (eta2 + sum(
alpha_array[n] * math.cos(2 * xi2 * n) * math.sinh(2 * eta2 * n)
for n in range(1, 6)))