DCIC2021_task1

赛题理解

赛题目的

通过对车辆数据的综合分析，对厦门岛内早高峰阶段潮汐点进行有效定位，进一步设计高峰期群智优化方案，缓解潮汐点供需问题。

赛题任务

• **任务一：**为更好地掌握早高峰潮汐现象的变化规律与趋势，参赛者需基于主办方提供的数据进行数据分析和计算模型构建等工作，识别出工作日早高峰07:00-09:00潮汐现象最突出的40个区域，列出各区域所包含的共享单车停车点位编号名称，并提供计算方法说明及计算模型，为下一步优化措施提供辅助支撑。

• **任务二：**参赛者根据任务一Top40区域计算结果进一步设计高峰期共享单车潮汐点优化方案，通过主动引导停车用户到邻近停车点位停车，进行削峰填谷，缓解潮汐点停车位（如地铁口）的拥堵问题。
  什么是潮汐？

  城市公共自行车从业者将发生在早晚高峰时段共享单车“借不到、还不进”的问题称之为“潮汐”现象。本题涉及的“潮汐现象”聚焦“还不进”的问题，识别出早高峰共享单车最淤积的40个区域

Baseline中 一些特殊的库的了解

%pylab inline

• % 表示这行代码是魔法命令magic
• inline 表示将图表内嵌到Nb中
• pylab 表示matplotlib与Ipython联合开发的使用matlab的套件

本质上，matplotlib与pylab这些模块其实功能都相同，程序运行的时候都在运行相同的code，不同的是导入模块的方式不同。

需要注意的是，matplotlib有两个使用接口，一种是状态机( state-machine )层的接口，通过pyplot模块来进行管理，这种是我们平时熟悉的调用方式，直接对数据进行plot。

另一种是面向对象的接口,这种是每一个画图的图层都是一个subplot，然后通过对子图对象的控制与操作，实现画图，并且可以实现多图的堆叠。

pylab将所有的功能函数（pyplot状态机函数，大部分时numpy里面的函数）全部导入其单独的命名空间内。为什么要这样做，是因为这样可以很好地与ipython（或者类似的IDE，比如pycharm）实现很好的交互模式，这个就和MATLAB差不多。

而pylab和pyplot的区别是，前者将numpy导入了其命名空间中。这样会使pylab表现的和matlab更加相似。现在来说我们经常使用pyplot，因为pyplot相比pylab更加纯粹。

Geohash

Geohash本质上是一种算法，用来将一个经纬度坐标，转换成一个字符串类型的值，这样的好处是，字符串可以被理解为是一个向量，而这样的坐标字符串就可以被用来比较，而单纯的经纬度坐标没办法比较大小。

事实上，尽管我们平时用一个二维数组[经度,纬度]来描述一个具体的位置。我们本意是通过这个数组去描述一个点，但我们这样做只能获取到这个点附近的一整块区域，这片区域有多大，取决于经纬度的取值精度。

所以Geohash的想法就是：将经纬度的二进制数进行组合，以奇数为纬度，偶数为经度组合，将获取到的经纬度二进制数以每5个数为一组，将每一组都进行转换成十进制数字。然后采用Base32对应编码进行转换可得到编码 wx4g0e这样的可比较的字符串，比如我们的经纬度都分了10次，那么最后生成的字符串的长度就是4，范围是20km，如果我们经纬度都分20次，那么最后生成的字符串的长度就是8，范围可以精确到19m。

Geohash其实就是将整个地图或者某个分割所得的区域进行一次划分，将每一个区域画成一块块矩形块，每个矩形块使用一个字符串表示，当我们需要查询附近的点时，通过自己的坐标计算出一个字符串，根据这个字符串定位到我们所在的矩形块，然后返回这个矩形块中的点。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-FW0CbBIO-1613737350706)(C:\Users\melon\AppData\Roaming\Typora\typora-user-images\image-20210219190655284.png)]

hnswlib

这个库我还没完全的安装上，还报着这样的错误[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-N8RuERWa-1613737350708)(C:\Users\melon\AppData\Roaming\Typora\typora-user-images\image-20210219201317269.png)]

但我还是找到了它的具体原理，其实HNSW是一个搜索算法，全名叫Hierarchical Navigable Small World，他的想法是：是ANN搜索领域基于图的算法，我们要做的是把D维空间中所有的向量构建成一张相互联通的图，并基于这张图搜索某个顶点的K个最近邻。

具体的流程如下：

1. 随机选择一个点作为查询起始点entry_point，把该点加入candidates中，同时加入visitedset
2. 遍历candidates，从candidates中选择距离查询点最近的点c，和results中距离查询点最远的点d进行比较，如果c和查询点q的距离大于d和查询点q的距离，则结束查询，说明当前图中所有距离查询点最近的点已经都找到了，或者candidates为空
3. 从candidates中删除点c
4. 查询c的所有邻居e，如果e已经在visitedset中存在则跳过，不存在则加入visitedset
5. 把比d和q距离更近的e加入candidates、results中，如果results未满，则把所有的e都加入candidates、results。如果results已满，则弹出和q距离最远的点
6. 循环2-5

库的导入 与 数据读取

常见库的导入

import os, codecs
import pandas as pd
import numpy as np

%pylab inline 
'''
% 表示这行代码是魔法命令magic
inline 表示将图表内嵌到Nb中
pylab 表示matplotlib与Ipython联合开发的使用matlab的套件
'''
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('svg') # 使用matplotlib 批量显示svg图片

from matplotlib import font_manager as fm, rcParams
import matplotlib.pyplot as plt

Populating the interactive namespace from numpy and matplotlib

读取数据集

共享的单车订单数据的聚合

#批量读取数据集并聚合
PATH = 'G:/Data Minning/DCIC2021/' 
bike_track = pd.concat([
    pd.read_csv(PATH + 'gxdc_gj20201221.csv'),
    pd.read_csv(PATH + 'gxdc_gj20201222.csv'),
    pd.read_csv(PATH + 'gxdc_gj20201223.csv'),
    pd.read_csv(PATH + 'gxdc_gj20201224.csv'),
    pd.read_csv(PATH + 'gxdc_gj20201225.csv')
])

#排序显示
# 按照单车ID和时间进行排序
bike_track = bike_track.sort_values(['BICYCLE_ID', 'LOCATING_TIME'])

bike_track.describe(include = 'all')

	BICYCLE_ID	LOCATING_TIME	LATITUDE	LONGITUDE	source	LOCATING_TIME1	date	day	time	hour
count	11920885	11920885	1.192088e+07	1.192088e+07	11920885	11920885	11920885	1.192088e+07	11920885	1.192088e+07
unique	103085	72000	NaN	NaN	1	72000	5	NaN	14400	NaN
top	1f14a3dfdcd045f8d340482751c9bfba	2020-12-24 08:22:42	NaN	NaN	CSXZGLZFJ	2020-12-24 08:22:42	2020-12-24	NaN	08:21:07	NaN
freq	3158	483	NaN	NaN	11920885	483	3096835	NaN	1641	NaN
mean	NaN	NaN	2.450373e+01	1.181120e+02	NaN	NaN	NaN	2.305376e+01	NaN	7.742647e+00
std	NaN	NaN	3.822377e-02	5.177484e-02	NaN	NaN	NaN	1.475251e+00	NaN	9.090383e-01
min	NaN	NaN	2.442406e+01	1.179086e+02	NaN	NaN	NaN	2.100000e+01	NaN	6.000000e+00
25%	NaN	NaN	2.447858e+01	1.180859e+02	NaN	NaN	NaN	2.200000e+01	NaN	7.000000e+00
50%	NaN	NaN	2.449319e+01	1.181127e+02	NaN	NaN	NaN	2.300000e+01	NaN	8.000000e+00
75%	NaN	NaN	2.452017e+01	1.181468e+02	NaN	NaN	NaN	2.400000e+01	NaN	8.000000e+00
max	NaN	NaN	2.610928e+01	1.193015e+02	NaN	NaN	NaN	2.500000e+01	NaN	9.000000e+00

bike_track['BICYCLE_ID'].nunique()

103085

地图可视化

import folium

m = folium.Map(location=[24.482426, 118.157606], zoom_start=12)
my_PolyLine=folium.PolyLine(locations=bike_track[bike_track['BICYCLE_ID'] == '000152773681a23a7f2d9af8e8902703'][['LATITUDE', 'LONGITUDE']].values,weight=5)
m.add_children(my_PolyLine)

<ipython-input-9-047b473c8b51>:3: FutureWarning: Method `add_children` is deprecated. Please use `add_child` instead.
  m.add_children(my_PolyLine)

Make this Notebook Trusted to load map: File -> Trust Notebook

def bike_fence_format(s):
    s = s.replace('[', '').replace(']', '').split(',')
    s = np.array(s).astype(float).reshape(5, -1)
    return s

# 共享单车停车点位（电子围栏）数据
bike_fence = pd.read_csv(PATH + 'gxdc_tcd.csv')
bike_fence['FENCE_LOC'] = bike_fence['FENCE_LOC'].apply(bike_fence_format)

# 共享单车订单数据
bike_order = pd.read_csv(PATH + 'gxdc_dd.csv')
bike_order = bike_order.sort_values(['BICYCLE_ID', 'UPDATE_TIME'])

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

from geopy.distance import geodesic

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]))

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]))

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)
)

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)
)

# bike_order
geohash.encode(24.521156, 118.140385, precision=6), \
geohash.encode(24.521156, 118.140325, precision=6)

('wsk52r', 'wsk52r')

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['wsk52r'].plot()
bike_outflow.loc['wsk52r'].plot()
plt.xticks(list(range(bike_inflow.shape[1])), bike_inflow.columns, rotation=40)
plt.legend(['Inflow', 'OutFlow'])

<matplotlib.legend.Legend at 0x2b0475be670>

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Aq4aIZx9-1613737350710)(output_21_1.svg)]

bike_inflow.loc['wsk596'].plot()
bike_outflow.loc['wsk596'].plot()
plt.xticks(list(range(bike_inflow.shape[1])), bike_inflow.columns, rotation=40)
plt.legend(['Inflow', 'OutFlow'])

<matplotlib.legend.Legend at 0x2b03b3241c0>

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-n9pCXa5V-1613737350712)(output_22_1.svg)]

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)
bike_fence['DENSITY'] = bike_fence['geohash'].map(bike_remain).fillna(0)

按照最近的经纬度

思路： 按照订单计算与停车点的距离计算潮汐点；

潮汐统计

方法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()

bike_density

	STREET	0
0	新丰路(火炬路至火炬北路)	264.109027
1	金榜西路0	261.504821
2	市政19	196.202314
3	市政2	192.822917
4	市政1	192.085638
...	...	...
988	莲前西路(金尚路至梧村隧道段)	0.000000
989	莲前西路(金尚路至东浦路)	0.000000
990	围里路	0.000000
991	围里路0	0.000000
992	龙虎西二里	0.000000

993 rows × 2 columns

方法2：基于KNN密度匹配

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

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

from sklearn.neighbors import NearestNeighbors

knn = NearestNeighbors(metric = "haversine", n_jobs=-1, algorithm='brute')
knn.fit(np.stack(bike_fence['FENCE_CENTER'].values))

NearestNeighbors(algorithm='brute', metric='haversine', n_jobs=-1)

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

dist, index = knn.kneighbors(bike_order[['LATITUDE','LONGITUDE']].values[:], n_neighbors=1)

也可以使用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)

bike_density['label'] = '0'
bike_density.iloc[:100, -1] = '1'

bike_density['BELONG_AREA'] ='厦门'
bike_density = bike_density.drop(0, axis=1)

bike_density.columns = ['FENCE_ID', 'FENCE_TYPE', 'BELONG_AREA']
bike_density.to_csv('result.txt', index=None, sep='|')