In [1]:
import collections
import functools
from imposm.parser import OSMParser
import itertools
from IPython.display import HTML
import json
from matplotlib import collections as mc
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.cm as cmx
from numpy import nan
import numpy as np
import pandas as pd
import pyproj
import requests
import scipy as sp
import rtree
import seaborn as sb
from scipy import signal
# import shapely
import shapely.geometry
%pylab inline
import data_munging
In [14]:
HTML('<iframe src="http://chihacknight.org/" width=800 height=400></iframe>')
Out[14]:
In [17]:
HTML('<iframe src="http://project-pavement.herokuapp.com/" width=800 height=400></iframe>')
Out[17]:
In [55]:
rides, readings = data_munging.read_raw_data()
raw_reading = readings.loc[0, :]
# This is here temporarily both to create a static dataset and because
# some of the newer android readings are bogus and need to be removed.
rides = rides.loc[rides['created_at'] < '2016-03-01', :]
readings = readings.loc[readings['created_at'] < '2016-03-01', :]
n, p = readings.shape
In [11]:
raw_reading
Out[11]:
In [57]:
# chi_readings = data_munging.filter_readings_to_chicago(readings)
# chi_rides = list(set(chi_readings.ride_id))
In [41]:
readings = readings.loc[np.logical_and(readings.gps_mph > 5, readings.gps_mph < 30), :]
In [182]:
# Visualization of Accelerometer Time Series
for i in random_indices:
colors.append(np.random.random(3))
sb.tsplot(readings['num_accel_' + axis][i][0:100], alpha=0.8, color=colors[-1])
fig = plt.gcf()
fig.set_size_inches(18.5, 10.5)
plt.xlabel('Accelerometer Values')
plt.ylabel('Force (Gravities)')
plt.title('Random sample of ' + str(sample_size) + ' ' + axis + ' Accelerometer Time Series')
plt.savefig('../dat/plots_for_presentation/z_accel.png')
In [56]:
readings = data_munging.clean_readings(readings)
readings = data_munging.add_proj_to_readings(readings, data_munging.NAD83)
In [177]:
axis = 'z'
sample_size = 10
random_indices = np.random.choice(readings.index, sample_size)
colors = []
I couldn't find anything practical to do with these.
In [42]:
readings_idx = data_munging.insert_readings_rtree(readings)
In [ ]:
total_intersections = []
In [175]:
for sample_in in range(100):
random_idx, random_point = data_munging.select_random_point(readings)
random_bb = data_munging.point_to_bb(*random_point, side_length=0.2)
intersecting_segments = list(readings_idx.intersection(random_bb))
intersecting_segments = [i for i in intersecting_segments if data_munging.calc_reading_diffs(readings.ix[i, :], readings.ix[random_idx, :]) < 0.4]
if len(intersecting_segments) > 1:
for i in intersecting_segments:
plt.plot([readings.ix[i, 'start_x'], readings.ix[i, 'end_x']], [readings.ix[i, 'start_y'], readings.ix[i, 'end_y']], alpha=0.5)
fig = plt.gcf()
fig.set_size_inches(18.5, 10.5)
plt.title('Intersections of ' + str(intersecting_segments))
plt.show()
In [44]:
# Use a threshold for 0.4
def calc_reading_diffs(reading0, reading1):
start0 = reading0[['start_x', 'start_y']].values
start1 = reading1[['start_x', 'start_y']].values
end0 = reading0[['end_x', 'end_y']].values
end1 = reading1[['end_x', 'end_y']].values
diff0 = np.linalg.norm(start0 - start1) + np.linalg.norm(end0 - end1)
diff1 = np.linalg.norm(start0 - end1) + np.linalg.norm(end0 - start1)
diff = min(diff0, diff1)
dist0 = np.linalg.norm(start0 - end0)
dist1 = np.linalg.norm(start1 - end1)
if dist0 == 0 or dist1 == 0:
return np.inf
return diff / (dist0 + dist1)
In [46]:
pairwise_comps = pd.DataFrame({'id0': ids0, 'id1': ids1})
for i in ('0', '1'):
for var in ('std_z', 'gps_speed', 'abs_mean_total'):
pairwise_comps[var + i] = readings[var][pairwise_comps['id' + i]].values
for i in ('0', '1'):
pairwise_comps['abs_mean_over_speed' + i] = pairwise_comps['abs_mean_total' + i] / pairwise_comps['gps_speed' + i]
In [190]:
pairwise_comps.plot(x='abs_mean_total0', y='abs_mean_total1', kind='scatter', alpha=0.5, title='Comparing Mean Absolute Acceleration')
fig = plt.gcf()
fig.set_size_inches(12.5, 12.5)
plt.show()
In [191]:
pairwise_comps.plot(x='abs_mean_over_speed0', y='abs_mean_over_speed1', kind='scatter', alpha=0.5, title='Comparing Absolute Acceleration Over Speed')
fig = plt.gcf()
fig.set_size_inches(12.5, 12.5)
plt.show()
fig.savefig('../dat/plots_for_presentation/mean_abs_accel_over_speed.png')
In [73]:
clean_chi_readings = pd.read_csv(data_munging.data_dir + 'clean_chi_readings.csv')
In [66]:
ride_ids = set(clean_chi_readings.ride_id)
print(ride_ids)
In [60]:
ride_id = 35
In [68]:
for ride_id in ride_ids:
try:
ax = clean_chi_readings.loc[readings['ride_id'] == ride_id, :].plot(x='start_lon', y='start_lat', title='Snapping GPS Trails Using OSRM')
clean_chi_readings.loc[readings['ride_id'] == ride_id, :].plot(x='snapped_lon', y='snapped_lat', ax=ax)
fig = plt.gcf()
fig.set_size_inches(18.5, 10.5)
plt.show()
except:
'oh well'
In [115]:
road_bumpiness = collections.defaultdict(list)
for index, reading in clean_chi_readings.iterrows():
if reading['gps_mph'] < 30 and reading['gps_mph'] > 5:
osm_segment = [(reading['snapped_lat'], reading['snapped_lon']),
(reading['next_snapped_lat'], reading['next_snapped_lon'])]
osm_segment = sorted(osm_segment)
if all([(not pd.isnull(lat_lon[0])) and (not pd.isnull(lat_lon[0])) for lat_lon in osm_segment]):
road_bumpiness[tuple(osm_segment)].append(reading['abs_mean_over_speed'])
In [117]:
total_road_readings = dict((osm_segment, len(road_bumpiness[osm_segment])) for osm_segment in road_bumpiness)
In [118]:
agg_road_bumpiness = dict((osm_segment, np.mean(road_bumpiness[osm_segment])) for osm_segment in road_bumpiness)
In [105]:
def find_seg_dist(lat_lon):
return data_munging.calc_dist(lat_lon[0][1], lat_lon[0][0], lat_lon[1][1], lat_lon[1][0])
In [135]:
dist_of_road = dict((osm_segment, find_seg_dist(osm_segment)) for osm_segment in road_bumpiness if find_seg_dist(osm_segment) < 200)
In [144]:
thickness_of_road = dict()
for osm_segment in dist_of_road:
thickness_of_road[osm_segment] = n(total_road_readings[osm_segment] / (dist_of_road[osm_segment] + 1)) ** 0.5
In [146]:
np.percentile(thickness_of_road.values(), np.arange(0, 100, 11))
Out[146]:
In [127]:
import matplotlib.colors as colors
plasma = cm = plt.get_cmap('plasma')
cNorm = colors.Normalize(vmin=0, vmax=1.0)
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=plasma)
In [187]:
plt.hist(total_road_readings.values(), bins= 40)
fig = plt.gcf()
fig.set_size_inches(12.5, 10.5)
In [158]:
std_road_bumpiness = dict((osm_segment, np.std(road_bumpiness[osm_segment])) for osm_segment in road_bumpiness)
In [192]:
for osm_segment in dist_of_road:
plt.plot([osm_segment[0][1], osm_segment[1][1]], [osm_segment[0][0], osm_segment[1][0]],
color=scalarMap.to_rgba(agg_road_bumpiness[osm_segment]), linewidth=2*thickness_of_road[osm_segment] + 4)
fig = plt.gcf()
fig.set_size_inches(12, 16)
plt.savefig('../dat/plots_for_presentation/final_agg_readings_per_osm_segment.png')
| Name | Photo | Role | Website |
|---|---|---|---|
| Michael Hassin | Rails and iOS | https://github.com/strangerloops | |
| Joshua Shin | Android | https://github.com/jjshin85 | |
| Nate Hutcheson | Group Dad | https://www.linkedin.com/in/nate-hutcheson-2946444 | |
| Stephen Vance | Data and Bike Policy Guru, Photos | http://www.stevevance.net/ |
| Name | Photo | Role | Website |
|---|---|---|---|
| Stephen Liu | UX and Project Management | http://www.iamstephenliu.com/ | |
| William Henderson | Backend (Ride Report) | https://twitter.com/quicklywilliam | |
| Evan Heidtmann | Backend (Ride Report) | https://twitter.com/evanheidtmann |
In [45]:
total_data_points = 0
total_tries = 0
examined_ids = set()
ids0 = list()
ids1 = list()
while total_data_points < 5000 and total_tries < 100000:
total_tries += 1
random_idx, random_point = data_munging.select_random_point(readings)
random_bb = data_munging.point_to_bb(*random_point, side_length=0.2)
intersecting_segments = list(readings_idx.intersection(random_bb))
intersecting_segments = set([i for i in intersecting_segments if data_munging.calc_reading_diffs(readings.ix[i, :], readings.ix[random_idx, :]) < 0.4])
intersecting_segments.difference_update(examined_ids)
examined_ids.update(intersecting_segments)
if intersecting_segments > 1:
for idx_pair in itertools.combinations(intersecting_segments, 2):
ids0.append(idx_pair[0])
ids1.append(idx_pair[1])
total_data_points += 1
suffices to look at OSM level since segments always take into account intersections. Negligible to consider differing pavement quality for different parts of a road segment for start and destination
In [7]:
for color, i in zip(colors, indices):
f, Pxx_den = signal.periodogram(readings['num_accel_' + axis][i][0:100])
plt.plot(f, Pxx_den, color=color, alpha=0.8)
plt.title('Power Spectrum for ' + axis + ' axis')
plt.xlabel('frequency [Hz]')
plt.ylabel('Power Spectrum Density')
fig = plt.gcf()
fig.set_size_inches(18.5, 10.5)
plt.show()
In [8]:
sample_size = 1000
indices = np.random.choice(n, sample_size)
for axis in ['z']:
for i in indices:
f, Pxx_den = signal.periodogram(readings['num_accel_' + axis][i][0:100])b
plt.plot(f, Pxx_den, alpha=0.5)
plt.title('Power Spectrum for ' + axis + ' axis')
plt.xlabel('frequency [Hz]')
plt.ylabel('Power Spectrum Density')
fig = plt.gcf()
fig.set_size_inches(18.5, 10.5)
plt.show()
Justification for Cycling
What is Chi Hack Night Ethos: empower cyclists and provide open data and tools to improve the urban cycling experience Collective of civically-minded hackers who meet at the Merch Mart
Obstacles to Cycling Congested roads No bicycle lanes Icy bicycle lanes Potholes/divots in road/cracked pavement Road closings
Pain Points Western and Ashland versus Damen Open Grate Bridges
Personalized Bicycle Routing Mark up OSM with relevant information Collect Surface Quality Data (Project Pavement) Collect Subjective Experience (Ride-Report and Stress Map) Integrate all of these data sources with OSRM and Leaflet-Routing-Machine
Project Pavement
Ride-Report
OSRM
Policy Implications Show bike equity Talk about how usability score could determine bicycle deserts
Future Directions
Ride Report Evan and William
Michael Hassin- iOS and rails app Steven Vance - data guy and bicycle outreach expert Nate Hutcheson - project management Joshua Shin - Android App (give him a job!)
Appendix and Cutting Room Floor