任务2：共享单车潮汐点分析

发布时间：2021-01-25，作者：阿水, 鱼遇雨欲语与余

## 学习目标

• 共享单车订单与停车点匹配
• 统计并可视化停车点潮汐情况
• 计算工作日早高峰07:00-09:00潮汐现象最突出的40个区域

## 经纬度匹配

在本赛题中如果想要完成具体潮汐计算，则需要将订单或轨迹与具体的停车点进行匹配，需要计算在不同时间下每个停车点：

• 停了多少车（订单终点为停车点）；
• 骑走多少车（订单起点为停车点）；

因此我们需要根据订单定位到具体的停车点：

### 停车点处理

首先我们对停车点进行处理，需要计算得出每个停车点的面积和中心经纬度：

# 得出停车点 LATITUDE 范围
bike_fence['MIN_LATITUDE'] = bike_fence['FENCE_LOC'].apply(lambda x: np.min(x[:, 1]))
bike_fence['MAX_LATITUDE'] = bike_fence['FENCE_LOC'].apply(lambda x: np.max(x[:, 1]))

# 得到停车点 LONGITUDE 范围
bike_fence['MIN_LONGITUDE'] = bike_fence['FENCE_LOC'].apply(lambda x: np.min(x[:, 0]))
bike_fence['MAX_LONGITUDE'] = bike_fence['FENCE_LOC'].apply(lambda x: np.max(x[:, 0]))

from geopy.distance import geodesic
# 根据停车点 范围 计算具体的面积
bike_fence['FENCE_AREA'] = bike_fence.apply(lambda x: geodesic(
    (x['MIN_LATITUDE'], x['MIN_LONGITUDE']), (x['MAX_LATITUDE'], x['MAX_LONGITUDE'])
).meters, axis=1)

# 根据停车点 计算中心经纬度
bike_fence['FENCE_CENTER'] = bike_fence['FENCE_LOC'].apply(
    lambda x: np.mean(x[:-1, ::-1], 0)
)

中心经纬度可以用来具体的距离匹配，当然如果用四个点的坐标会更加准确，但计算时间更大。

### Geohash经纬度匹配

在赛题中计算单车与停车点的关系主要的就在经纬度计算，如何将单车维度与停车点进行匹配？这个问题会影响具体的匹配精度，同时也会影响匹配的速度。在此我们使用geohash库来对经纬度进行编码，这样可以增加统计的速度，也可以动态调整经纬匹配的范围。

如果想深入了解geohash，可以参考以下链接：

• geohash介绍
• DCIC2020

通过geohash可以将一个经纬度编码为一个字符串，需要注意的是precision控制了编码精度，数值越大覆盖范围越小，数值越小覆盖范围越大。

需要注意：这里的precision与具体的定位点范围强相关，需要根据具体的赛题要求完成设置。

import geohash
bike_order['geohash'] = bike_order.apply(
    lambda x: geohash.encode(x['LATITUDE'], x['LONGITUDE'], precision=6), 
axis=1)

bike_fence['geohash'] = bike_fence['FENCE_CENTER'].apply(
    lambda x: geohash.encode(x[0], x[1], precision=6)
)

使用ws7gx9对单车订单数据进行索引，可以得到具体订单数据：

# bike_order[bike_order['geohash'] == 'ws7gx9']

## 区域流量与潮汐统计

在完成具体的经纬度匹配后，接下来就需要完成具体的区域流量统计，即统计某一范围内的不同时间的流量（入流量和出流量）。

首先对订单数据进行时间提取：

bike_order['UPDATE_TIME'] = pd.to_datetime(bike_order['UPDATE_TIME'])
bike_order['DAY'] = bike_order['UPDATE_TIME'].dt.day.astype(object)
bike_order['DAY'] = bike_order['DAY'].apply(str)

bike_order['HOUR'] = bike_order['UPDATE_TIME'].dt.hour.astype(object)
bike_order['HOUR'] = bike_order['HOUR'].apply(str)
bike_order['HOUR'] = bike_order['HOUR'].str.pad(width=2,side='left',fillchar='0')

# 日期和时间进行拼接
bike_order['DAY_HOUR'] = bike_order['DAY'] + bike_order['HOUR']

使用透视表统计每个区域在不同时间的入流量和出流量：

bike_inflow = pd.pivot_table(bike_order[bike_order['LOCK_STATUS'] == 1], 
                   values='LOCK_STATUS', index=['geohash'],
                    columns=['DAY_HOUR'], aggfunc='count', fill_value=0
)

bike_outflow = pd.pivot_table(bike_order[bike_order['LOCK_STATUS'] == 0], 
                   values='LOCK_STATUS', index=['geohash'],
                    columns=['DAY_HOUR'], aggfunc='count', fill_value=0
)

通过以上代码就可以统计完成区域流量，接下来可以做简单的绘图：

bike_inflow.loc['wsk593'].plot()
bike_outflow.loc['wsk593'].plot()
plt.xticks(list(range(bike_inflow.shape[1])), bike_inflow.columns, rotation=40)
plt.legend(['入流量', '出流量'])

bike_inflow.loc['wsk52r'].plot()
bike_outflow.loc['wsk52r'].plot()
plt.xticks(list(range(bike_inflow.shape[1])), bike_inflow.columns, rotation=40)
plt.legend(['入流量', '出流量'], prop = prop)

## 方法1：Geohash匹配计算潮汐

由于赛题需要统计工作日早高峰期间的潮汐现象，所以我们可以按照天进行单车流量统计：

bike_inflow = pd.pivot_table(bike_order[bike_order['LOCK_STATUS'] == 1], 
                   values='LOCK_STATUS', index=['geohash'],
                    columns=['DAY'], aggfunc='count', fill_value=0
)

bike_outflow = pd.pivot_table(bike_order[bike_order['LOCK_STATUS'] == 0], 
                   values='LOCK_STATUS', index=['geohash'],
                    columns=['DAY'], aggfunc='count', fill_value=0
)

根据入流量和出流量，可以计算得到每个位置的留存流量：

bike_remain = (bike_inflow - bike_outflow).fillna(0)

# 存在骑走的车数量 大于 进来的车数量
bike_remain[bike_remain < 0] = 0  

# 按照天求平均
bike_remain = bike_remain.sum(1)

这里假设我们需要统计街道维度的潮汐情况，我们可以先把街道信息提取，然后计算密度。这里我们需要计算每个街道不同停车点的留存车辆，所以不能重复统计。

# 总共有993条街
bike_fence['STREET'] = bike_fence['FENCE_ID'].apply(lambda x: x.split('_')[0])

# 留存车辆 / 街道停车位总面积，计算得到密度
bike_density = bike_fence.groupby(['STREET'])['geohash'].unique().apply(
    lambda hs: np.sum([bike_remain[x] for x in hs])
) / bike_fence.groupby(['STREET'])['FENCE_AREA'].sum()

# 按照密度倒序
bike_density = bike_density.sort_values(ascending=False).reset_index()

就可以得到潮汐情况最严重的道路：

## 方法2：距离匹配计算潮汐

如果使用Geohash来统计会存在一个问题，统计的方法会不准确，导致只能精确到街道信息。本节将使用经纬度距离匹配的方法来进行尝试，具体的思路为计算订单最近的停车点，进而计算具体的潮汐情况。

对于经纬度距离计算，可以直接使用sklearn中的NearestNeighbors，通过设置haversine距离可以很方便的完成最近停车点的计算。

from sklearn.neighbors import NearestNeighbors

# https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.DistanceMetric.html
knn = NearestNeighbors(metric = "haversine", n_jobs=-1, algorithm='brute')
knn.fit(np.stack(bike_fence['FENCE_CENTER'].values))

计算订单中对应的停车点位置：

# 需要11s左右
dist, index = knn.kneighbors(bike_order[['LATITUDE','LONGITUDE']].values[:20000], n_neighbors=1)

但是如果直接使用NearestNeighbors计算速度会非常慢，如果是全量定量订单数据可能需要较长时间。因此可以用hnsw做近似搜索，速度较快但精度差一点。

import hnswlib
import numpy as np

p = hnswlib.Index(space='l2', dim=2)
p.init_index(max_elements=300000, ef_construction=1000, M=32)
p.set_ef(1024)
p.set_num_threads(14)

p.add_items(np.stack(bike_fence['FENCE_CENTER'].values))

计算所有订单的停车位置：

index, dist = p.knn_query(bike_order[['LATITUDE','LONGITUDE']].values[:], k=1)

然后计算所有停车点的潮汐流量：

bike_order['fence'] = bike_fence.iloc[index.flatten()]['FENCE_ID'].values

bike_inflow = pd.pivot_table(bike_order[bike_order['LOCK_STATUS'] == 1], 
                   values='LOCK_STATUS', index=['fence'],
                    columns=['DAY'], aggfunc='count', fill_value=0
)

bike_outflow = pd.pivot_table(bike_order[bike_order['LOCK_STATUS'] == 0], 
                   values='LOCK_STATUS', index=['fence'],
                    columns=['DAY'], aggfunc='count', fill_value=0
)

bike_remain = (bike_inflow - bike_outflow).fillna(0)
bike_remain[bike_remain < 0] = 0  
bike_remain = bike_remain.sum(1)

计算停车点的密度：

bike_density = bike_remain / bike_fence.set_index('FENCE_ID')['FENCE_AREA']
bike_density = bike_density.sort_values(ascending=False).reset_index()
bike_density = bike_density.fillna(0)

## 改进方向

• 根据赛题提交的要求，可能与需要提交的区域信息存在差异，所以后续可能要根据赛题要求完成具体的计算，具体要求不同（从街道维度到停车点）经纬度编码的precision不同。
• 具体密度计算方法可能需要根据具体的流量和停车点面积完成计算，具体逻辑计算可能需要调整；