有很多工具可以用来分析这样的数据，但我选择的工具是（显然）Python。在这篇文章中，我想说明你如何开始分析这些数据，并使用PyData堆栈，即NumPy、Pandas、Matplotlib和Seaborn，将其与其他可用的数据源连接起来。在这里，我将看一下你可以用这些数据回答的一些基本问题。以后我希望能找到时间深入挖掘，提出一些更有趣、更有创意的问题--敬请关注
对于那些不熟悉的人来说，这篇文章是以Jupyter笔记本的形式组成的，它是一种开放的文档格式，结合了文本、代码、数据和图形，可以通过网络浏览器查看--如果你以前没有使用过它，我鼓励你尝试一下你可以在这里下载包含这篇文章的笔记本，用Jupyter打开它，然后开始对数据提出你自己的问题。
下载Pronto的数据
我们首先要下载数据（可在Pronto的网站上找到），你可以通过取消注释以下shell命令来完成（这里的感叹号是运行shell命令的特殊IPython语法）。总下载量约为70MB，解压后的文件约为900MB。
在[1]中：
# !curl -O https://s3.amazonaws.com/pronto-data/open_data_year_one.zip
# !unzip open_data_year_one.zip
接下来我们需要一些标准的Python包导入。
在[2]中：
%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import seaborn as sns; sns.set()
而现在我们用Pandas加载行程数据。
在[3]中：
trips = pd.read_csv('2015_trip_data.csv',
parse_dates=['starttime', 'stoptime'],
infer_datetime_format=True)
trips.head()
Out[3]:
| 行程_id | 开始时间 | 停止时间 | 自行车号 | 旅行时间 | from_station_name | to_station_name | from_station_id | to_station_id | 用户类型 | 性别 | 出生年月 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 431 | 2014-10-13 10:31:00 | 2014-10-13 10:48:00 | SEA00298 | 985.935 | 第2大道与Spring街 | 西方公园/西方大道S & S Washing... | CBD-06 | PS-04 | 年度会员 | 男性 | 1960 |
| 1 | 432 | 2014-10-13 10:32:00 | 2014-10-13 10:48:00 | SEA00195 | 926.375 | 第2大道与Spring街 | 西方公园/西方大道S & S Washing... | CBD-06 | PS-04 | 年度会员 | 男性 | 1970 |
| 2 | 433 | 2014-10-13 10:33:00 | 2014-10-13 10:48:00 | SEA00486 | 883.831 | 第2大道与Spring街 | 西方公园/西方大道S & S Washing... | CBD-06 | PS-04 | 年度会员 | 女性 | 1988 |
| 3 | 434 | 2014-10-13 10:34:00 | 2014-10-13 10:48:00 | SEA00333 | 865.937 | 第2大道与Spring街 | 西方公园/西方大道S & S Washing... | CBD-06 | PS-04 | 年度会员 | 女性 | 1977 |
| 4 | 435 | 2014-10-13 10:34:00 | 2014-10-13 10:49:00 | SEA00202 | 923.923 | 第2大道与Spring街 | 西方公园/西方大道S & S Washing... | CBD-06 | PS-04 | 年度会员 | 男性 | 1971 |
这个行程数据集的每一行都是一个人的一次骑行，数据包含14万多行!
探索不同时期的出行次数
让我们先来看看一年中每日出行次数的趋势
在[4]中：
# Find the start date
ind = pd.DatetimeIndex(trips.starttime)
trips['date'] = ind.date.astype('datetime64')
trips['hour'] = ind.hour
在[5]中：
# Count trips by date
by_date = trips.pivot_table('trip_id', aggfunc='count',
index='date',
columns='usertype', )
在[6]中：
fig, ax = plt.subplots(2, figsize=(16, 8))
fig.subplots_adjust(hspace=0.4)
by_date.iloc[:, 0].plot(ax=ax[0], title='Annual Members');
by_date.iloc[:, 1].plot(ax=ax[1], title='Day-Pass Users');
这张图显示了每天的趋势，按年度会员（顶部）和日票用户（底部）分开。有几个观察结果。
- 4月份短期通票使用量的大幅飙升可能是由于美国规划协会的全国会议，该周在西雅图市中心举行。唯一接近的其他时间是7月4日的周末。
- 日票用户似乎随着季节的变化呈现出稳定的起伏；年票用户的使用率没有随着秋季的到来而明显减弱。
- 年费会员和日票用户似乎都显示出一个明显的每周趋势。
让我们放大这个每周的趋势，通过对每周一天的所有乘坐次数进行平均。由于2015年1月前后模式的变化，我们将按年份和星期来划分数据。
在[7]中：
by_weekday = by_date.groupby([by_date.index.year,
by_date.index.dayofweek]).mean()
by_weekday.columns.name = None # remove label for plot
fig, ax = plt.subplots(1, 2, figsize=(16, 6), sharey=True)
by_weekday.loc[2014].plot(title='Average Use by Day of Week (2014)', ax=ax[0]);
by_weekday.loc[2015].plot(title='Average Use by Day of Week (2015)', ax=ax[1]);
for axi in ax:
axi.set_xticklabels(['Mon', 'Tues', 'Wed', 'Thurs', 'Fri', 'Sat', 'Sun'])
我们看到一个互补的模式：年度用户倾向于在周一至周五期间使用自行车（即作为通勤的一部分），而日票用户倾向于在周末使用自行车。然而，这种模式直到2015年初才完全形成，特别是对于年费会员：似乎在头几个月，用户还没有调整他们的通勤习惯以利用Pronto!
按工作日和周末查看平均每小时的行程也相当有趣。这需要进行一些操作。
在[8]中：
# count trips by date and by hour
by_hour = trips.pivot_table('trip_id', aggfunc='count',
index=['date', 'hour'],
columns='usertype').fillna(0).reset_index('hour')
# average these counts by weekend
by_hour['weekend'] = (by_hour.index.dayofweek >= 5)
by_hour = by_hour.groupby(['weekend', 'hour']).mean()
by_hour.index.set_levels([['weekday', 'weekend'],
["{0}:00".format(i) for i in range(24)]],
inplace=True);
by_hour.columns.name = None
现在我们可以绘制结果，看看每小时的趋势。
在[9]中：
fig, ax = plt.subplots(1, 2, figsize=(16, 6), sharey=True)
by_hour.loc['weekday'].plot(title='Average Hourly Use (Mon-Fri)', ax=ax[0])
by_hour.loc['weekend'].plot(title='Average Hourly Use (Sat-Sun)', ax=ax[1])
ax[0].set_ylabel('Average Trips per Hour');
我们看到 "通勤 "模式和 "休闲 "模式之间有明显的区别，前者在早上和晚上达到高峰（例如，工作日的年费会员），后者在下午早些时候有一个广泛的高峰（例如，周末的年费会员，以及一直以来的短期用户）。有趣的是，年票持有者在周末的平均行为似乎与日票使用者在工作日的平均行为几乎完全相同
对于那些读过我以前的文章的人，你可能会认识到这与我在弗里蒙特桥自行车计数中发现的模式非常相似。
旅行时间
接下来让我们来看看行程的持续时间。Pronto骑行的时间最长为30分钟；任何超过这个时间的单次使用都会产生使用费，前半小时为几美元，此后每小时约10美元。
让我们来看看年度会员和短期通行证持有者的出行时间的分布。
在[10]中：
trips['minutes'] = trips.tripduration / 60
trips.groupby('usertype')['minutes'].hist(bins=np.arange(61), alpha=0.5, normed=True);
plt.xlabel('Duration (minutes)')
plt.ylabel('relative frequency')
plt.title('Trip Durations')
plt.text(34, 0.09, "Free Trips\n\nAdditional Fee", ha='right',
size=18, rotation=90, alpha=0.5, color='red')
plt.legend(['Annual Members', 'Short-term Pass'])
plt.axvline(30, linestyle='--', color='red', alpha=0.3);
这里我加了一条红色的虚线，把免费乘车（左边）和产生使用费的乘车（右边）分开。看来，年卡用户对系统规则更为了解：只有一小部分行程分布的尾部超过30分钟。另一方面，约有四分之一的日票乘车时间超过半小时限制，并产生额外费用。我的直觉是，这些日票用户并不完全了解这种定价结构（"我付了一天的钱，对吗？"），而且很可能对这种体验不满意而离开。
估算行程距离
看一下行程的距离也很有意思。车站之间的距离没有包括在Pronto的数据发布中，所以我们需要通过其他来源找到它们。让我们从加载车站数据开始--由于一些旅行的起点和终点是Pronto的商店，我们将把它作为另一个 "车站 "加入。
在[11]中：
stations = pd.read_csv('2015_station_data.csv')
pronto_shop = dict(id=54, name="Pronto shop",
terminal="Pronto shop",
lat=47.6173156, long=-122.3414776,
dockcount=100, online='10/13/2014')
stations = stations.append(pronto_shop, ignore_index=True)
现在，我们需要找到一对纬度/伦度坐标之间的骑车距离。幸运的是，谷歌地图有一个距离的API，我们可以免费使用。
看了一下小字，免费使用限制在每天2500个距离，每10秒100个距离。有55个站点，我们有55美元/次54/2=1485美元的独特非零距离，所以我们可以在一天的几分钟内免费查询所有的距离（如果我们仔细做的话）。
为了做到这一点，我们将一次查询一个（部分）行，在两次查询之间等待10多秒（注意：我们也可以用googlemaps Python包来代替，但它需要获得一个API密钥）。
在[12]中：
from time import sleep
def query_distances(stations=stations):
"""Query the Google API for bicycling distances"""
latlon_list = ['{0},{1}'.format(lat, long)
for (lat, long) in zip(stations.lat, stations.long)]
def create_url(i):
URL = ('https://maps.googleapis.com/maps/api/distancematrix/json?'
'origins={origins}&destinations={destinations}&mode=bicycling')
return URL.format(origins=latlon_list[i],
destinations='|'.join(latlon_list[i + 1:]))
for i in range(len(latlon_list) - 1):
url = create_url(i)
filename = "distances_{0}.json".format(stations.terminal.iloc[i])
print(i, filename)
!curl "{url}" -o {filename}
sleep(11) # only one query per 10 seconds!
def build_distance_matrix(stations=stations):
"""Build a matrix from the Google API results"""
dist = np.zeros((len(stations), len(stations)), dtype=float)
for i, term in enumerate(stations.terminal[:-1]):
filename = 'queried_distances/distances_{0}.json'.format(term)
row = json.load(open(filename))
dist[i, i + 1:] = [el['distance']['value'] for el in row['rows'][0]['elements']]
dist += dist.T
distances = pd.DataFrame(dist, index=stations.terminal,
columns=stations.terminal)
distances.to_csv('station_distances.csv')
return distances
# only call this the first time
import os
if not os.path.exists('station_distances.csv'):
# Note: you can call this function at most ~twice per day!
query_distances()
# Move all the queried files into a directory
# so we don't accidentally overwrite them
if not os.path.exists('queried_distances'):
os.makedirs('queried_distances')
!mv distances_*.json queried_distances
# Build distance matrix and save to CSV
distances = build_distance_matrix()
这里是距离矩阵的第一个5x5部分的样子。
在[13]中：
distances = pd.read_csv('station_distances.csv', index_col='terminal')
distances.iloc[:5, :5]
Out[13]:
| BT-01 | BT-03 | BT-04 | BT-05 | CBD-13 | |
|---|---|---|---|---|---|
| 终端 | |||||
| --- | |||||
| BT-01 | 0 | 422 | 1067 | 867 | 1342 |
| BT-03 | 422 | 0 | 838 | 445 | 920 |
| BT-04 | 1067 | 838 | 0 | 1094 | 1121 |
| BT-05 | 867 | 445 | 1094 | 0 | 475 |
| CBD-13 | 1342 | 920 | 1121 | 475 | 0 |
让我们把这些距离转换为英里，并把它们加入到我们的行程数据中。
在[14]中：
stacked = distances.stack() / 1609.34 # convert meters to miles
stacked.name = 'distance'
trips = trips.join(stacked, on=['from_station_id', 'to_station_id'])
现在我们可以绘制出行距离的分布。
在[15]中：
fig, ax = plt.subplots(figsize=(12, 4))
trips.groupby('usertype')['distance'].hist(bins=np.linspace(0, 6.99, 50),
alpha=0.5, ax=ax);
plt.xlabel('Distance between start & end (miles)')
plt.ylabel('relative frequency')
plt.title('Minimum Distance of Trip')
plt.legend(['Annual Members', 'Short-term Pass']);
请记住，这显示的是车站之间可能的最短距离，因此是每次旅行的实际骑行距离的下限。许多旅行（尤其是日票用户）在几个街区内开始和结束。除此以外，行程的峰值在1英里左右，尽管一些极端的用户将他们的行程推到了4英里或更远。
估算乘客速度
鉴于这些距离，我们也可以计算出估计的骑行速度的下限。让我们来做这件事，然后看一下年度和短期用户的速度分布。
在[16]中：
trips['speed'] = trips.distance * 60 / trips.minutes
trips.groupby('usertype')['speed'].hist(bins=np.linspace(0, 15, 50), alpha=0.5, normed=True);
plt.xlabel('lower bound riding speed (MPH)')
plt.ylabel('relative frequency')
plt.title('Rider Speed Lower Bound (MPH)')
plt.legend(['Annual Members', 'Short-term Pass']);
有趣的是，这些分布是相当不同的，年度骑行者的平均推断速度更高。你可能会想在这里得出结论，年票会员比日票用户骑得更快，但单凭这些数据并不足以支持这一结论。如果年票用户倾向于通过最直接的路线从A点到B点，而日票用户倾向于蜿蜒曲折，间接到达目的地，也可以解释这个数据。我怀疑实际情况是这两种效应的某种混合。
看一下距离和速度之间的关系也是很有意义的。
在[17]中：
g = sns.FacetGrid(trips, col="usertype", hue='usertype', size=6)
g.map(plt.scatter, "distance", "speed", s=4, alpha=0.2)
# Add lines and labels
x = np.linspace(0, 10)
g.axes[0, 0].set_ylabel('Lower Bound Speed')
for ax in g.axes.flat:
ax.set_xlabel('Lower Bound Distance')
ax.plot(x, 2 * x, '--r', alpha=0.3)
ax.text(9.8, 16.5, "Free Trips\n\nAdditional Fee", ha='right',
size=18, rotation=39, alpha=0.5, color='red')
ax.axis([0, 10, 0, 25])
总的来说，我们看到较长的骑行往往更快--尽管这也受到上述相同的下限警告的限制。如上所述，作为参考，我绘制了免费旅行（红线以上）和需要额外费用的旅行（红线以下）的分离线。我们再次看到，年费会员在不超过半小时的限制方面要比日票用户精明得多--积分分布中的尖锐分界线表明用户在密切跟踪他们的时间，以避免额外的费用。
提升的趋势
在西雅图，人们经常提到的对共享单车可行性的担忧是，这是一个非常多山的城市--在启动之前，坐在椅子上的分析家们预测，将有一个稳定的自行车从上坡流向下坡，这将与其他挑战加在一起，拼出系统的灭亡（"当然，共享单车在其他地方是可行的，但它不可能在这里工作。西雅图很特别！我们很特别！"。我们就是这么特别！")
海拔数据不包括在数据发布中，但我们可以再次求助于谷歌地图API来获得我们所需要的数据；关于海拔API的描述，请看这个网站。
在这种情况下，免费使用的限制是每天2500个请求和每个请求512个海拔。由于我们只需要55个高程，我们可以通过一次查询完成。
在[18]中：
def get_station_elevations(stations):
"""Get station elevations via Google Maps API"""
URL = "https://maps.googleapis.com/maps/api/elevation/json?locations="
locs = '|'.join(['{0},{1}'.format(lat, long)
for (lat, long) in zip(stations.lat, stations.long)])
URL += locs
!curl "{URL}" -o elevations.json
def process_station_elevations():
"""Convert Elevations JSON output to CSV"""
import json
D = json.load(open('elevations.json'))
def unnest(D):
loc = D.pop('location')
loc.update(D)
return loc
elevs = pd.DataFrame([unnest(item) for item in D['results']])
elevs.to_csv('station_elevations.csv')
return elevs
# only run this the first time:
import os
if not os.path.exists('station_elevations.csv'):
get_station_elevations(stations)
process_station_elevations()
现在我们来读入高程数据。
In [19]:
elevs = pd.read_csv('station_elevations.csv', index_col=0)
elevs.head()
Out[19]:
| 高程 | 纬度 | 长度 | 分辨率 | |
|---|---|---|---|---|
| 0 | 37.351780 | 47.618418 | -122.350964 | 76.351616 |
| 1 | 33.815830 | 47.615829 | -122.348564 | 76.351616 |
| 2 | 34.274055 | 47.616094 | -122.341102 | 76.351616 |
| 3 | 44.283257 | 47.613110 | -122.344208 | 76.351616 |
| 4 | 42.460381 | 47.610185 | -122.339641 | 76.351616 |
为了让我们自己感觉好些，我们将仔细检查纬度和经度是否匹配。
在[20]中：
# double check that locations match
print(np.allclose(stations.long, elevs.lng))
print(np.allclose(stations.lat, elevs.lat))
True
True
现在我们可以通过站点数据将海拔与行程数据连接起来。
在[21]中：
stations['elevation'] = elevs['elevation']
elevs.index = stations['terminal']
trips['elevation_start'] = trips.join(elevs, on='from_station_id')['elevation']
trips['elevation_end'] = trips.join(elevs, on='to_station_id')['elevation']
trips['elevation_gain'] = trips['elevation_end'] - trips['elevation_start']
让我们看一下按骑行者类型划分的海拔高度分布。
在[22]中：
g = sns.FacetGrid(trips, col="usertype", hue='usertype')
g.map(plt.hist, "elevation_gain", bins=np.arange(-145, 150, 10))
g.fig.set_figheight(6)
g.fig.set_figwidth(16);
# plot some lines to guide the eye
for lim in range(60, 150, 20):
x = np.linspace(-lim, lim, 3)
for ax in g.axes.flat:
ax.fill(x, 100 * (lim - abs(x)),
color='gray', alpha=0.1, zorder=0)
我们在背景中绘制了一些阴影以帮助引导视线。同样，年度会员和短期用户之间有很大的区别：年度用户肯定表现出对下坡旅行的偏好（分布的左边），而每日用户没有表现得那么强烈，而是表现出对起点和终点在同一海拔高度（即同一车站）的骑行的偏好。
为了使海拔变化的影响更加量化，让我们来计算一下这些数字。
在[23]中：
print
print("total downhill trips:", (trips.elevation_gain < 0).sum()) print("total uphill trips: ", (trips.elevation_gain > 0).sum())
total downhill trips: 80532
total uphill trips: 50493
我们看到，第一年的下坡行程比上坡行程多了30,000次--这大约是60%。考虑到目前的使用水平，这意味着Pronto的工作人员必须每天平均将大约100辆自行车从低洼的站点运送到高处的站点。
天气
另一个常见的 "西雅图是特殊的 "反对自行车共享的可行性的论点是天气问题。让我们来看看骑行次数是如何随温度和降水变化的。
幸运的是，数据发布包括广泛的天气数据。
In[24]:
weather = pd.read_csv('2015_weather_data.csv', index_col='Date', parse_dates=True)
weather.columns
出[24]：
Index(['Max_Temperature_F', 'Mean_Temperature_F', 'Min_TemperatureF', 'Max_Dew_Point_F', 'MeanDew_Point_F', 'Min_Dewpoint_F', 'Max_Humidity', 'Mean_Humidity ', 'Min_Humidity ', 'Max_Sea_Level_Pressure_In ', 'Mean_Sea_Level_Pressure_In ', 'Min_Sea_Level_Pressure_In ', 'Max_Visibility_Miles ', 'Mean_Visibility_Miles ', 'Min_Visibility_Miles ', 'Max_Wind_Speed_MPH ', 'Mean_Wind_Speed_MPH ', 'Max_Gust_Speed_MPH', 'Precipitation_In ', 'Events'],
dtype='object')
让我们把这些天气数据与行程数据结合起来。
In [25]:
by_date = trips.groupby(['date', 'usertype'])['trip_id'].count()
by_date.name = 'count'
by_date = by_date.reset_index('usertype').join(weather)
现在我们可以看一下乘坐次数与温度和降水的关系，将数据按工作日和周末分开。
在[26]中：
# add a flag indicating weekend
by_date['weekend'] = (by_date.index.dayofweek >= 5)
#----------------------------------------------------------------
# Plot Temperature Trend
g = sns.FacetGrid(by_date, col="weekend", hue='usertype', size=6)
g.map(sns.regplot, "Mean_Temperature_F", "count")
g.add_legend();
# do some formatting
g.axes[0, 0].set_title('')
g.axes[0, 1].set_title('')
g.axes[0, 0].text(0.05, 0.95, 'Monday - Friday', va='top', size=14,
transform=g.axes[0, 0].transAxes)
g.axes[0, 1].text(0.05, 0.95, 'Saturday - Sunday', va='top', size=14,
transform=g.axes[0, 1].transAxes)
g.fig.text(0.45, 1, "Trend With Temperature", ha='center', va='top', size=16);
#----------------------------------------------------------------
# Plot Precipitation
g = sns.FacetGrid(by_date, col="weekend", hue='usertype', size=6)
g.map(sns.regplot, "Precipitation_In ", "count")
g.add_legend();
# do some formatting
g.axes[0, 0].set_ylim(-50, 600);
g.axes[0, 0].set_title('')
g.axes[0, 1].set_title('')
g.axes[0, 0].text(0.95, 0.95, 'Monday - Friday', ha='right', va='top', size=14,
transform=g.axes[0, 0].transAxes)
g.axes[0, 1].text(0.95, 0.95, 'Saturday - Sunday', ha='right', va='top', size=14,
transform=g.axes[0, 1].transAxes)
g.fig.text(0.45, 1, "Trend With Precipitation", ha='center', va='top', size=16);
对于温度和降水，我们看到的趋势与我们预期的方向一致（人们在温暖、晴朗的日子里骑行更多）。但也有一些有趣的细节：在工作周，年度会员和短期用户同样受到天气的影响。然而，在周末，年度会员似乎不太受天气的影响，而短期用户似乎更受影响。对于为什么会出现这种趋势，我没有任何好的理论--如果你有任何想法，请让我知道!
总结
从上面的讨论中可以看出很多东西，但我认为有几个主要的启发点，我们可以从这个数据中学习。
- 年度会员和日票用户在总体上表现出明显不同的行为：年度会员似乎主要使用Pronto在周一至周五从A点到B点的通勤，而短期用户主要在周末使用Pronto来探索城市的特定区域。
- 虽然年度会员似乎对定价结构很了解，但每四个短期通行证中就有一个超过了半小时的限制，会产生额外的使用费。为了客户的利益，Pronto也许应该努力让短期用户更好地了解这种定价结构。
- 海拔和天气对使用的影响正如你所期望的那样：下坡行程比上坡行程多60%，寒冷和下雨会大大减少一天的骑行次数。
