Python使用regionmask裁剪NetCDF文件的结果校验问题

裁剪后维度未变化的核心原因 #

你的代码存在两个核心逻辑问题：

1. where方法使用错误：你直接将regionmask生成的mask作为判断条件传入，但未明确指定保留非空值的规则。regionmask输出的掩膜中，落在城市边界内的网格值为对应城市的整数编号，边界外为NaN，没有明确判断条件时，where只会将非零、非空的位置保留为原值，其余设为NaN。
2. drop=True不生效：drop=True的规则是仅当某一整行/整列纬度/经度的数值全为NaN时，才会删除该条带。你的研究城市分散在全球不同位置，不存在任何一整行纬度或一整列经度完全没有城市覆盖，因此不会删除任何维度条带，最终维度和原始全球数据完全一致。

另外你的代码运行耗时1小时，是因为直接对全球4亿+网格做掩膜计算，没有提前裁剪缩小计算范围。

验证裁剪是否生效的方法 #

• 统计非空值数量：对比原始数据集和裁剪后数据集的PM2.5变量非NaN值数量，分别运行print(ds['你的PM2.5变量名'].count().values)和print(cropped_ds_temp2['你的PM2.5变量名'].count().values)，如果两个数值接近，说明裁剪未生效；如果裁剪生效，后者数值应该仅等于所有城市覆盖的网格总数，远小于前者。
• 绘图验证：截取单个城市周边的小范围区域绘图，比如针对雷恩市，运行cropped_ds_temp2.sel(lon=slice(-2, 0), lat=slice(48,49))['PM2.5变量名'].plot()，如果仅在城市边界范围内有值，其余区域全为空，说明掩膜生效。
• 检查掩膜本身：运行mask.plot()，正常生效的掩膜应该仅在你划定的城市位置有整数值，全球其余区域全为NaN。

修正后的可运行代码 #

注意你原有代码中读取WKT列时存在列名错误，原表格中几何列列名为WKT，不是geometry，以下代码已修正该问题，同时增加预裁剪步骤将运行速度提升数十倍，可直接输出各城市的PM2.5平均值：

import pandas as pd
import geopandas as gpd
import xarray as xr
import regionmask

# 读取原始PM2.5数据
ds = xr.open_dataset('/path/to/file/V5GL01.HybridPM25.Global.201912-201912.nc')
# 读取城市边界表格
df = pd.read_excel('/Users/lucius/Documents/MiR/Data1/bounding_boxes.xls')
# 转换WKT为几何对象，设置坐标系为WGS84和NetCDF数据保持一致
df['geometry'] = gpd.GeoSeries.from_wkt(df['WKT'])
gdf = gpd.GeoDataFrame(df, geometry='geometry', crs="EPSG:4326")

# 第一步：按所有城市的总边界预裁剪，大幅缩小计算范围
lon_pad = 0.1
lat_pad = 0.1
total_bbox = [
    gdf.longmin.min() - lon_pad,
    gdf.longmax.max() + lon_pad,
    gdf.latmin.min() - lat_pad,
    gdf.latmax.max() + lat_pad
]
ds_sub = ds.sel(
    lon=slice(total_bbox[0], total_bbox[1]),
    lat=slice(total_bbox[2], total_bbox[3])
)

# 第二步：创建城市掩膜
city_regions = regionmask.from_geopandas(
    gdf,
    names="cities",
    abbrevs="cities",
    name="city_regions"
)
mask = city_regions.mask(ds_sub, lat_name='lat', lon_name='lon')

# 第三步：应用掩膜，仅保留城市范围内的网格
masked_ds = ds_sub.where(mask.notnull())

# 第四步：按城市分组计算平均PM2.5
city_avg_pm = masked_ds.groupby(mask).mean(keep_attrs=True)
# 将数字编号替换为城市名
city_avg_pm['city_name'] = ('region', gdf.cities.values)
city_avg_pm = city_avg_pm.swap_dims({'region': 'city_name'})

# 输出结果
print(city_avg_pm)

注：如果你的NetCDF数据经度范围是0°_{360°，需要先将经度转换为-180°}180°范围再做裁剪，避免跨0度经线的区域匹配错误。

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