如何高效计算两个DataFrame中点对的距离（非O(n*m)算法）

针对15000个建筑、50000个医生的经纬度点对距离计算需求，直接用笛卡尔积会产生7.5亿次距离计算，效率极低。结合空间索引+批量大地距离计算的方法，可大幅减少计算量，同时支持可变半径筛选。

核心思路 #

1. 空间索引前置筛选：利用GeoDataFrame的R-tree空间索引，快速排除与目标建筑距离远超阈值的医生点，仅保留候选点对
2. 精确大地距离计算：对候选点对使用WGS84椭球参数计算精确 geodesic 距离，避免平面投影的误差
3. 批量处理优化：用批量查询替代逐行遍历，进一步提升大数据集下的处理速度

实现代码 #

方法1：批量空间索引查询+精确距离计算 #

import geopandas as gpd
import pandas as pd
from pyproj import Geod

# 测试数据
data_pdv = {'pdv': range(1, 6001),
            'latitude': [48.8566] * 3000 + [30.7128] * 3000,
            'longitude': [2.3522] * 6000}

data_pds = {'pds': range(1, 201),
            'latitude': [48.8588] * 200,
            'longitude': [2.2944] * 200}

# 转换为GeoDataFrame并设置WGS84坐标系
gdf_pdv = gpd.GeoDataFrame(
    data_pdv,
    geometry=gpd.points_from_xy(data_pdv['longitude'], data_pdv['latitude']),
    crs="EPSG:4326"
)
gdf_pds = gpd.GeoDataFrame(
    data_pds,
    geometry=gpd.points_from_xy(data_pds['longitude'], data_pds['latitude']),
    crs="EPSG:4326"
)

# 定义搜索半径（单位：米）
search_radius_m = 10000

# 初始化大地距离计算对象（基于WGS84椭球）
geod = Geod(ellps='WGS84')

# 创建医生GeoDataFrame的空间索引
pds_sindex = gdf_pds.sindex

# 计算地理坐标系下的近似搜索边界（1度≈111320米）
approx_degree_radius = search_radius_m / 111320

# 生成所有建筑点的搜索边界（minx, miny, maxx, maxy）
bboxes = gdf_pdv.geometry.apply(lambda geom: (
    geom.x - approx_degree_radius,
    geom.y - approx_degree_radius,
    geom.x + approx_degree_radius,
    geom.y + approx_degree_radius
)).tolist()

# 批量查询所有建筑点的候选医生索引
pdv_indices, pds_indices = pds_sindex.query_bulk(bboxes, predicate='intersects')

# 构建候选点对数据集
candidate_pairs = pd.DataFrame({
    'pdv_idx': pdv_indices,
    'pds_idx': pds_indices
})

# 关联建筑和医生的基础信息
candidate_pairs = candidate_pairs.merge(
    gdf_pdv[['pdv', 'latitude', 'longitude']],
    left_on='pdv_idx',
    right_index=True
).merge(
    gdf_pds[['pds', 'latitude', 'longitude']],
    left_on='pds_idx',
    right_index=True,
    suffixes=('_pdv', '_pds')
)

# 批量计算精确大地距离（单位：米）
_, _, distances = geod.inv(
    candidate_pairs['longitude_pdv'].tolist(),
    candidate_pairs['latitude_pdv'].tolist(),
    candidate_pairs['longitude_pds'].tolist(),
    candidate_pairs['latitude_pds'].tolist()
)
candidate_pairs['distance'] = distances

# 筛选符合距离要求的结果
output_df = candidate_pairs[candidate_pairs['distance'] < search_radius_m][['pdv', 'pds', 'distance']].reset_index(drop=True)

方法2：逐行遍历（适合小数据集调试） #

如果需要逐行验证逻辑，可使用以下更直观的版本：

import geopandas as gpd
import pandas as pd
from pyproj import Geod

# 数据初始化同方法1...

results = []
for idx, pdv_row in gdf_pdv.iterrows():
    # 生成当前建筑的搜索边界
    bbox = (
        pdv_row.geometry.x - approx_degree_radius,
        pdv_row.geometry.y - approx_degree_radius,
        pdv_row.geometry.x + approx_degree_radius,
        pdv_row.geometry.y + approx_degree_radius
    )
    # 获取候选医生索引
    candidate_pds_idx = list(pds_sindex.intersection(bbox))
    if not candidate_pds_idx:
        continue
    # 获取候选医生数据
    candidate_pds = gdf_pds.iloc[candidate_pds_idx]
    # 计算精确距离
    _, _, distances = geod.inv(
        [pdv_row.geometry.x]*len(candidate_pds),
        [pdv_row.geometry.y]*len(candidate_pds),
        candidate_pds.geometry.x.tolist(),
        candidate_pds.geometry.y.tolist()
    )
    # 筛选有效数据并添加建筑ID
    valid_pds = candidate_pds.copy()
    valid_pds['distance'] = distances
    valid_pds = valid_pds[valid_pds['distance'] < search_radius_m]
    valid_pds['pdv'] = pdv_row['pdv']
    results.append(valid_pds[['pdv', 'pds', 'distance']])

output_df = pd.concat(results, ignore_index=True)

效率对比 #

• 笛卡尔积方法：需计算6000200=120万次距离（测试数据），实际1500050000=7.5亿次，完全不可行
• 空间索引方法：测试数据中每个建筑的候选点仅约200个（因测试数据集中医生点集中），实际场景中候选点数量会远少于总医生数，计算量可降至百万级以内，效率提升数十倍甚至上百倍

注意事项 #

• 确保GeoDataFrame设置正确的CRS（EPSG:4326，即WGS84经纬度），否则空间索引和距离计算会出错
• 近似度数半径的计算是为了快速筛选，若需更精确的边界，可先将数据转换为以米为单位的投影坐标系（如EPSG:3857），再用buffer生成精确边界后做空间查询，但会增加投影转换的开销
• 可变半径需求下，仅需修改search_radius_m参数即可，无需重新构建GeoDataFrame，完美适配你的需求

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