BestHomeLocationFinder.py

#!/usr/bin/python3

# BestHomeLocationFinder by Hans Broeckx <hans.broeckx@gmail.com>

# DEPENDANCIES: bokeh, numpy, googlemaps

# USAGE: Fill in some points of interest latitude and longitude coordinates, max driving time and max car speed. Run the script.
# A few more settings are available in the Settings section.

# WARNING: This script uses the Google Maps Distance Matrix API (free API keys get 2500 request / day).
# Depending on the accuracy and the number of starting POIs, your API key might hit its daily limit.
# Free API keys get 2500 request / day. You can monitor your usage at this URL:
# https://console.developers.google.com/apis/api/distance_matrix_backend/quotas

from bokeh.io import output_file, show
from bokeh.models import GMapPlot, GMapOptions, ColumnDataSource, Circle, Patch, DataRange1d, PanTool, WheelZoomTool, BoxSelectTool
import numpy
import googlemaps

#### Settings

pois = dict( # points of interest coordinates around which to focus
    work=(51.123,5.789),
    parents=(51.456,5.456),
    favorite_bar=(51.789,5.123),
)
max_drive_time = 30 #in minutes
max_car_speed = 120 #in km/h
accuracy = 3 # 0 = 4 points in all directions // 1 = 8 points // 2 = 12 points // ...
google_api_key = 'YOUR_SERVER_API_KEY' # your server API key with access to Google Maps Distance Matrix API
zoom = 10 # initial map zoom level
color = 'blue' # POI and region marker color

##### Prep

gmaps = googlemaps.Client(key=google_api_key)

pois_data = ColumnDataSource(data=dict(
    lat=[i[0] for i in pois.values()],
    lon=[i[1] for i in pois.values()],
    names=list(pois.keys())
))

avg_latlon = numpy.mean(pois_data.data['lat']),numpy.mean(pois_data.data['lon'])

print('Centering map around coordinate {}'.format(avg_latlon))
max_drive_time = max_drive_time*60 #in seconds
max_car_speed = max_car_speed * 1000 / 3600 #in meters/second
max_distance = max_drive_time * max_car_speed
print('At {}m/s, the max distance we can reach is {}m over {}s.'.format(max_car_speed,max_distance,max_drive_time))

map_options = GMapOptions(lat=avg_latlon[0], lng=avg_latlon[1], map_type='roadmap', zoom=zoom)
plot = GMapPlot(x_range=DataRange1d(),y_range=DataRange1d(), map_options=map_options, title='Best Home Locations Finder')

##### Helper functions
def deg2rad(deg):
    return deg * numpy.pi / 180
def rad2deg(rad):
    return rad * 180 / numpy.pi
def latLonAtDistanceAndBearingFromLatLon(d,b,latlon):
    ad = d/6378140 #angular distance d/R (d=distance traveled, R=earth's radius)
    lat1 = deg2rad(latlon[0])
    lon1 = deg2rad(latlon[1])
    lat2 = numpy.arcsin( numpy.sin(lat1)*numpy.cos(ad) + numpy.cos(lat1)*numpy.sin(ad)*numpy.cos(b) )
    lon2 = lon1 + numpy.arctan2( numpy.sin(b)*numpy.sin(ad)*numpy.cos(lat1), numpy.cos(ad)-numpy.sin(lat1)*numpy.sin(lat2) )
    return rad2deg(lat2),rad2deg(lon2)

##### Number crunching
range_points = dict()
# Go over each starter POI
for poi_name, poi_latlon in pois.items():
    n_parts = 8+4*(accuracy-1)
    range_points[poi_name] = []
    for part in range(0,n_parts):
        rad = 2 * numpy.pi * part / n_parts
        print('Now analysing {} part {} of {} at {} rad from {}'.format(poi_name,part+1,n_parts,rad,poi_latlon))
        drive_time = max_drive_time
        distance = max_distance
        while drive_time >= max_drive_time:
            latlon = latLonAtDistanceAndBearingFromLatLon(distance,rad,poi_latlon)
            distance_matrix = gmaps.distance_matrix(poi_latlon,latlon)
            if distance_matrix['status'] == 'OK':
                drive_time = distance_matrix['rows'][0]['elements'][0]['duration']['value']
                print('Target point {} is a {}s drive away from POI.'.format(latlon, drive_time))
                if drive_time >= max_drive_time:
                    print('Too far away, trying 10% closer')
                    distance -= max_distance*0.1
                else:
                    print('Close enough!')
                    range_points[poi_name].append(latlon)

    init_range_data = ColumnDataSource(data=dict(
        lat=[i[0] for i in range_points[poi_name]],
        lon=[i[1] for i in range_points[poi_name]]
    ))
    init_rangearea_marker = Patch(x='lon',y='lat',fill_color=color,fill_alpha=0.1,line_color=None)
    plot.add_glyph(init_range_data,init_rangearea_marker)
    # Add some tiny markers to show the points we used?
    #init_rangepoint_marker = Circle(x='lon',y='lat',size=3,fill_color=color,fill_alpha=1.0,line_color=None)
    #plot.add_glyph(init_range_data,init_rangepoint_marker)

print('Adding starter coordinates to the map')
pois_marker = Circle(x='lon',y='lat',size=10,fill_color=color,fill_alpha=1.0,line_color=None)
plot.add_glyph(pois_data,pois_marker)

##### Render

plot.add_tools(PanTool(), WheelZoomTool(), BoxSelectTool())
output_file("gmap_plot.html")
show(plot)