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

Question

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

我一直在用 rasterstats.zonal_stats 做这件事，它工作得很好，但对于大型栅格来说速度很慢。最近我发现 rioxarray's example 不是将栅格裁剪到多边形，然后计算平均值。从我的实验来看，裁剪方法比 zonal_stats 快很多，结果是一样的。

因此我想知道除了时差之外是否还有其他原因会导致偏爱一种方法而不是另一种方法？为什么剪辑比 zonal_stats 快得多？

下面是两种方法的输出和时序，以及一段代码。完整的代码可以在 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 ---")

Answer 1

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

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

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

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

Compute polygon mean from raster: rasterstats.zonal_stats versus rio.clip

python

mask

raster

shapefile

python-xarray