Geopandas缓冲区异常：订阅区域内事件无法被正确检测

问题根源 #

你的问题出在使用了EPSG:3857（Web墨卡托投影）进行距离相关计算。Web墨卡托是为Web地图显示设计的投影，它在非赤道区域会产生明显的距离拉伸失真——纬度越高，失真越严重。你计算的两点实际地面距离分别是48662米和37417米，但用EPSG:3857生成的50000米、40000米缓冲区，实际覆盖的地面距离会小于预期值，导致事件点没有被包含在缓冲区内，最终返回空DataFrame。

解决方案 #

有两种可靠的方法解决这个问题，优先推荐第一种（更简洁）：

方案1：直接在地理坐标系（EPSG:4326）下生成米级缓冲区 #

GeoPandas 0.10+版本支持直接在WGS84（EPSG:4326）地理坐标系下，通过指定crs参数生成基于实际地面距离的缓冲区，无需转换到投影坐标系：

import pandas as pd
import geopandas as gpd

# 创建事件点GeoDataFrame
df_strike = pd.DataFrame(
    {'Latitude': [27.0779, 31.9974],
     'Longitude': [51.5144, 38.7078]})
gdf_events = gpd.GeoDataFrame(
    df_strike, 
    geometry=gpd.points_from_xy(df_strike.Longitude, df_strike.Latitude),
    crs='EPSG:4326'
)

# 创建订阅位置GeoDataFrame
SUB_LOCATION = pd.DataFrame(
    {'perimeter_id': [1370, 13858],
     'distance': [40.0, 50.0],
     'custom_lat': [31.6661, 26.6500],
     'custom_lon': [38.6635, 51.5700]})

gdf_locations = gpd.GeoDataFrame(
    SUB_LOCATION, 
    geometry=gpd.points_from_xy(SUB_LOCATION.custom_lon, SUB_LOCATION.custom_lat),
    crs='EPSG:4326'
)

# 直接在EPSG:4326下生成基于米的缓冲区
gdf_locations['geometry'] = gdf_locations.geometry.buffer(
    gdf_locations['distance'] * 1000, 
    crs='EPSG:4326'
)

# 执行空间连接
matching_entln = gpd.sjoin(gdf_locations, gdf_events, how='inner')
print(matching_entln)

方案2：使用UTM等距投影生成缓冲区 #

UTM（通用横轴墨卡托）是专为小范围区域设计的等距投影，能保证距离计算的准确性。需要为每个位置匹配对应的UTM分带：

import pandas as pd
import geopandas as gpd

def get_utm_crs(lon, lat):
    # 计算UTM带号
    utm_band = int((lon + 180) / 6) + 1
    # 根据纬度判断南北半球，返回对应EPSG代码
    epsg_code = f'EPSG:326{utm_band:02d}' if lat >= 0 else f'EPSG:327{utm_band:02d}'
    return epsg_code

# 创建事件点和订阅位置GeoDataFrame（同方案1）
df_strike = pd.DataFrame(
    {'Latitude': [27.0779, 31.9974],
     'Longitude': [51.5144, 38.7078]})
gdf_events = gpd.GeoDataFrame(
    df_strike, 
    geometry=gpd.points_from_xy(df_strike.Longitude, df_strike.Latitude),
    crs='EPSG:4326'
)

SUB_LOCATION = pd.DataFrame(
    {'perimeter_id': [1370, 13858],
     'distance': [40.0, 50.0],
     'custom_lat': [31.6661, 26.6500],
     'custom_lon': [38.6635, 51.5700]})

gdf_locations = gpd.GeoDataFrame(
    SUB_LOCATION, 
    geometry=gpd.points_from_xy(SUB_LOCATION.custom_lon, SUB_LOCATION.custom_lat),
    crs='EPSG:4326'
)

# 为每个位置分配对应UTM坐标系
gdf_locations['utm_crs'] = gdf_locations.apply(
    lambda row: get_utm_crs(row.custom_lon, row.custom_lat), axis=1
)

# 逐个生成准确缓冲区
buffers = []
for _, row in gdf_locations.iterrows():
    # 转换到UTM投影
    loc_utm = gpd.GeoSeries([row.geometry], crs='EPSG:4326').to_crs(row['utm_crs'])
    # 生成米级缓冲区
    buffer_utm = loc_utm.buffer(row['distance'] * 1000)
    # 转回EPSG:4326以便后续空间连接
    buffers.append(buffer_utm.to_crs('EPSG:4326').iloc[0])

gdf_locations['geometry'] = buffers

# 执行空间连接
matching_entln = gpd.sjoin(gdf_locations, gdf_events, how='inner')
print(matching_entln)

验证结果 #

两种方案都能生成符合预期的缓冲区，空间连接后会返回你期望的结果：

perimeter_id  distance  custom_lat  custom_lon                                           geometry  index_right  Latitude  Longitude
0          1370      40.0     31.6661     38.6635  POLYGON ((38.71080 31.66610, 38.71076 31.6623...            1   31.9974    38.7078
1         13858      50.0     26.6500     51.5700  POLYGON ((51.62573 26.65000, 51.62566 26.6455...            0   27.0779    51.5144

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