#!/usr/bin/env python import ogr import math import osr def lon_to_UTMzone(lon): return(int(lon + 180) / 6 + 1) def Terraserver_resolution_to_scale(resolution): return(int(math.log(resolution, 2) + 10)) def Terraserver_scale_to_resolution(scale): return(int(math.pow(2, scale - 10))) def Terraserver_resolution_to_multiplier(resolution): return(resolution * 200) def Terraserver_scale_to_multiplier(scale): return(Terraserver_resolution_to_multiplier(Terraserver_scale_to_resolution(scale))) zones = [ (-180, 1), (-174, 2), (-168.0001, 2), (-168, 3), ( -6, 30), ( 0, 31), ( 12, 33), ] scales = [ (10, 1, 200), (11, 2, 400), (12, 4, 800), (13, 8, 1600), (17, 128, 25600), (19, 512, 102400) ] #for (lon, zone) in zones: # print lon, zone, lon_to_UTMzone(lon) #for (scale, resolution, multiplier) in scales: # print scale, resolution, multiplier, Terraserver_scale_to_multiplier(scale), Terraserver_resolution_to_scale(resolution) # reference_srs = osr.SpatialReference() # reference_srs.ImportFromWkt(reference_projection) # reference_srsLatLong = reference_srs.CloneGeogCS() # ct_to_latlon = osr.CoordinateTransformation(reference_srs, reference_srsLatLong) # ct_from_latlon = osr.CoordinateTransformation(reference_srsLatLong, reference_srs) # return(ct_from_latlon.TransformPoint(lon, lat)) def Terraserver_lonlat_to_tileURL(lon, lat, resolution): multiplier = Terraserver_resolution_to_multiplier(resolution) scale = Terraserver_resolution_to_scale(resolution) zone = lon_to_UTMzone(lon) #print 'resolution=%d, multiplier=%d, scale=%d' % (resolution, multiplier, scale) #print 'lon=%f, lat=%f, zone=%d' % (lon, lat, zone) UTM_srs = osr.SpatialReference() UTM_srs.SetUTM(zone) lonlat_srs = osr.SpatialReference() lonlat_srs.SetWellKnownGeogCS('NAD83') UTM_to_lonlat = osr.CoordinateTransformation(UTM_srs, lonlat_srs) lonlat_to_UTM = osr.CoordinateTransformation(lonlat_srs, UTM_srs) (easting, northing, altitude) = lonlat_to_UTM.TransformPoint(lon, lat) X = easting / multiplier Y = northing / multiplier URL = 'http://terraserver-usa.com/tile.ashx?S=%d&T=1&X=%d&Y=%d&Z=%d' % ( scale, easting / multiplier, northing / multiplier, zone ) return URL # N 40.442769 W 86.937551, 12, X=631, Y=5596, Z=16 print Terraserver_lonlat_to_tileURL(-86.937551, 40.442769, 8)