緯度経度と平面直角座標の相互変換をPythonで実装する

# 1. はじめに
ここでは緯度経度 <-> 平面直角座標の実装（Python）のみ記す。

変換式の詳細や平面直角座標系の解説は、長くなるので別途ブログに書いた。
[緯度経度と平面直角座標の相互変換を実装するための数式 - sw1227's diary](https://sw1227.hatenablog.com/entry/2018/11/30/200702)

以下ではNumpyをimport済みとする。

```python
import numpy as np
```

# 2. 平面直角座標から経緯度への変換
中間的な変数が多数登場してややこしいが、それらの依存関係は下図のように整理できる。
[数式および図の詳細はこちらを参照](https://sw1227.hatenablog.com/entry/2018/11/30/200702#2-%E5%B9%B3%E9%9D%A2%E7%9B%B4%E8%A7%92%E5%BA%A7%E6%A8%99%E3%81%8B%E3%82%89%E7%B5%8C%E7%B7%AF%E5%BA%A6%E3%81%B8%E3%81%AE%E5%A4%89%E6%8F%9B%E5%BC%8F)のこと。

本記事の実装における変数名や式の順番はそちらの数式・図に準拠している。

![20181128210915.jpg](https://qiita-image-store.s3.amazonaws.com/0/58850/cf0c4e81-0894-949f-9b66-e8bba1546006.jpeg)

## 2.1. 実装

```python
def calc_lat_lon(x, y, phi0_deg, lambda0_deg):
    """ 平面直角座標を緯度経度に変換する
    - input:
        (x, y): 変換したいx, y座標[m]
        (phi0_deg, lambda0_deg): 平面直角座標系原点の緯度・経度[度]（分・秒でなく小数であることに注意）
    - output:
        latitude:  緯度[度]
        longitude: 経度[度]
        * 小数点以下は分・秒ではないことに注意
    """
    # 平面直角座標系原点をラジアンに直す
    phi0_rad = np.deg2rad(phi0_deg)
    lambda0_rad = np.deg2rad(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 np.array([A0, A1, A2, A3, A4, A5])
  
    def beta_array(n):
        b0 = np.nan # dummy
        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 np.array([b0, b1, b2, b3, b4, b5])
    
    def delta_array(n):
        d0 = np.nan # dummy
        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 np.array([d0, d1, d2, d3, d4, d5, d6])
    
    # 定数 (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 + np.dot(A_array[1:], np.sin(2*phi0_rad*np.arange(1,6))) )
    
    # (3) xi, etaの計算
    xi = (x + S_) / A_
    eta = y / A_

    # (4) xi', eta'の計算
    xi2 = xi - np.sum(np.multiply(beta_array[1:], 
                                  np.multiply(np.sin(2*xi*np.arange(1,6)),
                                              np.cosh(2*eta*np.arange(1,6)))))
    eta2 = eta - np.sum(np.multiply(beta_array[1:],
                                   np.multiply(np.cos(2*xi*np.arange(1,6)),
                                               np.sinh(2*eta*np.arange(1,6)))))
    
    # (5) chiの計算
    chi = np.arcsin( np.sin(xi2)/np.cosh(eta2) ) # [rad]
    latitude = chi + np.dot(delta_array[1:], np.sin(2*chi*np.arange(1, 7))) # [rad]

    # (6) 緯度(latitude), 経度(longitude)の計算
    longitude = lambda0_rad + np.arctan( np.sinh(eta2)/np.cos(xi2) ) # [rad]
    
    # ラジアンを度になおしてreturn
    return np.rad2deg(latitude), np.rad2deg(longitude) # [deg]
```

## 2.2. 簡易テスト
[国土地理院の換算サービス](http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/xy2blf.html)で変換したものと比較する。

計算結果を度,分,秒表記で比較できるよう、以下のような関数を用意しておく。

```python
def show_angle(deg):
    """ 小数点の角度[deg]を度,分,秒で表記 """
    d = int(np.floor(deg))
    m = int(np.floor((deg%1) * 60))
    s = ( ((deg%1)*60) % 1 ) * 60
    return """ {0}°{1:02d}'{2}" """.format(d, m, s) # 分は10の位を0埋めする
```

以下のように変換してみると、数値計算上のわずかな誤差を除き、国土地理院の換算サービスの結果と一致するはず。

```python
lat, lon = calc_lat_lon(11573.375, 22694.980, 33., 131.)
print(" latitude: {0}".format(show_angle(lat)))
print("longitude: {0}".format(show_angle(lon)))
# <<実行結果>>>
#  latitude:  33°06'14.856642798" 
# longitude:  131°14'35.3709252452" 
```

# 3. 経緯度から平面直角座標への変換
こちらの依存関係も下図のように整理できる。
[数式および図の詳細はこちらを参照](https://sw1227.hatenablog.com/entry/2018/11/30/200702#3-%E7%B5%8C%E7%B7%AF%E5%BA%A6%E3%81%8B%E3%82%89%E5%B9%B3%E9%9D%A2%E7%9B%B4%E8%A7%92%E5%BA%A7%E6%A8%99%E3%81%B8%E3%81%AE%E5%A4%89%E6%8F%9B%E5%BC%8F)のこと。

![20181129233652.jpg](https://qiita-image-store.s3.amazonaws.com/0/58850/ef4eb0cf-5e5c-f1ff-4797-78d7dfd15b78.jpeg)

## 3.1. 実装
```python
def calc_xy(phi_deg, lambda_deg, phi0_deg, lambda0_deg):
    """ 緯度経度を平面直角座標に変換する
    - input:
        (phi_deg, lambda_deg): 変換したい緯度・経度[度]（分・秒でなく小数であることに注意）
        (phi0_deg, lambda0_deg): 平面直角座標系原点の緯度・経度[度]（分・秒でなく小数であることに注意）
    - output:
        x: 変換後の平面直角座標[m]
        y: 変換後の平面直角座標[m]
    """
    # 緯度経度・平面直角座標系原点をラジアンに直す
    phi_rad = np.deg2rad(phi_deg)
    lambda_rad = np.deg2rad(lambda_deg)
    phi0_rad = np.deg2rad(phi0_deg)
    lambda0_rad = np.deg2rad(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 np.array([A0, A1, A2, A3, A4, A5])
  
    def alpha_array(n):
        a0 = np.nan # 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 np.array([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] # [m]
    S_ = ( (m0*a)/(1.+n) )*( A_array[0]*phi0_rad + np.dot(A_array[1:], np.sin(2*phi0_rad*np.arange(1,6))) ) # [m]
    
    # (3) lambda_c, lambda_sの計算
    lambda_c = np.cos(lambda_rad - lambda0_rad)
    lambda_s = np.sin(lambda_rad - lambda0_rad)
    
    # (4) t, t_の計算
    t = np.sinh( np.arctanh(np.sin(phi_rad)) - ((2*np.sqrt(n)) / (1+n))*np.arctanh(((2*np.sqrt(n)) / (1+n)) * np.sin(phi_rad)) )
    t_ = np.sqrt(1 + t*t)

    # (5) xi', eta'の計算
    xi2  = np.arctan(t / lambda_c) # [rad]
    eta2 = np.arctanh(lambda_s / t_)
    
    # (6) x, yの計算
    x = A_ * (xi2 + np.sum(np.multiply(alpha_array[1:],
                                       np.multiply(np.sin(2*xi2*np.arange(1,6)),
                                                   np.cosh(2*eta2*np.arange(1,6)))))) - S_ # [m]
    y = A_ * (eta2 + np.sum(np.multiply(alpha_array[1:],
                                        np.multiply(np.cos(2*xi2*np.arange(1,6)),
                                                    np.sinh(2*eta2*np.arange(1,6)))))) # [m]
    # return
    return x, y # [m]
```

## 3.2. 簡易テスト
[国土地理院の換算サービス](http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/bl2xyf.html)と比較し、一致することを確認した。

```python
x, y = calc_xy(36.103774791666666, 140.08785504166664, 36., 139+50./60)
print("x, y = ({0}, {1})".format(x, y))
# <<実行結果>>
# x, y = (11543.6883215, 22916.2435543)
```

# 4. 参考・リンク
- [緯度経度と平面直角座標の相互変換を実装するための数式 - sw1227's diary](https://sw1227.hatenablog.com/entry/2018/11/30/200702)
- [この記事の関数を実装したJupyter Notebook](http://nbviewer.jupyter.org/github/nabewari/PRCS2WGS/blob/master/coordinate_transformation.ipynb)
- [平面直角座標を換算して経緯度、子午線収差角及び縮尺係数を求める計算 - 国土地理院](http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/algorithm/xy2bl/xy2bl.htm)
- [経緯度を換算して平面直角座標、子午線収差角及び縮尺係数を求める計算 - 国土地理院](https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/algorithm/bl2xy/bl2xy.htm)
- [わかりやすい平面直角座標系](http://www.gsi.go.jp/sokuchikijun/jpc.html)
- [平面直角座標系（平成十四年国土交通省告示第九号）](http://www.gsi.go.jp/LAW/heimencho.html)
- [日本の測地系](http://www.gsi.go.jp/sokuchikijun/datum-main.html)