如何用Python基于shapefile掩膜计算栅格数据的统计特征？

方案说明 #

最优实现采用rioxarray+geopandas的技术栈，底层基于GDAL做空间运算，比手动逐栅格判断掩膜的效率高1-2个数量级，尤其适合大尺寸栅格场景。

依赖安装 #

先安装需要的第三方库：

pip install rioxarray geopandas xarray numpy

完整实现代码 #

import numpy as np
import geopandas as gpd
import xarray as xr
import rioxarray
from pyproj import CRS

# 1. 读取矢量文件（和你原有代码逻辑一致）
plume = gpd.read_file(shapefile)

# 2. 将现有numpy栅格数据封装为带空间信息的xarray对象
# 注意维度顺序要对应：假设你的data形状是 (lat维度长度, lon维度长度)
da = xr.DataArray(
    data,
    dims = ["y", "x"],
    coords = {
        "y": lat[:, 0] if lat.ndim == 2 else lat, # 兼容2D/1D纬度数组
        "x": lon[0, :] if lon.ndim == 2 else lon  # 兼容2D/1D经度数组
    }
)

# 3. 给栅格指定坐标系，这里默认是WGS84（和你原有可视化用的PlateCarree投影对应），如果是其他投影修改对应EPSG编码即可
da = da.rio.write_crs(CRS.from_epsg(4326), inplace=False)

# 4. 坐标系校验与转换：保证矢量和栅格坐标系一致
if plume.crs != da.rio.crs:
    plume = plume.to_crs(da.rio.crs)

# 5. 用shapefile掩膜裁剪栅格，drop=True会直接删掉掩膜外的栅格点
clipped_da = da.rio.clip(plume.geometry.values, plume.crs, drop=True, invert=False)

# 6. 计算所需统计值，skipna=True会自动跳过无效值/NaN
mean_val = clipped_da.mean(skipna=True).item()
median_val = clipped_da.median(skipna=True).item()
std_val = clipped_da.std(skipna=True).item()

# 输出结果
print(f"掩膜内均值：{mean_val}")
print(f"掩膜内中位数：{median_val}")
print(f"掩膜内标准差：{std_val}")

可选效果校验 #

你可以用裁剪后的栅格绘图验证掩膜效果，和你原有可视化逻辑对齐：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

fig = plt.figure()
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
clipped_da.plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree(), shading='auto')
plume.boundary.plot(ax=ax, color='red')
plt.show()

注意事项 #

• 如果你的栅格是规则格网（绝大多数NetCDF栅格都满足），上述代码可以直接使用；如果是不规则格网，可改用geopandas.sjoin做空间关联后再统计
• 掩膜裁剪时如果需要保留原栅格的形状、仅把掩膜外的值设为NaN，把drop=True改为drop=False即可

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