处理千万级GPS出租车数据时Python代码内存崩溃问题求助

针对5000万行GPS出租车数据匹配双向OSM路网时的内存崩溃问题，从匹配逻辑、数据处理流程、计算效率三个维度给出具体优化方案：

1. 缩小候选匹配范围，减少冗余数据 #

sjoin_nearest默认返回所有满足距离条件的链路，会导致中间结果数翻倍甚至更多。通过参数限制直接减少后续处理的数据量：

• 设置max_distance：根据实际场景（比如城市道路取50米）过滤过远链路，避免无效匹配。
• 限制k值：仅返回前2个最近候选（足够对比同链路双向或邻近少数链路）。

示例代码：

import geopandas as gpd

# 加载并预处理OSM路网（仅执行一次）
osm_roads = gpd.read_file("osm_bidirectional_roads.shp").to_crs("EPSG:3857")
# 优化数据类型：用更小的数值类型降低内存占用
osm_roads = osm_roads.astype({"osm_id": "int32", "direction": "int8"})

# 执行带限制的最近邻匹配
matched = gpd.sjoin_nearest(
    gps_points,
    osm_roads,
    max_distance=50,  # 仅匹配50米内的链路
    k=2,              # 仅保留前2个最近候选
    how="left",
    distance_col="match_dist"
)

2. 分块处理+即时释放内存，避免数据累积 #

采用“读取-处理-写入-释放”的流水线模式，每处理完一个块立即清理中间变量，防止内存逐步累积：

示例代码：

import pandas as pd
import numpy as np
import gc

# 路网仅加载一次
osm_roads = gpd.read_file("osm_bidirectional_roads.shp").to_crs("EPSG:3857")
osm_roads = osm_roads.astype({"osm_id": "int32", "direction": "int8"})

# 预计算路网链路的方向角（避免重复计算）
def calc_road_angle(line):
    start, end = line.coords[0], line.coords[-1]
    return np.arctan2(end[1] - start[1], end[0] - start[0])

osm_roads["road_angle"] = osm_roads["geometry"].apply(calc_road_angle).astype("float32")

# 分块读取GPS数据（每块10万行，可根据内存调整）
chunk_size = 100000
output_path = "final_matched_data.csv"

# 写入表头（仅第一次执行）
pd.DataFrame(columns=["gps_id", "lon", "lat", "speed", "final_osm_id", "final_direction"]).to_csv(
    output_path, index=False
)

for chunk_idx, chunk in enumerate(pd.read_csv("taxi_gps_data.csv", chunksize=chunk_size)):
    # 转换为GeoDataFrame并统一投影
    gps_chunk = gpd.GeoDataFrame(
        chunk,
        geometry=gpd.points_from_xy(chunk.lon, chunk.lat),
        crs="EPSG:4326"
    ).to_crs(osm_roads.crs)
    
    # 优化GPS数据类型
    gps_chunk = gps_chunk.astype({
        "lon": "float32",
        "lat": "float32",
        "speed": "float32",
        "gps_id": "int32"
    })
    
    # 计算GPS点的行驶方向角（假设数据已按taxi_id和timestamp排序）
    gps_chunk["prev_lon"] = gps_chunk.groupby("taxi_id")["lon"].shift(1)
    gps_chunk["prev_lat"] = gps_chunk.groupby("taxi_id")["lat"].shift(1)
    # 仅对有前序点的记录计算角度
    mask = gps_chunk["prev_lon"].notna()
    gps_chunk.loc[mask, "gps_angle"] = np.arctan2(
        gps_chunk.loc[mask, "lat"] - gps_chunk.loc[mask, "prev_lat"],
        gps_chunk.loc[mask, "lon"] - gps_chunk.loc[mask, "prev_lon"]
    ).astype("float32")
    
    # 执行匹配
    matched_chunk = gpd.sjoin_nearest(
        gps_chunk,
        osm_roads[["osm_id", "direction", "road_angle", "geometry"]],
        max_distance=50,
        k=2,
        how="left"
    )
    
    # 向量化筛选最优匹配（避免逐行循环）
    def get_best_match(group):
        if len(group) == 1:
            group["final_osm_id"] = group["osm_id"]
            group["final_direction"] = group["direction"]
            return group.drop_duplicates("gps_id")
        # 计算角度差（考虑圆周特性，取最小绝对值）
        angle_diff = np.abs(np.arctan2(
            np.sin(group["gps_angle"] - group["road_angle"]),
            np.cos(group["gps_angle"] - group["road_angle"])
        ))
        best_idx = angle_diff.argmin()
        best_row = group.iloc[best_idx]
        group["final_osm_id"] = best_row["osm_id"]
        group["final_direction"] = best_row["direction"]
        return group.drop_duplicates("gps_id")
    
    # 按GPS点分组处理
    final_chunk = matched_chunk.groupby("gps_id").apply(get_best_match).reset_index(drop=True)
    
    # 仅保留需要的列并写入文件
    final_chunk[["gps_id", "lon", "lat", "speed", "final_osm_id", "final_direction"]].to_csv(
        output_path, mode="a", header=False, index=False
    )
    
    # 强制释放内存
    del gps_chunk, matched_chunk, final_chunk
    gc.collect()
    print(f"完成第 {chunk_idx+1} 块处理")

3. 向量化替代循环，提升效率并降低内存 #

原角度差计算如果用逐行循环，会产生大量内存碎片且效率低下。改用numpy向量化操作，直接对分组后的数组计算，避免临时变量堆积：

核心角度差计算逻辑：

# 向量化计算角度差（处理圆周循环问题）
angle_diff = np.abs(np.arctan2(np.sin(gps_angle - road_angle), np.cos(gps_angle - road_angle)))
best_idx = angle_diff.argmin()

4. 用Dask-Geopandas实现自动并行分块 #

如果手动分块仍有压力，可使用Dask-Geopandas自动处理分块和并行计算，无需手动管理内存：

示例代码：

import dask_geopandas as dg

# 用Dask加载GPS数据（自动分块）
gps_dask = dg.read_file("taxi_gps_data.shp", chunksize=100000).to_crs("EPSG:3857")
osm_roads = gpd.read_file("osm_bidirectional_roads.shp").to_crs("EPSG:3857")
osm_roads["road_angle"] = osm_roads["geometry"].apply(calc_road_angle).astype("float32")

# 执行Dask版最近邻匹配
matched_dask = dg.sjoin_nearest(
    gps_dask,
    osm_roads,
    max_distance=50,
    k=2,
    how="left"
)

# 定义每个分区的处理函数
def process_partition(partition):
    # 计算GPS行驶角度
    partition["prev_lon"] = partition.groupby("taxi_id")["lon"].shift(1)
    partition["prev_lat"] = partition.groupby("taxi_id")["lat"].shift(1)
    mask = partition["prev_lon"].notna()
    partition.loc[mask, "gps_angle"] = np.arctan2(
        partition.loc[mask, "lat"] - partition.loc[mask, "prev_lat"],
        partition.loc[mask, "lon"] - partition.loc[mask, "prev_lon"]
    )
    # 筛选最优匹配
    def best_match(group):
        if len(group) == 1:
            return group.assign(final_osm_id=group["osm_id"], final_direction=group["direction"])
        angle_diff = np.abs(np.arctan2(np.sin(group["gps_angle"] - group["road_angle"]), np.cos(group["gps_angle"] - group["road_angle"])))
        best_idx = angle_diff.argmin()
        best_row = group.iloc[best_idx]
        return group.assign(final_osm_id=best_row["osm_id"], final_direction=best_row["direction"]).drop_duplicates("gps_id")
    
    return partition.groupby("gps_id").apply(best_match).reset_index(drop=True)[
        ["gps_id", "lon", "lat", "speed", "final_osm_id", "final_direction"]
    ]

# 应用函数并保存结果
result_dask = matched_dask.map_partitions(process_partition)
result_dask.to_csv("dask_matched_result.csv", single_file=True)

额外优化点 #

• 裁剪数据：用OSM路网的边界裁剪GPS点，仅处理路网覆盖范围内的点，减少无效计算。
• 精简列：读取GPS数据时仅加载gps_id、lon、lat、speed、taxi_id、timestamp等必需列，降低内存占用。
• 替换几何库：用PyGeos替代Shapely处理几何计算，可大幅降低内存占用并提升计算速度。

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