GIS Coordinate transformation

GIS 座標変換の計算ではまる

Unfortunately this text is in Japanese only.
座標からタイルのx,y座標を求めるがハマる。
OpenStreetMapの [Wiki : Slippy map tilenames] を参照して真球メルカトル投影の座標変換を行ったつもりが。。。

OpenStreetMap Wikiから転写。

単純な三角関数と対数なのに．．．






各言語のソースコードも載っている。なんて親切なページ(^_^;;
VB.NETのソースコード

Private Function CalcTileXY(ByVal lat As Single, ByVal lon As Single, ByVal zoom As Long) As Point CalcTileXY.X = CLng(Math.Floor((lon + 180) / 360 * 2 ^ zoom)) CalcTileXY.Y = CLng(Math.Floor((1 - Math.Log(Math.Tan(lat * Math.PI / 180) + 1 / Math.Cos(lat * Math.PI / 180)) / Math.PI) / 2 * 2 ^ zoom)) End Function

以下は表計算ソフトで試算した結果。

サンプルにGoogle Mapで調べた福智山山頂の座標を使用。座標系は多分WGS84。

ズームレベル５，６，７で実際に表示させてみるが正解のタイルＹ座標より南に計算されてしまう。














因みにズームレベル５の表計算の計算式は
=(2^ (B9-1) * (1 - (LOG(TAN(lat_rad) + SEC(lat_rad)) / PI() ) ))
何故だ！何故だ！何故だぁ！寝られん！

[memo]
真球メルカトル投影：EPSG:3857(旧：900913)　(電子国土WebシステムVer.4も準拠)
等経緯度投影：EPSG:4326

2014/01/05 追記
原因判明。
表計算ソフトのLOG()は基数省略すると底は１０だと．．．
e を底とする自然対数はLN()
うまく計算できました。