#!/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)