GeoPandas sjoin_nearest计算结果异常及匹配逻辑疑问求助

你遇到的这两个问题其实挺典型的，咱们一步步来拆解解决：

一、距离计算差异的核心原因：EPSG:3857的坐标失真 #

你把坐标转成EPSG:3857（Web墨卡托）来计算米级距离，这在中高纬度地区会产生严重的距离拉伸失真！Web墨卡托投影是为网页地图显示设计的，纬度越高，横向距离的拉伸误差就越大，你用到的北纬43度区域，这种失真已经非常明显了——这就是为什么GeoPandas算出187m，和haversine/geodesic的830+m结果差这么多的根本原因。

正确的距离计算方式： #

别再用3857做距离计算了，推荐两种靠谱的方案：

1. 直接基于WGS84（EPSG:4326）计算球面距离
   GeoPandas在pyproj≥2.6.0的支持下，能直接处理地理坐标系的球面距离计算，再通过pyproj的Geod工具转成米：

   # 保留原始EPSG:4326的GeoDataFrame，不用转3857
   from pyproj import Geod
   geod = Geod(ellps="WGS84")  # 用WGS84椭球计算
   
   # 以你的测试点为例
   lon1, lat1 = -80.242342, 43.139322
   lon2, lat2 = -80.2440, 43.1391
   _, _, distance_m = geod.inv(lon1, lat1, lon2, lat2)
   print(f"准确球面距离：{distance_m:.3f}米")
   

2. 使用对应区域的UTM投影
   你数据所在的北纬43度、西经80度属于UTM 17N带（EPSG:32617），转成这个等距投影后，计算的距离就是准确的米级数值，不会有3857的拉伸问题：

   # 转成UTM 17N
   D1df_utm = D1df.to_crs("EPSG:32617")
   D2df_utm = D2df.to_crs("EPSG:32617")
   # 此时用distance计算的就是真实米级距离
   

二、sjoin_nearest返回多个匹配结果的问题 #

在geopandas 0.12.2版本中，sjoin_nearest默认k=1（只取1个最近邻），但如果存在多个点和目标点的距离完全相等，就会返回所有这些匹配行。

你可以通过以下方式确保每个D1的行只返回一个最近点：

1. 显式指定k=1强制只取单个最近邻：

   onetwo = gp.sjoin_nearest(D1df_utm, D2df_utm, distance_col="distancemeters", how="left", k=1)
   

2. 如果还是有重复结果，说明存在等距点，可按距离排序后去重：

   # 按距离升序排序，保留每组D1行的第一个结果
   onetwo = onetwo.sort_values("distancemeters").drop_duplicates(subset=D1df.index.name, keep="first")
   

快速验证方法 #

你可以用单个测试点验证正确的距离计算：

# 创建测试点的GeoDataFrame
test_d1 = gp.GeoDataFrame(
    data={"Latitude": [43.139322], "Longitude": [-80.242342]},
    geometry=gp.points_from_xy([-80.242342], [43.139322], crs="EPSG:4326")
).to_crs("EPSG:32617")

test_d2 = gp.GeoDataFrame(
    data={"Latitude": [43.1391], "Longitude": [-80.2440]},
    geometry=gp.points_from_xy([-80.2440], [43.1391], crs="EPSG:4326")
).to_crs("EPSG:32617")

# 计算距离
distance = test_d1.geometry.distance(test_d2.geometry).iloc[0]
print(f"UTM投影下的准确距离：{distance:.3f}米")

这个结果会和haversine/geodesic的计算值几乎一致。

备注：内容来源于stack exchange，提问作者Bunny