python で 緯度経度と世界座標を相互変換してみる

やりたかったこと

• 緯度経度の情報を、メッシュ的に区画分割した単位で集計等を行いたい
• そのため、緯度経度を解像度を落した区画xyに配置したかった
• Google Maps APIの fromLatLngToPoint 同等のことをPythonでオフラインでやりたかった

やったこと

TRAIL NOTEさんの http://www.trail-note.net/tech/coordinate/ にあった情報を元に、pythonで、緯度経度と世界座標の変換を行う関数にしてみました。

実装

ひたすら↑の記載内容を写経しました。

from math import pi
from math import tanh
from math import sin
from math import asin
from numpy import arctanh

# refer from http://www.trail-note.net/tech/coordinate/

def tile2latlon(x, y, z):
    L = 85.05112878
    lon = ((x / 2.0**(z+7) )-1) * 180
    lat = 180/pi * (asin(tanh(-pi/2**(z+7)*y + arctanh(sin(pi/180*L)))))
    return [lat, lon]

def latlon2tile2(lat, lon, z):
    L = 85.05112878
    x = int((lon/180 + 1) * 2**(z+7))
    y = int( (2**(z+7) / pi * ( -arctanh(sin(pi*lat/180)) + arctanh(sin(pi*L/180)) ) ))
    return [x,y]

確認

既知の緯度経度を入力して、国土地理院の地図 で確認してみます。

#simbashi
lat = 35.666280
lon = 139.758375
# lat = 30.335927
# lon = 130.504283

pix = 256
# lat/lon to tile
for z in range(15,19):
    a,b = latlon2tile2(lat,lon,z)
    print [a,b]
    print 'http://cyberjapandata.gsi.go.jp/xyz/std/{0}/{1:d}/{2:d}.png'.format(z,a/pix,b/pix)
    a,b = tile2latlon(a,b,z)
    print [a,b]

実行するとこんな感じで出力されますので、元のlat/longとあっていること、zoomレベルが上がると精度が上がること。実際の地図での場所の確認を行いました。（国土地理院の地図はズームレベル18まで）

[7450910, 3303678]
http://cyberjapandata.gsi.go.jp/xyz/std/15/29105/12904.png
[35.66629207402928, 139.75836753845215]
[14901820, 6607356]
http://cyberjapandata.gsi.go.jp/xyz/std/16/58210/25809.png
[35.66629207402928, 139.75836753845215]
[29803640, 13214713]
http://cyberjapandata.gsi.go.jp/xyz/std/17/116420/51619.png
[35.66628335763601, 139.75836753845215]
[59607281, 26429426]
http://cyberjapandata.gsi.go.jp/xyz/std/18/232840/103239.png
[35.66628335763601, 139.75837290287018]

要確認事項

• 測地系とかよく分かっていません
• pyproj といったライブラリがあるようなのでこちらを使う方がいいのかもしれません。使い方をご存知の方は教えてください！