/calculate_path.py

import sqlalchemy as sa
import psycopg2
import pandas.io.sql as sqlio
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import folium
import json
import requests
from IPython import embed
import datetime
from folium import plugins
from passwords import postgres_password

plt.style.use('bmh')

MAP_CONN = psycopg2.connect(host='localhost', dbname='nyc_routing', user='postgres', password=postgres_password)

day_names = ['sun', 'mon', 'tues', 'wed', 'thurs', 'fri', 'sat']
DEGREES_PER_KM = 1 / 111.

def get_gps_from_address(address):
	response = requests.get('https://maps.googleapis.com/maps/api/geocode/json?address=' + str(address))
	resp_json_payload = response.json()
	lat = resp_json_payload['results'][0]['geometry']['location']['lat']
	lon = resp_json_payload['results'][0]['geometry']['location']['lng']

	return lon, lat

def calc_busiest_path_in_radius(source_id, target_id, source_lon, source_lat, r, day, hour, weight):
	"""Input source node id, target node id, source gps, day of week, hour
	of day returns busiest path within r kms between source and target.
	"""

	r_degrees = r * DEGREES_PER_KM
	sql = ("DROP TABLE IF EXISTS d_path; DROP TABLE IF EXISTS d_path2;" +
		"CREATE TABLE d_path AS SELECT seq, id1 AS node, id2 AS edge, cost FROM pgr_dijkstra('SELECT ways.gid as id, source::int4, target::int4, ways.cost / %(weight)s / (taxi_pickup_counts.green_%(day)s_%(hour)s + taxi_pickup_counts.yellow_%(day)s_%(hour)s)::float8 AS cost, ways.reverse_cost / %(weight)s / (taxi_pickup_counts.green_%(day)s_%(hour)s + taxi_pickup_counts.yellow_%(day)s_%(hour)s) AS reverse_cost FROM ways INNER JOIN taxi_pickup_counts ON ways.gid = taxi_pickup_counts.gid  WHERE ((ways.x1 - %(source_lon)s)^2 + (ways.y1 - %(source_lat)s)^2 ) < %(r_degrees)s^2', %(source_id)s, %(target_id)s, True, True);" +
	"CREATE TABLE d_path2 AS SELECT node, lon, lat, edge FROM ways_vertices_pgr, d_path where ways_vertices_pgr.id = d_path.node;" + 
	"SELECT node, lon, lat, edge, ST_AsGeoJSON(the_geom) as geojson FROM d_path2, ways_costs WHERE ways_costs.gid=d_path2.edge;")

	df = sqlio.read_sql(sql % {'source_id':source_id, 'target_id': target_id, 
		'source_lon': source_lon, 'source_lat': source_lat, 'r_degrees':r_degrees, 
		'day':day_names[day], 'hour':hour, 'weight': weight}, MAP_CONN)

	df.geojson = [json.loads(gj.replace("'", "\"")) for gj in df.geojson]
	df.geojson = [gj for gj in df.geojson]

	return df

def calc_astar_path(source_id, target_id, source_lon, source_lat, r, day, hour, weight):
	sql = "SELECT seq, id1 AS node, id2 AS edge, cost FROM pgr_AStar('SELECT ways.gid::INTEGER as id, source::INTEGER, target::INTEGER, ways.cost / %(weight)s / (taxi_pickup_counts.green_%(day)s_%(hour)s + taxi_pickup_counts.yellow_%(day)s_%(hour)s)::float8 AS cost, x1, y1, x2, y2, ways.reverse_cost / %(weight)s / (taxi_pickup_counts.green_%(day)s_%(hour)s + taxi_pickup_counts.yellow_%(day)s_%(hour)s) AS reverse_cost FROM ways INNER JOIN taxi_pickup_counts ON ways.gid = taxi_pickup_counts.gid  WHERE ((ways.x1 - %(source_lon)s)^2 + (ways.y1 - %(source_lat)s)^2 ) < %(r_degrees)s^2', %(source_id)s, %(target_id)s, True, True);"
	
	r_degrees = r * DEGREES_PER_KM
	df = sqlio.read_sql(sql % {'source_id':source_id, 'target_id': target_id, 
		'source_lon': source_lon, 'source_lat': source_lat, 'r_degrees':r_degrees, 
		'day':day_names[day], 'hour':hour, 'weight': weight}, MAP_CONN)

	df['geojson'] = pd.Series([json.loads(gj.replace("'", "\"")) for gj in df['geojson']], index=df.index)
	df['geojson'] = pd.Series([gj for gj in df.geojson], index=df.index)

	return path

def get_pickups_in_radius(lon, lat, r, hour, day):
	r_degrees = r * DEGREES_PER_KM
	sql = "SELECT pickup_longitude, pickup_latitude FROM yellow_trip_data WHERE pickup_hour=%(pickup_hour)s AND day_of_week=%(day_of_week)s AND ((pickup_longitude - %(lon)s)^2 + (pickup_latitude - %(lat)s)^2 ) < %(r_degrees)s^2;"
	df = sqlio.read_sql(sql % {'pickup_hour':hour, 'day_of_week':1, 'lon':lon, 'lat':lat, 'r_degrees':r_degrees}, MAP_CONN)
	df['latlon'] = zip(df['pickup_latitude'], df['pickup_longitude'])
	#df['geometry'] = df[['latlon']].applymap(lambda x:geometry.Point(x))

	return df

def get_ways_pickups_in_radius(lon, lat, r, hour, day):
	"""
	Takes lon, lat and returns all ways and pickup counts for 
	given hour of the day and day of the week.
	"""

	r_degrees = r * DEGREES_PER_KM
	sql = "SELECT ways.*, taxi_pickup_counts.%(day)s_%(hour)s_counts  FROM ways INNER JOIN taxi_pickup_counts ON ways.gid=taxi_pickup_counts.gid  WHERE ((ways.x1 - %(lon)s)^2 + (ways.y1 - %(lat)s)^2 ) < %(r_degrees)s^2;"
	df = sqlio.read_sql(sql % {'lon': lon, 'lat': lat, 'r_degrees': r_degrees, 'day': day_names[day], 'hour':hour}, MAP_CONN)

	return df

def get_busiest_node_id_in_radius(lon, lat, r, hour, day):
	
	r_degrees = r * DEGREES_PER_KM

	sql = ("SELECT ways.* FROM " + 
		"ways INNER JOIN taxi_pickup_counts ON ways.gid=taxi_pickup_counts.gid " + 
		"WHERE taxi_pickup_counts.yellow_%(day)s_%(hour)s + taxi_pickup_counts.green_%(day)s_%(hour)s = "  +
		"( SELECT max(taxi_pickup_counts.yellow_%(day)s_%(hour)s+taxi_pickup_counts.green_%(day)s_%(hour)s) FROM ways INNER JOIN taxi_pickup_counts ON ways.gid=taxi_pickup_counts.gid  WHERE ((ways.x1 - %(lon)s)^2 + (ways.y1 - %(lat)s)^2 ) < %(r_degrees)s^2 );")

	busiest_ways = sqlio.read_sql(sql % {'lon': lon, 'lat': lat, 'r_degrees': r_degrees, 'day': day_names[day], 'hour':hour}, MAP_CONN)
	d2 = ((busiest_ways.x1 - lon)**2 + (busiest_ways.y1 - lat)**2)
	#print(busiest_ways[d2 == d2.min()].source)

	return busiest_ways[d2 == d2.min()]['source'].iloc[0]

def get_nearest_node(lon, lat):
	"""Input lon, lat returns nearest osm_node"""
	sql = "SELECT * FROM ways_vertices_pgr ORDER BY the_geom <-> ST_GeometryFromText('POINT(%(lon)s %(lat)s )',4326) LIMIT 1;"

	df = sqlio.read_sql(sql % {'lon':lon, 'lat': lat}, MAP_CONN)

	return df['id'][0], df

def get_way(gid):
	sql = "SELECT * FROM ways WHERE gid=%(gid)s;"
	df = sqlio.read_sql(sql % {'gid':gid}, MAP_CONN)

	return df

def get_node(id):
	sql = "SELECT * FROM ways_vertices_pgr WHERE id=%(id)s;"
	df = sqlio.read_sql(sql % {'id':id}, MAP_CONN)

	return df

def get_nearest_way(lon, lat):
	"""Input lon, lat returns nearest osm_way"""
	sql = "SELECT * FROM ways ORDER BY the_geom <-> ST_GeometryFromText('POINT(%(lon)s %(lat)s )',4326) LIMIT 1;"

	df = sqlio.read_sql(sql % {'lon':lon, 'lat': lat}, MAP_CONN)

	return df['gid'][0], df

def generate_map(source_lon, source_lat, target_lon, target_lat, df_path, df_pickups):
	nyc_map = folium.Map(location=[source_lat, source_lon], zoom_start=15)
	grad={0: 'yellow', 0.5: 'blue', 1: 'green'}
	#pickups_map = plugins.HeatMap(df_pickups.latlon, gradient=grad, min_opacity=0.5, radius=5, blur=10)
	#nyc_map.add_children(pickups_map)
	[folium.PolyLine(np.array(gj['coordinates'])[:,[1,0]]).add_to(nyc_map) for gj in df_path.geojson]
	dt_string = str(datetime.datetime.today()).replace(' ', '_') 
	route_template = './taxi_path_app/templates/tmp/route_' + dt_string + '.html'
	nyc_map.simple_marker([source_lat, source_lon], marker_color='green')
	nyc_map.simple_marker([target_lat, target_lon], marker_color='red')
	nyc_map.save(route_template)

	return route_template

def main(source_lon, source_lat, day, hour, weight, target_lon=None, target_lat=None, r=2.0):
	source_id, df_source = get_nearest_node(source_lon, source_lat)
	#target_id, df_target = get_nearest_node(target_lon, target_lat)
	target_id = get_busiest_node_id_in_radius(source_lon, source_lat, r, hour, day)
	target_df = get_node(target_id)
	target_lon = target_df['lon'][0]; target_lat = target_df['lat'][0]
	df_path = calc_busiest_path_in_radius(source_id, target_id, source_lon, source_lat, r, day, hour, weight)
	#df_pickups = get_pickups_in_radius(source_lon, source_lat, r, hour, day)
	df_pickups = []
	route_template = generate_map(source_lon, source_lat, target_lon, target_lat, df_path, df_pickups)

	return route_template
	
if __name__ == '__main__':
	lon_insight = -73.9861860
	lat_insight = 40.7419650
	lon_home = -73.9576270
	lat_home = 40.7607880
	lon_world_trade = -74.0131
	lat_world_trade = 40.7130
	lon_bedford = -73.957504
	lat_bedford = 40.717883

	hour = 9
	day = 1
	r = 1.0
	source_id = 44670
	target_id = 29971
	weight = 1
	#df_path = calc_shortest_path(source_id, target_id, day=day, hour=hour, weight=weight)
	source_lon = -73.9861860
	source_lat = 40.7419650
	#df_pickups=get_pickups_in_radius(lon, lat, r, hour, day)
	main(source_lon, source_lat, day, hour, weight, target_lon=None, target_lat=None, r=2.0)