#!/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 .
#
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[+-])(?P\d{2}):(?P\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\d{4})- # Year
(?P\d{2})- # month
(?P\d{2})T # day
(?P\d{2}): # hour
(?P\d{2}): # minute
(?P\d{2}) # second
(?:(?:\.(?P\d{6}))?) # microsecond
(?PZ|[+-]\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 ."""
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()