🗺️
平面直角座標から緯度経度への変換
平面直角座標から緯度経度へ変換するための同じ数学的アルゴリズム（ガウス・クリューゲル逆変換と測地系変換）を実装しています。
import numpy as np
def calc_lat_lon(x, y, kei):
"""
平面直角座標系から緯度経度（世界測地系）に変換
Parameters
----------
x : float
平面直角座標のX座標 [m] (北方向)
y : float
平面直角座標のY座標 [m] (東方向)
kei : int
平面直角座標系の系番号 (1-8, 15-17)
例: 1=長崎・鹿児島, 2=福岡・佐賀等, 6=近畿地方, 8=関東・中部地方
Returns
-------
latitude : float
緯度 [度] (世界測地系)
longitude : float
経度 [度] (世界測地系)
Notes
-----
- 日本測地系（ベッセル楕円体）で計算後、世界測地系に変換
- kei=0の場合は (0.0, 0.0) を返す
"""
# =====================================================
# 1. 平面直角座標系の原点定義（日本測地系）
# =====================================================
# 各系番号に対応する原点の経度と緯度を定義
coordinate_origins = {
0: [0, 0], # ダミー（無効な系番号）
1: [129.5, 33.0], # I系（長崎・鹿児島）
2: [131.0, 33.0], # II系（福岡・佐賀・熊本・大分・宮崎）
3: [132 + 1.0/6.0, 36.0], # III系（山口・島根・広島）
4: [133.5, 33.0], # IV系（香川・愛媛・徳島・高知）
5: [134 + 2.0/6.0, 36.0], # V系（兵庫・鳥取・岡山）
6: [136.0, 36.0], # VI系（京都・大阪・福井・滋賀・三重・奈良・和歌山）
7: [137 + 1.0/6.0, 36.0], # VII系（石川・富山・岐阜・愛知）
8: [138.5, 36.0], # VIII系（新潟・長野・山梨・静岡）
15: [127.5, 26.0], # XV系（沖縄本島）
16: [124.0, 26.0], # XVI系（先島諸島）
17: [131.0, 26.0] # XVII系（大東諸島）
}
# 指定された系番号の原点を取得
lambda0_deg, phi0_deg = coordinate_origins[int(kei)]
# kei=0の場合は計算せずに返す
if kei == 0:
return 0.0, 0.0
# 原点座標をラジアンに変換
phi0_rad = np.deg2rad(phi0_deg)
lambda0_rad = np.deg2rad(lambda0_deg)
# =====================================================
# 2. 楕円体パラメータの定義（日本測地系：ベッセル楕円体）
# =====================================================
m0 = 0.9999 # 縮尺係数
a = 6377397.155 # 長半径 [m]
F = 299.152813 # 扁平率の逆数
n = 1.0 / (2*F - 1) # 第三扁平率
# =====================================================
# 3. 補助関数の係数配列を計算
# =====================================================
# Gauss-Kruger投影の展開式に使用する係数
# A係数（子午線弧長の計算用）
A0 = 1 + (n**2)/4.0 + (n**4)/64.0
A1 = -(3.0/2.0) * (n - (n**3)/8.0 - (n**5)/64.0)
A2 = (15.0/16.0) * (n**2 - (n**4)/4.0)
A3 = -(35.0/48.0) * (n**3 - (5.0/16.0)*(n**5))
A4 = (315.0/512.0) * (n**4)
A5 = -(693.0/1280.0) * (n**5)
A_array = np.array([A0, A1, A2, A3, A4, A5])
# β係数（正規緯度から等角緯度への変換用）
b1 = (1.0/2.0)*n - (2.0/3.0)*(n**2) + (37.0/96.0)*(n**3) - (1.0/360.0)*(n**4) - (81.0/512.0)*(n**5)
b2 = (1.0/48.0)*(n**2) + (1.0/15.0)*(n**3) - (437.0/1440.0)*(n**4) + (46.0/105.0)*(n**5)
b3 = (17.0/480.0)*(n**3) - (37.0/840.0)*(n**4) - (209.0/4480.0)*(n**5)
b4 = (4397.0/161280.0)*(n**4) - (11.0/504.0)*(n**5)
b5 = (4583.0/161280.0)*(n**5)
beta_array = np.array([np.nan, b1, b2, b3, b4, b5]) # b0はダミー
# δ係数（等角緯度から測地緯度への変換用）
d1 = 2.0*n - (2.0/3.0)*(n**2) - 2.0*(n**3) + (116.0/45.0)*(n**4) + (26.0/45.0)*(n**5) - (2854.0/675.0)*(n**6)
d2 = (7.0/3.0)*(n**2) - (8.0/5.0)*(n**3) - (227.0/45.0)*(n**4) + (2704.0/315.0)*(n**5) + (2323.0/945.0)*(n**6)
d3 = (56.0/15.0)*(n**3) - (136.0/35.0)*(n**4) - (1262.0/105.0)*(n**5) + (73814.0/2835.0)*(n**6)
d4 = (4279.0/630.0)*(n**4) - (332.0/35.0)*(n**5) - (399572.0/14175.0)*(n**6)
d5 = (4174.0/315.0)*(n**5) - (144838.0/6237.0)*(n**6)
d6 = (601676.0/22275.0)*(n**6)
delta_array = np.array([np.nan, d1, d2, d3, d4, d5, d6]) # d0はダミー
# =====================================================
# 4. 子午線弧長の計算
# =====================================================
# A_: 子午線弧長の係数
A_ = (m0 * a / (1.0 + n)) * A0
# S_: 原点における子午線弧長
S_ = (m0 * a / (1.0 + n)) * (
A0 * phi0_rad +
np.dot(A_array[1:], np.sin(2 * phi0_rad * np.arange(1, 6)))
)
# =====================================================
# 5. 正規化座標（ξ, η）の計算
# =====================================================
try:
xi = (x + S_) / A_ # 正規化されたX座標
eta = y / A_ # 正規化されたY座標
except Exception:
# エラー時は原点を返す
xi = 0
eta = 0
# =====================================================
# 6. 等角緯度座標（ξ', η'）の計算
# =====================================================
# β係数を使った逆変換
xi_prime = xi - np.sum(
beta_array[1:] *
np.sin(2 * xi * np.arange(1, 6)) *
np.cosh(2 * eta * np.arange(1, 6))
)
eta_prime = eta - np.sum(
beta_array[1:] *
np.cos(2 * xi * np.arange(1, 6)) *
np.sinh(2 * eta * np.arange(1, 6))
)
# =====================================================
# 7. 等角緯度χの計算
# =====================================================
chi = np.arcsin(np.sin(xi_prime) / np.cosh(eta_prime)) # [rad]
# =====================================================
# 8. 測地緯度・経度の計算（日本測地系）
# =====================================================
# δ係数を使って等角緯度から測地緯度に変換
latitude_jgs = chi + np.dot(delta_array[1:], np.sin(2 * chi * np.arange(1, 7))) # [rad]
# 経度の計算
longitude_jgs = lambda0_rad + np.arctan(np.sinh(eta_prime) / np.cos(xi_prime)) # [rad]
# ラジアンを度に変換
lat_jgs_deg = np.rad2deg(latitude_jgs)
lon_jgs_deg = np.rad2deg(longitude_jgs)
# =====================================================
# 9. 日本測地系から世界測地系への変換
# =====================================================
# TKY2JGD（測地系変換パラメータ）を使用した変換式
lat_wgs_deg = (
lat_jgs_deg
- 0.00010695 * lat_jgs_deg
+ 0.000017464 * lon_jgs_deg
+ 0.0046017
)
lon_wgs_deg = (
lon_jgs_deg
- 0.000046038 * lat_jgs_deg
- 0.000083043 * lon_jgs_deg
+ 0.010040
)
return lat_wgs_deg, lon_wgs_deg
また、以下を使って検証できます。
# =====================================================
# テスト・検証コード
# =====================================================
if __name__ == "__main__":
print("=" * 60)
print("平面直角座標→緯度経度変換のテスト")
print("=" * 60)
# テストケース1: VI系（近畿地方）の原点付近
print("\n【テスト1】VI系（近畿地方）の原点付近")
lat1, lon1 = calc_lat_lon(0, 0, 6)
print(f"入力: x=0m, y=0m, 系=6")
print(f"出力: 緯度={lat1:.8f}°, 経度={lon1:.8f}°")
print(f"期待値: 約36°N, 136°E付近")
# テストケース2: VIII系（関東地方）の東京付近の座標例
print("\n【テスト2】VIII系（関東地方）")
lat2, lon2 = calc_lat_lon(-50000, -30000, 8)
print(f"入力: x=-50000m, y=-30000m, 系=8")
print(f"出力: 緯度={lat2:.8f}°, 経度={lon2:.8f}°")
# テストケース3: 無効な系番号
print("\n【テスト3】無効な系番号（kei=0）")
lat3, lon3 = calc_lat_lon(1000, 2000, 0)
print(f"入力: x=1000m, y=2000m, 系=0")
print(f"出力: 緯度={lat3:.8f}°, 経度={lon3:.8f}°")
print(f"期待値: (0.0, 0.0)")
# テストケース4: XV系（沖縄）
print("\n【テスト4】XV系（沖縄本島）")
lat4, lon4 = calc_lat_lon(10000, 5000, 15)
print(f"入力: x=10000m, y=5000m, 系=15")
print(f"出力: 緯度={lat4:.8f}°, 経度={lon4:.8f}°")
print(f"期待値: 約26°N, 127.5°E付近")
print("\n" + "=" * 60)
print("検証完了")
print("=" * 60)
系が与えられていない場合は、以下のコードを加工して対応できます。
kei1 = ['長崎','鹿児島']
kei2 = ['福岡','佐賀','熊本','大分','宮崎'] # ,'鹿児島'
kei3 = ['山口','島根','広島']
kei4 = ['香川','愛媛','徳島','高知']
kei5 = ['兵庫','鳥取','岡山']
kei6 = ['京都','大阪','福井','滋賀','三重','奈良','和歌山']
kei7 = ['石川','富山','岐阜','愛知']
kei8 = ['新潟','長野','山梨','静岡']
kei15 = ['沖縄']
if ken in kei1:
lambda0_deg= 129.5; phi0_deg = 33.0
elif ken in kei2:
lambda0_deg= 131.0 ; phi0_deg = 33.0
elif ken in kei3:
lambda0_deg= 132 + 1.0 / 6.0 ; phi0_deg = 36.0
elif ken in kei4:
lambda0_deg= 133.5 ; phi0_deg = 33.0
elif ken in kei5:
lambda0_deg= 134 + 2.0 / 6.0 ; phi0_deg = 36.0
elif ken in kei6:
lambda0_deg= 136.0; phi0_deg = 36.0
elif ken in kei7:
lambda0_deg= 137 + 1.0 / 6.0 ; phi0_deg = 36.0
elif ken in kei8:
lambda0_deg= 138.5 ; phi0_deg = 36.0
elif ken in kei15:
lambda0_deg= 127.5 ; phi0_deg = 26.0
Discussion