网格中心点与复杂多边形/多多边形的包含检测性能优化求助

针对你遇到的性能瓶颈，我有几个基于shapely生态的优化方案，能大幅缩短运行时间，不用换库：

1. 用空间索引过滤候选多边形（最立竿见影的优化） #

你现在的问题是每个点都要遍历所有150个多边形，但绝大多数多边形和点根本不相交。用空间索引可以快速筛选出可能包含当前点的多边形，把每个点需要检查的多边形数量从150降到几个甚至1个。

Shapely自带STRtree（排序树）实现空间索引，代码示例：

from shapely.strtree import STRtree
from shapely.geometry import Point
import numpy as np

# 先给所有多边形构建空间索引
tree = STRtree(shapes)

# 把网格点批量转成Point对象（比嵌套循环逐个创建更高效）
# 假设x和y是二维数组，先扁平化
x_flat = x.flatten()
y_flat = y.flatten()
points = [Point(xi, yj) for xi, yj in zip(x_flat, y_flat)]

# 遍历每个点，先查候选多边形再做包含检测
for idx, p in enumerate(points):
    # 快速获取可能包含该点的多边形
    candidate_polygons = tree.query(p)
    for s in candidate_polygons:
        if s.contains(p):
            # 计算对应的网格i,j索引
            j = idx // nx
            i = idx % nx
            # 执行移除网格的逻辑
            break  # 找到匹配的多边形就停止，避免多余检查

2. 用GeoPandas做向量化批量运算 #

纯Python循环的开销很大，GeoPandas底层基于numpy和C扩展，能把空间运算批量处理，速度提升几十倍甚至上百倍。

代码示例：

import geopandas as gpd
from shapely.geometry import Point
import numpy as np

# 1. 把网格点转换成GeoSeries
x_flat = x.flatten()
y_flat = y.flatten()
points_geo = [Point(xi, yj) for xi, yj in zip(x_flat, y_flat)]
points_gs = gpd.GeoSeries(points_geo, crs="你的坐标系")  # 比如EPSG:4326

# 2. 把多边形转换成GeoDataFrame
polygons_gdf = gpd.GeoDataFrame({"geometry": shapes}, crs="你的坐标系")

# 3. 用空间连接快速找到所有在多边形内的点
# predicate='within' 表示点在多边形内部
joined = gpd.sjoin(points_gs.to_frame(name="geometry"), polygons_gdf, predicate="within")

# 4. 把结果映射回网格的i,j索引
# joined的index是原points_gs的索引，对应扁平化后的点顺序
for idx in joined.index:
    j = idx // nx
    i = idx % nx
    # 执行移除网格的逻辑

3. 简化复杂多边形 #

那个有64万个顶点的多边形肯定是性能杀手——顶点越多，contains检查的计算量越大。用shapely的simplify方法对多边形做简化，去掉冗余顶点，同时保留拓扑结构（比如孔洞）。

代码示例：

# tolerance是简化阈值，单位和你的地图坐标系一致，比如0.1米/度
# preserve_topology=True 确保孔洞、多边形连接关系不被破坏
simplified_shapes = [
    s.simplify(tolerance=0.1, preserve_topology=True) 
    for s in shapes
]
# 之后用简化后的simplified_shapes代替原shapes做检测

注意：需要根据你的业务精度要求调整tolerance，可以先测试不同值，找到速度和精度的平衡点。

4. 提前用边界框过滤网格 #

每个多边形都有自己的边界框（bbox），你可以先计算每个多边形的bbox，然后只检查落在这个bbox范围内的网格点，直接排除大量无关的网格。

比如：

for s in shapes:
    minx, miny, maxx, maxy = s.bounds
    # 找到x在[minx, maxx]，y在[miny, maxy]的网格索引
    # 假设x和y是网格中心点的数组，用numpy的布尔索引快速筛选
    mask_x = (x >= minx) & (x <= maxx)
    mask_y = (y >= miny) & (y <= maxy)
    # 遍历符合条件的网格点
    for j in np.where(mask_y)[0]:
        for i in np.where(mask_x)[0]:
            p = Point(x[i], y[j])
            if s.contains(p):
                # 执行移除逻辑

最优组合建议 #

优先尝试空间索引+GeoPandas向量化的组合，既能减少候选多边形数量，又能避免Python循环的开销，应该能把运行时间从几小时压缩到几分钟以内。如果还有超大多边形，再配合简化操作进一步提速。

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