[5315] | 1 | #!/usr/bin/env python |
---|
| 2 | #------------------------------------------------------- |
---|
| 3 | # Translates between lat/long and the slippy-map tile |
---|
| 4 | # numbering scheme |
---|
| 5 | # |
---|
| 6 | # http://wiki.openstreetmap.org/index.php/Slippy_map_tilenames |
---|
| 7 | # |
---|
| 8 | # Written by Oliver White, 2007 |
---|
| 9 | # This file is public-domain |
---|
| 10 | #------------------------------------------------------- |
---|
| 11 | from math import * |
---|
| 12 | |
---|
| 13 | def numTiles(z): |
---|
| 14 | return(pow(2,z)) |
---|
| 15 | |
---|
| 16 | def sec(x): |
---|
| 17 | return(1/cos(x)) |
---|
| 18 | |
---|
| 19 | def latlon2relativeXY(lat,lon): |
---|
| 20 | x = (lon + 180) / 360 |
---|
| 21 | y = (1 - log(tan(radians(lat)) + sec(radians(lat))) / pi) / 2 |
---|
| 22 | return(x,y) |
---|
| 23 | |
---|
| 24 | def latlon2xy(lat,lon,z): |
---|
| 25 | n = numTiles(z) |
---|
| 26 | x,y = latlon2relativeXY(lat,lon) |
---|
| 27 | return(n*x, n*y) |
---|
| 28 | |
---|
| 29 | def tileXY(lat, lon, z): |
---|
| 30 | x,y = latlon2xy(lat,lon,z) |
---|
| 31 | return(int(x),int(y)) |
---|
| 32 | |
---|
[5850] | 33 | def xy2latlon(x,y,z): |
---|
| 34 | n = numTiles(z) |
---|
| 35 | relY = y / n |
---|
| 36 | lat = mercatorToLat(pi * (1 - 2 * relY)) |
---|
| 37 | lon = -180.0 + 360.0 * x / n |
---|
| 38 | return(lat,lon) |
---|
| 39 | |
---|
[5315] | 40 | def latEdges(y,z): |
---|
| 41 | n = numTiles(z) |
---|
| 42 | unit = 1 / n |
---|
| 43 | relY1 = y * unit |
---|
| 44 | relY2 = relY1 + unit |
---|
| 45 | lat1 = mercatorToLat(pi * (1 - 2 * relY1)) |
---|
| 46 | lat2 = mercatorToLat(pi * (1 - 2 * relY2)) |
---|
| 47 | return(lat1,lat2) |
---|
| 48 | |
---|
| 49 | def lonEdges(x,z): |
---|
| 50 | n = numTiles(z) |
---|
| 51 | unit = 360 / n |
---|
| 52 | lon1 = -180 + x * unit |
---|
| 53 | lon2 = lon1 + unit |
---|
| 54 | return(lon1,lon2) |
---|
| 55 | |
---|
| 56 | def tileEdges(x,y,z): |
---|
| 57 | lat1,lat2 = latEdges(y,z) |
---|
| 58 | lon1,lon2 = lonEdges(x,z) |
---|
| 59 | return((lat2, lon1, lat1, lon2)) # S,W,N,E |
---|
| 60 | |
---|
| 61 | def mercatorToLat(mercatorY): |
---|
| 62 | return(degrees(atan(sinh(mercatorY)))) |
---|
| 63 | |
---|
[5850] | 64 | def tileSizePixels(): |
---|
| 65 | return(256) |
---|
| 66 | |
---|
[5767] | 67 | def tileLayerExt(layer): |
---|
| 68 | if(layer in ('oam')): |
---|
| 69 | return('jpg') |
---|
| 70 | return('png') |
---|
[5315] | 71 | |
---|
[5767] | 72 | def tileLayerBase(layer): |
---|
| 73 | layers = { \ |
---|
[18811] | 74 | "tah": "http://cassini.toolserver.org:8080/http://a.tile.openstreetmap.org/+http://toolserver.org/~cmarqu/hill/", |
---|
| 75 | #"tah": "http://tah.openstreetmap.org/Tiles/tile/", |
---|
[5767] | 76 | "oam": "http://oam1.hypercube.telascience.org/tiles/1.0.0/openaerialmap-900913/", |
---|
| 77 | "mapnik": "http://tile.openstreetmap.org/mapnik/" |
---|
| 78 | } |
---|
| 79 | return(layers[layer]) |
---|
| 80 | |
---|
| 81 | def tileURL(x,y,z,layer): |
---|
| 82 | return "%s%d/%d/%d.%s" % (tileLayerBase(layer),z,x,y,tileLayerExt(layer)) |
---|
| 83 | |
---|
[5315] | 84 | if __name__ == "__main__": |
---|
| 85 | for z in range(0,17): |
---|
| 86 | x,y = tileXY(51.50610, -0.119888, z) |
---|
| 87 | |
---|
| 88 | s,w,n,e = tileEdges(x,y,z) |
---|
| 89 | print "%d: %d,%d --> %1.3f :: %1.3f, %1.3f :: %1.3f" % (z,x,y,s,n,w,e) |
---|
| 90 | #print "<img src='%s'><br>" % tileURL(x,y,z) |
---|