タイルマップで緯度経度からタイルを求める
Googleマップで採用されている地図画像の形式であるタイルマップは、地理院地図でも使われるなど地図をネットで公開する際に便利な形式です。
これをgoogle Maps APIなどから使う際には、自動的に表示範囲の地図画像が読み込まれて表示されるので、詳しい仕組みを知らなくてもいいのですが、スタンドアローンのソフトで表示する場合には、自作する必要があります。
そこで、必要そうな式を作ってみました。
まず、メルカトル図法に関して、その投影法の式は次のようになっています。
メルカトル図法の式
緯度lat 経度lon からメルカトル図法の座標に変換(Rは地球の半径)
X=π*R*(lon/180)
Y=-R*Log(Tan((45+lat/2)*π/180))
lat=90になるとタンジェントの値が無限大となるので、メルカトル図法の地図では極点を示せません。
タイルマップの表示範囲と分割
タイルマップは、このメルカトル図法を東西方向と南北方向を正方形で切り出して、内部を分割し（ズームレベル）、小縮尺から大縮尺までの地図を表示します。
地理院地図の説明 http://portal.cyberjapan.jp/help/development.html#siyou-zm
ズームレベルを1つあげると、ファイルの数は4倍になります。
タイルマップの南北端
メルカトル図法では、南北端は無限大なりますが、タイルマップでは南北に限界があります。
メルカトル投影後に南北方向の長さを東西方向に合わせ、正方形に切り出します。
東西方向の幅は、半球なら
X = π*R*(180/180) = π*R
となります。
これが赤道から南北端までの距離に等しくなる緯度を次式から求めていくと、
R*Log(Tan((45+lat/2)*π/180))=R*180*π/180
Log(Tan((45+lat/2)*π/180))=π
Tan((45+lat/2)*π/180)=e^π
Tan(π/4+π/360*lat)=e^π
lat=2 * atn(e^π) * 180 /π - 90
= 85.0511287798067
と求まりました。
緯度経度とズームレベルを指定して、当該タイルのXY値を取得
さて、いよいよ本題です。表示されている地図にタイルマップのタイルを表示するには、緯度経度とズームレベルを指定してタイルを特定する必要があります。
ズームレベル(z)ごとの縦横の分割数は、2^zとなります。
そのため、メルカトル投影後のタイルの幅・高さ(W(z))は、
W(z)=(2*π*R)/2^z
となります。
そこで、経度lonのX値は
X(lon,z)=int((π*R*(lon/180)+π*R))/W(z))
=int(π*R*(lon/180+1))/W(z))
=int(π*R*(lon/180+1))/(2*π*R)/2^z))
=int((lon/180+1)*2^z/2)
緯度latのY値は、
Y(lat,z)=int((-R*Log(Tan((45+lat/2)*π/180))+π*R)/W(z))
=int((R*(-Log(Tan((45+lat/2)*π/180))+π)/W(z))
=int((R*(-Log(Tan((45+lat/2)*π/180))+π)/(2*π*R)/2^z)
=int(((-Log(Tan((45+lat/2)*π/180))+π)*2^z/(2*π))
となります。地球の半径Rは計算では使わなくてもすみます。
指定したタイル(x,y,z)の北端緯度、西端経度を求める
指定した位置のタイルが特定されたら、当該タイルの四隅の緯度経度を求めることで、地図上に表示できます。
西端経度は、
lon(x,z)=((x*W(z)/(2*π*R))*360-180
=((x*(2*π*R)/2^z/(2*π*R))*360-180
=(x/2^z)*360-180
北端緯度は、まず投影後の座標mapyを求め、
mapy=(y / 2 ^ z) * 2 * π- π
メルカトル投影の逆変換を行って求めます。
lat(y,z)=2 * atn(e^(-mapy)) * 180 / π - 90
これをプログラムに組み込めば、自作ソフトにタイルマップを表示できるようになります。
これをgoogle Maps APIなどから使う際には、自動的に表示範囲の地図画像が読み込まれて表示されるので、詳しい仕組みを知らなくてもいいのですが、スタンドアローンのソフトで表示する場合には、自作する必要があります。
そこで、必要そうな式を作ってみました。
まず、メルカトル図法に関して、その投影法の式は次のようになっています。
メルカトル図法の式
緯度lat 経度lon からメルカトル図法の座標に変換(Rは地球の半径)
X=π*R*(lon/180)
Y=-R*Log(Tan((45+lat/2)*π/180))
lat=90になるとタンジェントの値が無限大となるので、メルカトル図法の地図では極点を示せません。
タイルマップの表示範囲と分割
タイルマップは、このメルカトル図法を東西方向と南北方向を正方形で切り出して、内部を分割し（ズームレベル）、小縮尺から大縮尺までの地図を表示します。
地理院地図の説明 http://portal.cyberjapan.jp/help/development.html#siyou-zm
ズームレベルを1つあげると、ファイルの数は4倍になります。
タイルマップの南北端
メルカトル図法では、南北端は無限大なりますが、タイルマップでは南北に限界があります。
メルカトル投影後に南北方向の長さを東西方向に合わせ、正方形に切り出します。
東西方向の幅は、半球なら
X = π*R*(180/180) = π*R
となります。
これが赤道から南北端までの距離に等しくなる緯度を次式から求めていくと、
R*Log(Tan((45+lat/2)*π/180))=R*180*π/180
Log(Tan((45+lat/2)*π/180))=π
Tan((45+lat/2)*π/180)=e^π
Tan(π/4+π/360*lat)=e^π
lat=2 * atn(e^π) * 180 /π - 90
= 85.0511287798067
と求まりました。
緯度経度とズームレベルを指定して、当該タイルのXY値を取得
さて、いよいよ本題です。表示されている地図にタイルマップのタイルを表示するには、緯度経度とズームレベルを指定してタイルを特定する必要があります。
ズームレベル(z)ごとの縦横の分割数は、2^zとなります。
そのため、メルカトル投影後のタイルの幅・高さ(W(z))は、
W(z)=(2*π*R)/2^z
となります。
そこで、経度lonのX値は
X(lon,z)=int((π*R*(lon/180)+π*R))/W(z))
=int(π*R*(lon/180+1))/W(z))
=int(π*R*(lon/180+1))/(2*π*R)/2^z))
=int((lon/180+1)*2^z/2)
緯度latのY値は、
Y(lat,z)=int((-R*Log(Tan((45+lat/2)*π/180))+π*R)/W(z))
=int((R*(-Log(Tan((45+lat/2)*π/180))+π)/W(z))
=int((R*(-Log(Tan((45+lat/2)*π/180))+π)/(2*π*R)/2^z)
=int(((-Log(Tan((45+lat/2)*π/180))+π)*2^z/(2*π))
となります。地球の半径Rは計算では使わなくてもすみます。
指定したタイル(x,y,z)の北端緯度、西端経度を求める
指定した位置のタイルが特定されたら、当該タイルの四隅の緯度経度を求めることで、地図上に表示できます。
西端経度は、
lon(x,z)=((x*W(z)/(2*π*R))*360-180
=((x*(2*π*R)/2^z/(2*π*R))*360-180
=(x/2^z)*360-180
北端緯度は、まず投影後の座標mapyを求め、
mapy=(y / 2 ^ z) * 2 * π- π
メルカトル投影の逆変換を行って求めます。
lat(y,z)=2 * atn(e^(-mapy)) * 180 / π - 90
これをプログラムに組み込めば、自作ソフトにタイルマップを表示できるようになります。
この記事へのコメント