Python实现：将GeoDataFrame转为0.5°分辨率均值聚合二维网格

数据加载与样例 #

加载GeoDataFrame代码 #

gdf = gpd.GeoDataFrame(
    ds.to_pandas(),
    geometry=gpd.points_from_xy(ds["CENLON"], ds["CENLAT"]),
    crs="EPSG:4326",
)

数据样例 #

CENLON   CENLAT O1REGION O2REGION   AREA  ...  ZMAX  ZMED  SLOPE  
index                                               ...                       
0      -146.8230  63.6890        1        2  0.360  ...  2725  2385   42.0   
1      -146.6680  63.4040        1        2  0.558  ...  2144  2005   16.0   
2      -146.0800  63.3760        1        2  1.685  ...  2182  1868   18.0   
3      -146.1200  63.3810        1        2  3.681  ...  2317  1944   19.0   
4      -147.0570  63.5510        1        2  2.573  ...  2317  1914   16.0   
...          ...      ...      ...      ...    ...  ...   ...   ...    ...   
216424  -37.7325 -53.9860       19        3  0.042  ...   510  -999   29.9   
216425  -36.1361 -54.8310       19        3  0.567  ...   830  -999   23.6   
216426  -37.3018 -54.1884       19        3  4.118  ...  1110  -999   16.8   
216427  -90.4266 -68.8656       19        1  0.011  ...   270  -999    0.4   
216428   37.7140 -46.8972       19        4  0.528  ...  1170  -999    9.6   

实现步骤 #

1. 导入依赖库 #

import geopandas as gpd
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

2. 给每个点分配网格索引 #

0.5度分辨率的全球网格包含720个经度单元（-180°到180°）和360个纬度单元（-90°到90°）。计算每个点对应的网格索引：

# 计算经度索引：范围0-719
gdf['lon_idx'] = ((gdf['CENLON'] + 180) / 0.5).astype(int)
# 计算纬度索引：范围0-359
gdf['lat_idx'] = ((gdf['CENLAT'] + 90) / 0.5).astype(int)

3. 按网格聚合计算均值 #

以网格索引为分组键，计算每个网格内O1REGION的均值：

# 分组聚合
grid_mean = gdf.groupby(['lon_idx', 'lat_idx'])['O1REGION'].mean().reset_index()

如果需要按AREA列加权计算均值，可替换为：

grid_mean = gdf.groupby(['lon_idx', 'lat_idx']).apply(
    lambda x: np.average(x['O1REGION'], weights=x['AREA'])
).reset_index(name='O1REGION_mean')

4. 生成720x360的二维矩阵 #

初始化空矩阵并填充聚合结果，无数据的网格单元会保留为NaN：

# 创建360行（纬度）×720列（经度）的空矩阵
grid_matrix = np.full((360, 720), np.nan)

# 将聚合值填入对应网格位置
for _, row in grid_mean.iterrows():
    grid_matrix[int(row['lat_idx']), int(row['lon_idx'])] = row['O1REGION']

5. 用Cartopy可视化验证 #

通过绘图确认矩阵的空间分布：

# 创建PlateCarree投影的绘图对象
fig, ax = plt.subplots(figsize=(12, 6), subplot_kw={'projection': ccrs.PlateCarree()})

# 绘制矩阵，origin='lower'对应纬度从南到北
im = ax.imshow(
    grid_matrix,
    origin='lower',
    extent=[-180, 180, -90, 90],
    transform=ccrs.PlateCarree(),
    cmap='viridis'
)

# 添加海岸线参考
ax.coastlines()

# 添加颜色条与标题
plt.colorbar(im, ax=ax, label='O1REGION 均值')
plt.title('全球0.5度分辨率O1REGION均值矩阵')
plt.show()

注意事项 #

• 确保CENLON和CENLAT为EPSG:4326坐标系的经纬度，若数据使用其他CRS需先转换
• 无数据的网格单元值为NaN，可根据需求用np.nan_to_num()或插值方法填充
• 网格索引计算时，若经纬度刚好落在边界上，astype(int)会向下取整，可根据需求调整索引规则

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