-
Notifications
You must be signed in to change notification settings - Fork 0
/
4上下车点的可视化+聚类结果.py
207 lines (185 loc) · 9.19 KB
/
4上下车点的可视化+聚类结果.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
"""***python调用包的引用部分***"""
# 科学计算基础库，用于存储和处理大型矩阵
import numpy as np
# 基于 NumPy 的数据操作库，提供了高效地操作大型数据集所需的工具和方法
import pandas as pd
# 绘图库，支持多种可视化功能
import matplotlib as mpl
import matplotlib.pyplot as plt
# 基于Pandas和Matplotlib的地理数据处理、绘图库，以及GeoPandas的依赖项
import geopandas as gpd
# 用于查找符合特定规则的文件路径名
import glob
# 处理日期和时间的标准库
import time,datetime
# 生成随机数的模块
import random as rd
# 数学库
import math
from math import tan, atan, acos, sin, cos, asin, sqrt, radians
# 平面特征进行集合理论分析和操作,几何对象的基本类型是点、曲线和曲面
from shapely.geometry import LineString,Point
from shapely.ops import nearest_points
# 经纬坐标模块
import pyproj
from pyproj import Transformer
# 基于Python语言的机器学习工具
from sklearn import metrics
from sklearn.cluster import KMeans
from sklearn.cluster import DBSCAN
# 高级Python地理空间绘图库
import geoplot as gplt
# 动态时间规整（DTW）算法的python实现
from fastdtw import fastdtw
# tslearn库中dtw的计算模块
from tslearn.metrics import dtw
# 可以处理插值、积分、优化、图像处理、常微分方程数值解的求解、信号处理等问题的软件包
import scipy
import scipy.cluster.hierarchy as sch
from scipy.spatial.distance import pdist, squareform, euclidean
# 用于创建、操作和研究复杂网络的结构、动态和功能
import networkx as nx
from networkx.algorithms import community
# HTTP库
import requests
"""***代码函数定义部分***"""
# 对某辆出租车轨迹进行分段，提取上下车点
def get_sub_trajectory(df):
loads=[]
no_loads=[]
pick_up=[]
drop_off=[]
# 记录每一段轨迹的开始
idx1 = -1
# 记录每一段轨迹的结束
idx2 = -1
# 记录原始的出租车状态：空车/重车
old_status = ''
for index, row in df.iterrows():
status = row['空车/重车']
# 初始化，当对索引号为0的时
if index==0:
idx1 = index
old_status = status
# 判断状态是否转变，当发生了改变（status != old_status）时
if status != old_status:
# 将两次分界点之间的状态作为一段轨迹，取df的子集df[轨迹开始index1:轨迹结束idx2 + 1]，此处+1是*****
sub_df = df[idx1:idx2 + 1]
# 当状态为'重车'时，说明此时出租车是载客的状态，记录在loads轨迹集中
if old_status=='重车':
loads.append(sub_df)
# 对于下车点记录经纬度(row['经度'],row['纬度'])，记录在drop_off点集中
drop_off.append((row['经度'],row['纬度']))
else:
# 当状态不为‘重车’，即为‘空车’时，说明此时出租车是空车的状态，记录在no_loads数据集中
no_loads.append(sub_df)
# 对于上车点记录经纬度(row['经度'],row['纬度'])，记录在pick_up点集中
pick_up.append((row['经度'],row['纬度']))
# 当前时刻即为下一个时段的开始
idx1 = index
idx2 = index
# 状态更新为下一个时段的状态
old_status = status
else:
# 当状态没有发生改变时，轨迹结束index往后推一个
idx2 = index
# 对最后一个时段的数据进行处理并归档
sub_df = df[idx1:idx2 + 1]
if old_status=='重车':
loads.append(sub_df)
else:
no_loads.append(sub_df)
return loads, no_loads, pick_up, drop_off
"""***代码正文部分***"""
# 记录开始执行代码的时间
time_maincode_start = time.localtime()
print("程序“4上下车点的可视化+聚类结果”的主体代码开始执行的时间为：%s\n" % time.strftime("%Y-%m-%d %H:%M:%S", time_maincode_start))
# 读取经过排序处理后的得到的轨迹数据csv
data_track = pd.read_csv('data/processed/output_data/经过预处理后的轨迹数据201811.csv')
# 对出租车轨迹进行分段，提取上下车点
x, xx, pick_up, drop_off = get_sub_trajectory(data_track)
# 对上下车点构造dataframe便于操作
pick_df = pd.DataFrame(pick_up)
pick_df.columns = ['经度','纬度']
pick_df['类型'] = '上车点'
drop_df = pd.DataFrame(drop_off)
drop_df.columns = ['经度','纬度']
drop_df['类型'] = '下车点'
# 将上车点和下车点的dataframe合并，并重置index，防止索引重复
points_df = pick_df.append(drop_df).reset_index(drop=True)
# 基于上下车点dataframe构造GeoDataFrame数据集
gdf = gpd.GeoDataFrame(
points_df, geometry=gpd.points_from_xy(points_df['经度'], points_df['纬度']),crs=4326)
# 使matplotlib支持中文字体
plt.rcParams['font.family']=['SimHei']
# 绘图，分别指定渲染字段、颜色表、显示图例、点大小、图片大小
ax = gdf.plot(column='类型',cmap='coolwarm',legend=True,markersize=2,figsize=(10,10))
# 设置标题
plt.title("201811 所有出租车的上下车点经纬度显示")
# 叠加武汉市路网
wuhan_road = gpd.GeoDataFrame.from_file('data/road/WuhanAllroad/osmWHmainRoad.shp')
wuhan_road.plot(ax=ax, linewidth=0.5, alpha=0.5, color='grey')
# 设置子图的显示经纬度范围
ax.set_ylim([30,31.5])
ax.set_xlim([113.5,115])
# 保存图片在路径下：data/processed/output_picture/
plt.savefig('data/processed/output_picture/所有出租车的上下车点经纬度显示201811')
# 提取上下车点的地理坐标集合，为聚类做准备
location_points_df = np.array(points_df.loc[:, ['经度', '纬度']])
# 分割字图并设置整个图幅的大小为2000*1000
fig, axes = plt.subplots(1,2,figsize=(20,10))
# 使用k均值聚类kmeans,指定10个类
# sklearn.cluster.KMeans()包含fit(x[,y])方法，用来执行k均值聚类
# 推荐使用fit_predict(x[,y])计算簇的中心并且预测每个样本对应的簇类别，相当于先调用fit(X)再调用predict(X)
kmeans_value = KMeans(n_clusters=10).fit(location_points_df)
# 将得到的K均值聚类结果作为一个新的字段加入上下车点集中
points_df['K'] = kmeans_value.labels_
# 基于加入聚类后的上下车点dataframe构造GeoDataFrame数据集并绘图显示
gdf_kmeans = gpd.GeoDataFrame(
points_df, geometry=gpd.points_from_xy(points_df['经度'], points_df['纬度']),crs=4326)
plt.rcParams['font.family']=['SimHei']
# 生成第一幅子图
ax1 = gdf_kmeans.plot(ax=axes[0],column='K',cmap='coolwarm',markersize=2)
# DBSCAN聚类
# Density-Based Spatial Clustering of Applications with Noise（具有噪声的基于密度的聚类方法）
# 是一种基于密度的空间聚类算法，将具有足够密度的区域划分为簇，并在具有噪声的空间数据库中发现任意形状的簇，它将簇定义为密度相连的点的最大集合。
# 传统的DBSCAN密度聚类算法，需要：邻域阈值(Eps)和点数阈值(minPts)2个参数来对数据集进行聚类
DBSCAN_value = DBSCAN(eps=0.01).fit(location_points_df)
# 将得到的DBSCAN聚类结果作为一个新的字段加入上下车点集中
points_df['DB'] = DBSCAN_value.labels_
# 核心点
core = points_df[points_df['DB']>-1]
# 边界点
border = points_df[points_df['DB']==-1]
# 核心点构造GeoDataFrame并进行可视化
gdf_dbscan = gpd.GeoDataFrame(
core, geometry=gpd.points_from_xy(core['经度'], core['纬度']),crs=4326)
plt.rcParams['font.family']=['SimHei']
# 生成第二幅子图
ax2 = gdf_dbscan.plot(ax=axes[1],column='DB', cmap='coolwarm', markersize=3)
# 边界点构造GeoDataFrame并进行可视化
gdf_dbscan0 = gpd.GeoDataFrame(
border, geometry=gpd.points_from_xy(border['经度'], border['纬度']),crs=4326)
# 核心点的基础上叠加边界点
gdf_dbscan0.plot(ax=axes[1], column='DB', color='black', markersize=1)
# 设置主标题和子标题
plt.suptitle('201811 K均值和DBSCAN聚类后的上下车点经纬度显示')
ax1.set_title('K均值聚类后的上下车点地图可视化')
ax2.set_title('DBSCAN聚类后的上下车点地图可视化')
# 叠加武汉市路网
wuhan_road.plot(ax=ax1, linewidth=0.5, alpha=0.5, color='grey')
wuhan_road.plot(ax=ax2, linewidth=0.5, alpha=0.5, color='grey')
# 设置子图范围
ax1.set_ylim([30,31.5])
ax1.set_xlim([113.5,115])
ax2.set_ylim([30,31.5])
ax2.set_xlim([113.5,115])
plt.savefig('data/processed/output_picture/K均值聚类和DBSCAN聚类后的上下车点经纬度显示201811')
# 记录结束执行代码的时间
time_maincode_end = time.localtime()
print("\n程序“4上下车点的可视化+聚类结果”的主体代码结束执行的时间为：%s" % time.strftime("%Y-%m-%d %H:%M:%S", time_maincode_end))
time_maincode_start = datetime.datetime.strptime(time.strftime("%Y-%m-%d %H:%M:%S", time_maincode_start),"%Y-%m-%d %H:%M:%S")
time_maincode_end = datetime.datetime.strptime(time.strftime("%Y-%m-%d %H:%M:%S", time_maincode_end),"%Y-%m-%d %H:%M:%S")
# 结束时间减去开始时间得到运行的时间
maincode_seconds = (time_maincode_end - time_maincode_start).seconds
print("程序“4上下车点的可视化+聚类结果”的主体代码执行的时间为：%d 秒" % maincode_seconds)