Python中基于Shapely的海量点路匹配向量化高效实现咨询

我在Python中有两个DataFrame：一个约15万条呼叫数据，每条都对应一个地理位置，另一个约5万条街道数据，每条都对应一条地理路径。目标是根据每个呼叫的位置，将最近街道的头节点、尾节点ID追加到呼叫DataFrame中。

我已完成基础数据转换，对应以下两个算法，同时自行实现了最近街道匹配逻辑，但性能达不到要求：

算法一：呼叫点坐标转Shapely Point #

% Algorithm One: given two columns of latitude & longitude, create a new Point 
def call_iter():
    points = []
    for index, row in calls.iterrows():
        points.append(Point(row['Incident Latitude'], row['Incident Longitude']))
    
    return points % appended to the call dataframe

算法二：街道字符串路径转Shapely LineString #

% Algorithm Two: given a string column containing coordinate data, construct a LineString 
def street_iter():     
    paths = []
    for geo in streets.geometry:
        l = []
        for t in geo.split():
            try:
                t = t.strip('(,)')
                l.append(float(t))
            except ValueError:
                pass
    
        p = []
        for i in range(0, len(l), 2):
            p.append(Point(l[i], l[i+1]))
    
        paths.append(LineString(p))
    return paths  % appended to the street dataframe

算法三：最近街道匹配 #

依托Shapely的line.distance(point)方法实现匹配，功能正常但单条呼叫处理耗时1-2秒，全量15万条数据处理耗时过高，无法满足批量处理需求：

% Algorithm Three: find the closest street (head 'u' and tail 'v' nodes) to each call
def build_matrix():
    heads = []
    tails = []

    for i_c, r_c in calls.iterrows():
        print(i)
        p = r_c[4]
        head_min = -1
        tail_min = -1
        dist_min = float('inf')
        min_group = []
    
        for i_s, r_s in streets.iterrows():
            l = r_s[5].distance(p)
            if dist_min > l:
                head_min = r_s['u'] % head node
                tail_min = r_s['v'] % tail node
                dist_min = l
                min_group = []
                min_group.append(r_s)
            if dist_min == l:
                min_group.append(r_s)
            
        if len(min_group) > 1:
            choice = secrets.choice(min_group) % randomly selects an arc
            head_min = choice['u']
            tail_min = choice['v']
        
        heads.append(head_min)
        tails.append(tail_min)

    return (heads, tails) % both appended to the calls dataframe

核心优化思路是用向量化操作替换Python层级的循环，用空间索引将匹配复杂度从O(MN)降到O(MlogN)，全流程如下：

1. 前置准备 #

提前安装geopandas、rtree两个库，用于空间操作和空间索引。

2. 基础数据转换优化（替换算法一、算法二） #

import geopandas as gpd
import re
from shapely.geometry import LineString
import numpy as np

# 算法一优化：批量生成呼叫点，无需逐行迭代
calls['geometry'] = gpd.points_from_xy(calls['Incident Longitude'], calls['Incident Latitude'])
# 转为GeoDataFrame，注意crs参数替换为你的实际坐标参考系，WGS84经纬度填EPSG:4326
calls_gdf = gpd.GeoDataFrame(calls, crs="EPSG:4326")

# 算法二优化：批量解析街道路径，替换多层循环
def parse_linestring(geo_str):
    # 正则批量提取所有数值，比逐字符处理效率高
    coords = list(map(float, re.findall(r'-?\d+\.?\d*', geo_str)))
    # 每两个数值组成一个坐标对
    point_pairs = list(zip(coords[::2], coords[1::2]))
    return LineString(point_pairs)

streets['geometry'] = streets['geometry'].apply(parse_linestring)
streets_gdf = gpd.GeoDataFrame(streets, crs=calls_gdf.crs)

3. 最近街道匹配优化（替换算法三） #

用geopandas内置的最近邻空间连接方法，底层自动使用空间索引筛选候选，无需手动嵌套循环：

# 给街道加随机排序字段，处理多条街道距离完全相同时随机选择的需求
streets_gdf['random_sort'] = np.random.rand(len(streets_gdf))

# 最近邻空间连接，每个呼叫点匹配最近的街道
matched = gpd.sjoin_nearest(
    calls_gdf, 
    streets_gdf[['u', 'v', 'geometry', 'random_sort']],
    how='left',
    distance_col='distance_to_street' # 可保留匹配距离用于后续校验
)

# 处理同距离的情况：按距离、随机序排序后取第一个
matched = matched.sort_values(['distance_to_street', 'random_sort']).groupby(level=0).first()

# 把匹配到的头、尾节点追加到原呼叫DataFrame
calls[['street_head_u', 'street_tail_v']] = matched[['u', 'v']]

性能说明 #

原嵌套循环处理15万条呼叫需要几十甚至上百小时，优化后全流程仅需几秒到几十秒，可满足多批量数据处理需求。
如果需要计算米级的物理距离，可将两个GeoDataFrame先投影到对应地区的平面坐标系后再执行匹配。

内容的提问来源于stack exchange，提问作者Jordan MacLachlan