使用Rasterio Mask提取单多边形栅格时的重叠错误问题

嘿，这个场景我之前处理过！批量mask能成功但单个多边形报错，核心原因大概率是你选中的那个多边形其实和栅格范围没有空间交集——rasterio.mask在批量处理时会自动跳过不重叠的几何（只要有部分几何重叠就能生成结果），但单个几何完全不重叠的话就会直接抛出Input shapes do not overlap raster错误。当然也可能是几何本身有问题（比如自相交、空几何），所以我们需要先做预处理，再优雅地逐个提取。

下面是一套完整的解决方案，包含几何校验、空间过滤、错误处理，以及灵活的输出选项：

第一步：预处理几何列表 #

先过滤空几何，同时修复无效的多边形（比如自相交的情况），避免后续操作踩坑：

import rasterio
import fiona
import numpy as np
from shapely.geometry import shape, mapping
from shapely.validation import make_valid

# 读取Shapefile并预处理几何
with fiona.open("test.shp", "r") as shp_file:
    shp_crs = shp_file.crs
    cleaned_features = []
    for feat in shp_file:
        geom = feat["geometry"]
        # 跳过空几何
        if not geom:
            continue
        # 转换为shapely对象，修复无效几何
        shapely_geom = shape(geom)
        if not shapely_geom.is_valid:
            shapely_geom = make_valid(shapely_geom)
        # 转换回GeoJSON格式，方便rasterio使用
        cleaned_features.append(mapping(shapely_geom))

第二步：逐个提取多边形对应的栅格数据 #

我们写一个封装函数，先判断多边形和栅格是否相交，再执行mask操作，同时加入错误处理，确保整个流程不会因为单个多边形失败而中断：

def extract_single_polygon_raster(raster_path, polygon_geom):
    """提取单个多边形对应的栅格数组和空间变换参数"""
    with rasterio.open(raster_path) as src:
        # 获取栅格的边界范围，转为shapely多边形用于空间判断
        raster_bounds = src.bounds
        raster_polygon = shape({
            "type": "Polygon",
            "coordinates": [[
                (raster_bounds.left, raster_bounds.bottom),
                (raster_bounds.right, raster_bounds.bottom),
                (raster_bounds.right, raster_bounds.top),
                (raster_bounds.left, raster_bounds.top),
                (raster_bounds.left, raster_bounds.bottom)
            ]]
        })
        
        # 判断当前多边形是否和栅格有交集
        poly_shapely = shape(polygon_geom)
        if not poly_shapely.intersects(raster_polygon):
            print(f"⚠️ 多边形与栅格无重叠，跳过")
            return None, None
        
        # 执行mask操作，处理nodata值
        out_img, out_transform = rasterio.mask.mask(
            src,
            [polygon_geom],
            crop=True,
            nodata=src.nodata if src.nodata is not None else np.nan
        )
        # 单波段栅格可以移除多余的维度（从(1, h, w)转为(h, w)）
        if out_img.shape[0] == 1:
            out_img = out_img.squeeze()
        return out_img, out_transform

# 遍历所有预处理后的几何，逐个提取
raster_file = "sat_img_B01.jp2"
extracted_results = []
for idx, poly in enumerate(cleaned_features):
    print(f"正在处理第 {idx+1} 个多边形...")
    img_array, transform = extract_single_polygon_raster(raster_file, poly)
    if img_array is not None:
        extracted_results.append((img_array, transform))
        # 可选：保存为Numpy数组文件（方便后续数据分析）
        np.save(f"polygon_{idx}_raster.npy", img_array)
        # 可选：保存为GeoTIFF文件（保留空间坐标信息）
        # with rasterio.open(
        #     f"polygon_{idx}_output.tif",
        #     "w",
        #     driver="GTiff",
        #     height=img_array.shape[0],
        #     width=img_array.shape[1],
        #     count=1,
        #     dtype=img_array.dtype,
        #     crs=shp_crs,
        #     transform=transform,
        #     nodata=np.nan
        # ) as dst:
        #     dst.write(img_array, 1)

关键细节说明 #

• 几何有效性修复：用shapely.make_valid处理自相交、无效的多边形，避免mask操作时的隐性错误。
• 空间交集预判：提前过滤掉和栅格不重叠的多边形，直接跳过，从根源避免报错。
• 容错性处理：即使个别多边形处理失败，整个循环也会继续运行，保证批量任务的完整性。
• 灵活输出：既可以保存为Numpy数组（适合后续机器学习、数值分析），也可以保存为GeoTIFF（适合GIS软件查看）。

如果还是有个别多边形报错，可以单独打印该多边形的坐标，检查是否确实在栅格范围内，或者是否存在CRS的隐性问题（比如看起来CRS一致，但实际坐标单位不同？不过你说批量能成功，这个概率很低）。

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