最近新幹線に乗って地方へ出かける機会が多く、 GPS トラッキングとかして遊んでるのですが、 GPS ユニットの出力ログが日本測地系なので Google Maps 上で表示させるには測地系を変換する必要が出てきました。
なので、測地系パラメタを使う方法 を参考に、三次元直交座標を経由して変換する方法で実装してみました。ページにあったのは Perl だったので、今回は Python で書き直してみた。ていうか Python 3.0 系って 2.x 系から文法まで変わってるのな。びっくらこいた。
# coding: utf-8 import math # Radian RAD = math.pi / 180 # 日本測地系における赤道半径、扁平率、第一離心率 TKY = 6377397.155 tireT = 1/299.152813 ecT = tireT * (2 - tireT) # 世界測地系における赤道半径、扁平率、第一離心率 WGS = 6378137 tireW = 1/298.257223 ecW = tireW * (2 - tireW) def E2R(B, L, H, Eq, Ec): # 楕円体(B, L, H)から直行座標系(X, Y, Z)へ # Ellipoid to Rectangular Coordinate System Brad = B * RAD Lrad = L * RAD rN = Eq / math.sqrt(1 - Ec * math.pow(math.sin(Brad), 2)) x = (rN + H) * math.cos(Brad) * math.cos(Lrad) y = (rN + H) * math.cos(Brad) * math.sin(Lrad) z = (rN * (1 - Ec) + H) * math.sin(Brad) return (x, y, z) def R2E(X, Y, Z, Eq, Ec): # 直行座標系(X, Y, Z)から楕円体(B, L, H)へ # Rectangular Coordinate System to Ellipoid BdA = math.sqrt(1 - Ec) P = math.sqrt(X * X + Y * Y) T = math.atan2(Z, P * BdA) B = math.atan2(Z + Ec * Eq / BdA * math.pow(math.sin(T), 3), P - Ec * Eq * math.pow(math.cos(T), 3)) L = math.atan2(Y, X) H = P / math.cos(B) - Eq / math.sqrt(1 - Ec * math.pow(math.sin(B), 2)) return (B / RAD, L / RAD, H) def TKY2WGS(lat, lng, hgt): # 日本測地系から世界測地系へ # 引数は 緯度, 経度, 高度 dx = -148 dy = 507 dz = 681 RCS = E2R(lat, lng, hgt, TKY, ecT) Ellips = R2E(RCS[0] + dx, RCS[1] + dy, RCS[2] + dz, WGS, ecW) return(Ellips) def WGS2TKY(lat, lng, hgt): # 世界測地系から日本測地系へ # 引数は 緯度, 経度, 高度 dx = +128 dy = -481 dz = -664 RCS = E2R(lat, lng, hgt, WGS, ecW) Ellips = R2E(RCS[0] + dx, RCS[1] + dy, RCS[2] + dz, TKY, ecT) return(Ellips) # 緯度、経度は度数を浮動小数点形式で指定。 # ラジアンや度/分/秒(Degree/Min/Sec: DMS)形式の場合は事前に変換のこと。 # 緯度 latitude = 35 + 20/60 + 39.984328/3600 # 経度 longitude = 138 + 35/60 + 8.086122/3600 # 高度 height = 0 print(TKY2WGS(latitude, longitude, height)) print(WGS2TKY(latitude, longitude, height))
Python は初めてなのでこういう書き方でいいんだろうか。いまいちよく解らない。クラス分けとかするほどでもないしなぁと思いつつ、関数ベタ書きで。
一応ついでなので世界測地系から日本測地系への変換関数も実装してみた。日本→世界→日本、とやっても、まぁ当たり前っちゃ当たり前だけどやっぱり少し誤差が出る。
Python の利点は｛カッコ｝がないから処理に集中できること、とか言われるけど、逆に制御構造を一つ外すと制御構造中にある文のインデントを全て戻さないと即エラーで止まってしまうので、デバッグ作業が非常にやり辛い。いかにも取って付けた感のある宣伝文句だなぁと。 C 系の書き方でも｛カッコ｝の処理はプログラマ側の問題じゃねーのかと思うんだが…。
あと書き慣れてないせいなのかもしれないけど、なんか頭にヘッドギア付けられてるような感じ。さーっと書くんじゃなくて、こういうしきたり通りに書きなさいと半ば強制されているような雰囲気がある。まぁ教育用途にはもってこいなんだろうけど…慣れないとしんどそう。
これを使って処理してできた KML は、 http://beta.feedvalidator.org/ で Validate できる。あと Google Maps って拡張子見てないので、中身が KML であれば、拡張子が XML でも問題ない。これだと Google Sites に投げたファイルを見ることができる。
