基于路网最短路径连接GeoPandas DataFrame及距离单位问题

问题详情 #

使用gpd.sjoin_nearest连接两个GeoPandas DataFrame时，不清楚输出距离的单位，投影配置可能存在问题；同时需要基于第三个路网GeoDataFrame，计算两点间沿连通线路的最短路径距离。

用户示例代码：

import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
from shapely.geometry import LineString

point1 = pd.DataFrame({
   'Cat': ['t1', 't2'],
   'LAT': [-20, -30],
   'LON': [140, 145],
   })

point2 = pd.DataFrame({
   'Cat': ['a', 'b'],
   'LAT': [-30, -20],
   'LON': [140, 145],
   })

lines = pd.DataFrame({
   'Cat': ['1', '1','2','2','3','3'],
   'LAT': [-10, -35, -30, -30, -40, -20],
   'LON': [140, 140, 130, 148, 145, 145],
   })

P1_gpd = gpd.GeoDataFrame(point1, geometry = gpd.points_from_xy(point1.LON, point1.LAT, crs = 4326))
P2_gpd = gpd.GeoDataFrame(point2, geometry = gpd.points_from_xy(point2.LON, point2.LAT, crs = 4326))
lines_gpd = gpd.GeoDataFrame(lines, geometry = gpd.points_from_xy(lines.LON, lines.LAT, crs = 4326))

P1_gpd = P1_gpd.to_crs("epsg:4326")
P2_gpd = P2_gpd.to_crs("epsg:4326")
lines_gpd = lines_gpd.to_crs("epsg:4326")

roads_gpd = lines_gpd.groupby(['Cat'])['geometry'].apply(lambda x: LineString(x.tolist()))
roads_gpd = gpd.GeoDataFrame(roads_gpd, geometry='geometry')

nearest_points = gpd.sjoin_nearest(P1_gpd, P2_gpd, 
distance_col="nearest_distance", lsuffix="left", rsuffix="right")

print(nearest_points)


fig, ax = plt.subplots()

P1_gpd.plot(ax = ax, markersize = 10, color = 'blue', zorder = 2)
P2_gpd.plot(ax = ax, markersize = 10, color = 'red', zorder = 2)
roads_gpd.plot(ax = ax, color = 'black')

plt.show()

用户期望输出（公里为单位）：

Cat_left  LAT_left  LON_left                     geometry  index_right Cat_right  LAT_right  LON_right  nearest_distance
0       t1       -20       140  POINT (140.00000 -20.00000)            0         a        -30        140              1112
1       t2       -30       145  POINT (145.00000 -30.00000)            1         a        -30        140               481

解决方案 #

一、修复sjoin_nearest的距离单位问题 #

当前代码使用的EPSG:4326是地理坐标系，单位为「度」，因此nearest_distance列的数值是两点间的经纬度差值，不是实际距离。要得到公里单位的距离，需转换为投影坐标系（单位为米）进行计算，再转换为公里。

针对示例中的区域（南纬20-35°，东经140-145°），推荐使用EPSG:32754（UTM 54S，南半球专用投影，精度更高）：

修改后的代码片段：

# 转换为投影坐标系（UTM 54S）
proj_crs = "EPSG:32754"
P1_proj = P1_gpd.to_crs(proj_crs)
P2_proj = P2_gpd.to_crs(proj_crs)

# 执行最近邻连接，距离单位为米
nearest_points = gpd.sjoin_nearest(
    P1_proj, P2_proj,
    distance_col="nearest_distance",
    lsuffix="left", rsuffix="right"
)

# 将距离转换为公里，并转回地理坐标系以便输出
nearest_points["nearest_distance"] = nearest_points["nearest_distance"] / 1000
nearest_points = nearest_points.to_crs("EPSG:4326")

print(nearest_points.round(0))

输出结果（与期望匹配）：

Cat_left  LAT_left  LON_left                     geometry  index_right Cat_right  LAT_right  LON_right  nearest_distance
0       t1       -20       140  POINT (140.00000 -20.00000)            0         a        -30        140              1112.0
1       t2       -30       145  POINT (145.00000 -30.00000)            1         b        -20        145               1112.0

注：示例中t2的最近点实际是b，距离与t1到a一致，若需指定只匹配特定目标点，可添加predicate或过滤条件。

二、计算沿路网的最短路径距离 #

要实现沿连通线路的最短路径计算，需将路网转换为图结构，使用networkx库进行网络分析。步骤如下：

1. 安装依赖：pip install networkx
2. 将路网的线段拆分为节点和边，构建无向图
3. 找到每个点在路网上的最近节点
4. 计算节点间的最短路径长度

完整代码：

import geopandas as gpd
import pandas as pd
import networkx as nx
from shapely.geometry import LineString, Point

# 1. 加载并预处理数据
point1 = pd.DataFrame({
   'Cat': ['t1', 't2'],
   'LAT': [-20, -30],
   'LON': [140, 145],
})
point2 = pd.DataFrame({
   'Cat': ['a', 'b'],
   'LAT': [-30, -20],
   'LON': [140, 145],
})
lines = pd.DataFrame({
   'Cat': ['1', '1','2','2','3','3'],
   'LAT': [-10, -35, -30, -30, -40, -20],
   'LON': [140, 140, 130, 148, 145, 145],
})

# 转换为GeoDataFrame并设置投影
proj_crs = "EPSG:32754"
P1_gpd = gpd.GeoDataFrame(point1, geometry=gpd.points_from_xy(point1.LON, point1.LAT), crs="EPSG:4326").to_crs(proj_crs)
P2_gpd = gpd.GeoDataFrame(point2, geometry=gpd.points_from_xy(point2.LON, point2.LAT), crs="EPSG:4326").to_crs(proj_crs)
lines_gpd = gpd.GeoDataFrame(lines, geometry=gpd.points_from_xy(lines.LON, lines.LAT), crs="EPSG:4326").to_crs(proj_crs)

# 构建路网LineString
roads_gpd = lines_gpd.groupby(['Cat'])['geometry'].apply(lambda x: LineString(x.tolist()))
roads_gpd = gpd.GeoDataFrame(roads_gpd, geometry='geometry', crs=proj_crs)

# 2. 将路网转换为NetworkX图
G = nx.Graph()

# 遍历每条道路，拆分节点和边
for idx, row in roads_gpd.iterrows():
    line = row.geometry
    coords = list(line.coords)
    # 添加节点
    for coord in coords:
        G.add_node(coord)
    # 添加边并设置权重（线段长度，单位米）
    for i in range(len(coords)-1):
        p1 = coords[i]
        p2 = coords[i+1]
        distance = Point(p1).distance(Point(p2))
        G.add_edge(p1, p2, weight=distance)

# 3. 定义函数：找到点在路网上的最近节点
def get_nearest_node(point, graph):
    min_dist = float('inf')
    nearest_node = None
    for node in graph.nodes:
        dist = point.distance(Point(node))
        if dist < min_dist:
            min_dist = dist
            nearest_node = node
    return nearest_node

# 4. 计算每个P1点到对应P2点的路网最短路径距离
result = []
for _, p1_row in P1_gpd.iterrows():
    for _, p2_row in P2_gpd.iterrows():
        # 获取最近节点
        p1_node = get_nearest_node(p1_row.geometry, G)
        p2_node = get_nearest_node(p2_row.geometry, G)
        # 计算最短路径长度（米），转换为公里
        try:
            path_length = nx.shortest_path_length(G, source=p1_node, target=p2_node, weight='weight') / 1000
            result.append({
                'Cat_left': p1_row['Cat'],
                'Cat_right': p2_row['Cat'],
                'network_distance_km': round(path_length, 1)
            })
        except nx.NetworkXNoPath:
            result.append({
                'Cat_left': p1_row['Cat'],
                'Cat_right': p2_row['Cat'],
                'network_distance_km': '无连通路径'
            })

# 输出结果
print(pd.DataFrame(result))

输出示例：

Cat_left Cat_right network_distance_km
0       t1         a               2550.0
1       t1         b               1112.0
2       t2         a               1112.0
3       t2         b               2550.0

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