如何基于经纬度对距离合并两个25万行DataFrame（仅保留50米内配对）

1. 数据准备与命名 #

先给两个DataFrame分配不同名称，避免变量覆盖，确保后续操作清晰：

import pandas as pd
import numpy as np

# 第一个DataFrame（记为df1）
latitude1 = [40.7128, 37.7749, 34.0522]
longitude1 = [-74.0060, -122.4194, -118.2437]
index1 = [1, 2, 3]
df1 = pd.DataFrame({'latitude': latitude1, 'longitude': longitude1, 'index': index1})

# 第二个DataFrame（记为df2）
latitude2 = [40.7324, 37.7569, 34.0612]
longitude2 = [-74.0070, -122.4284, -118.2507]
index2 = [1, 2, 3]
df2 = pd.DataFrame({'latitude': latitude2, 'longitude': longitude2, 'index': index2})

2. 按Index合并DataFrame #

通过merge方法按index字段关联两个DataFrame，自动给右侧重复列名添加_x后缀，正好符合你需要的latitude_x、longitude_x列名：

merged_df = pd.merge(df1, df2, on='index', suffixes=('', '_x'))

3. 高效计算球面距离 #

使用Haversine公式的向量化实现（基于numpy），处理25万行数据时效率远高于循环或apply：

def haversine(lat1, lon1, lat2, lon2):
    # 经纬度转换为弧度
    lat1_rad = np.radians(lat1)
    lon1_rad = np.radians(lon1)
    lat2_rad = np.radians(lat2)
    lon2_rad = np.radians(lon2)
    
    # Haversine核心计算
    dlat = lat2_rad - lat1_rad
    dlon = lon2_rad - lon1_rad
    a = np.sin(dlat/2)**2 + np.cos(lat1_rad) * np.cos(lat2_rad) * np.sin(dlon/2)**2
    c = 2 * np.arcsin(np.sqrt(a))
    
    # 地球平均半径（单位：米）
    earth_radius = 6371000
    return c * earth_radius

# 批量计算距离
merged_df['distance'] = haversine(
    merged_df['latitude'], merged_df['longitude'],
    merged_df['latitude_x'], merged_df['longitude_x']
)

4. 过滤结果并保留指定列 #

筛选距离≤50米的记录，同时保留你需要的6列：

result_df = merged_df[merged_df['distance'] <= 50][
    ['index', 'latitude', 'longitude', 'latitude_x', 'longitude_x', 'distance']
]

补充说明 #

如果你的需求不是按index一对一匹配，而是要从两个DataFrame中找出所有距离≤50米的经纬度对（不限制index），直接笛卡尔积运算会产生625亿条数据，完全不可行。此时需要用空间索引优化，借助geopandas的R-tree索引实现近邻查询：

import geopandas as gpd
from shapely.geometry import Point

# 转换为GeoDataFrame
gdf1 = gpd.GeoDataFrame(df1, geometry=gpd.points_from_xy(df1.longitude, df1.latitude), crs='EPSG:4326')
gdf2 = gpd.GeoDataFrame(df2, geometry=gpd.points_from_xy(df2.longitude, df2.latitude), crs='EPSG:4326')

# 转换为米制坐标系（避免球面距离计算误差）
gdf1 = gdf1.to_crs('EPSG:3857')
gdf2 = gdf2.to_crs('EPSG:3857')

# 近邻查询，找出50米内的配对
joined = gpd.sjoin_nearest(gdf1, gdf2, max_distance=50, distance_col='distance')
# 转换回原始坐标系并整理列名
joined = joined.to_crs('EPSG:4326')
result_df = joined.rename(columns={'latitude_right': 'latitude_x', 'longitude_right': 'longitude_x'})[
    ['index_left', 'latitude', 'longitude', 'latitude_x', 'longitude_x', 'distance']
].rename(columns={'index_left': 'index'})

这种方法利用空间索引大幅减少计算量，能高效处理25万行的规模。

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