如何用Python地理空间包模拟ArcGIS Intersect工具处理海量线要素GeoDataFrame

问题根源 #

ArcGIS的Intersect工具是提取两个线图层的空间交集片段，同时保留双方属性，且不会产出重叠要素。你用how='union'会保留两个图层的所有要素（包括不相交部分），这直接导致了重叠；而属性缺失是因为union后部分要素只来自单个图层，导出shapefile时NaN被自动转成0是shapefile的固有特性（不支持空值，用默认值填充）。

正确操作步骤 #

1. 拆分MULTILINESTRING #

先把多线拆成单线，保留原要素索引方便后续关联：

model_layer_explode = model_layer.explode(index_parts=True)

2. 用交集模式替代并集 #

ArcGIS的Intersect本质是求空间交集，所以用how='intersection'调用overlay：

info_layer_gdf = gpd.overlay(orig_layer, model_layer_explode, how='intersection', keep_geom_type=True)

keep_geom_type=True确保输出始终是线要素，避免意外生成点或面。

3. 解决属性缺失问题 #

如果还是有部分要素没拿到双方属性，大概率是因为线段完全重合但拓扑上没有交点（比如同一条线但分割点不同）。这时候先按model_layer的节点拆分orig_layer，再空间关联属性：

from shapely.ops import unary_union

# 提取model_layer所有节点作为拆分点
split_points = unary_union([geom for line in model_layer_explode.geometry for geom in line.coords])

# 拆分orig_layer的线段
split_lines = orig_layer.geometry.apply(lambda x: x.split(split_points))
split_orig = split_lines.explode(index_parts=True).reset_index(drop=True)
split_orig = gpd.GeoDataFrame(split_orig, columns=['geometry'], crs=orig_layer.crs)

# 空间关联双方属性
info_layer_gdf = gpd.sjoin(split_orig, orig_layer, how='left', predicate='within')
info_layer_gdf = gpd.sjoin(info_layer_gdf, model_layer_explode, how='left', predicate='within')

这样拆分后的每一段线都会匹配到包含它的orig_layer和model_layer要素属性，不会出现NaN。

4. 清理重叠要素 #

如果还有重叠，直接按几何和属性去重：

# 按几何及所有属性字段去重
info_layer_gdf = info_layer_gdf.drop_duplicates(
    subset=['geometry'] + list(orig_layer.columns) + list(model_layer.columns)
)

也可以先修复拓扑错误再去重：

info_layer_gdf['geometry'] = info_layer_gdf.geometry.buffer(0)

额外提示 #

• Shapefile不支持空值，导出时NaN会被转成0或空字符串，要保留空值建议用GeoJSON格式
• 提前用buffer(0)修复两个输入图层的拓扑错误（比如自相交），能减少GeoPandas处理时的异常

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