Python中快速求取两个多线要素图层相交点的高性能方法

问题背景 #

现有两个线要素图层：layer A（黄色，约23000条线要素）、layer B（蓝色，约50000条线要素），需要提取两个图层线要素的所有交点（即图示红色节点），已掌握两个图层所有节点坐标：

当前已知的逐对调用shapely相交方法的实现计算量过高，示例代码如下：

from shapely.geometry import LineString

line1 = LineString([(x1,y1), (x2,y2), (x3,y3)])
line2 = LineString([(x4,y4), (x5,y5)])

output = line1.intersection(line2)

核心问题是暴力两两遍历的时间复杂度达O(n*m)，总计需要执行11.5亿次相交判断，绝大多数计算是无意义的。

可落地的提速方法 #

• 优先用空间索引做预筛选，过滤99%以上的无效计算
  直接使用shapely内置的STRtree空间索引，先给要素量更大的图层构建边界框索引，对每一条待计算的线，只从索引中取出边界框和它重叠的候选要素，再对候选要素做精确相交判断即可，不需要和另一图层的所有线做计算。
  参考实现：

  from shapely.geometry import LineString
  from shapely.strtree import STRtree
  
  # lines_a、lines_b分别为layer A、layer B的LineString对象列表
  # 给要素更多的layer B构建空间索引
  index_b = STRtree(lines_b)
  intersect_nodes = []
  
  for line_a in lines_a:
      # 仅查询空间上可能相交的候选线
      candidates = index_b.query(line_a)
      for line_b in candidates:
          res = line_a.intersection(line_b)
          # 筛选点类型的交点结果
          if res.geom_type == "Point":
              intersect_nodes.append(res)
          elif res.geom_type == "MultiPoint":
              intersect_nodes.extend(list(res.geoms))
  

• 升级依赖库版本，调用矢量化接口
  如果你用的是2.0以下版本的shapely，直接升级到shapely 2.0+，底层整合了pygeos的矢量化计算能力，单条相交判断的速度能提升5-10倍；也可以直接用geopandas的空间连接接口sjoin，底层已经封装好了空间索引和批量运算逻辑，不需要手动写循环，代码更简洁，速度也更快。
• 提前做数据预处理
  计算前先去除图层内的重复线、冗余节点，对精度要求不高的场景可以适当简化线要素的节点数，也能进一步压缩计算量。

不要尝试手写纯Python的线段相交判断逻辑替换shapely的intersection方法：shapely的相交计算是调用底层GEOS的C语言实现，运算效率远高于纯Python代码，速度慢的核心原因是对大量空间上完全不相邻的线对做了无意义计算，而非intersection方法本身性能差。

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