python gps 道路_ARCGIS+PYTHON——GPS轨迹点的路网投影与流量计算

一分钟读完全文

通过本文你可以了解到，如何将CSV格式的经纬度点位投影到最近的路段上，并生成新的投影点，以及利用投影后的点位计算路段流量。在操作方法上，同时介绍了手动操作与ARCPY操作。在大批量点位的情况下，建议不使用手动操作，而纯用ARCPY实现。写在前面

拿到轨迹点后，很多时候都会有投影到路网上的需求，这是很多研究或可视化的第一步。利用ARCGIS可以很快将点位与路段进行空间邻近连接。而在此过程中，也往往存在许多问题，如轨迹点往往量非常大，以百度GPS为例，其采集精度为1s，深圳短短15min内就会产生两百多万条记录(还只是我拿到的，不知道有没有做抽样)，因此处理速度与效率往往很关键。另一方面，本文介绍的临近匹配只是非常粗暴的一种匹配，无法实现一些精细化的需求，比如交叉口处的方位偏移，高架和下部道路的区分等等，如果需要考虑这些问题，还需要进一步的在现有基础上加入规则。在ARCMAP中如何操作

很简单，一个NEAR工具箱解决所有问题。

这个工具箱可以返回每个点最近的路段ID，并设定搜索半径(超过该半径则不匹配最近路段)，同时，可以返回往最近路段上的投影点以及最近距离。

一些小TIP：搜索半径有什么意义：理论上，每一个点都可以找到一个最邻近的路段，但实际上，很多噪声点位或跨界点位是不应该被分配路段的。因此，需要预先定义规则，对合理点位进行框选。这时搜索半径便有用了，如设置为50米，作为最大偏移容许值。

返回的最近距离没有单位怎么办：需要设置投影坐标，如深圳则设置UTM ZONE 49投影

以及，在使用他之前你需要做的数据准备是：将CSV中的GPS点位以“添加XY数据”的方式加入，然后导出为SHP文件

加载已有的路网SHP文件，并保证路网与GPS点位在同一坐标系内(如84)

最终返回的数据存在在点的属性表里，打开可以看到表的最后多了几列。利用这些信息，即可生成新的投影点了。ARCPY实现

200万点位，勉强还能手动操作一下，如果直接每天2个亿的点位下来，手动基本是不可能了——还好还有ARCPY。

PS：惊喜的发现10.5的ARCPY里面装了好多外部包，熊猫包也有了，不需要再手动切换环境重新PIP一堆包了，所以软件还是新的好~~

好了，直接上代码。

# 将数据导入

import arcpy

from arcpy import env

import pandas as pd

import os

import glob

from simpledbf import Dbf5

import numpy

import time

env.workspace = r"C:UsersSonghua HuDesktopGPS to RoadTEST"

os.chdir(r"C:UsersSonghua HuDesktopGPS to RoadTEST")

# 导入XY坐标数据

def createxyfromtable(in_Table, x_coords, y_coords, out_Layer, point_shp):

# Make the XY event layer...

arcpy.MakeXYEventLayer_management(in_Table, x_coords, y_coords, out_Layer)

# Process: Copy Features

arcpy.CopyFeatures_management(out_Layer, point_shp, "", "0", "0", "0")

# Print total rows

print(arcpy.GetCount_management(out_Layer))

start = time.time()

# CSV文件转点SHP

createxyfromtable("GPS_WGS_1.csv", "LON", "LAT", "GPS_POINT", "GPS_POINT.shp")

# 调用NEAR工具箱

arcpy.Near_analysis("GPS_POINT.shp", "road.shp", "50 Meters", "LOCATION")

# 生成新的GPS点

createxyfromtable("GPS_POINT.dbf", "NEAR_X", "NEAR_Y", "GPS_POINT_NEW", "GPS_POINT_NEW.shp")

# 报时间

end = time.time()

print(end - start)

# 计算流量

RAW_DBF = Dbf5('GPS_POINT.dbf')

RAW_DF = RAW_DBF.to_dataframe()

volume = RAW_DF.groupby(['NEAR_FID', 'VID']).count().reset_index()[['NEAR_FID', 'VID']].groupby(

'NEAR_FID').count().reset_index()

volume = volume[volume['NEAR_FID'] > 0]

volume.columns = ['NEAR_FID', 'volume_t']

# 数据表单的链接

tuples = [tuple(x) for x in volume.values]

array = numpy.array(tuples, numpy.dtype([('NEAR_FID', numpy.int32), ('volume_t', numpy.int32)]))

arcpy.DeleteField_management("road.dbf", ['volume_t']) # 防止重名插入不了

arcpy.da.ExtendTable("road.dbf", "FID", array, "NEAR_FID")

# 报时间

end = time.time()

print(end - start)

50几行代码，基本把需要的功能都实现了。理一理其中的几个关键点：createxyfromtable函数可以把CSV/DBF格式的经纬度点位直接转为SHP，默认为84坐标系。

Near_analysis就是在调用NEAR工具箱

两次调用createxyfromtable是分别将CSV的原始文件和NEAR处理后的新点位数据转化为SHP。

Dbf5可以快速将DBF转化为熊猫里面的DF格式(大内存会报错，后面修复再更)。

两个groupby是在计算流量：同一辆车只算一次流量，轨迹点数量不等于流量！！

array处理：为了之后与DBF进行连接，DF貌似不能直接做JOIN。

DeleteField_management是防止非KEY重名列，重名就没法做表的连接了。

ExtendTable是将计算出的流量连接到道路属性表中。

——————200万点位做下来，大约是900秒——————

最后上一张投影效果图(直接拿汇报的PPT截图了)：

以及流量图：

基本还是符合现状的，可以看到一些干道流量明显要粗，支路则普遍很细。