使用沿线网络的最短点连接两个geopandas df - python

Question

我正在使用

gpd.sjoin_nearest

加入两个 geopandas 数据框。这会返回最近的点，但我正在尝试计算输出的测量单位是什么？我已经完成了以公里为单位的近似计算，但我认为投影不起作用。

我还打算合并第三个 geopandas 数据框（

lines

），其中点之间的最近距离必须沿着这些线行进。所以要相互联系。

是否可以导入连接线网络（

roads

）并测量相同数据帧之间的最短距离，但距离必须通过连接线？

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, -20, -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()

从

t1

到

t2

和

P1_gpd

最近的点是

a

中的点

P2_gpd

。

我的目标是将距离转换为公里。我已经做了原始计算。第一个点的距离为1107公里，第二个点的距离为482公里。

预期输出：

  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

Answer 1

实际上，这不仅仅是测量单位的问题，这更是一个图表问题。您想要计算 P2（a 或 b）中最近点到 P1（t1 和 t2）沿道路（即不在自由笛卡尔平面中）的最近点之间的距离（单位为公里）。

如果是这样，我个人会使用momepy/networkx的原始方法。

我们首先使用

gdf_to_nx

构建图表：

PX_gpd = pd.concat([P1_gpd, P2_gpd])

spl_lines = set(
    split(l, p)
    for l in roads_gpd.clip(box(*PX_gpd.total_bounds)).geometry
    for p in PX_gpd.geometry
)

G = momepy.gdf_to_nx(
    gpd.GeoDataFrame(geometry=list(spl_lines), crs="EPSG:4326")
        .explode()
        .to_crs("EPSG:20354"),
    approach="primal",
    multigraph=False,
    length="nearest_distance",
)

cat_names = (
    PX_gpd.to_crs("EPSG:20354").pipe(
        lambda x: x.set_index(
            x.get_coordinates().agg(tuple, axis=1)
        )["Cat"].to_dict()
    )
)

nx.set_node_attributes(G, {n: {"name": a} for n, a in cat_names.items()})

然后，我们带回一个 GeoDataFrame 并仅保留 (

P1<=>P2

) 之间的最小距离：

ps, ls = momepy.nx_to_gdf(G)
node_cols = ["node_start", "node_end"]
mapper_ids = ps.set_index("nodeID")["name"]

ndis = (
    ls.drop(columns=node_cols)
    .join(
        ls[node_cols]
        .replace(mapper_ids)
        .apply(lambda r: pd.Series(sorted(r)), axis=1)
        .set_axis(node_cols, axis=1)
    )
    .eval("nearest_distance = nearest_distance / 1e3")
    .loc[ # comment this chain if you need to preserve self-P1-paths
        lambda x: ~x[node_cols].agg(
            lambda r: set(r).issubset(P1_gpd["Cat"]), axis=1
        )
    ]
    .pipe(
        lambda x: x.loc[
            x.groupby(
                x[node_cols].isin(P1_gpd["Cat"].tolist())
                    .all().idxmax()
            )["nearest_distance"].idxmin()
        ]
    )
)

输出：

                             geometry  nearest_distance node_start node_end
0  LINESTRING (395267.479 7788031....       1107.445824          a       t1
1  LINESTRING (403427.466 6680615....        482.440476          a       t2

NB：由于您的输入采用经度/纬度（即 EPSG:4326），因此您需要一个预计的 CRS。因此，请随意识别正确的投影 CRS，并记住，当平铺到地图上时，不同的投影 CRS 会以不同的方式扭曲地球表面。这会导致距离计算略有变化。例如，对于澳大利亚的特定区域（在您的 MRE 中），EPSG:20354 似乎提供了与您的期望非常匹配的距离计算。

绘图（参见完整代码）：