使用Python和Pandas分析Pronto CycleShare數據
這是一篇非常不錯的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]:
| trip id | starttime | stoptime | bikeid | tripduration | fromstation name | tostation name | fromstation id | tostation_id | usertype | gender | birthyear | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 431 | 2014-10-13 10:31:00 | 2014-10-13 10:48:00 | SEA00298 | 985.935 | 2nd Ave & Spring St | Occidental Park / Occidental Ave S & S Washing... | CBD-06 | PS-04 | Annual Member | Male | 1960 |
| 1 | 432 | 2014-10-13 10:32:00 | 2014-10-13 10:48:00 | SEA00195 | 926.375 | 2nd Ave & Spring St | Occidental Park / Occidental Ave S & S Washing... | CBD-06 | PS-04 | Annual Member | Male | 1970 |
| 2 | 433 | 2014-10-13 10:33:00 | 2014-10-13 10:48:00 | SEA00486 | 883.831 | 2nd Ave & Spring St | Occidental Park / Occidental Ave S & S Washing... | CBD-06 | PS-04 | Annual Member | Female | 1988 |
| 3 | 434 | 2014-10-13 10:34:00 | 2014-10-13 10:48:00 | SEA00333 | 865.937 | 2nd Ave & Spring St | Occidental Park / Occidental Ave S & S Washing... | CBD-06 | PS-04 | Annual Member | Female | 1977 |
| 4 | 435 | 2014-10-13 10:34:00 | 2014-10-13 10:49:00 | SEA00202 | 923.923 | 2nd Ave & Spring St | Occidental Park / Occidental Ave S & S Washing... | CBD-06 | PS-04 | Annual Member | Male | 1971 |
這個行程數據集的每一行是由一個人單獨騎行，共包含超過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]:
| BT-01 | BT-03 | BT-04 | BT-05 | CBD-13 | |
|---|---|---|---|---|---|
| terminal | |||||
| 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 |
讓我們將這些距離轉換為英里，并將它們加入我們的行程數據：
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 獲取我們需要的數據;。
在這種情況下，自由使用限制為每天 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]:
| elevation | lat | lng | resolution | |
|---|---|---|---|---|
| 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 |
為了驗證結果，我們需要仔細檢查緯度和經度是否匹配：
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%，寒冷和下雨也會明顯減少當天的騎行數量。
來自：http://ipfans.github.io/2017/02/analyzing-pronto-cycleshare-data-with-python-and-pandas/