Python实现：线段最近气象站匹配及网络图扩缩调整

一、高性能线段-最近气象站点匹配实现 #

原有嵌套循环方案时间复杂度为O(n*m)，数据量上涨后性能会快速下降，替换为空间索引+粗筛+精算的方案可将复杂度降至O((n+m)logm)，万级边+千级站点场景下性能可提升数十至上百倍。
核心依赖为geopandas、shapely、scipy、pandas、numpy，实现逻辑如下：

• 先将所有数据从WGS84经纬度坐标系转换为研究区对应的等距投影坐标系（单位为米，直接支持米级距离计算，避免球面距离计算的性能损耗）
• 基于站点坐标构建KDTree空间索引，对每条线段先用外接矩形做范围粗筛，仅对落在候选范围内的少量站点计算精确的点到线段距离，取1km邻域内距离最近的站点作为匹配结果

import geopandas as gpd
import pandas as pd
from shapely.geometry import Point, LineString
from scipy.spatial import KDTree
import numpy as np

# ---------------------- 1. 数据加载与投影转换 ----------------------
# 加载气象站点数据，替换为你的本地数据路径
stations_df = pd.read_csv("met_stations.csv")
stations_gdf = gpd.GeoDataFrame(
    stations_df,
    geometry=gpd.points_from_xy(stations_df.StationLongitude, stations_df.StationLatitude),
    crs="EPSG:4326" # WGS84经纬度坐标系
)
# 替换为研究区对应的UTM投影EPSG编码，例如中国中东部常用EPSG:32650，投影后坐标单位为米
proj_crs = "EPSG:32650"
stations_proj = stations_gdf.to_crs(proj_crs)

# 加载网络边数据，替换为你的本地数据路径
edges_df = pd.read_csv("network_edges.csv")
edges_gdf = gpd.GeoDataFrame(
    edges_df,
    geometry=edges_df.apply(lambda row: LineString([
        (row.StartLongitude, row.StartLatitude),
        (row.EndLongitude, row.EndLatitude)
    ]), axis=1),
    crs="EPSG:4326"
)
edges_proj = edges_gdf.to_crs(proj_crs)

# ---------------------- 2. KDTree构建与快速匹配 ----------------------
# 提取站点投影坐标构建空间索引
station_coords = np.array([(geom.x, geom.y) for geom in stations_proj.geometry])
station_tree = KDTree(station_coords)
SEARCH_RADIUS = 1000 # 1km邻域阈值，单位米

match_records = []
for edge_idx, edge_row in edges_proj.iterrows():
    line = edge_row.geometry
    # 用线段外接矩形做粗筛，快速圈定候选站点范围
    minx, miny, maxx, maxy = line.bounds
    search_centroid = [(minx+maxx)/2, (miny+maxy)/2]
    search_radius = SEARCH_RADIUS + max(maxx-minx, maxy-miny)/2
    candidate_idx = station_tree.query_ball_point(search_centroid, r=search_radius)
    
    nearest_station = None
    min_dist = float("inf")
    # 仅对候选站点计算精确点到线距离
    for s_idx in candidate_idx:
        station_geom = stations_proj.iloc[s_idx].geometry
        dist = line.distance(station_geom)
        if dist <= SEARCH_RADIUS and dist < min_dist:
            min_dist = dist
            nearest_station = stations_proj.iloc[s_idx]
    
    match_records.append({
        "edge_idx": edge_idx,
        "StartCity": edge_row.StartCity,
        "EndCity": edge_row.EndCity,
        "nearest_station_id": nearest_station.StationID if nearest_station is not None else None,
        "nearest_station_name": nearest_station.StationName if nearest_station is not None else None,
        "dist_to_station": min_dist if nearest_station is not None else None
    })
edge_station_match = pd.DataFrame(match_records)
initial_match_rate = edge_station_match["nearest_station_id"].notna().sum() / len(edge_station_match)

注意：禁止直接在WGS84经纬度坐标系下做米级距离计算，经纬度单位为角度，计算结果和实际距离偏差可达数倍到数十倍。如果研究区跨多个UTM带，可临时用EPSG:3857（Web墨卡托）做近似计算，精度足以满足1km邻域判断需求。

二、拓扑保持的线段调整方案（满足70%匹配率要求） #

调整核心原则：不改变边连接的起止城市对应关系（即网络图拓扑结构完全不变），仅沿线段原有走向等比例延长线段，最小化几何调整幅度。
逻辑说明：只要气象站点到线段所在直线的垂直距离≤1km，将线段沿走向延长至覆盖站点到直线的垂足位置，站点就会落入线段1km邻域范围内，且线段始终在原两城市连线上，不会改变网络拓扑。
实现时优先调整仅需最小延长幅度即可匹配的线段，避免无意义的大尺度改动：

# ---------------------- 3. 线段调整满足70%匹配率 ----------------------
TARGET_MATCH_RATE = 0.7
required_match_count = int(np.ceil(len(edges_proj) * TARGET_MATCH_RATE))
current_match_count = edge_station_match["nearest_station_id"].notna().sum()
final_edges = edges_proj.copy()

if current_match_count < required_match_count:
    unmatched_edges = edge_station_match[edge_station_match["nearest_station_id"].isna()]
    adjust_plan = []
    for _, match_row in unmatched_edges.iterrows():
        line = edges_proj.loc[match_row["edge_idx"], "geometry"]
        start_point = Point(line.coords[0])
        end_point = Point(line.coords[-1])
        original_len = line.length
        
        # 找离线段所在直线最近的站点
        min_perp_dist = float("inf")
        nearest_foot = None
        for _, s_row in stations_proj.iterrows():
            s_geom = s_row.geometry
            foot = line.line_interpolate_point(line.project(s_geom))
            perp_dist = s_geom.distance(foot)
            if perp_dist < min_perp_dist:
                min_perp_dist = perp_dist
                nearest_foot = foot
        
        # 仅保留垂直距离在1km内、可通过延长匹配的边
        if min_perp_dist <= SEARCH_RADIUS:
            dist_foot_start = nearest_foot.distance(start_point)
            dist_foot_end = nearest_foot.distance(end_point)
            extend_start = max(0, dist_foot_start - original_len/2)
            extend_end = max(0, dist_foot_end - original_len/2)
            # 计算最小延长比例，留1%余量避免边界误差，设置最多延长2倍的阈值
            extend_ratio = min(max(extend_start, extend_end)/(original_len/2) + 0.01, 2)
            adjust_plan.append({
                "edge_idx": match_row["edge_idx"],
                "extend_ratio": extend_ratio
            })
    
    # 按延长幅度从小到大排序，优先调整改动最小的边
    adjust_plan = sorted(adjust_plan, key=lambda x:x["extend_ratio"])
    for plan in adjust_plan:
        if current_match_count >= required_match_count:
            break
        idx = plan["edge_idx"]
        original_line = edges_proj.loc[idx, "geometry"]
        start = np.array(original_line.coords[0])
        end = np.array(original_line.coords[-1])
        direction = end - start
        dir_norm = direction / np.linalg.norm(direction)
        new_half_len = original_line.length/2 * (1 + plan["extend_ratio"])
        mid_point = (start + end)/2
        # 生成沿原方向延长的新线段
        new_start = mid_point - dir_norm * new_half_len
        new_end = mid_point + dir_norm * new_half_len
        final_edges.loc[idx, "geometry"] = LineString([new_start, new_end])
        current_match_count += 1

# 转回WGS84经纬度坐标系输出结果
final_edges_wgs84 = final_edges.to_crs("EPSG:4326")
edge_station_match.to_csv("edge_station_relation.csv", index=False, encoding="utf-8-sig")
final_edges_wgs84.to_file("adjusted_network.gpkg", driver="GPKG", encoding="utf-8")

如果数据量超过10万条边，可将KDTree替换为BallTree、或用dask_geopandas做并行空间计算进一步压缩运行时间。

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