測地座標をメルカトル図法に変換
目次
- あらすじ
- メルカトル図法
- 換算式
- ミラー図法
- 換算式
あらすじ
どうせ簡単だろうと高を括って測地座標を直交座標に変換する関数の作成に挑戦。 思いの外難しい用語と複数の規格の存在で苦戦するもなんとか完成させた。 ようやく完成させたはいいが使い所がないことに気づき、結局メルカトル図法への変換を実装することに！！
事前に十分計画してから行動しましょう。
メルカトル図法
正角円筒図法ともいう。等角航路（地球上の2点を結ぶ最短航路）が直線で表される。 赤道付近の陸地は正確な大きさだが、高緯度ほど面積が大きく描かれる。 一方で角度は正しいので、十分狭い範囲だけを見ると形が正しく表示される。そのため昔のGoogle Mapでも採用されていた時期もある。
詳しくはWikipediaへ
https://ja.wikipedia.org/wiki/メルカトル図法
換算式
Wikipediaに計算式が載っています。
https://ja.wikipedia.org/wiki/メルカトル図法#投影法の表式
ありがたや。
経度：λ
緯度：φ
地図の中央経度：λ0
地図上の点：x, y
経緯度→メルカトル図法
x = λ-λ0
y = gd^-1(φ)
メルカトル図法→経緯度
λ = x + λ0
φ = gd(y)
まって。グーデルマン関数って何。思ってたより難しそうなんですけど…。 お、Wikipediaに計算式が載っています。
gd(x) = 2 arctan(e^x) - 1/2 * π
HSP3で実装するとこんな感じ。
#defcfunc gd_func double p_phi
; 2.7182818284590452354 ; ネイピア数 e
return 2. * atan( powf(2.7182818284590452354, p_phi) ) - M_PI / 2.
逆グーデルマン関数も載ってますね。
gd^-1(x)
= arsinh( tan(x) )
= ln( tan(1/4π+1/2π) )
#defcfunc gd_i_func double p1
;logf = ln = 自然対数 = 底がeのlog
return logf(absf( tan( M_PI / 4. + p1 / 2. ) ))
難しそうな式ですが、そのまま書き写せばいいだけなので実装は簡単です。グーデルマン関数おそるるに足りず。（意味はさっぱりわかりませんが。）
あとは緯度をグーデルマン関数に突っ込むだけでメルカトル図法に変換できるので実装例は書かなくてもいいですよね。
注意点として、HSP3で角度の計算をする場合、ラジアンじゃないといけないのでdeg2rad関数とrad2deg関数使うくらいですか。
地球1週＝2π の大きさの平面地図座標が相互変換できるようになります。
λ0が0なら東経0度（グリニッジ子午線）が地図の左右中心になります。
ミラー図法
メルカトル図法は有名ですが、北極点や南極点が描けないことや、北極・南極近くが大きく表示されすぎます。 ミラー図法というのがちょうど良さそうなので調べてみました。
- ミラー図法
- 円筒図法の一種。主に世界地図に用いられる。メルカトル図法の、南北両極が無限遠点になってしまうという問題を改善した図法。
角度も直線航路も面積も正しさを持たない。
正しい情報は、Wikipedia参照。
なるほど、大まかに世界を見渡せる地図としては便利そうです。ゲーム向きかもしれません。
換算式
Wikipediaに計算式が載っていました。 緯度を4/5倍してグーデルマン関数に入れるだけですね。
経緯度→ミラー図法
x = λ-λ0
y = 5/4 * gd^-1(4/5*φ)
逆変換は載っていませんが、簡単な計算式なのですぐに導出できます。
ミラー図法→経緯度
λ = x + λ0
φ = 5/4 * gd(4/5*y)
メルカトル図法とほとんど同じ式ですね。
いろいろな地図
Wikipediaにはいろいろな地図（地図投影法）が紹介されています。 変換式も載っているようです。専門職の人には馴染み深い地図ばかりかもしれませんが、私にはさっぱりです。
モジュール
ということで、前回の直交座標系への相互変換モジュールにメルカトル図法、ミラー図法への相互変換を実装したモジュールです。 HSP3で使用することができます。