把地震目录按县界摊开：GeoPandas 空间连接的一次踩坑

手上有一份 4.2 万条的区域地震目录（CSV，字段是 time、lat、lon、mag、depth），想按县界统计每个县的历史最大震级和事件数，用来判断哪几个县在应急演练里该重点排。环境是 macOS + Python 3.11、pandas 2.2、geopandas 0.14、shapely 2.0，县界是省测绘部门公开的 GeoJSON，约 2800 个多边形。

第一步直接 gpd.sjoin，结果跑了 6 分钟还没完，而且内存冲到 3G 多。查下来是坐标参考不一致：地震目录是 WGS84 经纬度（EPSG:4326），县界 GeoJSON 没写 CRS，geopandas 默认当成无投影。先统一：

import geopandas as gpd
import pandas as pd

pts = pd.read_csv("catalog.csv")
gdf = gpd.GeoDataFrame(
    pts,
    geometry=gpd.points_from_xy(pts.lon, pts.lat),
    crs="EPSG:4326",
)

counties = gpd.read_file("counties.geojson")
counties = counties.set_crs("EPSG:4326", allow_override=True)
counties = counties.to_crs("EPSG:3857")
gdf = gdf.to_crs("EPSG:3857")

counties["geometry"] = counties.buffer(0)  # 修自相交
joined = gpd.sjoin(gdf, counties[["name", "geometry"]], predicate="within", how="left")

关键改动有三个：投影到 3857 让空间索引有意义；buffer(0) 修掉公开数据里常见的自相交多边形，否则 within 会直接抛 TopologyException；how="left" 保留落在边界外或海上的点，不悄悄丢数据。

换完之后 4.2 万点 × 2800 面大约 8 秒。代价是 3857 在高纬度形变明显，做距离量算不能用它；我这里只做「点在哪个面内」，形变不影响归属。

会失效的地方：一是跨县界的大地震，震中点在哪个县不代表影响范围，统计最大震级时最好用震中距做缓冲而不是行政区；二是县界数据如果本身有缝隙，within 会把点判成无归属，how="left" 之后要单独数一下 name 为空的行，我这边有 37 条，基本都在海上；三是 shapely 2.0 下 buffer(0) 对某些极细长多边形会返回空几何，得先 is_valid 再筛。

另外提醒一句，sjoin 结果里 index_right 是县界行号，聚合时别直接 groupby("name")，同名县在不同省很常见，我一开始就合错了两个。

你们做行政区统计时，边界上的点按 within 还是 intersects 处理？我倾向 within 但会漏掉正好落在边界线上的事件。