如何在 Python 中高效遍历矢量要素并逐个裁剪栅格影像

本文介绍使用 `rioxarray` 和 `geopandas` 对多边形矢量要素进行循环迭代，逐个裁剪栅格影像并提取统计信息的正确方法，解决因误用 `.loc` 导致的 `indexingerror: too many indexers` 错误。

在遥感与地理空间分析中，常需对栅格影像按多个矢量面（如行政区、地块、样方）分别裁剪并计算统计值（如均值、标准差、像元数量等）。虽然可将每个面单独保存为独立 shapefile 再调用 raster.rio.clip()，但面对成千上万的要素时，这种“文件拆分+重复读写”的方式效率极低且易出错。理想做法是直接在内存中遍历 GeoDataFrame 的每一行，并构造符合 rioxarray.clip() 要求的几何输入。

关键问题在于：raster.rio.clip() 接收的是 一个 GeoDataFrame（含 geometry 列）或包含 {‘type’: …, ‘coordinates’: …} 字典的列表，而非单个 shapely.geometry 对象。原代码中：

shp = row.loc[index, 'geometry']  # ❌ 返回单个 Polygon/MultiPolygon 对象 raster_clipped = raster.rio.clip(shp.geometry.apply(mapping))  # ❌ shp 无 .geometry 属性，且 mapping() 作用对象错误

不仅逻辑混乱（row.loc[index, …] 在 iterrows() 中本身已按行索引，无需再用 loc 取索引），还误将单个几何体当作 GeoDataFrame 处理，导致 IndexingError 或 AttributeError。

✅ 正确做法是：将每行构造成单行 GeoDataFrame，再调用 .geometry.apply(mapping) 得到 GeojsON-like 字典列表（即使只有一个要素，也必须是 list[dict]）：

立即学习“Python免费学习笔记（深入）”；

import geopandas as gpd import rioxarray as rxr from shapely.geometry import mapping  # 读取数据 shpfile = gpd.read_file('shapefile.shp') raster = rxr.open_rasterio('raster.tif')  # 遍历每个多边形要素并裁剪 for idx, row in shpfile.iterrows():     # ✅ 构造仅含当前要素的 GeoDataFrame     single_gdf = gpd.GeoDataFrame([row], geometry='geometry', crs=shpfile.crs)      # ✅ 传入 clip：要求 geometry 列存在，且返回 list of geojson dicts     try:         clipped = raster.rio.clip(single_gdf.geometry.apply(mapping), drop=True)          # 示例：计算该区域栅格均值（忽略 nodata）         mean_val = float(clipped.mean(dim=['x', 'y'], skipna=True).item())         print(f"Feature {idx}: mean = {mean_val:.4f}")      except ValueError as e:         print(f"Feature {idx} skipped — {e}")  # 如无重叠、全 nodata 等情况

⚠️ 注意事项：

• drop=True（默认）会移除被裁剪后全为 nodata 的像元，推荐保留；
• 若矢量与栅格坐标系不一致，务必先统一 CRS：shpfile = shpfile.to_crs(raster.rio.crs)；
• 对于超大栅格或海量要素，建议添加 raster.rio.nodata 显式指定无效值，并考虑使用 dask 延迟计算提升性能；
• rioxarray.clip() 要求输入几何必须与栅格 CRS 一致，否则会报 CRSError。

最终，该方法避免了磁盘 I/O 开销，内存可控，且结构清晰，适用于从几十到数十万个面要素的批量处理任务。