如何用shapefile掩膜NetCDF并计算各多边形的均值与距平

以下是可直接运行的完整代码和流程说明：

完整可运行代码 #

# Load libraries
%matplotlib inline

import regionmask
import numpy as np
import xarray as xr
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt

# ---------------------- 1. 数据读取 ----------------------
# 读取美国州级shapefile
us = gpd.read_file('./shp/state_cus.shp')
# 统一转为WGS84投影，确保和NetCDF格点数据投影匹配
us = us.to_crs(epsg=4326)

# 读取逐月气温格点数据
ds = xr.open_mfdataset('./temp/monthly_mean_t2m_*.nc')
# 若原始气温单位为开尔文，执行以下代码转为摄氏度，变量名根据实际数据调整
# ds['t2m'] = ds['t2m'] - 273.15

# ---------------------- 2. 生成州级掩膜 ----------------------
# 生成每个格点对应的州标签，state_abbr为shapefile中存储州缩写的字段，根据实际属性调整
state_mask = regionmask.mask_geopandas(us, ds.lon, ds.lat, wrap_lon=False)
ds['state'] = state_mask

# ---------------------- 3. 计算2000-2019年各州平均气温 ----------------------
# 先对时间维度求平均得到20年格点均值，再按州分组求区域平均
ds_20y_avg = ds.t2m.mean(dim='time').groupby(ds.state).mean()
# 生成初始结果表
result_df = pd.DataFrame({
    '州': us['state_abbr'].values,
    '平均气温（avg_temp）': ds_20y_avg.values
})

# ---------------------- 4. 计算2019年气温距平 ----------------------
# 计算2000-2010年基准期各州平均气温
baseline_avg = ds.sel(time=slice('2000-01-01', '2010-12-31')).t2m.mean(dim='time').groupby(ds.state).mean()
# 计算2019年各州平均气温
temp_2019 = ds.sel(time=slice('2019-01-01', '2019-12-31')).t2m.mean(dim='time').groupby(ds.state).mean()
# 距平=2019年均温 - 基准期均温
result_df['2019年气温距平（anom_temp2019）'] = (temp_2019 - baseline_avg).values

# ---------------------- 5. 结果输出 ----------------------
# 打印结果，也可通过to_csv方法导出为本地文件
print(result_df)

关键注意事项 #

• 投影对齐是掩膜生效的前提，若出现全空值先检查shapefile和NetCDF的坐标系是否一致
• 若格点数据的经纬度字段名为longitude/latitude，需要将代码中对应的ds.lon/ds.lat同步修改
• shapefile的州缩写字段名需和实际属性匹配，若字段名为其他值，替换代码中的state_abbr即可
• 时间筛选的范围需要和数据的时间格式匹配，部分仅标注年月的数据可省略日维度，直接写slice('2000-01', '2010-12')

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