#!/usr/bin/python
# A script to make a mashup of Garmin Training Center histories.
# By Stephen Skory
# stephenskory (aatt) yahoo.com
# last updated 26 June 2008

# ---- DEFAULT VALUES ----- #
# distance between radial bins in meters
radial_binsize = 200
# distance between farthest points in meters
dist_max = 1000000
# distance between unique points in meters
dist_unique = 100
# Set your own central point
# [SetLat,SetLon] = [32,-117] # N.B. This is roughly San Diego, CA...
# read in from XML file or Training History.gtc?
# put in 'XML' or 'gtc'
# Don't worry, on windows it should default to XML on its own.
read_in_type = 'gtc'
# mac os x only, set a specific gtc file, comment to use default location
#gtc_file = '/your/random/path/Training Center.gtc'

# ---- YOU SHOULDN'T NEED TO EDIT ANYTHING BELOW ---- #

import math,string,re,os,sqlite3 as sql,sys

pi = float(math.atan(1) * 4.0) # 3.141...
r = float(6378135.0) # radius of earth in meters
e = math.exp(1) # 2.71727...
gps_conversion =  11930464.71 # this is my best guess as how to convert
                              # the units in the SQLite file into real units

# convert decimal to radians
def DtoRad(deg):
    rad = float(deg) * float(pi) / 180.0
    return rad

# a lat1
# b lon1
# c lat2
# d lon2
# calculate distance between two points
def D(a,b,c,d):
    (a,b,c,d) = (DtoRad(a),DtoRad(b),DtoRad(c),DtoRad(d))
    temp = r * math.acos(math.sin(a)*math.sin(c)+math.cos(a)*math.cos(c)*math.cos(b-d))
    return temp

# with two lists, the IDs of points in a bin, and the list of unique points in that bin
# add to uniqueBins
def addToUniqueBins(bin):
    tempPoints = PointBins[bin]
    tempUnique = uniqueBins[bin]
    for point in tempPoints:
        min = 2 * dist_unique # some big number larger than distance between points
        step = 0
        for unique in tempUnique:
            # in the unlikely ocurrance of a repeated point, zero separation gives D() problems.
            if ((points[point].Lat == unique[0]) and (points[point].Lon == unique[1])):
                dist = 0
            else:
                dist = D(points[point].Lat,points[point].Lon,unique[0],unique[1])
            if (dist<min):
                min = dist
                if (dist<dist_unique):
                    mark = step
            step += 1
        
        #either add the point to an existing point, or make a new one
        if (min < dist_unique):
            # increase the count if we're close to an existing point
            tempUnique[mark] = (tempUnique[mark][0],tempUnique[mark][1],tempUnique[mark][2]+1)
        if (min >= dist_unique):
            # add a new entry to tempUnique if this point is dist_unique or farther from any existing points
            tempUnique.append((points[point].Lat,points[point].Lon,1))
        
    # put in the new tempUnique into UniqueBins
    uniqueBins[bin] = tempUnique
    print 'There are %d unique points in bin %d' % (len(tempUnique),bin)


# search for garmin tcx files, add their names to the list
def Get_TCXs(arg, dirname, names):
    for name in names:
        if (name.endswith('tcx')):
            string = str(name)
            #print string
            TCXnames.append(string)


class Points:

    def __init__(self, Lat, Lon):
        self.Lat = Lat # the latitude of the point
        self.Lon = Lon # the longitude of the point

def openDB(filename):
    global cursor, conn
    conn = sql.connect(filename)
    cursor = conn.cursor()

def findLatLonAvgDB():
    global cursor, conn
    cursor.execute('SELECT AVG(ZLATITUDE),AVG(ZLONGITUDE) FROM ZCDTRACKPOINT')
    line = cursor.fetchone()
    temp = list(line)
    return temp


def selectLatLon():
    global cursor,conn,points,PointBins,totalpoints
    ID = 0
    cursor.execute('SELECT ZLATITUDE,ZLONGITUDE FROM ZCDTRACKPOINT')
    line = cursor.fetchone()
    while line:
        # get the Lat,Lon pair
        [Lat,Lon] = list(line)
        # sometimes it's empty, so we skip it
        if ((Lat) and (Lon)):
            # convert to real Lat,Lon units
            Lat = float(Lat) / gps_conversion
            Lon = float(Lon) / gps_conversion
            #print 'SetLat,Lon = %f %f' %(SetLat,SetLon)
            # calculate distance to the set point
            dist = D(Lat,Lon,SetLat,SetLon)
            # calculate the bin
            bin = int(dist)/radial_binsize
            #print 'bin = %d' % bin
            # only if this point is inside our search radius do we add it
            if (bin <= numberofbins-1):
                # record the point
                points.append(Points(Lat,Lon))
                # bin the point
                PointBins[bin].append(ID)
                ID += 1
        line = cursor.fetchone()
    return ID

# ---- PROCEEDURAL STUFF ----- #

# mac os x only
if ((sys.platform == 'darwin') and (read_in_type == 'gtc')):
    home_dir = os.environ["HOME"]
    try:
        print 'Reading in %s' % gtc_file
    except NameError:
        gtc_file = home_dir + "/Library/Application Support/Garmin/Training Center/Training Center.gtc"
        print 'Reading in %s' % gtc_file

# windows hard set to XML
if (sys.platform == 'win32'):
    read_in_type = 'XML'

# make the list of points
points = []
# contains the list of garmin history files
TCXnames = []
# the lists of unique points in that bin, with Lat,Lon,count
uniqueBins = []
# a list of bins, each bin contains a list of particle IDs in that bin.
PointBins = []

# make bins based on maximum distance and reasonable annulus size. There is probably a more
# elegant way to do this.
numberofbins = int(dist_max/radial_binsize)
for i in range(numberofbins):
    PointBins.append([])
    uniqueBins.append([])

# ----- BEGIN XML ----- #
if (read_in_type == 'XML'):
    # search for tcx files in the current directory
    os.path.walk(os.getcwd(), Get_TCXs, None)
    
    
    # find an intial starting point for the center of our annulus
    fp = open(TCXnames[0],'r')
    line = fp.readline()
    while line:
        if (line.count('Latitude')>0):
            line = re.split('>|<',line)
            SetLat = float(line[2])
            # Longitude is always just after Latitude
            line = fp.readline()
            line = re.split('>|<',line)
            SetLon = float(line[2])
            break;
        line = fp.readline()
    
    fp.close()
    
    # now read in all the files, binning the points by distance from SetLat,SetLon
    ID = 0
    for TCXfile in TCXnames:
        
        fp = open(TCXfile,'r')
        print 'Reading in %s' % TCXfile
        
        line = fp.readline()
        while line:
            if (line.count('Latitude')>0):
                line = re.split('>|<',line)
                Lat = float(line[2])
                # Longitude is always just after Latitude
                line = fp.readline()
                line = re.split('>|<',line)
                Lon = float(line[2])
                dist = D(Lat,Lon,SetLat,SetLon)
                # calculate the bin
                bin = int(dist)/radial_binsize
                # only if this point is inside our search radius do we add it
                if (bin <= numberofbins-1):
                    # record the point
                    points.append(Points(Lat,Lon))
                    # bin the point
                    PointBins[bin].append(ID)
                    ID += 1
            line = fp.readline()
        
        
        fp.close()
    
    
    totalpoints = ID
    print 'read in %d points' % totalpoints

# ----- END XML ------- #

# ----- BEGIN GTC ----- #
if (read_in_type == 'gtc'):
    # open connection to file
    openDB(gtc_file)
    # find the Lat,Lon center point if not already set
    try:
        print 'SetLat,SetLon = %f, %f' % (SetLat,SetLon)
    except NameError:
        [SetLat,SetLon] = findLatLonAvgDB()
        # fix SetLat,Lon
        SetLat = float(SetLat) / gps_conversion
        SetLon = float(SetLon) / gps_conversion
    totalpoints = selectLatLon();
    print 'read in %d points' % totalpoints


# ----- END GTC ------- #

# put the first point in each non-zero bin in PointBins as the first entry of 
# the corresponding member of uniqueBins
# also remove it from the bin so it's not counted twice
for bin in range(numberofbins):
    if (len(PointBins[bin]) > 0):
        uniqueBins[bin].append((points[PointBins[bin][0]].Lat,points[PointBins[bin][0]].Lon,1))
        PointBins[bin].pop(0)


# using the points in PointBins, and the points in uniqueBins, we'll add points to the
# lists in uniqueBins
totalsofar = 0
for bin in range(numberofbins):
    if (len(PointBins[bin]) > 0):
        print 'uniquing bin %d, member count %d, distance from center is %.2f KM' % (bin,len(PointBins[bin]),float(bin)*float(radial_binsize)/1000.0)
        addToUniqueBins(bin)
        totalsofar += len(PointBins[bin])
        print 'Done %d/%d = %f percent\n\n' % (totalsofar,totalpoints,float(totalsofar)/float(totalpoints))


# print out the contents of uniqueBins
fp = open('UniquePoints.txt','w')
fp.write('latitude, longitude, n\n')
flag = 1
outline = 0
gaps = []
gaplines = []
last = 0
for bin in range(numberofbins):
    if (len(uniqueBins[bin]) > 0):
        last = bin
        for i in range(len(uniqueBins[bin])):
            line = '%s, %s, %1.5f\n' % (str(uniqueBins[bin][i][0]),str(uniqueBins[bin][i][1]),math.log(uniqueBins[bin][i][2]+2,e))
            fp.write(line)
            outline += 1
    if ((bin != last) and (gaps.count(last)==0)):
        gaps.append(last)
        gaplines.append(outline)
        flag = 1
    if ((flag == 1) and (bin == last)):
        gaps.append(last)
        gaplines.append(outline)
        flag = 0

fp.close()

# uncomment below if you're interested in splitting the output of your file into distinct geographical
# areas. The numbers printed below represent where the output in FinalPoints.txt gives a jump
# larger than radial_binsize between tracks. If you can't figure out what I mean by studying what's going
# on above with gaps and gaplines, don't worry about it.
#print gaps
#print gaplines