bus_router.py

import csv, os, subprocess, json, urllib, re, gpolyencode
from pprint import pprint
from geojson import LineString, Point, FeatureCollection, Feature, GeometryCollection
import geojson
import sys
import argparse

class AutoVivification(dict):
	"""Implementation of perl's autovivification feature."""
	def __getitem__(self, item):
		try:
			return dict.__getitem__(self, item)
		except KeyError:
			value = self[item] = type(self)()
			return value		


def decode(point_str):
	'''Decodes a polyline that has been encoded using Google's algorithm
	http://code.google.com/apis/maps/documentation/polylinealgorithm.html
	
	This is a generic method that returns a list of (latitude, longitude) 
	tuples.
	
	:param point_str: Encoded polyline string.
	:type point_str: string
	:returns: List of 2-tuples where each tuple is (latitude, longitude)
	:rtype: list
	
	'''
			
	# sone coordinate offset is represented by 4 to 5 binary chunks
	coord_chunks = [[]]
	for char in point_str:
		
		# convert each character to decimal from ascii
		value = ord(char) - 63
		
		# values that have a chunk following have an extra 1 on the left
		split_after = not (value & 0x20)		 
		value &= 0x1F
		
		coord_chunks[-1].append(value)
		
		if split_after:
				coord_chunks.append([])
		
	del coord_chunks[-1]
	
	coords = []
	
	for coord_chunk in coord_chunks:
		coord = 0
		
		for i, chunk in enumerate(coord_chunk):					
			coord |= chunk << (i * 5) 
		
		#there is a 1 on the right if the coord is negative
		if coord & 0x1:
			coord = ~coord #invert
		coord >>= 1
		coord /= 100000.0
					
		coords.append(coord)
	
	# convert the 1 dimensional list to a 2 dimensional list and offsets to 
	# actual values
	points = []
	prev_x = 0
	prev_y = 0
	for i in xrange(0, len(coords) - 1, 2):
		if coords[i] == 0 and coords[i + 1] == 0:
			continue
		
		prev_x += coords[i + 1]
		prev_y += coords[i]
		# a round to 6 digits ensures that the floats are the same as when 
		# they were encoded
		points.append((round(prev_x, 6), round(prev_y, 6)))
	
	return points

def processPolylines():
	encoder = gpolyencode.GPolyEncoder()
	json_data=open('data.txt')
	datadir = os.path.join(os.getcwd(), 'data')
	gtfsdir = os.path.join(datadir, 'gtfs')
	geojsondir = os.path.join(datadir, 'geojson')
	polydir = os.path.join(datadir, 'polylines')
	data = json.load(json_data, object_hook=_decode_dict)
	# pprint(data)
	json_data.close()
	with open(gtfsdir + "/shapes.txt", 'wb') as shapesfile:
		shapeswriter = csv.writer(shapesfile)
		shapeswriter.writerow(["shape_id","shape_pt_sequence","shape_dist_traveled","shape_pt_lon","shape_pt_lat"])
		for trip, stops in data.items():
			print trip
			
			legpoints = []
			jsonpoints = []

			for i in range(20):
				filepath = os.path.join(polydir, trip + "_" + str(i) + ".json")
				if(os.path.exists(filepath)):
					gmaps=open(filepath)
					
					linedata = json.load(gmaps)
					print trip + "_" + str(i)
					if args.dir == 'goog':
						if linedata['status'] != "OK":
							continue
						print linedata['routes'][0]['overview_polyline']['points']
						points = decode(linedata['routes'][0]['overview_polyline']['points'])
						# print points
						for point in points:
							dictpoint = {'x': point[0], 'y': point[1]}
							
							legpoints.append(dictpoint)
						gmaps.close()
					elif args.dir == 'osrm':
						if linedata['status'] != 0:
							continue
						print linedata['route_geometry']
						points = decode(linedata['route_geometry'])
						# print points
						for point in points:
							dictpoint = {'x': point[0]/10, 'y': point[1]/10}
							
							legpoints.append(dictpoint)

						gmaps.close()

			# print legpoints
			
			# print ls.geojson
			if not legpoints:
				continue
			else:
				simplified = simplify(legpoints, .0002, True)
				# print "new:" + str(simplified)
				count = 0
				for point in simplified:
					jsonpoints.append((point['x'], point['y']))
					shppoint = [point['x'], point['y']]
					shppoint.insert(0, trip)
					shppoint.insert(1, count)
					shppoint.insert(2, "")
					shapeswriter.writerow(shppoint)
					count += 1
				ls = LineString(jsonpoints)
				gc = GeometryCollection([ls])

				gtfsfile = os.path.join(geojsondir, trip + '.geojson')

				with open(gtfsfile, 'wb') as tripgeo:
					geojson.dump(gc, tripgeo)

def modifyTrips():
	datadir = os.path.join(os.getcwd(), 'data')
	gtfsdir = os.path.join(datadir, 'gtfs')
	
	with open(gtfsdir + "/trips.txt", 'r+') as tripsfile:
		tripsreader = csv.DictReader(tripsfile)
		keys = tripsreader.fieldnames
		with open(gtfsdir + "/trips_new.txt", 'wb') as tripsnew:
			tripswriter = csv.DictWriter(tripsnew, keys)
			tripswriter.writeheader()
			for row in tripsreader:
				
				# print row
				newtrip = row['trip_headsign'].replace(" ","").replace("/","") + "_" + row['route_id']
				row['shape_id'] = newtrip
				# print row
				#  remove that pesky last comma
				# print row[0]
				tripswriter.writerow(row)
				# tripsnew.write(string)
				

def _decode_list(data):
	rv = []
	for item in data:
		if isinstance(item, unicode):
			item = item.encode('utf-8')
		elif isinstance(item, list):
			item = _decode_list(item)
		elif isinstance(item, dict):
			item = _decode_dict(item)
		rv.append(item)
	return rv

def _decode_dict(data):
	rv = {}
	for key, value in data.iteritems():
		if isinstance(key, unicode):
			key = key.encode('utf-8')
		if isinstance(value, unicode):
			value = value.encode('utf-8')
		elif isinstance(value, list):
			value = _decode_list(value)
		elif isinstance(value, dict):
			value = _decode_dict(value)
		rv[key] = value
	return rv

def getDirections():
	json_data=open('data.txt')
	with open("env.json", 'rb') as envfile:
		env = json.load(envfile)
		google_key = env["google_key"]
	datadir = os.path.join(os.getcwd(), 'data')
	polydir = os.path.join(datadir, 'polylines')
	data = json.load(json_data, object_hook=_decode_dict)
	# pprint(data)
	json_data.close()
	
	for trip, stops in data.items():
		with open("log.txt", 'a') as log:
			log.write(trip + '\n')
		print trip
		# print stops
		stopcount = 1
		segmentcount = 0
		origin = ""
		dest = ""
		points = ""
		previous = {}
		waypoints = ""
		osrmpoints = []
		last = stops['stops'][-1]
		
		for index, stop  in enumerate(stops['stops']):
			fname = polydir + "/" + trip + "_" + str(segmentcount) + '.json'
			lastCheck = cmp(stop, last)
			# print lastCheck
			if stopcount == 1:
				# print "first stop"
				origin = stop['lat'] + "," + stop['lon']

				if segmentcount > 0:
					origin = stops['stops'][index - 1]['lat'] + "," + stops['stops'][index - 1]['lon']
					waypoints += stop['lat'] + "," + stop['lon'] + "|"
					osrmpoints.append(stop['lat'] + "," + stop['lon'])

				if lastCheck == 0:
					dest = stop['lat'] + "," + stop['lon']
					if args.dir == 'goog':
						directionscall(google_key, stop, origin, dest, waypoints, fname)
					elif args.dir == 'osrm':
						osrmDirectionsCall(stop, origin, dest, osrmpoints, fname)
					waypoints = ""
					osrmpoints = []
					segmentcount += 1
					continue

				stopcount += 1

			elif stopcount == 9 or lastCheck == 0:
				dest = stop['lat'] + "," + stop['lon'] 
				if args.dir == 'goog':
					directionscall(google_key, stop, origin, dest, waypoints, fname)
				elif args.dir == 'osrm':
					osrmDirectionsCall(stop, origin, dest, osrmpoints, fname)
				stopcount = 1
				waypoints = ""
				osrmpoints = []
				segmentcount += 1
				continue
			else:
				waypoints += stop['lat'] + "," + stop['lon'] + "|"
				osrmpoints.append(stop['lat'] + "," + stop['lon'])
				stopcount += 1


def directionscall(google_key, stop, origin, dest, waypoints, fname):
	print "getting dirs..."
	base = 'https://maps.googleapis.com/maps/api/directions/json?'
	params = urllib.urlencode({'origin': origin, 'destination': dest, 'waypoints': waypoints, 'sensor': 'false','key': google_key})
	# print params
	# if waypoints == "":
	with open("log.txt", 'a') as log:
		log.write(base + params + '\n')
	response = urllib.urlopen(base + params)
	data = json.load(response)
	with open(fname, 'w') as outfile:
		json.dump(data, outfile)

def osrmDirectionsCall(stop, origin, dest, osrmpoints, fname):
	print "getting dirs..."
	base = 'http://router.project-osrm.org/viaroute?'
	viastring = ""
	for point in osrmpoints:
		viastring += 'loc=' + point + '&'

	params = 'loc=' + origin + '&' + viastring + 'loc=' + dest
	# params = urllib.urlencode({'loc': origin, 'loc': dest, 'waypoints': waypoints, 'sensor': 'false','key': google_key})
	print params
	# if waypoints == "":
	with open("log.txt", 'a') as log:
		log.write(base + params + '\n')
	response = urllib.urlopen(base + params)
	data = json.load(response)
	with open(fname, 'w') as outfile:
		json.dump(data, outfile)

def shapesToGeojson():
	json_data=open('data.txt')
	datadir = os.path.join(os.getcwd(), 'data')
	gtfsdir = os.path.join(datadir, 'gtfs')
	geojsondir = os.path.join(datadir, 'geojson')
	data = json.load(json_data, object_hook=_decode_dict)
	json_data.close()
	with open(gtfsdir + "/shapes.txt", 'rb') as shapesfile:
		shapesreader = csv.DictReader(shapesfile)
		keys = shapesreader.fieldnames

		jsonpoints = []
		features = []
		currentTrip = ''
		for i, point in enumerate(shapesreader):
			if point['shape_pt_sequence'] == '0':
				print 'creating trip'
				currentTrip = point['shape_id']
				if i > 0:
					ls = LineString(jsonpoints)
					feature = Feature(geometry=ls, properties={"shape_id": currentTrip})
					# print feature
					features.append(feature)
					jsonpoints = []
			else:
				pnt = (float(point['shape_pt_lon']), float(point['shape_pt_lat']))
				# print pnt
				jsonpoints.append(pnt)

		# write linestring for last shape
		ls = LineString(jsonpoints)
		feature = Feature(geometry=ls, properties={"shape_id": currentTrip})
		print feature
		features.append(feature)
		jsonpoints = []		
		fc = FeatureCollection(features)
		print fc

		geojsonfile = os.path.join(geojsondir, 'shapes.geojson')

		with open(geojsonfile, 'wb') as tripgeo:
			geojson.dump(fc, tripgeo)


def geojsonToShapes():
	datadir = os.path.join(os.getcwd(), 'data')
	gtfsdir = os.path.join(datadir, 'gtfs')
	geojsondir = os.path.join(datadir, 'geojson')

	with open(gtfsdir + "/shapes_new.txt", 'wb') as shapesfile:
		shapeswriter = csv.writer(shapesfile)
		shapeswriter.writerow(["shape_id","shape_pt_sequence","shape_dist_traveled","shape_pt_lon","shape_pt_lat"])

		geojsonfile = os.path.join(geojsondir, 'shapes.geojson')

		with open(geojsonfile, 'rb') as fc:
			geo_fc = geojson.load(fc)
			# print geo_fc
			for feature in geo_fc['features']:
				for i, coord in enumerate(feature['geometry']['coordinates']):
					shapeswriter.writerow([feature['properties']['shape_id'],i,'',coord[0],coord[1]])


def processGtfs():
	datadir = os.path.join(os.getcwd(), 'data/gtfs')
	stopspath = os.path.join(datadir, 'stops.txt')
	timespath = os.path.join(datadir, 'stop_times.txt')
	tripspath = os.path.join(datadir, 'trips.txt')

	trips = AutoVivification()
	names = {}
	with open(tripspath, 'rb') as tripsfile:
		tripsreader = csv.DictReader(tripsfile)
		for trip in tripsreader:
			
			tripname = trip['trip_headsign'].replace(" ","").replace("/","") + "_" + trip['route_id']
			if tripname in names.values():
				continue
			else:
				names[trip['trip_id']] = tripname
				# print tripname
	with open('names.json', 'wb') as namesfile:
		json.dump(names, namesfile)

	# read times into dict
	with open(timespath, 'rb') as timesfile:
		timesreader = csv.DictReader(timesfile)
		tripscount = 0
		stopscount = 0
		curr_trip = ""
		for time in timesreader:
			# print time['trip_id']
			# print names[time['trip_id']]
			
			tripname = ""
			if time['trip_id'] in names.keys():
				tripname = names[time['trip_id']]
			else:
				continue

			if not (tripname in trips):
				trips[tripname]['stops'] = []
				nextStop = {'id': time["stop_id"]}
				trips[tripname]['stops'].append(nextStop)
				curr_trip = tripname
				tripscount += 1

			elif tripname == curr_trip:
				nextStop = {'id': time["stop_id"]}
				trips[tripname]['stops'].append(nextStop)
				stopscount += 1
				curr_trip = tripname
			else:
				continue
	# print trips
	# read stops into dict
	with open(stopspath, 'rb') as stopsfile:
		stopsreader = csv.DictReader(stopsfile)
		for stop in stopsreader:
			# print stop['stop_id']
			for trip in trips:
				print trip
				for i in trips[trip]['stops']:
					# print i
					# print trips[trip]['stops']
					if stop['stop_id'] == i['id']:
						i['lat'] = stop['stop_lat']
						i['lon'] = stop['stop_lon']
	with open('data.txt', 'w') as outfile:
		json.dump(trips, outfile)

	# read latlng for stop_ids from stops.txt
	
	# urlib for https://maps.googleapis.com/maps/api/directions/json?origin=33.983518,-84.084661&destination=33.762069,-84.384326&waypoints=33.765978,-84.387597&sensor=false&key=AIzaSyD2KTHZHT8Bl-JzgF3yI1t7Ln05udSu318

	# get json back

	# decode routes[0].legs[0].steps.0.polyline
	# iterate on legs and steps...
	# add points to array (each with their own number)
	# 




def getSquareDistance(p1, p2):
	"""
	Square distance between two points
	"""
	dx = p1['x'] - p2['x']
	dy = p1['y'] - p2['y']

	return dx * dx + dy * dy


def getSquareSegmentDistance(p, p1, p2):
	"""
	Square distance between point and a segment
	"""
	x = p1['x']
	y = p1['y']

	dx = p2['x'] - x
	dy = p2['y'] - y

	if dx != 0 or dy != 0:
		t = ((p['x'] - x) * dx + (p['y'] - y) * dy) / (dx * dx + dy * dy)

		if t > 1:
			x = p2['x']
			y = p2['y']
		elif t > 0:
			x += dx * t
			y += dy * t

	dx = p['x'] - x
	dy = p['y'] - y

	return dx * dx + dy * dy


def simplifyRadialDistance(points, tolerance):
	length = len(points)
	prev_point = points[0]
	new_points = [prev_point]

	for i in range(length):
		point = points[i]

		if getSquareDistance(point, prev_point) > tolerance:
			new_points.append(point)
			prev_point = point

	if prev_point != point:
		new_points.append(point)

	return new_points


def simplifyDouglasPeucker(points, tolerance):
	length = len(points)
	markers = [0] * length  # Maybe not the most efficent way?

	first = 0
	last = length - 1

	first_stack = []
	last_stack = []

	new_points = []

	markers[first] = 1
	markers[last] = 1

	while last:
		max_sqdist = 0

		for i in range(first, last):
			sqdist = getSquareSegmentDistance(points[i], points[first], points[last])

			if sqdist > max_sqdist:
				index = i
				max_sqdist = sqdist

		if max_sqdist > tolerance:
			markers[index] = 1

			first_stack.append(first)
			last_stack.append(index)

			first_stack.append(index)
			last_stack.append(last)

		# Can pop an empty array in Javascript, but not Python, so check
		# the length of the list first
		if len(first_stack) == 0:
			first = None
		else:
			first = first_stack.pop()

		if len(last_stack) == 0:
			last = None
		else:
			last = last_stack.pop()

	for i in range(length):
		if markers[i]:
			new_points.append(points[i])

	return new_points


def simplify(points, tolerance=0.1, highestQuality=True):
	sqtolerance = tolerance * tolerance

	if not highestQuality:
		points = simplifyRadialDistance(points, sqtolerance)

	points = simplifyDouglasPeucker(points, sqtolerance)

	return points

if __name__ == '__main__':
	parser = argparse.ArgumentParser()
	parser.add_argument("-d", "--dir", type=str, metavar='osrm', help='specify the directions provider, either "goog" or "osrm"')
	parser.add_argument("-s", "--shapes", action='store_true', help='create shapes.txt from GeoJSON')
	parser.add_argument("-l", "--lines", action='store_true', help='process polylines if directions calls have already been made')
	parser.add_argument("-t", "--trips", action='store_true', help='modify trips file with new shape_ids (creates new file)')
	parser.add_argument("-g", "--geojson", action='store_true', help='create GeoJSON from shapes.txt')
	args = parser.parse_args()

	if args.dir:
		print "directions turned on using " + args.dir
	else:
		args.dir = 'osrm'
	if args.shapes:
		print "converting GeoJSON to shapes.txt"
		geojsonToShapes()
		sys.exit()
	if args.geojson:
		print "converting shapes.txt to GeoJSON"
		shapesToGeojson()
		sys.exit()
	if args.lines:
		print "processing polylines (no directions call or GTFS processing) using polylines from " + args.dir
		processPolylines()
		sys.exit()
	if args.trips:
		print 'modify trips file with new shape_ids (creates new file)'
		modifyTrips()
		sys.exit()

	processGtfs()
	getDirections()
	processPolylines()
	modifyTrips()