Python中基于GeoDataFrame裁剪NetCDF文件计算城市PM2.5均值

原代码无法提取所有城市数据的核心问题有4个：

• 范围取值错误：读取NetCDF边界时仅取了第一个城市（Rennes）的经纬度范围，没有计算所有城市的整体外接矩形，其余城市覆盖的网格根本没被纳入读取范围。
• 循环逻辑漏洞：月份循环写为range(12,1,-1)，永远不会读取1月数据；个位数月份未补前导零，会匹配错误文件名。
• 关键字段丢失：空间连接后输出时丢弃了城市名称字段，后续无法按城市维度统计平均值。
• 存在性能与兼容性问题：逐网格读取NetCDF变量效率极低，使用了已被新版GeoPandas弃用的参数和旧版CRS写法，容易触发警告或报错，且未关闭打开的NetCDF文件易造成句柄泄漏。

import os
import glob
import numpy as np
import pandas as pd
import geopandas as gpd
from netCDF4 import Dataset
from shapely.geometry import Point

# 读取城市边界数据
df = pd.read_excel('/Users/lucius/Documents/MiR/Data1/bounding_boxes.xls')
df['WKT'] = gpd.GeoSeries.from_wkt(df['WKT'])
gdf = gpd.GeoDataFrame(df, geometry='WKT', crs="EPSG:4326")

path = "/Users/lucius/Documents/MiR/Data1/"

for year in range(2020, 2000, -1):
    for month in range(12, 0, -1):
        month_str = f"{month:02d}"
        outputfile = os.path.join(path, f'Monthly_csv/Asia.{year}{month_str}.csv')
        if os.path.exists(outputfile):
            print(f"{year}-{month_str} already processed, skip")
            continue
        
        fn = os.path.join(path, f'Monthly/V5GL01.HybridPM25.Asia.{year}{month_str}-{year}{month_str}.nc')
        ds_temp = Dataset(fn, format="NETCDF4")
        lon_arr = ds_temp.variables['LON'][:]
        lat_arr = ds_temp.variables['LAT'][:]
        minlon, maxlon = lon_arr.min(), lon_arr.max()
        minlat, maxlat = lat_arr.min(), lat_arr.max()
        n_lon, n_lat = len(lon_arr), len(lat_arr)

        # 计算所有城市的总外接矩形范围
        total_minx, total_miny, total_maxx, total_maxy = gdf.total_bounds

        # 计算对应网格索引，增加越界保护
        i_min = int(np.floor((total_minx - minlon) / 0.01))
        i_max = int(np.ceil((total_maxx - minlon) / 0.01)) + 1
        i_min, i_max = max(0, i_min), min(n_lon, i_max)

        j_min = int(np.floor((maxlat - total_maxy) / 0.01))
        j_max = int(np.ceil((maxlat - total_miny) / 0.01)) + 1
        j_min, j_max = max(0, j_min), min(n_lat, j_max)

        # 一次性读取目标范围的PM2.5切片，提升读取效率
        pm25_slice = ds_temp.variables['PM25'][j_min:j_max, i_min:i_max]
        ds_temp.close()

        points, i_s, j_s, pol = [], [], [], []
        for idx_i, i in enumerate(range(i_min, i_max)):
            for idx_j, j in enumerate(range(j_min, j_max)):
                centro_x = minlon + 0.01 * i
                centro_y = maxlat - 0.01 * j
                points.append(Point(centro_x, centro_y))
                i_s.append(i)
                j_s.append(j)
                pol.append(pm25_slice[idx_j, idx_i].tolist())

        # 构建网格中心点GeoDataFrame
        pnts = gpd.GeoDataFrame(
            {'geometry': points, 'i': i_s, 'j': j_s, 'pol': pol},
            crs="EPSG:4326"
        )
        pnts['lon'] = pnts.geometry.x
        pnts['lat'] = pnts.geometry.y

        # 空间连接匹配网格与所属城市
        pnts_within = gpd.sjoin(pnts, gdf, how="inner", predicate='intersects')
        pnts_within['year'] = year
        pnts_within['month'] = month
        # 保留城市名称、时间、PM2.5值与经纬度字段
        pnts_within = pnts_within[['cities', 'year', 'month', 'pol', 'lon', 'lat']]
        pnts_within.to_csv(outputfile, mode='w', index=False)
        print(f"-> {year}-{month_str} processing completed")

所有月份处理完成后，运行以下代码即可直接得到每个城市逐月的PM2.5平均值：

# 批量读取所有处理完成的网格文件
all_csv = glob.glob(os.path.join(path, "Monthly_csv/*.csv"))
df_all = pd.concat((pd.read_csv(f) for f in all_csv), ignore_index=True)

# 按城市、年、月分组计算平均值
city_pm25_mean = df_all.groupby(['cities', 'year', 'month'], as_index=False)['pol'].mean()
city_pm25_mean.to_csv(os.path.join(path, 'city_monthly_pm25_average.csv'), index=False)

• 代码保留了原逻辑中的纬度索引计算方式，如果之前Rennes的提取结果经纬度匹配正确，无需调整；若出现纬度偏移，检查NetCDF文件中LAT数组的排序方向，对应调整j值计算逻辑即可。
• 增加了索引越界保护，如果当前读取的是亚洲区域NetCDF，不在亚洲范围内的城市（如Rennes、Copenhagen）会被自动过滤，不会触发索引报错。
• 对1-9月自动补前导零，避免文件名匹配错误。

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