gpxmatch.py

gpxmatch.py

#!/usr/bin/env python

# This hacky script takes a GPX file and compares it to OSM ways in
# its surroundings. It looks for points that are not close to an 
# OSM way. If it finds any, it will add them to an output GPX file.
# 
# Good:
#
# * It works
#
# Bad:
#
# * It is slow
# * It is ugly
# * It compares GPX track points to way-nodes, so it breaks on sparsely noded ways.
#
# (c) Martijn van Exel 2016. WTFPL.

import sys
import gpxpy
import overpass
from haversine import haversine

def usage():
	print('Usage: gpxmatch.py infile.gpx outfile.gpx threshold_in_m')

def get_bbox(points, buffer=0):
	west = 180
	south = 90
	east = -180
	north = -90
	for point in points:
		west = min(west, point.longitude)
		south = min(south, point.latitude)
		east = max(east, point.longitude)
		north = max(north, point.latitude)
	return (south-buffer, west-buffer, north+buffer, east+buffer)

def gpx_point_to_coord(point, lat_first=False):
	if not lat_first:
		return (point.longitude, point.latitude)
	return (point.latitude, point.longitude)

if __name__ == '__main__':
	gpx_points=[]
	osm_coords = []
	out_points = []
	if len(sys.argv) < 4:
		usage()
		sys.exit(-1)
	threshold = int(sys.argv[3])
	with open(sys.argv[1], 'r') as gpxfile:
		gpx = gpxpy.parse(gpxfile)
		for track in gpx.tracks:
			print('processing track...')
			for seg in track.segments:
				print('processing segment...')
				for point in seg.points:
					gpx_points.append(point)
	bbox = get_bbox(gpx_points)
	overpass_query = 'way["highway"]({}, {}, {}, {})'.format(*bbox)
	api = overpass.API()
	print('getting OSM data from Overpass...')
	ways = api.Get(overpass_query).features
	if not ways:
		print('no OSM ways in vicinity')
		sys.exit(0)
	for way in ways:
		osm_coords.extend(way.geometry.coordinates)
	print('calculating distances...')
	progress = 0
	for point in gpx_points:
		min_distance = 10000
		for osm_coordinate in osm_coords:
			min_distance = min(haversine(gpx_point_to_coord(point), osm_coordinate), min_distance)
		if min_distance * 1000 > threshold:
			out_points.append({'point': point, 'distance': min_distance})
		progress = progress + 1
		sys.stdout.write("progress: {}/{}   \r".format(progress, len(gpx_points)))
		sys.stdout.flush()
	print('creating output GPX...')
	gpx = gpxpy.gpx.GPX()
	gpx_track = gpxpy.gpx.GPXTrack()
	gpx.tracks.append(gpx_track)
	gpx_segment = gpxpy.gpx.GPXTrackSegment()
	gpx_track.segments.append(gpx_segment)
	for point in out_points:
		out_point = point.get('point')
		out_point.description = str(point.get('distance'))
		gpx_segment.points.append(out_point)
	with open(sys.argv[2], 'w') as out_file:
		out_file.write(gpx.to_xml())
	print('done')