本文介绍如何利用 rasterio.features.shapes 高效地将二值网格数组（如 2000×2000 的 0/1 矩阵）批量矢量化为 Shapely 几何对象，避免逐单元构造导致的性能瓶颈。
本文介绍如何利用 `rasterio.features.shapes` 高效地将二值网格数组（如 2000×2000 的 0/1 矩阵）批量矢量化为 shapely 几何对象，避免逐单元构造导致的性能瓶颈。
在地理空间数据处理中，常需将栅格化的占用地图（例如机器人导航中的障碍物栅格）转换为矢量几何以支持空间分析、可视化或与 GIS 工具交互。直接遍历每个值为 1 的像素并创建 shapely.geometry.Polygon 再调用 unary_union，时间复杂度高，对大型数组（如 2000×2000）极易成为性能瓶颈。
推荐方案是使用 rasterio.features.shapes —— 该函数基于底层 GDAL 实现，采用连通区域检测（4-或8-邻域）一次性提取所有同值区域的边界轮廓，效率远超纯 Python 循环。
✅ 基础实现（获取所有 1 区域的多边形）
import numpy as np
import rasterio.features
import shapely.geometry
# 示例：4×4 二值栅格
grid = np.array([
[0, 0, 0, 0],
[0, 1, 1, 0],
[0, 1, 1, 0],
[0, 0, 0, 0]
])
# 提取所有非零值（value=1）区域的几何形状
# mask 参数可指定仅处理 grid == 1 的位置（更精准）
shapes = list(rasterio.features.shapes(
grid,
mask=(grid == 1), # 关键：只对值为 1 的像元进行矢量化
connectivity=4 # 或设为 8，控制邻接方式（影响边界平滑度）
))
# 转换为 Shapely Polygon 对象列表（仅保留 value=1 的几何）
polygons = [
shapely.geometry.shape(geom)
for geom, value in shapes
if value == 1.0
]
# 若需合并为单一多边形（如所有障碍物视为一个整体）
from shapely.ops import unary_union
unioned_shape = unary_union(polygons)⚙️ 注意事项与优化建议
-
connectivity=4vs8：4仅考虑上下左右邻接，生成更“方正”的边界；8包含对角线邻接，边界更紧凑，但可能合并本应分离的小区域。 -
坐标系与分辨率：
rasterio.features.shapes输出的是像素坐标（行列索引）。若需真实地理坐标，须配合仿射变换（transform参数），常见于.tif文件读取场景；纯数组则默认(0,0)为左上角原点，单位为像素。 -
内存与精度：对超大数组（如 2000×2000），
shapes返回的 GeoJSON-like 几何可能包含大量顶点。可后续调用polygon.simplify(tolerance)进行 Douglas-Peucker 简化。 -
替代方案对比：
skimage.measure.find_contours适用于单连通轮廓提取，但不支持多区域自动分组；cv2.findContours功能强大但依赖 OpenCV，且输出格式需手动转 Shapely。
? 总结
rasterio.features.shapes 是将二值栅格批量转换为 Shapely 几何的最优实践——它底层高度优化、API 简洁，并天然支持掩膜过滤与邻域配置。结合 shapely.geometry.shape() 和 unary_union，即可在毫秒级完成数千像素级栅格的矢量化，显著优于手动构建 + 合并的低效模式。
