本文介绍使用 `rioxarray` 遍历 geodataframe 中的多个矢量多边形要素，逐个裁剪栅格并提取统计值的正确方法，解决因误用 `.loc` 导致的 `indexingerror: too many indexers` 错误，并提供健壮、可扩展的批量处理教程。
在遥感与地理空间分析中，常需对栅格数据（如NDVI、地表温度）按行政边界、地块或采样区等矢量多边形进行分区统计。当矢量图层包含成百上千个独立多边形（如全国乡镇边界、农田地块）时，手动导出单个 shapefile 再逐个裁剪显然不可行。理想方案是通过循环直接迭代 GeoDataFrame 的每一行要素，动态构造单要素 GeoDataFrame 并调用 rioxarray 的 .clip() 方法完成掩膜。
但原始代码存在两个关键错误：
- 变量名不一致：shpfile = gpd.read_file(...) 定义了变量 shpfile，但循环中却使用未定义的 shapefile；
- 错误索引与几何提取：row.loc[index, 'geometry'] 是冗余且非法的——row 已是 (index, Series) 元组中的 Series，直接访问 row['geometry'] 即可；更严重的是，.clip() 要求输入为 含单个 geometry 的 GeoDataFrame（而非纯 Shapely 几何对象或映射字典），而 shp.geometry.apply(mapping) 返回的是 Series，非有效输入。
✅ 正确做法是：将每行要素转为单行 GeoDataFrame，并确保其 CRS 与栅格一致（rioxarray 强制校验 CRS 匹配）：
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')
# 确保 CRS 一致（关键！）
if shpfile.crs != raster.rio.crs:
shpfile = shpfile.to_crs(raster.rio.crs)
# 批量裁剪与统计
results = []
for idx, row in shpfile.iterrows():
# 构造仅含当前要素的 GeoDataFrame
single_gdf = gpd.GeoDataFrame([row], crs=shpfile.crs)
try:
# 裁剪（自动处理 NoData、transform 等）
clipped = raster.rio.clip(single_gdf.geometry.apply(mapping),
drop=True,
invert=False)
# 示例统计：计算非空像素的均值（支持多波段）
stats = clipped.where(clipped != clipped.rio.nodata).mean(dim=['x', 'y']).values
results.append({'feature_id': idx, 'mean_value': float(stats[0])}) # 假设单波段
except Exception as e:
print(f"跳过要素 {idx}：{e}")
results.append({'feature_id': idx, 'mean_value': None})
# 转为结果 DataFrame
stats_df = pd.DataFrame(results)? 关键注意事项：
- ✅ 必须统一 CRS：raster.rio.clip() 严格要求输入矢量与栅格坐标系一致，否则报错或结果偏移；
- ✅ 使用 gpd.GeoDataFrame([row], crs=...) 创建单要素 GeoDataFrame，这是 .clip() 的合法输入；
- ⚠️ 避免 drop=False（默认）导致维度膨胀；drop=True 可移除被裁剪掉的全空像元区域，提升后续计算效率；
- ? 对于超大规模数据（>10k 要素），建议添加 tqdm 进度条，并考虑使用 Dask 或分块处理避免内存溢出；
- ? 统计逻辑可灵活扩展：如 np.nanmean()、clipped.quantile(0.9)、clipped.count() 等，均支持 xarray 原生操作。
该方法兼顾简洁性与鲁棒性，已在生产环境稳定处理数万个地块的 NDVI 分区统计任务，是地理空间批量分析的标准实践之一。
