#!/usr/bin/env python # coding: utf-8 # 平面直角座標と緯度経度の変換 # ===== # *** # ## 0. 参考 # # - [計算式 - 国土地理院](http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/algorithm/xy2bl/xy2bl.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) # ## 1. 計算式（平面直角座標 -> 緯度経度） # # ### 1.1. 注意 # [国土地理院の計算式](http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/algorithm/xy2bl/xy2bl.htm)はラジアンと度がごっちゃになっていて曖昧なので微調整した。 # 以下に出てくる三角関数はすべてラジアンを受け取り、逆三角関数はすべてラジアンを返すものとする。 # # # ### 1.2. 記号 # # - $x, y$ : 緯度経度に変換したい$x, y$座標$[m]$ # - $\phi_0, \lambda_0$：平面直角座標系原点の緯度経度$[rad]$（日本には19個あるらしい:[わかりやすい平面直角座標系](http://www.gsi.go.jp/sokuchikijun/jpc.html), [平面直角座標系（平成十四年国土交通省告示第九号）](http://www.gsi.go.jp/LAW/heimencho.html)） # - $a, F$：楕円体の長半径$[m]$及び逆扁平率（値は[日本の測地系](http://www.gsi.go.jp/sokuchikijun/datum-main.html)参照。準拠楕円体はITRF座標系GRS80楕円体で良いらしい） # > 2001年以前の測地基準点成果は、緯度・経度においては日本測地系に基づいた数値で準拠楕円体はベッセル楕円体でしたが、現在の測地測量成果は世界測地系（測地成果2011）と呼び、準拠楕円体はITRF座標系GRS80楕円体です。[日本測地系と世界測地系](http://www.gsi.go.jp/sokuchikijun/datum-main.html#p5) # - $m_0$：平面直角座標系の$x$軸上における縮尺係数(0.9999) # ### 1.3. 計算 # 1.2.に記した$x, y, \phi_0, \lambda_0, a, F, m_0$を既知として、平面直角座標$x, y$を変換した緯度$\phi$と経度$\lambda$を計算する # # #### (1) $n, A_i, \beta_i, \delta_i$ の計算 # 逆扁平率$F$をもとに以下のような$n$を定める。 # $$ # n = \frac{1}{2F-1} # $$ # この$n$から以下の $A_i(i=0, 1, ..., 5), \: \beta _i(i=1, 2, ..., 5), \: \delta_i(i=1, 2, ..., 6)$を定義する。 # # - $A_i(i=0, 1, ..., 5)$ # $$ # \begin{eqnarray} # \begin{cases} # A_0 = 1 + \frac{n^2}{4} + \frac{n^4}{64} & \\ # A_1 = -\frac{3}{2} \bigl( n-\frac{n^3}{8}-\frac{n^5}{64} \bigr) & \\ # A_2 = \frac{15}{16} \bigl( n^2-\frac{n^4}{4} \bigr) & \\ # A_3 = -\frac{35}{48} \bigl( n^3-\frac{5}{16}n^5 \bigr) & \\ # A_4 = \frac{315}{512} n^4 & \\ # A_5 = -\frac{693}{1280}n^5 & # \end{cases} # \end{eqnarray} # $$ # # - $\beta _i(i=1, 2, ..., 5)$ # $$ # \begin{eqnarray} # \begin{cases} # \beta_1 = \frac{1}{2}n - \frac{2}{3}n^2 + \frac{37}{96}n^3 -\frac{1}{360}n^4 - \frac{81}{512}n^5 & \\ # \beta_2 = \frac{1}{48}n^2 + \frac{1}{15}n^3 - \frac{437}{1440}n^4 + \frac{46}{105}n^5 & \\ # \beta_3 = \frac{17}{480}n^3 - \frac{37}{840}n^4 - \frac{209}{4480}n^5 & \\ # \beta_4 = \frac{4397}{161280}n^4 - \frac{11}{504}n^5 & \\ # \beta_5 = \frac{4583}{161280}n^5 & # \end{cases} # \end{eqnarray} # $$ # # - $\delta_i(i=1, 2, ..., 6)$ # $$ # \begin{eqnarray} # \begin{cases} # \delta_1 = 2n - \frac{2}{3}n^2 - 2n^3 + \frac{116}{45}n^4 + \frac{26}{45}n^5 - \frac{2854}{675}n^6 & \\ # \delta_2 = \frac{7}{3}n^2 - \frac{8}{5}n^3 - \frac{227}{45}n^4 + \frac{2704}{315}n^5 + \frac{2323}{945}n^6 & \\ # \delta_3 = \frac{56}{15}n^3 - \frac{136}{35}n^4 - \frac{1262}{105}n^5 + \frac{73814}{2835}n^6 & \\ # \delta_4 = \frac{4279}{630}n^4 - \frac{332}{35}n^5 - \frac{399572}{14175}n^6 & \\ # \delta_5 = \frac{4174}{315}n^5 - \frac{144838}{6237}n^6 & \\ # \delta_6 = \frac{601676}{22275}n^6 & # \end{cases} # \end{eqnarray} # $$ # # #### (2) $\overline{S}_{\phi_0}, \overline{A}$ の計算 # 今までで定義したものから以下のように $\overline{S}_{\phi_0}, \overline{A}$を計算する。 # # ここでは$\phi_0$をラジアンとしているため、[国土地理院の計算式](http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/algorithm/xy2bl/xy2bl.htm)の$\rho^{"}$を除いていることに注意。 # $$ # \overline{S}_{\phi_0} = \frac{m_0 a}{1+n} \biggl( A_0 \phi_0 + \sum_{j=1}^5 A_j sin(2j\phi_0) \biggr) # $$ # # $$ # \overline{A} = \frac{m_0 a}{1+n} A_0 # $$ # # #### (3) $\xi, \eta$ の計算 # (2)で計算した$\overline{S}_{\phi_0}, \overline{A}$と変換したい座標$x, y$から、以下のように$\xi, \eta$を計算する。 # \begin{eqnarray} # \begin{cases} # \xi = \frac{x + \overline{S}_{\phi_0}}{\overline{A}} & \\ # \eta = \frac{y}{\overline{A}} & # \end{cases} # \end{eqnarray} # # #### (4) $\xi^{'}, \eta^{'}$ の計算 # (3)で計算した$\xi, \eta$と(1)で計算した$\beta_i$から、以下のように$\xi^{'}, \eta^{'}$を計算する。 # $$ # \begin{eqnarray} # \begin{cases} # \xi^{'} = \xi - \sum_{j=1}^5 \beta_j sin(2j\xi) cosh(2j\eta) & \\ # \eta^{'} = \eta - \sum_{j=1}^5 \beta_j cos(2j\xi) sinh(2j\eta) & # \end{cases} # \end{eqnarray} # $$ # # #### (5) $\chi$ の計算 # (4)で計算した $\xi^{'}, \eta^{'}$ から、以下のように $\chi\: [rad]$を計算する。 # # $$ # \chi = sin^{-1} \biggl( \frac{sin \: \xi^{'}}{cosh \: \eta{'}} \biggr) \:\: [rad] # $$ # # #### (6) 緯度$\phi$, 経度$\lambda$の計算 # ここまでで求めたものから、以下のように緯度$\phi$, 経度$\lambda$ [rad]を計算する。 # ここでは$\phi, \lambda$ をラジアンとしているため、[国土地理院の計算式](http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/algorithm/xy2bl/xy2bl.htm)の$\rho^{"}$を除いていることに注意。 # # - 緯度 # $$ # \phi = \chi + \sum_{j=1}^6 \delta_j sin(2j\chi) \:\: [rad] # $$ # - 経度 # $$ # \lambda = \lambda_0 + tan^{-1} \biggl( \frac{sinh \: \eta^{'}}{cos \: \xi{'}} \biggr) \:\: [rad] # $$ # ## 2. 実装（平面直角座標 -> 緯度経度の関数） # ### 2.1. import # In[1]: import numpy as np # ### 2.2. 平面直角座標 -> 緯度経度の関数 # In[2]: 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.3. 実数の角度[度]を33°06′14.85664″みたいに度・分・秒表記のstringに変換する関数 # In[3]: 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埋めする # ### 2.4. テスト # [国土地理院の換算サービス](http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/xy2blf.html)と比較 # # -> OK # In[4]: 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)) # ## 3. 実装（緯度経度 -> 平面直角座標の関数） # ### 3.1. 緯度経度 -> 平面直角座標の関数 # In[5]: 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)と比較 # In[6]: x, y = calc_xy(36.103774791666666, 140.08785504166664, 36., 139+50./60) print "x, y = ({0}, {1})".format(x, y) # -> OK # - 日本の範囲（ざっくり）からランダムに緯度経度をとって来て、緯度経度 -> 平面直角座標 -> 緯度経度で元に戻ることを確認 # In[7]: for i in range(10000): # 緯度経度を適当に取ってくる lat = 20 + 25*np.random.rand() lon = 120 + 30*np.random.rand() # 平面直角座標に変換 平面直角座標系はとりあえず2系のみ x, y = calc_xy(lat, lon, 33., 131.) # 再び緯度経度に変換 lat2, lon2 = calc_lat_lon(x, y, 33., 131.) assert ((lat-lat2)**2 + (lon-lon2)**2) < 0.001, "lat, lon = ({0}, {1})".format(lat, lon) print "Test passed!" # ## 4. 鍋割山 # Folium (Leaflet.js)を使い、平面直角座標から鍋割山荘をプロットする。 # # [Folium documentaion](https://folium.readthedocs.io/en/latest/index.html) # ### 4.1. import # In[8]: import folium # ### 4.2. 鍋割山荘の位置 # [国土地理院の換算サービス](http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/xy2blf.html)で「地図上で選択」を行うと、鍋割山荘の平面直角座標は # $$ # (x, y) = (-61474,-62795) \:\: [m] # $$ # また、平面直角座標系の原点は9系（36度0分0秒, 139度50分0秒）らしい。 # # まずはこれを緯度経度に変換する。 # In[9]: nabewari_lat, nabewari_lon = calc_lat_lon(-61474, -62795, 36., 139+50./60) print " latitude: {0}".format(show_angle(nabewari_lat)) print "longitude: {0}".format(show_angle(nabewari_lon)) # ### 4.3. Foliumによるプロット # 計算された鍋割山荘の緯度経度を中心とし、国土地理院の[地理院タイル](https://maps.gsi.go.jp/development/ichiran.html)を利用して地図を描画する。 # # 鍋割山荘の位置にはマーカーを配置する。 # In[10]: # map nabewari_map = folium.Map(location=[nabewari_lat, nabewari_lon], zoom_start=15, tiles='http://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png', attr='© 地理院タイル contributors') # marker folium.Marker([nabewari_lat, nabewari_lon], popup=u'鍋割山荘').add_to(nabewari_map) # plot inline with jupyter nabewari_map # #### [色別標高図](https://maps.gsi.go.jp/development/ichiran.html#relief)を使ってみる # In[11]: # map nabewari_map = folium.Map(location=[nabewari_lat, nabewari_lon], zoom_start=15, tiles='https://cyberjapandata.gsi.go.jp/xyz/relief/{z}/{x}/{y}.png', attr='© 地理院タイル contributors') # marker folium.Marker([nabewari_lat, nabewari_lon], popup=u'鍋割山荘').add_to(nabewari_map) # plot inline with jupyter nabewari_map