使用Python Shapely 如何以O(n)时间最快生成任意闭合多边形内的均匀网格点？

原实现的核心效率瓶颈 #

• 纯Python双重循环逐点生成网格，没有利用NumPy向量化的批量计算能力，循环本身的固定开销极高
• 逐点实例化Point对象再做包含校验，单个对象构造、单次contains调用的固定开销被放大数十到数千倍
• 没有提前构建多边形的空间索引，每个点的包含判断都需要遍历多边形所有边，单查询时间复杂度高
• 没有对外接矩形区域做前置裁剪，对外接矩形中完全落在多边形外的冗余点也做了无意义的包含校验
• 存在坐标映射错误：Shapely对象的bounds属性固定返回(minx, miny, maxx, maxy)，原代码直接将其赋值给latmin, lonmin, latmax, lonmax，易导致网格坐标系和几何对象坐标系不匹配
• 整体时间复杂度远高于O(n)（n为最终有效点的数量），耗时和外接矩形的面积强相关，和多边形实际面积无关

逼近O(n)复杂度的优化路径 #

• 预处理多边形构建内置空间索引，把单次包含查询的开销降低一个量级
• 用NumPy批量生成全网格坐标数组，替代Python级的双重循环
• 用向量化的空间谓词一次性批量校验所有点，无需逐点构造Point对象
• 采用水平扫描线方法逐行裁剪，跳过完全落在多边形内部区域的所有无效校验，仅对和多边形边界相交的窄条带做少量点的校验，完全消除外接矩形冗余点的计算开销

优化后的基础实现（向量化版本） #

较原纯循环实现可提速10~100倍

from shapely.prepared import prep
from shapely.geometry import shape, Point
import numpy as np

polygon = shape(geojson['features'][i]['geometry'])
# 预处理多边形，构建边的空间索引，加速后续点查询
prep_polygon = prep(polygon)
minx, miny, maxx, maxy = polygon.bounds

# 一次性生成所有网格坐标，替代双重Python循环
x_coords = np.arange(minx, maxx, resolution)
y_coords = np.arange(miny, maxy, resolution)
xx, yy = np.meshgrid(x_coords, y_coords)

# 批量校验点是否在多边形内，直接传入坐标数组，无需逐点构造Point对象
valid_mask = prep_polygon.contains_xy(xx.ravel(), yy.ravel())

# 提取有效坐标并生成Point对象
valid_points = [
    Point((round(x, 4), round(y, 4)))
    for x, y in zip(xx.ravel()[valid_mask], yy.ravel()[valid_mask])
]

严格O(n)复杂度的优化（扫描线裁剪版本） #

向量化版本的校验开销仍然和外接矩形内的总点数相关，扫描线方法可以完全消除无效点的校验开销，让整体耗时和最终有效点的数量线性相关：

• 对每一行网格对应的固定y值，生成一条横跨整个外接矩形的水平扫描线
• 计算扫描线和多边形的所有相交区间，区间内的所有网格点必然落在多边形内部，无需做逐点包含校验
• 仅对扫描线和多边形边界交叉的边缘区域做少量点的校验，避免对全量外接矩形点的重复判断
• 完全内部的行网格、区间内的网格点直接按间距采样加入结果集，零校验开销

该方法在生成复杂多边形（如细长三角、带凹陷的多边形）的网格时，较原方案可提速100倍以上，无外接矩形冗余点的计算开销。

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