这是一篇非常不错的 pandas 分析入门文章,在此简单翻译摘录如下。
本周,西雅图的自行车共享系统 Pronto CycleShare 一周岁了。 为了庆祝这一点,Pronto 提供了从第一年的数据缓存,并宣布了 Pronto Cycle Share 的数据分析挑战。
你可以用很多工具分析这些数据,但我的选择工具是 Python。 在这篇文章中,我想展示如何开始分析这些数据,并使用 PyData 技术栈,即 NumPy,Pandas,Matplotlib 和 Seaborn 与其他可用的数据源。
这篇文章以 Jupyter Notebook 形式组织,它是一种开放的文档格式。结合了文本、代码、数据和图形,并且通过 Web 浏览器查看。本文中的内容可以下载 对应的 Notebook 文件,并通过 Jupyter 打开。
下载 Pronto 的数据
我们可以从 Pronto 官网 下载对应的 数据文件。总下载大约 70MB,解压缩的文件大约 900MB。
接下来我们需要导入一些 Python 包:
In [2]:
%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import seaborn as sns; sns.set()
现在我们使用 Pandas 加载所有的行程记录:
In [3]:
trips = pd.read_csv('2015_trip_data.csv',
parse_dates=['starttime', 'stoptime'],
infer_datetime_format=True)
trips.head()
Out[3]:
这个行程数据集的每一行是由一个人单独骑行,共包含超过 140,000 条数据。
探索时间与行程的关系
让我们先看看一年中每日行程次数的趋势。
In [4]:
# Find the start date
ind = pd.DatetimeIndex(trips.starttime)
trips['date'] = ind.date.astype('datetime64')
trips['hour'] = ind.hour
In [5]:
# Count trips by date
by_date = trips.pivot_table('trip_id', aggfunc='count',
index='date',
columns='usertype', )
In [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 月份左右模式的变化,我们按照年份和星期几进行拆分:
In [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 的通勤习惯。
查看平日和周末的平均每小时骑行也很有趣。这需要一些操作:
In [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
现在我们可以绘制结果来查看每小时的趋势:
In [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 分钟; 任何长于此的单次使用,在前半个小时都会产生几美元的使用费,此后每小时大约需要十美元。
让我们看看年费会员和临时使用者的旅行持续时间的分布:
In [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 的服务点之间开始和结束,我们将其添加为一个 “车站”:
In [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)
现在我们需要找到两对纬度 / 经度坐标之间的骑车距离。幸运的是,Google 地图有一个距离 API,我们可以免费使用。
从文档中知道,我们每天免费使用的限制为每天最多 2500 个距离,每 10 秒最多 100 个距离。现在有 55 个站,我们有(55 * 54/2) = 1485 个非零距离,所以我们可以在几天内免费查询所有车站之间的距离。
为此,我们一次查询一行,在查询之间等待 10 + 秒(注意:我们可能还会使用 googlemaps Python 包 ,但使用它需要获取 API 密钥)。
In [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 距离矩阵:
In [13]:
distances = pd.read_csv('station_distances.csv', index_col='terminal')
distances.iloc[:5, :5]
Out[13]:
让我们将这些距离转换为英里,并将它们加入我们的行程数据:
In [14]:
stacked = distances.stack() / 1609.34 # convert meters to miles
stacked.name = 'distance'
trips = trips.join(stacked, on=['from_station_id','to_station_id'])
现在我们可以绘制行程距离的分布:
In [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 英里左右,也有一些用户将他们的旅行距离扩展到四英里或更长的距离。
骑手速度
给定这些距离,我们还可以计算估计骑行速度的下限。 让我们这样做,然后看看年费用户和临时用户的速度分布:
In [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,那么这些数据也可以被解释,而临时用户倾向于绕行并间接到达他们的目的地。我怀疑现实是这两种效应的混合。
还要看看距离和速度之间的关系:
In [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])
总的来说,我们看到较长的路途速度更快 - 虽然这受到与上述相同的下限影响。如上所述,作为参考,我绘制了需要的红线用于区分额外费用(低于红线)和免费费用(红线以上)。我们再次看到,年度会员对于不超过半小时的限制比每天通过用户更加精明 - 点的分布的指向了用户注意了他们使用的时间,以避免额外的费用。
海拔高度
在西雅图自行车分享服务的可行性的一个焦点是,西雅图是一个丘陵城市。在服务发布之前,一些分析师预测,西雅图会有源源不断的自行车上坡下坡,所以并不适合分享单车系统的落地。
数据版本中不包含海拔高度数据,但我们可以转到 Google Maps API 获取我们需要的数据; 请参阅 此网站 了解海拔 API 的描述。
在这种情况下,自由使用限制为每天 2500 个请求,每次请求最多包含 512 个海拔高度。 由于我们只需要 55 个海拔高度,我们可以在单个查询中执行:
In [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]:
为了验证结果,我们需要仔细检查纬度和经度是否匹配:
In [20]:
# double check that locations match
print(np.allclose(stations.long, elevs.lng))
print(np.allclose(stations.lat, elevs.lat))
True
True
现在我们可以将海拔数据与行程数据整合:
In [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']
让我们来看看海拔数据和会员类型的分布关系:
In [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)
我们在背景中绘制了一些阴影以帮助引导分析。 年度会员和临时用户之间有很大的区别:年费用户非常明显的表示出偏好下坡行程(左侧的分布),而临时用户表现并不明显,而是表示喜欢骑乘开始并在相同高度结束。
为了使海拔数据变化的影响更加明显,我们做一些计算:
In [23]:
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
我们看到,第一年下坡比上坡多出了 3 万次 - 这是大约 60% 以上。 根据目前的使用水平,这意味着 Pronto 工作人员必须每天从海拔较低的服务点运送大约 100 辆自行车到高海拔服务点。
天气
另一个常见的反对循环共享的可行性的论点是天气。让我们来看看出行数量随着温度和降水量的变化。
幸运的是,数据包括了大范围的天气数据:
In [24]:
weather = pd.read_csv('2015_weather_data.csv', index_col='Date', parse_dates=True)
weather.columns
Out[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')
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)
现在我们可以看看按工作日和周末为纬度,查看出行数量随温度和降水量的变化:
In [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 进行工作日的通勤。临时用户则是周末使用 Pronto 探索城市的特定区域。
- 尽管年费用户对定价策略有所了解,但是四分之一的行程还是超过了半小时的限制,并产生了额外的费用。为了客户的权益,Pronto 应该更好告知用户这种定价策略。
- 海拔和天气会影响使用,正如你预计的一样:下坡比上坡多 60%,寒冷和下雨也会明显减少当天的骑行数量。