国土地理院の計算式はラジアンと度がごっちゃになっていて曖昧なので微調整した。 以下に出てくる三角関数はすべてラジアンを受け取り、逆三角関数はすべてラジアンを返すものとする。
2001年以前の測地基準点成果は、緯度・経度においては日本測地系に基づいた数値で準拠楕円体はベッセル楕円体でしたが、現在の測地測量成果は世界測地系(測地成果2011)と呼び、準拠楕円体はITRF座標系GRS80楕円体です。日本測地系と世界測地系
1.2.に記した$x, y, \phi_0, \lambda_0, a, F, m_0$を既知として、平面直角座標$x, y$を変換した緯度$\phi$と経度$\lambda$を計算する
逆扁平率$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)$を定義する。
今までで定義したものから以下のように $\overline{S}_{\phi_0}, \overline{A}$を計算する。
ここでは$\phi_0$をラジアンとしているため、国土地理院の計算式の$\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 $$(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}
(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} $$
(4)で計算した $\xi^{'}, \eta^{'}$ から、以下のように $\chi\: [rad]$を計算する。
$$ \chi = sin^{-1} \biggl( \frac{sin \: \xi^{'}}{cosh \: \eta{'}} \biggr) \:\: [rad] $$ここまでで求めたものから、以下のように緯度$\phi$, 経度$\lambda$ [rad]を計算する。 ここでは$\phi, \lambda$ をラジアンとしているため、国土地理院の計算式の$\rho^{"}$を除いていることに注意。
import numpy as np
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]
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埋めする
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"
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]
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)
-> OK
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!"
Test passed!
import folium
国土地理院の換算サービスで「地図上で選択」を行うと、鍋割山荘の平面直角座標は $$ (x, y) = (-61474,-62795) \:\: [m] $$ また、平面直角座標系の原点は9系(36度0分0秒, 139度50分0秒)らしい。
まずはこれを緯度経度に変換する。
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))
latitude: 35°26'38.087211585" longitude: 139°08'29.8934370248"
# 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='© <a href="https://maps.gsi.go.jp/development/ichiran.html">地理院タイル</a> contributors')
# marker
folium.Marker([nabewari_lat, nabewari_lon], popup=u'鍋割山荘').add_to(nabewari_map)
# plot inline with jupyter
nabewari_map
# 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='© <a href="https://maps.gsi.go.jp/development/ichiran.html">地理院タイル</a> contributors')
# marker
folium.Marker([nabewari_lat, nabewari_lon], popup=u'鍋割山荘').add_to(nabewari_map)
# plot inline with jupyter
nabewari_map