ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

transbigdata出租车GPS轨迹分析:清洗、OD矩阵与热力图实战

transbigdata出租车GPS轨迹分析:清洗、OD矩阵与热力图实战 简介基于Python和transbigdata的出租车轨迹数据可视化分析源码包面向数据分析师、交通领域研究者及需要开展移动轨迹学习的开发人员主要解决出租车GPS轨迹数据的清洗、起讫点分析、网格聚合至地图可视化等环节的快速复现问题。该方案充分利用transbigdata库在空间大数据处理上的优势将复杂运算封装成简洁函数降低使用门槛。整套资源共包含二十六个文件具体包括十五个脚本文件、两个交互式笔记、两个表格数据、两个结构数据、说明文档及开源许可等压缩包总大小约七十MB。脚本文件负责轨迹预处理、质量检查与地图绘制交互式笔记可在浏览器中分步运行并观察结果表格与结构化数据内置上海、深圳两地的样例轨迹说明文档对使用方式作出介绍开源许可则明确项目引用边界。项目采用数据、处理、展示三层分离设计结构清晰便于替换自身数据或调整参数可视化结果同时支持静态图表与动态地图能直观呈现出行时空规律。当前已有四百七十人学习借助该源码可快速跑通出租车轨迹分析整套流程并为城市热点识别、通行效率评估、智能调度等相关研究提供可扩展的实验基础。1. 出租车GPS轨迹分析为什么我推荐这套基于transbigdata的源码出租车GPS数据是我接触过最“脏”的城市数据之一定位漂移、重复上报、超速、空车绕路任何噪声都会在统计结果里被放大。用Python做轨迹可视化分析看起来只是把点画到地图上但真正产出可用的OD矩阵和热点图中间隔着坐标系对齐、轨迹切分、网格聚合三道工序。这套源码把工序拆进了15个Python脚本和2个Jupyter Notebook用transbigdata这个专门做城市数据空间分析的第三方库做底层输入上海和深圳的出租车GPS CSV就能得到上下车点分布、轨迹热力图和OD流向图。对想复现的人来说它最值钱的部分不在最终那张图而在taxigps.py、odprocess.py、grids.py这几个脚本对GPS脏数据的容忍方式以及Jupyter Notebook里一步一步交互调参的节奏。2. transbigdata的地理网格与preprocess.py的清洗逻辑2.1 出租车GPS数据集的字段与噪声打开shenzhen_taxi_gps.csv前几列一般是下面这组字段。不同城市命名略有差异但核心字段就那么几个。字段名类型示例值说明vehicle_idstrT1234车辆唯一标识lngfloat113.9821经度latfloat22.5472纬度timestr2019-07-01 12:00:01上报时间间隔约10秒stateint10为空车1为载客speedfloat45.2瞬时速度单位km/hdirectionint128行驶方向16方位取值拿到CSV直接画图会看到很多诡异现象。城市峡谷里同一辆车10秒内坐标可能跳出去500米等红灯时相同坐标连续上报好几次GPRS缓存还会让时间戳倒挂。这些噪声不清理画出来的轨迹线会横穿楼体OD统计也会多出一堆假行程。所以第一步永远是清洗而不是可视化。2.2 为什么transbigdata把经纬度提前切成网格transbigdata的核心抽象是把连续经纬度离散成网格编号。它不直接比较两个点是否相近而是先把点映射到网格列和网格行后续聚合、邻域搜索、热力统计都在这两个整数列上做。import transbigdata as tbd gridsize 0.01 # 网格边长单位度 df[[LONCOL, LATCOL]] tbd.GPS_to_grid( df[lng], df[lat], gridsize, gridsize )这段代码把原始经纬度变成网格编号。gridsize按度数传0.01度在中纬度地区大约是1公里0.002度大约是200米。transbigdata内部用向量化方式计算几百万行也不会慢。网格编号之间保留了空间邻近关系后面做邻域聚合时可以直接把相邻的LONCOL和LATCOL拼在一起算省掉空间索引的开销。网格粒度不能盲目选小点密度不够时热区会碎成芝麻。常用粒度如下。网格大小实际边长中纬度适合分析的问题注意点0.002约200m路口级热区点少时热区碎片化严重0.005约500m商圈、热点片区推荐起步值0.01约1km城市级OD和热点边界处信息会被抹平0.02约2km跨区通勤细节丢失明显我一般会先用0.01跑一遍看图是否过密或过疏再决定回到0.005还是跳到0.02。这个参数在Jupyter Notebook里改起来非常顺手改完重跑一格就行。2.3 preprocess.py的清洗流程preprocess.py把清洗步骤整理成了可复用的函数顺序一般是去空值、过滤越界、去重、时间排序、速度过滤。以下是一个简化实现代码逻辑与源码里保持一致函数名可能略有出入。def clean_gps_data(df, lonlng, latlat, timetime, vehiclevehicle_id, statestate, speedspeed, bbox(113.5, 21.5, 122.5, 32.5)): df df.dropna(subset[lon, lat, time]) df df[ (df[lon] bbox[0]) (df[lon] bbox[2]) (df[lat] bbox[1]) (df[lat] bbox[3]) ] df df.drop_duplicates(subset[vehicle, time, lon, lat]) df df.sort_values([vehicle, time]).reset_index(dropTrue) if speed is not None: df df[(df[speed] 0) (df[speed] 80)] return dfbbox传的是经纬度范围源码默认值覆盖了上海和深圳所在区域。换成其他城市前必须重设这个框否则会把市外点全部删掉。速度上限80km/h是出租车在市区道路的上限但如果数据里有跨城高速轨迹这个值要抬到140。时间排序必须放在过滤之后否则乱序点可能干扰去重逻辑。清洗后的结果可以用项目中的quality.py输出一行质量指标总点数、删除点数、去重比例、平均采样间隔快速判断清洗是否过激。3. 从轨迹到ODodprocess.py与grids.py的聚合实战3.1 用载客状态切分出行段出租车GPS轨迹里最有价值的信息是载客状态变化。state从0变到1的时刻记录的是上车位置从1变到0的时刻记录的是下车位置。odprocess.py的核心逻辑就是扫描每辆车的状态序列找到这些跳变点。def extract_od(df, vehiclevehicle_id, statestate, timetime, lonlng, latlat): od_list [] for veh_id, trace in df.groupby(vehicle, sortFalse): trace trace.sort_values(time).reset_index(dropTrue) prev_state trace[state].shift(1) start_mask (trace[state] 1) (prev_state ! 1) end_mask (trace[state] ! 1) (prev_state 1) for start_idx in start_mask.index[start_mask]: end_idx end_mask.index[end_mask] end_idx end_idx[end_idx start_idx] if len(end_idx) 0: continue start trace.loc[start_idx] end trace.loc[end_idx[0]] od_list.append({ O_lng: start[lon], O_lat: start[lat], D_lng: end[lon], D_lat: end[lat], start_time: start[time], end_time: end[time], duration_min: ( pd.to_datetime(end[time]) - pd.to_datetime(start[time])).seconds / 60 }) return pd.DataFrame(od_list)这里用shift(1)比较状态变化避免在行级别做循环速度会快很多。注意载客状态的取值在不同数据源里不统一有的是0空车1载客有的反过来。拿到新数据先value_counts看取值再跑extract_od。state字段缺失的数据不能用这个逻辑只能靠停留时间识别上下车点。3.2 剔除无效OD起点终点不能重合时长不能为负状态跳变检测出来之后还要做一层业务过滤。odprocess.py里主要做四件事删除O点和D点坐标相同的记录删除时长小于1分钟或大于180分钟的记录删除平均速度超过合理阈值的记录以及压制状态抖动。od od[(od[O_lng] ! od[D_lng]) | (od[O_lat] ! od[D_lat])] od od[(od[duration_min] 1) (od[duration_min] 180)]状态抖动是最隐蔽的问题。有些车载终端重启会产生一个瞬时错误状态把一条空车轨迹误判成载客。源码里的做法是设置一个最小持续时长参数两次状态跳变之间如果不足30秒就忽略这次跳变继续沿用原状态。这个参数在odprocess.py里通常叫min_interval按秒传。3.3 用grids.py做OD网格聚合拿到OD表之后需要把出发点和到达点分别落到网格上再统计每个网格的出发量和到达量。这就是grids.py的职责。import transbigdata as tbd gridsize 0.01 od[LONCOL], od[LATCOL] tbd.GPS_to_grid( od[O_lng], od[O_lat], gridsize, gridsize ) od[DLONCOL], od[DLATCOL] tbd.GPS_to_grid( od[D_lng], od[D_lat], gridsize, gridsize ) origin_stat od.groupby([LONCOL, LATCOL]).size().reset_index(namecount) origin_stat[lng], origin_stat[lat] tbd.grid_to_center( origin_stat[LONCOL], origin_stat[LATCOL], gridsize, gridsize )GPS_to_grid生成网格编号groupby统计频次最后grid_to_center把网格编号还原成中心点经纬度。画图时直接拿中心点出热力不需要再保留原始点坐标。如果想看网格与网格之间的OD流量可以用透视矩阵od_matrix od.pivot_table( index[LONCOL, LATCOL], columns[DLONCOL, DLATCOL], valuesduration_min, aggfunccount ).fillna(0)这个矩阵是后续计算网格间通勤量和平均时长的基础。网格粒度到0.005时矩阵会非常稀疏大部分格子是0这是正常现象不需要填成平均值或中位数。4. 可视化成图plotmap.py与visualizion.py的叠加思路4.1 GeoJSON地图文件怎么用项目里的shanghai.json和sz.json是市级边界GeoJSON。plotmap.py负责把GeoJSON里的多边形解析出来画到matplotlib坐标轴上作为轨迹图的地理底图。GeoJSON相对shapefile的优势是纯文本Jupyter里直接读不需要安装GDAL。import json import matplotlib.pyplot as plt from matplotlib.path import Path from matplotlib.patches import PathPatch def plot_map(ax, geo_json, lw0.6, facecolornone, edgecolor0.6): for feature in geo_json[features]: geom feature[geometry] if geom[type] Polygon: coords geom[coordinates] elif geom[type] MultiPolygon: coords [c for p in geom[coordinates] for c in p] else: continue for ring in coords: path Path(ring) ax.add_patch(PathPatch(path, lwlw, facecolorfacecolor, edgecoloredgecolor))这里把每一块行政区的边界都添加为PathPatch。facecolor默认透明只显示边界线如果想画区块底色可以传类似#f5f5f5的颜色。项目内的plotmap.py还自动计算了所有polygon的总体边界框并设置到ax.set_xlim和ax.set_ylim上这样底图不会因为单块行政区的坐标跳动而位移。要注意GeoJSON的坐标环是闭合的首位两个点相同Path会正常处理不需要额外去重。项目里与可视化直接相关的文件如下。文件类型在可视化中的作用shanghai.jsonGeoJSON上海市区边界底图sz.jsonGeoJSON深圳市区边界底图frame.pngPNG示例输出展示可视化期望达到的效果visualizion.pyPython封装热力图、OD图的绘制函数plotmap.pyPython封装底图边界绘制函数4.2 把网格统计画成热力图热力图有两条画法直接用hexbin画原始点密度或者先网格聚合再画色块。数据量到百万级时建议用后者hexbin对几百万点的计算开销明显偏高。visualizion.py里封装的热图函数核心逻辑等价于下面的代码。import numpy as np def plot_hotmap(ax, grid_stat, alpha0.75, cmapYlOrRd): lons grid_stat[lng].values lats grid_stat[lat].values cnt grid_stat[count].values img ax.scatter(lons, lats, s8, cnp.log1p(cnt), cmapcmap, alphaalpha, edgecolorsnone) plt.colorbar(img, axax, shrink0.6) return axs8是散点大小具体值依赖网格粒度。网格是0.01度时建议在6到10之间太大会糊成一片太小则热区不明显。透明度alpha在0.5到0.8之间比较合适太低底图透得太清楚热区比例感反而弱。颜色值我用了np.log1p(cnt)而不是原始cnt这一点很关键。出租车上下车点服从长尾分布核心商圈上千次郊区可能几次直接用原始值画色带超过90%的颜色都会集中在小数值上log变换能把中间层次拉开。4.3 OD箭头的批量绘制与视觉优化OD流向图里每一条线代表一个网格到另一个网格的行程。最忌讳的是把所有OD线都画出来几万条线叠在一起就是一坨黑色毛线团。我一般按流量排序后取前100条线线宽用归一化流量映射。from matplotlib.collections import LineCollection od_top od_flow.sort_values(count, ascendingFalse).head(100) lines [ [(row[O_lng], row[O_lat]), (row[D_lng], row[D_lat])] for _, row in od_top.iterrows() ] lc LineCollection(lines, linewidths2, cmapBlues, alpha0.6) ax.add_collection(lc)用LineCollection画线比循环plt.plot快画1000条时优势不明显几万条时能差两个量级。cmapBlues会把线条颜色映射到数值上但如果流量最大值和最小值差距太大颜色层次也不明显同样建议对count做log变换。如果OD图方向特别乱先回到第3章检查state字段是否被正确解析问题通常不在画图而在OD提取。5. 在Jupyter里跑通源码并验证结果5.1 运行环境准备依赖安装很直接transbigdata、pandas、matplotlib、jupyter就够了。geopandas不是必须的因为底图用的是GeoJSON加PathPatch不需要矢量化计算。pip install transbigdata pandas matplotlib jupyter两个Notebook分别对应上海和深圳按单元格顺序跑即可。最容易在这里踩坑的是CSV路径。代码默认读取相对路径shenzhen_taxi_gps.csv如果文件放在data目录下需要改成data/shenzhen_taxi_gps.csv。上海和深圳的CSV字段命名并不完全一致常见差异是时间列叫time或gps_time载客状态叫state或passenger_status建议在Notebook最前面加一个rename统一列名。5.2 两个让人卡住的问题第一是内存占用。几百万行DataFrame在Jupyter里同时保留原始点、OD表、网格统计三个副本16GB内存会告急。处理完一个阶段后及时用df.drop(columns[...])释放中间列或者用del df配合gc.collect()。第二是清洗速度慢。如果总执行时间超过5分钟先怀疑是不是有逐行apply。preprocess.py应该用向量化操作如果自己改了代码把apply换成shift、where、groupby组合速度能提升几十倍。Jupyter Notebook的好处是每个单元格的执行时间会显示可以用来定位具体卡在哪一步。5.3 验证网格聚合没有丢点一个实用的验证方法是把网格统计和原始点分组统计对在一起比确认聚合过程没有引入系统误差。raw df.groupby([LONCOL, LATCOL]).size().rename(raw_cnt).reset_index() grid origin_stat.rename(columns{count: grid_cnt}) merge raw.merge(grid, howouter, on[LONCOL, LATCOL], indicatorTrue) diff (merge[raw_cnt].fillna(0) - merge[grid_cnt].fillna(0)).abs() print(diff.max()) print(merge[_merge].value_counts())diff.max()如果为0说明每个网格的点数完全一致聚合没有丢点。如果出现大量left_only代表聚合时用了不同的gridsize回到GPS_to_grid调用处检查网格参数是否一致。把gridsize换成0.005后重新执行这一步左右两边比例会同步变化这时diff.max()仍然应为0。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进