Coordenadas a tiles

#!/usr/bin/python
"""Usage:

./bounds.py
"""

# o==object to search
# NW = north-west
# SE = south-east

import math

TILE_SIZE = 256
olat = 45.5
olng = -3.81

olat = 43.1742630005
olng = -2.2827100753

iNWlat = 85.0511287798066 
iNWlng = -180.0

iSElat = -85.0511287798066 
iSElng = 180.0

iX = 0
iY = 0

Y = {}
X = {}

zooms = 16

def tile_number (zoom):
    """ """
    return math.pow(2,zoom)

def circumference (zoom):
    """ """
    return TILE_SIZE*tile_number(zoom)
        
def radius (zoom):
    """ """
    return circumference(zoom) / 2 * math.pi

def latToY(latitudeDegrees,zoom):
    latitude = math.radians(latitudeDegrees);
    y = radius(zoom)/2.0 * math.log((1.0 + math.sin(latitude)) / (1.0 - math.sin(latitude)))
    return y
    
def yToLat(y,zoom):
    """ """
    latitude =  (math.pi/2) - (2 * math.atan(math.exp(-1.0 * y / radius(zoom))))
    return math.degrees(latitude);

def longToX(longitudeDegrees,zoom):
    """ """
    longitude = math.radians(longitudeDegrees)
    return (radius(zoom) * longitude)


def xToLong(x,zoom):
    """ """
    longRadians = x/radius(zoom);
    longDegrees = math.degrees(longRadians)
    rotations = math.floor((longDegrees + 180)/360)
    return longDegrees - (rotations * 360)

def main():
    """ """ 


    aNWlat = iNWlat
    aNWlng = iNWlng
    aSElat = iSElat
    aSElng = iSElng
    for i in range(0,zooms+1):

#
#        en cada nivel de zoom ...
#        los puntos aNW (aNWlat,aNWlng) y aSE (aSElat,aSElng) definen la tile (a) que contiene el punto o (olat,olng)
#
#       aNW-----------------------
#        |                       |
#        |                       |
#        |                       |
#        |       tile (a)        |
#        |                       |
#        |                       |
#        |                       |
#        |    o                  |
#        |                       |
#        -----------------------aSE
#
#
#        # partimos la tile (a) en cuatro tiles: tile (n0), tile (n1), tile (n2), tile (n3)
#
#        # nNW (nNWlat,nNWlng)
#
#        
#      nNW0---------nNW1----------
#        |           |           |
#        | tile (n0) | tile (n1) |
#        |           |           |
#        |          nNW3         |
#      nNW2---------nSE0--------nSE1
#        |           |           |
#        | tile (n2) | tile (n3) |
#        |           |           |
#        |    o      |           |
#        -----------nSE2--------nSE3
#

        nNWlat = [aNWlat,aNWlat,yToLat((latToY(aNWlat,i)+latToY(aSElat,i))/2,i),yToLat((latToY(aNWlat,i)+latToY(aSElat,i))/2,i)]
        nNWlng = [aNWlng,(aNWlng+aSElng)/2.,aNWlng,(aNWlng+aSElng)/2.]

        nSElat = [yToLat((latToY(aNWlat,i)+latToY(aSElat,i))/2,i),yToLat((latToY(aNWlat,i)+latToY(aSElat,i))/2,i),aSElat,aSElat]
        nSElng = [(aNWlng+aSElng)/2.,aSElng,(aNWlng+aSElng)/2.,aSElng]
        
        for j in range(0,4):
            if nSElat[j] < olat < nNWlat[j] and nNWlng[j] < olng < nSElng[j]:
                aNWlat = nNWlat[j]
                aNWlng = nNWlng[j]
                aSElat = nSElat[j]
                aSElng = nSElng[j]

#           y elegimos aquella que contiene al punto o (olat,olng)
#
#              -------------------------
#              |           |           |
#              |           |           |
#              |           |           |
#              |           |           |
#             aNW-----------------------
#              |           |           |
#              | tile (a)  |           |
#              |           |           |
#              |    o      |           |
#              -----------aSE-----------
#

                lngdif = (iSElng-iNWlng)/tile_number(i+1)
                X = int((olng-iNWlng)/lngdif)
                
                latdif = (iSElat-iNWlat)/tile_number(i+1)
                ydif = (latToY(iNWlat,i)-latToY(iSElat,i))
                Y = int((olat-iNWlat)/ydif)                
                print zooms-i,nNWlat[j],nNWlng[j],nSElat[j],nSElng[j],'-----  ',X
                """
                falseEasting = -1.0 * circumference(i)/2.0
                falseNorthing =  circumference(i)/2.0
                thisy = latToY(iNWlat,i)+falseEasting
                print thisy
                """


if __name__ == "__main__":
    main()