GhPythonで平面直角座標を緯度経度に変換する

#はじめに
[**「緯度経度と平面直角座標の相互変換をPythonで実装する」**](https://qiita.com/sw1227/items/e7a590994ad7dcd0e8ab)を参考に Grasshopper の GhPython で平面直角座標を緯度経度に変換してみた。**@sw1227** さん大変参考になりました！ありがとうございます！
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/433170/5ef6e946-8e0f-85ab-9d85-12eb605d0b6b.png)

#変えたところ
GhPython では numpy を使えないので・・・
　・math で使える関数に置き換えた
　・np.array をリストに置き換えた
　・np.dot を for文で書き直した
といったあたりを変更

#GhPythonで作成
Type hint はすべて float にする
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/433170/b3fb6bbc-f9a6-7b36-9a16-d3176d51cab2.png)

コードはこちら

```python
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))

```

#参考
[緯度経度と平面直角座標の相互変換をPythonで実装する](https://qiita.com/sw1227/items/e7a590994ad7dcd0e8ab)
[緯度経度と平面直角座標の相互変換を実装するための数式](https://qiita.com/sw1227/items/e7a590994ad7dcd0e8ab)