座標の変換（緯度経度から世界座標、ピクセル座標、タイル座標）

以下はTrailNoteさんの

TrailNote : 座標の変換（世界座標、ピクセル座標、タイル座標、緯度・経度）

http://www.trail-note.net/tech/coordinate/

からのひたすら写経ですが、

点群として与えられてる地形データーをクラスター状の格子データに変換して立体地図を作ってみたかったので、その前段として、緯度経度からGoogleマップや国土地理院のタイル地図などで使われているタイル地図内の各ピクセルの緯度経度が知りたく、世界座標、ピクセル座標、タイル座標に変換する関数をpythonで書いてみました。

ソースは以下

三角関数とか逆双曲線関数とか使います。

#!/usr/bin/env python
# -*- coding: utf-8 -*- 

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

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

import pandas as pd
import numpy as np



def fromLatLngToWordLatLon(lat, lon, z):
    L = 85.05112878
    wordLon = (2**(z+7)) * ((lon/180) + 1)
    wordLat = (2**(z+7)/pi) * (-1. * arctanh(sin(lat * pi/180)) + arctanh(sin(L * pi/180)))
    return wordLat, wordLon

def fromLatLngToPoint(lat,lon, z):
    L = 85.05112878
    pointX = int( (2**(z+7)) * ((lon/180) + 1) )
    pointY = int( (2**(z+7)/pi) * (-1 * arctanh(sin(lat * pi/180)) + arctanh(sin(L * pi/180))) )
    return pointX, pointY


def fromPointToLatLng(pixelLat,pixelLon, z):
    L = 85.05112878
    lon = 180 * ((pixelLon / 2.0**(z+7) ) - 1)
    lat = 180/pi * (asin(tanh(-pi/2**(z+7)*pixelLat + arctanh(sin(pi/180*L)))))
    return lat, lon

確認のため、TrailNoteさん例に載ってる利尻山（北海道）、雲取山（東京)、宮之浦岳（鹿児島）の山頂座標の緯度経度から世界座標、ピクセル座標、タイル座標を求めてみてるのに

上のソースを、こんな感じで実行すると

    # 1タイルのピクセル数
    pix = 256

    #利尻山
    lat=45.178506
    lon=141.242035

    for z in range(15,19):
        # 世界座標に変換
        wordLat, wordLon = fromLatLngToWordLatLon(lat, lon, z)
        print [wordLat, wordLon]

        # ピクセル座標に変換
        pointX, pointY = fromLatLngToPoint(lat, lon, z)
        print [pointX, pointY]
        tileURL = 'http://cyberjapandata.gsi.go.jp/xyz/std/{0}/{1:d}/{2:d}.png'.format(z, pointX/pix, pointY/pix)
        print tileURL
        
        # ピクセル座標を緯度経度に変換
        a, b = fromPointToLatLng(pointY, pointX, z)
        print [a, b]


    #雲取山
    lat=35.855499
    lon=138.943905
    
    for z in range(15,19):
        # 世界座標に変換
        wordLat, wordLon = fromLatLngToWordLatLon(lat, lon, z)
        print [wordLat, wordLon]

        # ピクセル座標に変換
        pointX, pointY = fromLatLngToPoint(lat, lon, z)
        print [pointX, pointY]
        tileURL = 'http://cyberjapandata.gsi.go.jp/xyz/std/{0}/{1:d}/{2:d}.png'.format(z, pointX/pix, pointY/pix)
        print tileURL
        
        # ピクセル座標を緯度経度に変換
        a, b = fromPointToLatLng(pointY, pointX, z)
        print [a, b]


    #宮之浦岳
    lat=30.335927
    lon=130.504283
    
    for z in range(15,19):
        # 世界座標に変換
        wordLat, wordLon = fromLatLngToWordLatLon(lat, lon, z)
        print [wordLat, wordLon]

        # ピクセル座標に変換
        pointX, pointY = fromLatLngToPoint(lat, lon, z)
        print [pointX, pointY]
        tileURL = 'http://cyberjapandata.gsi.go.jp/xyz/std/{0}/{1:d}/{2:d}.png'.format(z, pointX/pix, pointY/pix)
        print tileURL
        
        # ピクセル座標を緯度経度に変換
        a, b = fromPointToLatLng(pointY, pointX, z)
        print [a, b]

結果はこうなります。

なお、下の地図画像は国土地理院さんの地図タイルをタイルURLの下にそのまま貼り付けてみました。

[3011700.7215389027, 7485481.957603555]
[3011700, 7485481]
http://cyberjapandata.gsi.go.jp/xyz/std/15/29240/11764.png

[45.178527827298865, 141.24199390411377]
[6023401.443077805, 14970963.91520711]
[6023401, 14970963]
http://cyberjapandata.gsi.go.jp/xyz/std/16/58480/23528.png

[45.17851270178206, 141.2420153617859]
[12046802.88615561, 29941927.83041422]
[12046802, 29941927]
http://cyberjapandata.gsi.go.jp/xyz/std/17/116960/47057.png

[45.17851270178206, 141.24202609062195]
[24093605.77231122, 59883855.66082844]
[24093605, 59883855]
http://cyberjapandata.gsi.go.jp/xyz/std/18/233921/94115.png

[45.178508920402244, 141.24203145503998]
[3298244.7938754167, 7431931.647317333]
[3298244, 7431931]
http://cyberjapandata.gsi.go.jp/xyz/std/15/29030/12883.png

[35.855526613165544, 138.9438772201538]
[6596489.5877508335, 14863863.294634666]
[6596489, 14863863]
http://cyberjapandata.gsi.go.jp/xyz/std/16/58061/25767.png

[35.85550922179457, 138.94389867782593]
[13192979.175501667, 29727726.589269333]
[13192979, 29727726]
http://cyberjapandata.gsi.go.jp/xyz/std/17/116123/51535.png

[35.855500526107654, 138.94389867782593]
[26385958.351003334, 59455453.178538665]
[26385958, 59455453]
http://cyberjapandata.gsi.go.jp/xyz/std/18/232247/103070.png

[35.855500526107654, 138.94390404224396]
[3451877.727625328, 7235274.201133511]
[3451877, 7235274]
http://cyberjapandata.gsi.go.jp/xyz/std/15/28262/13483.png

[30.3359539507533, 130.50427436828613]
[6903755.455250656, 14470548.402267022]
[6903755, 14470548]
http://cyberjapandata.gsi.go.jp/xyz/std/16/56525/26967.png

[30.33593543109004, 130.50427436828613]
[13807510.910501312, 28941096.804534044]
[13807510, 28941096]
http://cyberjapandata.gsi.go.jp/xyz/std/17/113051/53935.png

[30.33593543109004, 130.50427436828613]
[27615021.821002625, 57882193.60906809]
[27615021, 57882193]
http://cyberjapandata.gsi.go.jp/xyz/std/18/226102/107871.png

[30.33593080117366, 130.50427973270416]