从栅格计算多边形平均值：rasterstats.zonal_stats 与 rio.clip答案

【问题标题】：Compute polygon mean from raster: rasterstats.zonal_stats versus rio.clip从栅格计算多边形平均值：rasterstats.zonal_stats 与 rio.clip
【发布时间】：2021-10-14 00:01:10
【问题描述】：

我有一个光栅文件和一个包含多边形的 shapefile。对于每个多边形，我想计算该多边形内所有栅格像元的平均值。

我一直在使用rasterstats.zonal_stats 执行此操作，它运行良好，但对于大型栅格来说非常慢。最近我遇到了rioxarray's example，而不是将栅格裁剪为多边形，然后计算平均值。从我的实验来看，裁剪方法比 zonal_stats 快很多，结果是一样的。

因此，我想知道除了时差之外是否还有其他因素会导致一种方法优于另一种方法？以及为什么剪辑比zonal_stats 快这么多？

下面是两个方法的输出和时序，以及一个sn-p的代码。完整代码可以在here找到。

很高兴能对此有所了解:)

--- Raster stats ---
    ADM1_EN   mean_adm
0   Central  14.624690
1  Northern   3.950312
2  Southern  20.649534
--- Raster stats: 15.52 seconds ---

--- Rio clip ---
           mean_adm
Central   14.624689
Northern   3.950313
Southern  20.649534
--- Rio clip: 0.28 seconds ---

#load the data
ds=rioxarray.open_rasterio(chirps_path,masked=True)
ds=ds.rio.write_crs("EPSG:4326")
ds_date=ds.sel(time="2020-01-01").squeeze()

#use rasterstats
start_time = time.time()
gdf = gpd.read_file(hdx_adm1_path)[["ADM1_EN", "geometry"]]
gdf["mean_adm"] = pd.DataFrame(
            zonal_stats(
                vectors=gdf,
                raster=ds_date.values,
                affine=ds_date.rio.transform(),
                nodata=np.nan,
                all_touched=False
            )
        )["mean"]
print("--- Raster stats ---")
display(gdf[["ADM1_EN","mean_adm"]])
print(f"--- Raster stats: {(time.time() - start_time):.2f} seconds ---")

#use clipping
start_time = time.time()
gdf = gpd.read_file(hdx_adm1_path)[["ADM1_EN", "geometry"]]
df_adm=pd.DataFrame(index=gdf.ADM1_EN.unique())
for a in gdf.ADM1_EN.unique():
    gdf_adm=gdf[gdf.ADM1_EN==a]

    da_clip = ds_date.rio.set_spatial_dims(
        x_dim="x", y_dim="y"
    ).rio.clip(
        gdf_adm["geometry"], all_touched=False
    )

    grid_mean = da_clip.mean(dim=["x", "y"],skipna=True).rename("mean_adm")
    df_adm.loc[a,"mean_adm"]=grid_mean.values
print("--- Rio clip ---")
display(df_adm)
print(f"--- Rio clip: {(time.time() - start_time):.2f} seconds ---")

【问题讨论】：

标签： python raster mask shapefile python-xarray

【解决方案1】：

从自学的角度，我意识到我在这里犯了一个错误，我认为报告是有益的。

这完全取决于您运行函数的顺序。当我使用rioxarray.open_rasterio 加载栅格数据时，此数据不会加载到内存中。此后剪辑和zonal_stats 仍然必须将其加载到内存中。当两者中的一个在另一个之前被调用时，后者利用了数据已经加载的事实。

您也可以选择在调用函数之前调用ds.load()。在我的测试中，这并没有改变总计算时间。

【讨论】：