Shapely中如何查找MultiPolygon内包含指定Point的Polygon

Shapely的MultiPolygon由多个互不相交的Polygon构成（仅允许单点接触的边界场景），直接调用polys.contains(pt)仅能返回点是否被包含的布尔值，该方法底层确实基于R树实现了高效查询，但不会直接返回命中的子Polygon。不需要全量遍历polys.geoms逐一检测的高效实现思路如下：

核心方案：复用Shapely内置STRtree空间索引 #

Shapely自带的排序瓦片R树（STRtree）和MultiPolygon.contains底层用的空间索引逻辑一致，提前构建一次索引即可支持多次高效点查询，时间复杂度为O(logn)，远优于全量遍历的O(n)，尤其适合子Polygon数量大的场景。
示例代码：

from shapely.strtree import STRtree

# 索引仅需构建1次，后续多次查询可直接复用，不要每次查询都重新建索引
poly_list = list(polys.geoms)
spatial_index = STRtree(poly_list)

# 第一步：用R树快速筛出候选多边形，过滤掉绝大多数不可能命中的要素
candidate_polys = spatial_index.query(pt)
# 第二步：对少量候选做精确包含判断，排除仅边界接触的情况（可按需选择contains/intersects/within等判断逻辑）
target_polygon = next((poly for poly in candidate_polys if poly.contains(pt)), None)

注意事项 #

• 由于MultiPolygon的子多边形互不重叠，点落在多边形内部时，精确判断后最多只会得到1个命中结果；如果点刚好落在两个子多边形的公共交点上，可能返回多个结果，可根据业务需求自行处理边界点归属
• Shapely 1.x版本中STRtree的query方法返回的是要素索引id，需要自行维护id到Polygon对象的映射表；2.0+版本直接返回几何对象，无需额外映射
• 当子Polygon总数小于100时，全量遍历和走索引的性能差异几乎可以忽略，不需要额外引入索引逻辑；如果是批量点查询、或者子Polygon数量过千，索引的性能优势会非常明显
• 不存在直接从MultiPolygon.contains的执行过程中直接取命中子Polygon的方法：Shapely是GEOS库的Python绑定，底层GEOS执行空间判断时不会把子要素的引用暴露到上层Python接口，无法直接截获中间计算结果

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