import csv
import math
import os.path
import urllib

### Function to Read in trace from csv ###
def traceImportCSV(fname):
    import csv
    trace = []
    for lat,lon in csv.reader(open(fname,'r')):
        trace.append([float(lat),float(lon)])

    return trace

### Function to Read in trace from GPX ###
def traceImportGPX(fname):
    import re
    trace = []
    for line in open(fname,'r'):
        # This will match "lat" and "lon" flags as long as they are on the same line.
        # Does not differentiate between trkpt and other points,so make sure that 
        # the file is just the track you are interested in.
        matchLat = re.search(r'.* lat=\"(\S*)\".*',line)
        matchLon = re.search(r'.* lon=\"(\S*)\".*',line)                
        if matchLat != None and matchLon != None:
            lat = matchLat.group(1)
            lon = matchLon.group(1)
            trace.append([float(lat),float(lon)])

    return trace

### Function to get trace boundries ###
def traceBoundaries(trace):

    lat = zip(*trace)[0]
    lon = zip(*trace)[1]

    return {"north":max(lat),"south":min(lat),"east":max(lon),"west":min(lon)}


### Function to convert lat,lon degrees to tile x/y number ### 
def deg2num(lat_deg, lon_deg, zoom):
  lat_rad = math.radians(lat_deg)
  n = 2.0 ** zoom
  xtile = int((lon_deg + 180.0) / 360.0 * n)
  ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
  return (xtile, ytile)

### Function to convert xy tile to NW corner of tile ###
def num2deg(xtile, ytile, zoom):
  n = 2.0 ** zoom
  lon_deg = xtile / n * 360.0 - 180.0
  lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
  lat_deg = math.degrees(lat_rad)
  return (lat_deg, lon_deg)

### Determine tile range given boundaries and zoom ###
def determineTileRange(boundaries,zoom):
    Xmax,Ymin = deg2num(boundaries["north"],boundaries["east"],zoom)
    Xmin,Ymax = deg2num(boundaries["south"],boundaries["west"],zoom)
    return {"xMin":Xmin,"xMax":Xmax,"yMin":Ymin,"yMax":Ymax}

### Take a tile range and download them (if not locally present) ###
def getTiles(xyRange,zoom):
    #set acive directory to that of the script
    currentdir = os.curdir
    tileDir = os.path.join(currentdir,"tiles")

    tileServerUrl = "http://tile.stamen.com/watercolor/"

    #create a list of all the x and y coordinates to download
    xRange = range(xyRange["xMin"],xyRange["xMax"]+1)
    yRange = range(xyRange["yMin"],xyRange["yMax"]+1)

    for x in xRange:
        for y in yRange:
            #define the file name
            tileFileName = str(y)+".jpg"

            #define the local path as well as the complete path to the local and remote files
            localPath = os.path.join(tileDir,str(zoom),str(x))
            localFile = os.path.join(localPath,tileFileName)
            remoteFile = tileServerUrl+str(zoom)+"/"+str(x)+"/"+str(y)+".jpg"

            #check if the file exists locally
            if not os.path.isfile(localFile):                
                print "retrieving "+remoteFile
                #if local directory doesn't yet exist, make it
                if not os.path.isdir(localPath):
                    os.makedirs(localPath)
                #retrieve the file from the server and save it    
                urllib.urlretrieve(remoteFile,localFile)

### Merge tiles into one image ###
def mergeTiles(xyRange,zoom,filename):
    import Image
    tileSize = 256
    tileDir = os.path.join(os.curdir,"tiles",str(zoom))

    out = Image.new( 'RGB', ((xyRange["xMax"]-xyRange["xMin"]+1) * tileSize, (xyRange["yMax"]-xyRange["yMin"]+1) * tileSize) ) 

    imx = 0;
    for x in range(xyRange["xMin"], xyRange["xMax"]+1):
        imy = 0
        for y in range(xyRange["yMin"], xyRange["yMax"]+1):
            tileFile = os.path.join(tileDir,str(x),str(y)+".jpg")
            tile = Image.open(tileFile)
            out.paste( tile, (imx, imy) )
            imy += tileSize
        imx += tileSize

    out.save(os.path.join(os.curdir,filename))

### Draw Path Image ###
def drawTraceMask(trace,xResolution,yResolution,traceBoundaries,zoom,filename):
    import Image
    import ImageDraw

    # Get XY number of NW and SE corner tiles
    xy_nw = deg2num(traceBoundaries["north"],traceBoundaries["west"],zoom)
    xy_se = deg2num(traceBoundaries["south"],traceBoundaries["east"],zoom)

    # get lat,lon of corners
    # (since the function returns the NW corner of a tile,
    # we need lat,lon of X+1,Y+1 for the SE corner)
    NW = num2deg(xy_nw[0],xy_nw[1],zoom)
    SE = num2deg(xy_se[0]+1,xy_se[1]+1,zoom)

    # The image boundaries are actually different, because
    # they are the boundaries of the tiles, not the trace
    # define the new boundaries
    mapBoundaries = {}
    mapBoundaries["north"] = NW[0]
    mapBoundaries["south"] = SE[0]
    mapBoundaries["west"] = NW[1]
    mapBoundaries["east"] = SE[1]

    # Offset to ensure that NW corner is 0,0
    latOffset = -(mapBoundaries["north"])
    latDivisor = mapBoundaries["north"]-mapBoundaries["south"]
    lonOffset = -(mapBoundaries["west"])
    lonDivisor = mapBoundaries["east"]-mapBoundaries["west"]

    out = Image.new( 'RGB', (xResolution, yResolution) )
    draw = ImageDraw.Draw(out)
    
    firstRun = True
    for lat,lon in trace:
        # Convert zeroed lat,lon into x,y coordinates
        # this will need correction for northern hemisphere        
        x = abs(int(xResolution*((lon + lonOffset)/lonDivisor)))
        y = abs(int(yResolution*((lat + latOffset)/latDivisor)))

        if firstRun:
            firstRun = False
        else:
            draw.line((x,y,xPrev,yPrev),fill="white",width=10)
        xPrev = x
        yPrev = y
    del draw
    out.save(os.path.join(os.curdir,filename))


# define parameters
zoom = 12
trace = traceImportCSV("test_trace_4.csv")
# determine the boundaries of the trace
boundaries_trace = traceBoundaries(trace)
# determine xy numbers of boundary tiles
tileRange = determineTileRange(boundaries_trace,zoom)
# count number of tiles in x and y direction
xTiles = tileRange["xMax"]-tileRange["xMin"]
yTiles = tileRange["yMax"]-tileRange["yMin"]
numTiles = xTiles*yTiles
# download tiles if needed
getTiles(tileRange,zoom)
# merge tiles into oneimage
mergeTiles(tileRange,zoom,"output-map.jpg")
# draw the path
# Note: the range in "tilerange" refers to the NW corner, but our image extends on block further
drawTraceMask(trace,256*(xTiles+1),256*(yTiles+1),boundaries_trace,zoom,"output-mask.jpg")

print "boundaries: ",boundaries_trace
print "tile X,Y range: ",tileRange
print "x tiles: ",xTiles
print "y tiles: ",yTiles
print "num tiles: ",numTiles