如何解决rasterstats zonal_stats报错TypeError: invalid path or file？

使用Python的rasterstats模块zonal_stats函数计算shapefile中每个地块内栅格数据的平均太阳辐照度值时，执行以下代码：

print(irradiance.read().shape)

得到输出：(1, 2852, 2425)

随后执行：

with rioxarray.open_rasterio("C:/Users/Daniele/OneDrive/PhD/GIS + ESOM/GetSolarIrr/Solar_Irradiance/Yearly/Int_AreaSol_3.tif") as src:
    stats = zonal_stats(particles,src)

触发TypeError报错：

TypeError: invalid path or file: <xarray.DataArray (band: 1, y: 3000, x: 2000)>
[6000000 values with dtype=int32]
Coordinates:

• band (band) int32 1
• x (x) float64 7.6e+05 7.6e+05 7.6e+05 ... 7.8e+05 7.8e+05 7.8e+05
• y (y) float64 4.09e+06 4.09e+06 4.09e+06 ... 4.06e+06 4.06e+06
  spatial_ref int32 0
  Attributes: (12/15)
  AREA_OR_POINT: Area
  BandName: Band_1
  RepresentationType: ATHEMATIC
  STATISTICS_COVARIANCES: 8498387208.982427
  STATISTICS_MAXIMUM: 1552768
  STATISTICS_MEAN: 1286479.1959047
  ... ...
  STATISTICS_SKIPFACTORY: 1
  STATISTICS_STDDEV: 92186.69757065
  _FillValue: -1
  scale_factor: 1.0
  add_offset: 0.0
  long_name: Band_1

原以为是数据类型（float64/int32）导致问题，但不知如何解决，需修复该错误并了解计算shapefile中多个多边形内像素均值的简便方法。

1. 修复TypeError错误 #

报错核心原因是zonal_stats不接受xarray.DataArray作为栅格输入，它仅支持栅格文件路径字符串、rasterio的DatasetReader对象，或是numpy数组+仿射变换参数。以下是三种可行修复方式：

方式一：直接传入栅格文件路径（最简便） #

无需用rioxarray提前打开栅格，直接将文件路径传给zonal_stats，函数会自动用rasterio加载文件：

stats = zonal_stats(particles, "C:/Users/Daniele/OneDrive/PhD/GIS + ESOM/GetSolarIrr/Solar_Irradiance/Yearly/Int_AreaSol_3.tif")

方式二：从rioxarray对象中提取rasterio Dataset #

如果必须用rioxarray打开栅格（比如需要提前预处理），可以提取其底层的rasterio Dataset对象传入：

with rioxarray.open_rasterio("C:/Users/Daniele/OneDrive/PhD/GIS + ESOM/GetSolarIrr/Solar_Irradiance/Yearly/Int_AreaSol_3.tif") as src:
    stats = zonal_stats(particles, src.rio.dataset)

方式三：传入numpy数组+仿射变换 #

提取栅格的numpy数据（注意单波段需挤压band维度）和仿射变换参数，一起传给函数：

with rioxarray.open_rasterio("C:/Users/Daniele/OneDrive/PhD/GIS + ESOM/GetSolarIrr/Solar_Irradiance/Yearly/Int_AreaSol_3.tif") as src:
    # 挤压单波段维度，得到2D numpy数组
    raster_data = src.squeeze().values
    # 获取栅格的仿射变换信息
    affine_transform = src.rio.transform()
    stats = zonal_stats(particles, raster_data, affine=affine_transform)

2. 计算多边形内像素均值的简便方法 #

除了rasterstats.zonal_stats，还有两种常用简便方案：

方案一：用rasterstats直接指定统计量 #

调用zonal_stats时通过stats参数指定只计算均值，结果会直接返回每个多边形的均值：

# 仅计算均值，结果列表每个元素是包含'mean'键的字典
stats = zonal_stats(particles, "C:/Users/Daniele/OneDrive/PhD/GIS + ESOM/GetSolarIrr/Solar_Irradiance/Yearly/Int_AreaSol_3.tif", stats='mean')

# 若用geopandas存储地块数据，可直接将均值合并到属性表
import geopandas as gpd
gdf = gpd.read_file(particles)
gdf['mean_irradiance'] = [stat['mean'] for stat in stats]

方案二：geopandas + rioxarray 裁剪计算 #

利用rioxarray的rio.clip方法对每个多边形裁剪栅格，再计算均值，适合需要可视化裁剪结果的场景：

import geopandas as gpd
import rioxarray

# 读取地块shapefile
gdf = gpd.read_file(particles)
# 读取栅格并挤压单波段维度
src = rioxarray.open_rasterio("C:/Users/Daniele/OneDrive/PhD/GIS + ESOM/GetSolarIrr/Solar_Irradiance/Yearly/Int_AreaSol_3.tif").squeeze()

# 遍历每个多边形，裁剪后计算均值并写入属性表
gdf['mean_irradiance'] = gdf.apply(lambda row: src.rio.clip([row.geometry]).mean().item(), axis=1)

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