#!/usr/bin/python # # A library for parsing, creating, and transforming GPX files, which # are an XML format for GPS logs. # # Seth Golub http://www.sethoscope.net/geophoto/ # # # 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 2 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. # # I don't bother to provide a copy of the GNU General Public License # along with this program, but you can get one from the Free Software # Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA import os import sys from time import mktime,strptime,tzset from xml.parsers.xmlproc import xmlproc from math import * class TrackPoint: def distance(self, other): pi = 3.14159265358979 lat1=self.lat * pi / 180 lat2=other.lat * pi / 180 lon1=self.lon * pi / 180 lon2=other.lon * pi / 180 # http://en.wikipedia.org/wiki/Great-circle_distance # Of the several formula listed, this one retains the most accuracy at short distances dist = atan(sqrt(pow(cos(lat2)*sin(abs(lon1-lon2)),2) + pow(cos(lat1)*sin(lat2) - sin(lat1)*cos(lat2)*cos(lon1-lon2),2)) / (sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon1-lon2))) dist *= 6372795 # convert from radians to meters return dist # TODO: annotate points that start or end a tracklog so geotagging can consider log boundaries class _GPXParser(xmlproc.Application): def doc_start(self): self.tag_name_stack = [] self.points = [] self.current_point = None def handle_start_tag(self, name, attrs): self.tag_name_stack.append(name) if name == 'trkpt': self.current_point = TrackPoint() self.current_point.lat = float(attrs['lat']) self.current_point.lon = float(attrs['lon']) def handle_end_tag(self, name): popped_name = self.tag_name_stack.pop() if ( name != popped_name ): print 'expected %s, but ended %s' % (popped_name, name) assert(0) if name == 'trkpt': self.points.append(self.current_point) self.current_point = None def handle_data(self, data, start, end): os.environ['TZ'] = 'GMT' # this is awful, but I can't seem to fix it otherwise tzset() if self.tag_name_stack[-1] == 'ele': self.current_point.elevation = float(data[start:end]) elif self.tag_name_stack[-1] == 'time' \ and self.tag_name_stack[-2] == 'trkpt': self.current_point.time_str = data[start:end] self.current_point.time = mktime(strptime(data[start:end], '%Y-%m-%dT%H:%M:%SZ')) # TODO: make the algorithm symmetric # TODO: discard outliers # TODO: scrap this entirely, replace with Kalman smoothing def _moving_average(data, window_size): window_sum = 0 ave = [] for i in range(min(len(data), window_size)): window_sum += data[i] ave.append(window_sum / window_size) for i in range(min(len(data), window_size) - 1, len(data)): window_sum += data[i] window_sum -= data[i-window_size] ave.append(window_sum / window_size) return ave class TrackLog: def __init__(self, filename): self.load_from_file(filename) def load_from_file(self, filename): gpx_parser = _GPXParser() p=xmlproc.XMLProcessor() p.set_application(gpx_parser) p.parse_resource(filename) self.points = gpx_parser.points def distance(self, start=0, end=-1): if len(self.points) <= start + 1: return 0.0 total = 0.0 prev = self.points[start] for point in self.points[start+1:end]: total += point.distance(prev) prev = point return total def find_nearest_in_time(self, target_time): # TODO: interpolate closest = None for point in self.points: if not closest or abs(point.time - target_time) < abs(closest.time - target_time): closest = point return closest # TODO: instead of a moving average, use Kalman smoothing def compute_updown(self, window_size=3): if len(self.points) <= 1: return 0 moving_ave = _moving_average([x.elevation for x in self.points], window_size) climb = 0 descent = 0 prev = moving_ave[0] for point in moving_ave[1:]: diff = point - prev if diff > 0: climb += diff else: descent -= diff prev = point return (climb, descent) def max_elevation(self): return max([x.elevation for x in self.points]) def min_elevation(self): return min([x.elevation for x in self.points]) def instantaneous_grades(self): grade = lambda x,y: (x.elevation - y.elevation) / distance(x,y) return map(grade, self.points[:-1], self.points[1:]) # Each trackpoint has some measurable value to the shape of the # track. To reduce the track, we remove the least valuable point, # one by one, greedily. First we compute all the values, then we # remove points until we feel like stopping. def _compute_point_value(self, i): if i == 0 or i == len(self.points)-1: return sys.maxint A = self.points[i-1] B = self.points[i] C = self.points[i+1] # We assume our coordinate system (Mercator) is square, # which is only true at the equator. Also the world is flat. ACxlen = A.lat - C.lat ACylen = A.lon - C.lon AClen = sqrt(ACxlen * ACxlen + ACylen * ACylen) ABCarea = abs((A.lat - B.lat) * (A.lon - C.lon) - (A.lon - B.lon) * (A.lat - C.lat)) return ABCarea / AClen # gives us the diversion distance def _store_point_value(self, i): self.points[i].value_to_track = self._compute_point_value(i) def prepare_to_reduce(self): for i in range(len(self.points)): self._store_point_value(i) # We could make this faster by keeping a list of points sorted by # value, perhaps in a heap. Optimization last! def reduce_by_one_point(self): min_i = None for i in range(1, len(self.points)-2): # skip first & last if not min_i or self.points[i].value_to_track < self.points[min_i].value_to_track: min_i = i if not min_i: return None point = self.points[min_i] del self.points[min_i] self._store_point_value(min_i) self._store_point_value(min_i-1) return point def save_to_file(self, filename, start=0, numpoints=-1): if numpoints < 0: numpoints = len(self.points) close_when_done = 0 if filename == '-': fh = sys.stdout else: fh = open(filename, 'w') close_when_done = 1 fh.write('\n' '\n' '\n' '\n' ) for point in self.points[start:start+numpoints]: fh.write('\n' % (point.lat, point.lon)) if point.elevation: fh.write('%f\n' % (point.elevation,)) if point.time_str: fh.write('\n' % (point.time_str,)) fh.write('\n') fh.write('\n' '\n' '\n') if close_when_done: fh.close()