python交通数据分析：车辆轨迹点进行路网匹配并更新轨迹点位置

Python交通数据分析：车辆轨迹点进行路网匹配并更新轨迹点位置

前言

对车辆轨迹点进行路网匹配并更新轨迹点的位置。原理很简单，就是搜索道路缓冲区内的轨迹点，识别到它的id和位置，对车辆轨迹点进行匹配，并将轨迹点更新在道路的位置。适用于一些轨迹点的数据预处理场景。

---

一、数据准备

1.1 环境和库

Arcgis 10.8
我用的10.8，任意版本应该都可以

Python2.7

Arcpy 库
ArcPy是一个 Python 包，可提供以实用高效的方式通过 Python 执行地理数据分析、数据转换、数据管理和地图自动化。
Arcpy要和Arcgis一起用，必须要有arcgis桌面版或者arcgis pro存在
ArcPy仅能依赖于ArcGIS/GeoScene平台存在，因为ArcPy的大部分功能接口，只是通过Python调度ArcGIS/GeoScene底层的核心接口，所以不能脱离ArcGIS/GeoScene平台存在
arcpy官网

1.2 路网数据

实验所使用的路网数据，是经过数据预处理后的道路中心线数据，也对立交桥、交叉路口等区域进行了简化。格式为gis数据中的线要素（polyline），本人用道路中心线提取的。此外，路网需要拓扑关系已经经过arcgis的拓扑检查。路网能够完全构成网络。

什么是gis中的拓扑检查
拓扑错误

能够看到路网数据线要，能够保证每一条路的两端，至少有一端有与它相连接的另一条道路

OBJECTID是程序识别道路id的所用字段，如果还有其他要导入的路网属性，可以添加在后面。也可以在完成路网匹配后，做一些join等操作。

1.3 轨迹点要素

轨迹点的数据为包含了经度、纬度和其他属性的csv格式的文件，命名为**“gps.csv”**
点要素包含了**‘CARNO’, ‘LONGITUDE’, ‘LATITUDE’,‘POSITION_TIME’, ‘newid’**五个属性
分别表示车牌号、经度、纬度、定位时间、轨迹点唯一的id
在实验里使用经度、纬度、唯一索引的id就可以了，其他的属性可以自己添加。要注意的是，定位时间（如1990-1-1 12:12:12）这类python里能够表示为time时间格式的字段，在arcpy里是没办法识别和导出的。

1.4 建立地理数据库

arcpy需要使用工作空间
什么是地理数据库gbd

在arcgis map中打开catalog

添加要创建地理数据库的文件夹

右键创建一个地理数据库，后缀名为gbd

将路网数据 road_network.shp 放在gbd里

二、

需要三个py的文件

1.导入数据，csv转换为shp

创建名为importData.py的py文件

代码如下：

# coding=utf-8
import arcpy

'''
做了csv转化为要素类的接口，让roadMatching.py调用的
方法需要返回转为POINT的文件列表
'''

# 获得数据的的坐标系
spatial_re = arcpy.Describe("D:/gisdata/test.gdb/road_network").spatialReference
arcpy.env.workspace = 'D:/gisdata/test.gdb'#工作空间为你创建的gbd地理数据库

def csvChange(lines, outName):
    arcpy.CreateFeatureclass_management('D:/gisdata/test.gdb', outName, 'POINT',
                                        spatial_reference=spatial_re)
    print ('create')
    
    #以下为点要素的属性
    arcpy.AddField_management(outName, 'CARNO', 'TEXT')
    arcpy.AddField_management(outName, 'LONGITUDE', 'TEXT')
    arcpy.AddField_management(outName, 'LATITUDE', 'TEXT')
    arcpy.AddField_management(outName, 'POSITION_TIME', 'TEXT')
    arcpy.AddField_management(outName, 'newid', 'TEXT')


    point = arcpy.Point()  # create an empty Point object
    rows = arcpy.InsertCursor(outName)
    # 表格数据第一行是表头，所以跳过表头
    for line in lines[1:]:
        row = rows.newRow()
        line = line.split(',')
     
        point.X = float(line[1][:])#x是经度，1是经度所在的列索引
        point.Y = float(line[2][:])#y是纬度，2是纬度所在的列索引
        pointGeometry = arcpy.PointGeometry(point)
        row.shape = pointGeometry
        
        row.CARNO = line[0][:]
        row.LONGITUDE = line[1][:]
        row.LATITUDE = line[2][:]
        row.POSITION_TIME = line[3][:]
        row.newid = line[4][:]

        rows.insertRow(row)
        print(line)


#csv格式的文件转换为点要素
def csvToPoint(inpath, outName):
    with open(inpath, 'r') as fr:
        lines = fr.readlines()
        csvChange(lines, outName)
    print ('csv to point over!')


if __name__ == '__main__':

    print (csvToPoint(r'轨迹点的源文件路径.csv', 'taxi'))

2.更新轨迹点位置

创建名为roadMatching.py的py文件
代码如下：

# coding=utf-8
import arcpy
//工作空间为你所创建的地理数据库位置
arcpy.env.workspace = 'D:/gisdata/test.gdb'

#这里的100米为搜索范围，可以自定义，搜索范围越大则成功匹配的点越多，因为有的轨迹点会搜不到指定半径内的道路
def updateGPSdata(firstInterOut, road_name):
    arcpy.Near_analysis(firstInterOut, road_name, "100 Meters", location='LOCATION')
    print ('near analysis')

    # 提取最近邻的坐标点，将最新的坐标点存入dic中
    dic = {}
    #OBJECTID为路网每一条路的唯一id
    cursor = arcpy.da.SearchCursor(firstInterOut, ['OBJECTID'])
    for row in cursor:
        dic[row[0]] = [row[1], row[2]]
    del cursor
    del row
    print ('search')

    with arcpy.da.UpdateCursor(firstInterOut, ['OBJECTID', 'SHAPE@XY','LONGITUDE','LATITUDE']) as cursor:
        for row in cursor:
            # 更改点的坐标
            row[1] = dic[row[0]]
            # 更新点的属性表中的经纬度数据，经纬度数据采用平面投影坐标
            row[2] = dic[row[0]][0]
            row[3] = dic[row[0]][1]
            cursor.updateRow(row)
    print ('finish')

---

3.路网匹配

创建名为roadMatching的py文件
代码如下：

# encoding=utf-8

import arcpy
from updateGPS import updateGPSdata
from importData import csvToPoint


arcpy.env.workspace = 'D:/gisdata/test.gdb'

def doIntersect(roadBuffer, points, intersectOut):
    print ('Intersect')
    arcpy.Intersect_analysis([roadBuffer, points], intersectOut, 'ALL', '#', 'INPUT')

def getProperty(intersectOut):
    cursor = arcpy.SearchCursor(intersectOut)
    for row in cursor:
        print (row.getValue('NAME'))

def deleteShapefile(shapefile):
    arcpy.Delete_management(shapefile)
    print ('delete shapeFile over')


if __name__ == '__main__':
    # 要转化为要素类的csv表
    inpath = r'轨迹点的源文件路径' + 'gps.csv'
    road = 'road_network'
    pointName = 'gps'#要生成的点要素的名称
    csvToPoint(inpath, pointName)
    updateGPSdata(pointName, road)
    # 相交分析的结果
    interOutName = 'gps_inner'
    doIntersect(road, pointName, interOutName)
    print('inter is done')
    arcpy.ExportXYv_stats(interOutName,['CARNO', 'LONGITUDE', 'LATITUDE', 'POSITION_TIME', 'newid', 'NEAR_FID'], "COMMA",'输出路网匹配后轨迹点的文件路径' + 'gps_100.csv', "ADD_FIELD_NAMES")#NEAR_FID为所在路网的道路id

总结

用了一个比较简单的方法对路网进行匹配，并且将点要素的位置更新在路网