#!/usr/bin/env python
# -*- coding: UTF-8 -*-
# 
#   SnelheidHistogram.py
#   Created by Maarten Sneep on 05-12-09
#
#   © Maarten Sneep, 2009
# 
#   This program is free software: you can redistribute it and/or modify
#   it under the terms of the GNU General Public License as published by
#   the Free Software Foundation, either version 3 of the License, or
#   (at your option) any later version.
#
#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#
#   You should have received a copy of the GNU General Public License
#   along with this program.  If not, see <http://www.gnu.org/licenses/>.
#

from __future__ import division, with_statement
import os, sys, os.path
import re
import random
import tempfile
import zipfile
from datetime import tzinfo, timedelta, datetime
import math
import numpy as np
from optparse import OptionParser
from xml.parsers import expat

__author__    = "Maarten Sneep"
__version__   = "1.0"
__date__      = "05-12-09"
__copyright__ = u"Copyright © 2009 Maarten Sneep"
__license__   = "GPL-3"
__program__   = "SnelheidHistogram"

"""This program SnelheidHistogram takes a GPX file, 
and creates a histogram of the speed of travel."""

class Histogram(object):
    '''Gather a histogram'''
    def __init__(self, speed=None, filename=None, comment=None, lang=None, verbose=None):
        speed = float(speed) if speed is not None else 1.0
        if speed <= 0: speed = 1.0
        self._lang = lang.lower() if lang else 'nl' 
        self._nbins = int(200/speed)
        self._maxspeed = -100.0
        self._speedbin = 200/self._nbins
        self._distances = np.zeros(self._nbins)
        self._times = np.zeros(self._nbins)
        self._totaldistance = 0.0
        self._totaltime = 0.0
        self._filename = filename
        self._comment = comment
        self._pointcount = 0
        self._verbose = True if verbose else False
    
    def addpoint(self, distance=None, timedifference=None):
        if timedifference <= timedelta(0):
            return
        # calculate time difference in seconds (float)
        t = timedifference.seconds + timedifference.microseconds/1000000
        speed = 3.6*(distance/t) # in km/hour
            
        # range checks
        speedbin = int(speed/self._speedbin)
        if speedbin >= self._nbins:
            if self._verbose:
                sys.stderr.write("Speed %f is out of range\nIgnoring datapoint.\n" % speed)
            return
        
        # add values to histogram
        self._distances[speedbin] += distance
        self._times[speedbin] += t
        
        # bookkeeping
        self._totaldistance += distance
        self._totaltime += t
        
        if self._maxspeed < speed:
            self._maxspeed = speed
        self._pointcount +=1
        
    def maxspeed(self):
        return self._maxspeed
    def distance(self):
        return self._totaldistance
    def time(self):
        return self._totaltime
    def average_speed(self):
        return 3.6*(self._totaldistance/self._totaltime)
    def average_moving_speed(self):
        return 3.6*(self._totaldistance/(self._totaltime-self._times[0]))
    def maxspeedbin(self, max=None):
        maxspd = min(max, self.maxspeed()) if max else self.maxspeed()
        return int(math.floor(maxspd/self._speedbin))
    def most_common(self, max=None):
        times = self.time_distribution(max)
        speed = self.speedaxis(max)
        idx = np.argmax(times[1:])
        return speed[idx+1]
    def numpoints(self):
        return self._pointcount
    def speedaxis(self, max=None):
        return np.arange(self._speedbin/2, 200.0, self._speedbin)[0:self.maxspeedbin(max)+1]
    def distance_distribution(self,max=None):
        return self._distances[0:self.maxspeedbin(max)+1]
    def time_distribution(self,max=None):
        return self._times[0:self.maxspeedbin(max)+1]
    def plot(self, filename, max=None, keep_zero=None, smallfig=None):
        '''Create a plot of the histograms'''
        from matplotlib import pyplot as plt
        from matplotlib.font_manager import FontProperties
        
        distances = self.distance_distribution(max)
        times = self.time_distribution(max)
        speed = self.speedaxis(max)
        
        if keep_zero:
            time_speed = speed
        else:
            times = times[1:]
            time_speed = speed[1:]
        
        fnt = FontProperties()
        fnt.set_family('sans-serif')
        if os.uname()[1] == 'lumpy-gravy.local':
            fnt.set_name('Myriad Pro')
        else:
            fnt.set_name('Verdana')
        
        if smallfig:
            fnt.set_size(9)
        else:
            fnt.set_size(11)
        notefnt = fnt.copy()
        
        if smallfig:
            notefnt.set_size(7)
        else:
            notefnt.set_size(9)
        
        if smallfig:
            fig = plt.figure(figsize=(4,2.7), dpi=96)
        else:
            fig = plt.figure(figsize=(6,4), dpi=96)
        ax1 = fig.add_subplot(111)
        
        dist_line = ax1.plot(speed, distances/1000, 'b')
        ax1.plot(speed, distances/1000, 'bo')
        
        for tick in ax1.xaxis.get_major_ticks():
            tick.label1.set_fontproperties(fnt)
        for tick in ax1.yaxis.get_major_ticks():
            tick.label1.set_fontproperties(fnt)
        
        ax2 = ax1.twinx()
        
        time_line = ax2.plot(time_speed, times/60, 'r')
        ax2.plot(time_speed, times/60, 'ro')
        
        for tick in ax2.yaxis.get_major_ticks():
            tick.label2.set_fontproperties(fnt)
        
        if self._comment:
            ax1.set_title(self._comment, fontproperties=fnt)
        
        if max: ax1.set_xlim(0,max)
        
        if self._lang == 'nl':
            ax1.set_xlabel('Snelheid in km/uur', fontproperties=fnt)
            ax1.set_ylabel('Afstand in km', fontproperties=fnt)
            ax2.set_ylabel('Tijd in minuten', rotation=-90, fontproperties=fnt)
            fig.legend((dist_line, time_line),
                ('Afstand (links)', 'Tijd (rechts)'), 
                (0.65, 0.7), 
                prop=notefnt, 
                fancybox=False, 
                frameon=False)
            plt.figtext(0.15, 0.7, 
                "Afstand: %.1f km\nTijd: %d min\nGemiddeld bew. snelheid: %.1f km/uur\nGemiddelde snelheid: %.1f km/uur" % (self.distance()/1000, 
                    self.time()/60, 
                    self.average_moving_speed(), 
                    self.average_speed()), 
                fontproperties=notefnt)
        else:
            ax1.set_xlabel('Speed in km/hour', fontproperties=fnt)
            ax1.set_ylabel('Distance in km', fontproperties=fnt)
            ax2.set_ylabel('Time in minutes', rotation=-90, fontproperties=fnt)
            fig.legend((dist_line, time_line),
                ('Distance (left)', 'Time (right)'), 
                (0.65, 0.7), 
                prop=notefnt, 
                fancybox=False, 
                frameon=False)
            plt.figtext(0.15, 0.7, 
                "Distance: %.1f km\nTime: %d min\nAverage moving speed: %.1f km/uur\nOverall average speed: %.1f km/uur" % (self.distance()/1000, 
                    self.time()/60, 
                    self.average_moving_speed(), 
                    self.average_speed()), 
                fontproperties=notefnt)
        
        l = plt.axvline(x=self.average_moving_speed(), linewidth=1, color='0.5')
        l = plt.axvline(x=self.average_speed(), linewidth=1, color='0.75')
        
        name,ext = os.path.splitext(filename)
        if ext[1:] in ["emf", "eps", "pdf", "png", "ps", "raw", "rgba", "svg", "svgz"]:
            plt.savefig(filename, 
                dpi=96, orientation='landscape', 
                format=ext[1:], bbox_inches='tight')
        else:
            plt.show()
        
    def __str__(self):
        '''Print a histogram'''
        distances = self.distance_distribution()
        times = self.time_distribution()
        speed = self.speedaxis()
        
        result = ["# Produced by "+__program__]
        result.append("# "+__copyright__)
        result.append("# Licensed under "+__license__)
        if self._filename: 
            result.append("# file: " + self._filename)
        if self._comment: 
            result.append("# " + self._comment)
        result.append("# speed in km/hour, distance in m, time in s")
        result.append("# total distance: %f km" % (self.distance()/1000))
        result.append("# total time: %f minutes" % (self.time()/60))
        # Moving average (leave out time standing still)
        result.append("# (Moving) average speed: %f km/hour" %
                       (self.average_moving_speed()))
        # overall average
        result.append("# (Overall) average speed: %f km/hour" %
                       (self.average_speed()))
        result.append("# number of data points %d" % (self.numpoints()))
        result.append("#")
        result.append("speed\ttimes\tdistances")
        
        for ii in range(len(times)):
            result.append("%.1f\t%.0f\t%.1f" %(speed[ii], times[ii], distances[ii]))
        
        return "\n".join(result)+"\n"
    
class TZ(tzinfo):
    """time zone object"""
    
    def __init__(self, stroffset):
        super(TZ,self).__init__()
        if stroffset.upper() == "Z":
            self._offset = timedelta(0)
            self._name = "UTC"
        else:
            m = re.search(r'''(?P<sign>[+-])(?P<hour>\d{2}):(?P<minute>\d{2})''', 
                          stroffset, re.I)
            sign = 1 if m.group('sign') == '+' else -1
            self._offset = timedelta(hour=int(m.group('hour'))*sign, 
                                     minute=int(m.group('minute'))*sign)
            self._name = "unknown"
    
    def utcoffset(self, dt):
        return self._offset

    def tzname(self, dt):
        return self._name

    def dst(self, dt):
        return timedelta(0)

class Parser(object):
    """This is a minimal GPX parser based on ExPat"""
    
    def __init__(self, histogram=None):
        self._parser = expat.ParserCreate()
        self._parser.StartElementHandler = self.start
        self._parser.EndElementHandler = self.end
        self._parser.CharacterDataHandler = self.data
        self._parser.buffer_text = True
        self._loc1 = None
        self._loc2 = None
        self._time1 = None
        self._time2 = None
        self._elevation1 = None
        self._elevation2 = None
        self.elevation = False
        self.time = False
        self.ontrack = False
        self._histogram = histogram
        self.re = re.compile(r'''(?P<year>\d{4})- # Year
                (?P<month>\d{2})- # month
                (?P<day>\d{2})T # day
                (?P<hour>\d{2}): # hour
                (?P<minute>\d{2}): # minute
                (?P<second>\d{2})  # second
                (?:(?:\.(?P<microsecond>\d{6}))?) # microsecond
                (?P<tzone>Z|[+-]\d{2}:\d{2}) # timezone (+ Zulu)''', 
                re.I|re.X)
    
    def setHistogram(self, speed=None, climbrate=None):
        if self._histogram is not None:
            del self._histogram
        self._histogram = Histogram(speed=speed, climbrate=climbrate)
    
    def histogram(self):
        return self._histogram
    
    def parse(self, file):
        '''Kick off the parser'''
        with open(file, 'r') as fp:
            self._parser.ParseFile(fp)
    
    def __del__(self):
        '''Kill the parser.'''
        del self._parser # get rid of circular references

    def start(self, tag, attrs):
        '''ExPat StartElementHandler.'''
        if tag == "trkpt" and self.ontrack:
            # put the lat and lon into the appropriate locations
            self.settrackpointlocation(attrs)
        elif tag == "ele" and self.ontrack:
            # The next data element will be an elevation
            self.elevation = True
        elif tag == "time" and self.ontrack:
            # The next data element will be a time
            self.time = True
        elif tag == "trkseg" and self.ontrack:
            # reset the location when starting a new track segment
            # This will avoid most jumps
            self._loc1 = None
        elif tag == "trk":
            self.ontrack = True

    def end(self, tag):
        '''ExPat EndElementHandler'''
        if tag == "trkpt" and self.ontrack:
            # A complete trackpoint is in: process it.
            self.processtrackpoint()
        elif tag == "ele" and self.ontrack:
            # We have just processed an elevation
            self.elevation = False
        elif tag == "time" and self.ontrack:
            # We have just processed a time
            self.time = False
        elif tag == "trk":
            self.ontrack = False

    def data(self, data):
        '''ExPat CharacterDataHandler'''
        if not data.strip():
            return
        if self.elevation and self.ontrack:
            # the current character data is an elevation
            self.settrackpointelevation(data)
        elif self.time and self.ontrack:
            # the current character data is a time
            self.settrackpointtime(data)
    
    def settrackpointelevation(self, data):
        '''Parse the elevation information'''
        self._elevation2 = float(data.strip())
    
    def settrackpointtime(self, data):
        '''Parse the time information 
        
in the character data of the time element. 
The string format is the ISO time format 
(including timezone and potentially microseconds).
Parsing is done using a regular expression (created at init time)'''
        m = self.re.search(data)
        tz = TZ(m.group('tzone'))
        micro = 0 if m.group('microsecond') is None else int(m.group('microsecond'))
        self._time2 = datetime(int(m.group('year')), 
                               int(m.group('month')), 
                               int(m.group('day')),
                               int(m.group('hour')),
                               int(m.group('minute')), 
                               int(m.group('second')), 
                               micro, tz)
    
    def settrackpointlocation(self, attrs):
        '''Set the second location'''
        self._loc2 = (float(attrs['lat']), float(attrs['lon']))
    
    def processtrackpoint(self):
        '''All information for a track-point is in.

Process the information that has been gathered. If the first
location is set, we have a pair to process. At the end of this
method we shift the data from the second to the first point.'''
        if self._loc1 is not None:
            d = distance(self._loc1, self._loc2)
            timediff = self._time2 - self._time1
            self._histogram.addpoint(distance=d, 
                timedifference=timediff)
        self._loc1 = self._loc2
        self._time1 = self._time2
        self._elevation1 = self._elevation2
    
def Main():
    """Main program "SnelheidHistogram". """
    (options, args) = cmdline()
    
    hist = Histogram(speed=options.binsize, 
                     filename="; ".join(args), 
                     comment=options.comment,
                     lang=options.lang,
                     verbose=options.verbose)
    for arg in args:
        if os.path.isfile(arg):
            p = Parser(hist)
            p.parse(arg)
            del p
    
    if options.outfile:
        if options.outfile == '-':
            sys.stdout.write(str(hist))
        else:
            with open(options.outfile, 'w') as fp:
                fp.write(str(hist))
    
    if options.plot:
        hist.plot(options.plot, 
            max=options.max, 
            keep_zero=options.keep_zero,
            smallfig=options.smallfig)
    del hist
    
def cmdline():
    """The command-line parser object"""
    usage=u"""usage: %prog [options] infile [infile] [...]
purpose: read the gpx-file, and produce a histogram of the speed you did

Copyright © 2009  Maarten Sneep

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>."""
    
    parser = OptionParser(usage=usage, 
                          version="%prog "+__version__+" ("+__date__+")")
    parser.add_option("-v", "--verbose", action="store_true", 
        dest="verbose", default=False, 
        help="Be chatty")
    parser.add_option("-b", "--binsize", dest="binsize", type="float", 
        default=1.0, 
        help="The bin size for the histogram, defaults to 1 km/hour.")
    parser.add_option("-m", "--max", dest="max", type="float",
        default=65.0, help="Give the expected maximum speed.")
    parser.add_option("-c", "--comment", dest="comment", default="", 
        metavar="STRING", help="A comment to put into the output")
    parser.add_option("-o", "--output", dest="outfile",
        default="", metavar="FILE", 
        help="Put the histogram textual output in FILE.")
    parser.add_option("-p", "--plot", dest="plot", default="", 
        metavar="FILE", help="Plot the graph to FILE. Format is determined from the extension.")
    parser.add_option("-k", "--keep-zero", action="store_true", 
        dest="keep_zero", default=False, 
        help="Include the time that you did not move in the plot")
    parser.add_option("-s", "--small", action="store_true", 
        dest="smallfig", default=False,
        help="Make a small figure")
    parser.add_option("-l", "--lang", dest="lang", default="en",
        metavar="LANGUAGE", help="The language of the output ('nl' or 'en')")
    return parser.parse_args()

def distance(loc1, loc2):
    """Calculate geodesic distance (in m) between two points specified by
latitude/longitude (in numeric degrees) using Vincenty inverse formula 
for ellipsoids. 

Arguments are tuples of (lat, lon) pairs. 

Code based on a Javascript version by © 2002-2008 Chris Veness. 

That code contains the following: "You are welcome to re-use these scripts 
[without any warranty express or implied] provided you retain my copyright notice 
and when possible a link to my website (under a LGPL license)."

The originating web page: http://www.movable-type.co.uk/scripts/latlong-vincenty.html
"""
    a = 6378137 # in m
    b = 6356752.3142 # in m
    f = 1/298.257223563
    deg2rad = math.pi / 180
    lat1, lon1 = (loc1[0] * deg2rad, loc1[1] * deg2rad)
    lat2, lon2 = (loc2[0] * deg2rad, loc2[1] * deg2rad)
    L = (lon2 - lon1)
    U1 = math.atan((1-f) * math.tan(lat1))
    U2 = math.atan((1-f) * math.tan(lat2))
    sinU1 = math.sin(U1)
    cosU1 = math.cos(U1)
    sinU2 = math.sin(U2)
    cosU2 = math.cos(U2)
    Lambda = L
    LambdaP = 2*math.pi
    iterLimit = 20
    
    while math.fabs(Lambda-LambdaP) > 1e-12 and iterLimit > 0:
        sinLambda = math.sin(Lambda)
        cosLambda = math.cos(Lambda)
        sinSigma = math.sqrt((cosU2*sinLambda)**2 + (cosU1*sinU2-sinU1*cosU2*cosLambda)**2)
        if (sinSigma==0): return 0 # co-incident points
        cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda
        sigma = math.atan2(sinSigma, cosSigma)
        # sinSigma ≠ 0, already caught by condition above.
        sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
        cosSqAlpha = 1 - sinAlpha**2
        try:
            cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha
        except:
            cos2SigmaM = 0 # equatorial line: cosSqAlpha=0
        C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
        LambdaP = Lambda;
        Lambda = L + (1-C) * f * sinAlpha * (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM**2)))
        iterLimit -= 1
    
    if iterLimit == 0:
        return None  # formula failed to converge
    
    uSq = cosSqAlpha * (a**2 - b**2) / (b**2)
    A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)))
    B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)))
    deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM**2) - B/6*cos2SigmaM*(-3+4*sinSigma**2)*(-3+4*cos2SigmaM**2)))
    return b*A*(sigma-deltaSigma);

if __name__ == "__main__":
    Main()