如何高效将0-1网格数组转换为Shapely对象？

嘿，我之前处理过超大网格转Shapely对象的需求，逐个单元格合并确实慢到让人头大！给你两个亲测高效的方案，对付2000x2000的网格完全没问题：

方案一：连通域分组+轮廓提取 #

先把所有连续的1区域分组，再为每组生成单个多边形，避免逐个单元格拼接的低效操作：

1. 用scipy.ndimage.label标记所有连通的1区域，把分散的1单元格按连通性归类：

import numpy as np
from scipy.ndimage import label
from shapely.geometry import Polygon, MultiPolygon
from skimage.measure import find_contours

# 假设你的0-1网格数组是grid，shape为(2000,2000)
# grid = 你的实际数组
structure = np.ones((3,3), dtype=int)  # 这里是8连通规则，改成(2,2)可以用4连通
labeled_grid, region_count = label(grid, structure=structure)

2. 提取每个区域的轮廓并转成Shapely多边形：

# 提取所有1区域的轮廓，level=0.5刚好区分0和1的边界
contours = find_contours(grid, level=0.5)

polygons = []
for contour in contours:
    # 注意：find_contours返回的坐标是(y, x)，如果你的坐标系是(x,y)，记得交换顺序
    poly = Polygon(contour)
    # 偶尔会生成自相交的无效多边形，用buffer(0)快速修复
    if not poly.is_valid:
        poly = poly.buffer(0)
    polygons.append(poly)

# 最后合并成单个Polygon或MultiPolygon
final_shape = polygons[0] if len(polygons) == 1 else MultiPolygon(polygons)

方案二：用rasterio直接生成矢量形状（最推荐！） #

这个方法底层基于GDAL实现，速度快到飞起，处理2000x2000的网格基本秒出结果：

import numpy as np
import rasterio.features
from shapely.geometry import shape, MultiPolygon

# 生成所有1区域的矢量形状生成器，mask参数指定只处理值为1的单元格
shape_gen = rasterio.features.shapes(grid, mask=grid == 1)

polygons = []
for geom, _ in shape_gen:
    shapely_poly = shape(geom)
    # 修复可能的无效多边形
    if not shapely_poly.is_valid:
        shapely_poly = shapely_poly.buffer(0)
    polygons.append(shapely_poly)

# 合并结果
final_shape = polygons[0] if len(polygons) == 1 else MultiPolygon(polygons)

额外提醒 #

• 坐标系适配：如果你的每个单元格对应实际空间中的正方形（比如边长为10米），可以用仿射变换调整坐标：

from affine import Affine
from shapely.ops import transform

# 假设单元格边长为cell_size，左下角起始坐标为(x0, y0)
cell_size = 1
x0, y0 = 0, 0
# 构建变换矩阵，把网格索引转成实际坐标
affine_transform = Affine.translation(x0 - cell_size/2, y0 - cell_size/2) * Affine.scale(cell_size, cell_size)
# 对最终形状应用变换
final_shape = transform(affine_transform, final_shape)

• 无效多边形修复：buffer(0)是Shapely里修复自相交、无效多边形的小技巧，几乎能解决90%的无效问题。

备注：内容来源于stack exchange，提问作者none none