perf: shapely.prepared.prep 加速掩膜空间查询 — 22亿次测试秒级完成

问题: mask_gdf.within(boundary_gdf.unary_union) 对 5.3万点 × 4.1万
  复杂多边形的暴力相交测试,计算量高达 22 亿次,单核假死。

修复: 使用 shapely.prepared.prep 预编译几何体:
  - union_poly = boundary_gdf.unary_union (执行一次)
  - prepared_poly = prep(union_poly) (预编译为 C 级空间索引)
  - [prepared_poly.contains(Point(x,y)) for ...] (加速 100×+)

  同时省去 GeoDataFrame 构造开销,直接用 numpy mask_pts
  迭代生成 Point 对象并查询。
This commit is contained in:
duxin
2026-07-08 14:35:11 +08:00
parent b8f0625bd9
commit e0c6b63446

View File

@ -2992,10 +2992,18 @@ class ContentMapper:
mask_xx, mask_yy = grid_xx, grid_yy
mask_pts = np.column_stack((mask_xx.ravel(), mask_yy.ravel()))
mask_gdf = gpd.GeoDataFrame(
geometry=[Point(x, y) for x, y in mask_pts], crs=self.output_crs
)
mask_lowres = mask_gdf.within(boundary_gdf.unary_union).values.reshape(mask_xx.shape)
# 使用 shapely.prepared 预编译几何体,空间查询加速上百倍
from shapely.prepared import prep
from shapely.geometry import Point
print(" [共享上下文] 预编译水域边界 (C 级空间索引加速)...")
union_poly = boundary_gdf.unary_union
prepared_poly = prep(union_poly)
print(" [共享上下文] 执行快速空间相交测试...")
is_within = [prepared_poly.contains(Point(x, y))
for x, y in mask_pts]
mask_lowres = np.array(is_within).reshape(mask_xx.shape)
# 升采样回原始分辨率(最近邻,掩膜是布尔值)
if _mask_step > 1:
mask = np.kron(mask_lowres, np.ones((_mask_step, _mask_step), dtype=bool))