使用rasterio与geopandas时如何避免多边形被错误填充

核心原因 #

你遇到的孔洞被填充的问题，根源是忽略了rasterio.features.shapes()返回的嵌套坐标结构：

• 当处理带孔洞的1值区域时，shapes()返回的每个shape对象中，shape[0]["coordinates"]是一个列表：第一个元素是区域的外边界环，后续元素都是被1值包围的0值孔洞的内边界环。
• 你当前的代码只取了coordinates[0]（外边界）来创建Polygon，相当于只保留了外框，完全丢弃了孔洞的内边界信息，最终生成的多边形自然会填充原本的0值孔洞区域。

即使你通过shape[1] == 1筛选了1值区域，也无法解决这个问题——因为带孔洞的1值区域本身就是一个属性为1的连通shape，只是它的坐标结构包含了内外多层环，你没有正确解析这部分结构。

修正代码 #

要正确保留孔洞，需要把外边界和所有内边界（孔洞）都传入Polygon的构造函数：

def flooded_to_shapefile(array, transform, filename):
    shapes = rasterio.features.shapes(array, transform=transform)
    
    polygons = []
    for shape in shapes:
        if shape[1] == 1:
            coords = shape[0]["coordinates"]
            # 提取外边界和孔洞内边界
            exterior_ring = coords[0]
            interior_rings = coords[1:] if len(coords) > 1 else []
            # 创建带孔洞的多边形
            poly = shapely.geometry.Polygon(exterior_ring, interior_rings)
            polygons.append(poly)
    
    # 创建GeoDataFrame并保存
    df = gpd.GeoDataFrame({"geometry": polygons})
    df.set_geometry("geometry", inplace=True)
    df.to_file(filename)

补充说明 #

• shapely.geometry.Polygon的第二个参数接受一个内边界列表，每个内边界是一个坐标环，用于定义多边形中的孔洞。
• 确保你的shapely版本支持这种构造方式（大多数现代版本都支持）。

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