为何GeoPandas中LINESTRING与MULTIPOLYGON的Overlay交集结果为空？

我想要按健康区域拆分不列颠哥伦比亚省的道路网，具体来说，为每个健康区域保存该区域内道路的Shapefile。道路的几何类型为LINESTRING，而健康区域的几何类型为MULTIPOLYGON。

我原以为使用GeoPandas的overlay方法并设置how=intersection可以实现该需求，但结果显示交集为空。然而将道路与健康区域叠加可视化时，两者明显存在重叠，因此我怀疑自己对GeoPandas处理不同几何类型交集的方式存在误解。

最小复现代码 #

import geopandas as gpd
import matplotlib.pyplot as plt

health_regions_df = gpd.read_file(
    '../data/canada.geojson'
    )

bc_roads_df = gpd.read_file(
    '../data/geobase_nrn_rrn/bc/nrn_rrn_bc_shp_en.zip',
    layer='NRN_bc_14_0_ROADSEG'
    )

bc_roads_df.to_crs(
    'EPSG:4326', 
    inplace=True
)

bc_hrs = health_regions_df[health_regions_df['province'] == 'British Columbia']

intersection = gpd.overlay(bc_hrs, bc_roads_df, how='intersection')

print(intersection)

运行输出 #

Empty GeoDataFrame
Columns: [fid, hr_uid, province, engname, frename, count, NID, ROADSEGID, ADRANGENID, DATASETNAM, SPECVERS, ACCURACY, ACQTECH, PROVIDER, CREDATE, REVDATE, METACOVER, ROADCLASS, RTNUMBER1, RTNUMBER2, RTNUMBER3, RTNUMBER4, RTNUMBER5, RTENAME1FR, RTENAME2FR, RTENAME3FR, RTENAME4FR, RTENAME1EN, RTENAME2EN, RTENAME3EN, RTENAME4EN, EXITNBR, NBRLANES, PAVSTATUS, PAVSURF, UNPAVSURF, STRUCTID, STRUCTTYPE, STRUNAMEEN, STRUNAMEFR, L_ADDDIRFG, L_HNUMF, L_HNUML, L_STNAME_C, L_PLACENAM, R_ADDDIRFG, R_HNUMF, R_HNUML, R_STNAME_C, R_PLACENAM, CLOSING, ROADJURIS, SPEED, TRAFFICDIR, geometry]
Index: []

1. overlay方法的局限性 #

GeoPandas的overlay方法核心设计用于多边形与多边形的空间操作，尽管文档提及支持其他几何类型，但在处理线（LINESTRING）与多边形（MULTIPOLYGON）的交集时，对数据精度、拓扑完整性的要求极高，细微的CRS偏差或拓扑错误都可能导致返回空结果。

2. 正确的处理方式：使用clip或空间连接 #

针对“提取多边形范围内线要素”的需求，推荐使用以下两种更可靠的方法：

方法一：clip函数（最直接）

clip是GeoPandas专门用于提取指定多边形范围内几何要素的函数，完美支持线与多边形的组合：

# 强制对齐两个数据集的CRS
bc_hrs = bc_hrs.to_crs(bc_roads_df.crs)

# 按健康区域拆分并保存
for hr_name, hr_poly in bc_hrs.groupby('engname'):
    # 裁剪当前健康区域内的道路
    hr_roads = gpd.clip(bc_roads_df, hr_poly)
    # 保存为Shapefile（替换空格避免命名问题）
    hr_roads.to_file(f'../data/bc_roads_{hr_name.replace(" ", "_")}.shp')

方法二：空间连接+裁剪

先通过空间连接筛选出与健康区域相交的道路，再精准裁剪到对应区域：

# 空间连接，筛选出与BC省健康区域相交的道路
joined_roads = gpd.sjoin(bc_roads_df, bc_hrs, how='inner', predicate='intersects')

# 按健康区域分组保存
for hr_uid, road_group in joined_roads.groupby('hr_uid'):
    # 获取对应健康区域的多边形
    target_poly = bc_hrs[bc_hrs['hr_uid'] == hr_uid].geometry.iloc[0]
    # 裁剪道路到该多边形范围
    clipped_group = road_group[road_group.geometry.intersects(target_poly)]
    # 保存文件
    clipped_group.to_file(f'../data/bc_roads_hr_{hr_uid}.shp')

3. 额外排查项 #

• CRS一致性验证：通过print(bc_roads_df.crs)和print(bc_hrs.crs)确认两者CRS完全匹配，即使都是EPSG:4326，也要检查是否存在轴顺序等细微差异。
• 拓扑错误修复：如果原始数据存在拓扑问题，可通过buffer(0)修复：

  bc_hrs['geometry'] = bc_hrs['geometry'].buffer(0)
  bc_roads_df['geometry'] = bc_roads_df['geometry'].buffer(0)
  

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