測地座標を直交座標

測地座標を直交座標に変換

2020/05/18

はじめに

　HSPTV!掲示板に「日本地図を表示したい」という質問が書き込まれていました。「拡大してもきれいな表示がいい」とのことなので、ベクタ形式のデータを表示したいんだろうなという感じの質問です。質問そのものは、緯度経度からメルカトル図法の座標を計算して出すだけで解決していました。

　そこでふと緯度経度を手軽に扱える形式に変換できる関数があると3D表示するには都合良さそう！極座標から直交座標に変換するだけだから簡単だろ。ちょっと真面目に高精度に実装してみよう！と考えてしまったのでHSP3で実装してみました。

地図データを入手

　この質問を見ていて、まず地図のデータはどこから持ってくるんだろうと思って少し探してみました。

国土交通省国土数値情報　ダウンロードサービス
http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v2_3.html

　詳細なデータなのですごくサイズが大きいですが、入手できるんですね。早速ダウンロードしてみました。HSP3からゲームなどで使うには大きすぎるので、データを間引く必要はありそうですが使えそうです。座標の数値データは…なるほど緯度経度です。地図データだからあたりまえですね。 3D座標に変換するには計算が必要そうです。

計算式を入手

　メルカトル図法は、緯度経度の数値をそのままY,X座標に置き換えればそれっぽく書けそうな気がするので今回は直交座標系（X,Y,Zの3次元座標）への変換を考えてみました。（※そんなことはなかった。）

　まずは実装する換算式と実装後の答え合わせが必要なので探してみると良いサイトがありました。

緯度・経度と地心直交座標の相互換算
https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/transf.html

　さすが国土地理院さんいい仕事してます。極座標から直交座標に変換するだけだから簡単だろと思ったのに意外とめんどくさそうです。いやになりますね。しかし足りないパラメータもあるようなのでもう少し調べてみました。というか楕円体の計算式とか覚えてないのでその辺りの式ですね。

師匠の散歩
http://tancro.e-central.tv/grandmaster/excel/radius.html

　とても見やすい式で綺麗にまとめてあります。ありがたや。
しかし…楕円体高？うさぎとり…半径？ジオゾイド高？離心率？聞いたこと無い言葉ばかり。ここからか。ここからなのかッ！

用語解説

　調べているとわからない言葉や聞いたことある言葉だけど、知ってた意味と違う言葉がたくさん出てきたのでまとめてみました。

測地座標系: 地球上の位置を緯度、経度、楕円体高で表す座標系。
地心直交座標系: 地球の重心を中心としたXYZの直交座標系。右手系。
赤道上の経度0度（グリニッジ子午線）方向がX軸方向。
X軸を東に90度回した方向がY方向。
中心から北極に向かう方向がZ方向。
ECEF 座標 (Earth-Centered Earth-Fixed)と呼ばれるものと同じようです。（但し、同じと書いてある資料はなぜか見つかりませんでした。厳密には違うのかもしれません。）
楕円体高: 地球を仮想的に表した楕円体（地球楕円体）の表面からの高さ。
楕円体高＝ジオイド高＋標高
ジオイド高: 地球表面の重力分布は実は均一ではありません。海洋より内陸の方が質量が大きいからそうなる理屈ですね。したがって平均海面高さ（標高0m）も地域によって異なります。地球楕円体から標高0m（ジオイド）までの高さを「ジオイド高」と呼びます。
日本のジオイド高さはおおよそ35m前後ぐらいです。日本は東京湾の海面水位が基準らしいのですが、数値わかりませんでした。
標高: 標高0m（ジオイド）からの高さ。
標高＝楕円体高－ジオイド高
卯酉線曲率半径: 卯酉線（ぼうゆうせん）は地球の赤道のことのようです。楕円体の長半径。算出するために、測地系、長半径、扁平率が必要そう。
離心率: よくわからないので知りたい人は Wikipedia先生と数学系サイト、数学の教科書見てください。しかし計算式はあるので万事解決です。算出するために、測地系、長半径、扁平率が必要そう。
測地系: 地球の楕円の形と座標系の規格。GRS80やらWGS84などいろいろあり、各国で採用している測地系は異なる事がある。測地系が決まると長半径、扁平率の値が決まる。
日本はGRS80を採用している。…こいう規格は世界統一してほしいですよね。勘弁してほしいです。
長半径: 楕円の長軸の半分。地球は球ではなく楕円なので必要な値。測地系が決まれば決まる値。
扁平率: 球に比べてどれくらい潰れているかを表す値。地球は球ではなく楕円なので必要な値。測地系が決まれば決まる値。

参考資料
ジオイドを知る - 国土地理院
https://www.gsi.go.jp/buturisokuchi/grageo_geoid.html
楕円体高、ジオイド高、標高についてはこのサイトがわかりやすいと思います。

ジオイド高計算
https://vldb.gsi.go.jp/sokuchi/surveycalc/geoid/calcgh/calcframe.html
日本限定で任意の位置のジオイド高を調べることができます。

測地座標系→地心直交座標系

　測地座標（緯度経度と楕円体高）から地心直交座標を算出してみます。

　まずは楕円関係の計算式から離心率と卯酉線曲率半径を計算します。測地系は GRS80 を使用し、緯度を変数 ido として与えます。
師匠の散歩
http://tancro.e-central.tv/grandmaster/excel/radius.html

ido : 緯度
kdo : 経度
he : 楕円体高

; 長半径 [m]
a = 6378137.0

; 扁平率
f = 1.0 / 298.257222101

; 短半径	Semi-minor axis
b = a * ( 1. - f )

; 離心率	Eccentricity
e    = sqrt( 2. * f - f * f )
e2   = e * e
sini = sin(ido)
W    = sqrt( 1. - e2 * sini*sini)

; 卯酉線曲率半径	Prime vertical radius of curvature
N = a / W

　換算式を使って地心直交座標を計算します。
緯度・経度と地心直交座標の相互換算
https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/transf.html
の「計算式」を参照。

lat : 緯度
lon : 経度
he : 楕円体高
e2 : 離心率の2乗
N : 卯酉線曲率半径

; 地心直交座標系
X = (N + he) * cos(lat) * cos(lon)
Y = (N + he) * cos(lat) * sin(lon)
Z = (N * (1. - e2) + he) * sin(lat)

入力例
緯度： 36.103774792 度
経度：140.087855042 度
楕円体高：65.840 m

計算結果
X座標：-3957314.622 m
Y座標： 3310254.134 m
Z座標： 3737540.044 m

大丈夫そうです。

地心直交座標系→測地座標系

　せっかくなので逆もやってみます。地心直交座標から測地座標（緯度経度と楕円体高）を算出してみます。計算式を見ると反復計算を何回か繰り返して結果を出すようです。直交座標系への変換より少し複雑ですね。

　まずは離心率を出します。

; 長半径 [m]
a = 6378137.0

; 扁平率
f = 1.0 / 298.257222101

;	短半径	Semi-minor axis
b = a * ( 1. - f )

;	離心率	Eccentricity
e    = sqrt( 2. * f - f * f )
e2   = e * e

経度→緯度→楕円体高の順に計算したほうが良さそうです。

in_X, in_Y, in_Z : 地心直交座標( x,y,z )
lat : 緯度
lon : 経度
he : 楕円体高
e2 : 離心率の2乗
N : 卯酉線曲率半径

P   = sqrt( in_X * in_X + in_Y * in_Y)

;	経度
lon = atan( in_Y, in_X)

;	緯度
lat0 = atan( in_Z, P * (1. - e2) )
repeat
	;	卯酉線曲率半径
	sini = sin(lat0)
	N    = a / sqrt( 1. - e2 * sini*sini)
	;	緯度
	lat = atan( in_Z, P - e2 * N * cos(lat0) )
	if absf(lat - lat0) <= 0.000000000001 {
		; 10e-12 以下
		break
	}
	lat0 = lat
loop

;	楕円体高
sini = sin(lat)
N    = a / sqrt( 1. - e2 * sini*sini)
he   = P / cos(lat) - N

　HSP3は整数と実数で関数を変えなければいけなかったり、計算の順番や変数代入時も注意しないと実数計算のつもりが整数になってしまっていることがあるので注意が必要です。 absf関数をabsと書いてしまったり、「1.」と書いて実数にすべきところを「1」と書いてしまって計算が合わなくて何度か悩みました。

入力例
in_X = -3957314.62177
in_Y = 3310254.13387
in_Z = 3737540.04441

計算結果
緯度： 36.103774792 度
経度：140.087855042 度
楕円体高：65.840 m

これも大丈夫そうです。

モジュール

　今回の計算はHSP3で使いやすいようにモジュール化しました。興味があれば使ってください。

mod経緯度

　これで緯度経度の数値データが有れば、プログラムで3D空間に正確な描画をすることができますね。 2点間の直線距離とか方角も正確に算出できます！…何に使うか知りませんが。何に使うんだ。日本からブラジルまでの直線距離ってもマントル突き抜ける直線ルートだぞ。何に使うんだ。もしかして普通に極座標として計算してもあまり変わらない計算結果だったのでは…ここまでなかなかに大変だったんだが…私は一体…。あれ…もしかしてこれメルカトル図法とかそいういう関数も実装しないと使いどころのないモジュールなのでは？

次回…メルカトル図法！