在Python中以0.5度分辨率将纬度、经度、数据点转换为矩阵（2D网格）

Question

我有一个地理数据框，加载如下：

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

看起来像：

print(gdf)
          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   

我想以 0.5 度分辨率（720x360 世界地图）创建“01REGION”列的 2D 矩阵（世界地图），并以平均值作为聚合方法。我该如何做到这一点（最好使用 cartopy？）

Answer 1

对于任何看起来相同的人：

## Create a 0.5 degree map of pixels with corresponding RGI regions
# Extract relevant data
rgi = gdf['01REGION'].values
lon = gdf['CENLON'].values
lat = gdf['CENLAT'].values
# Set resolution and create grid
resolution = 0.5
x_range = np.arange(lon.min(), lon.max()+resolution, resolution)
x_range = np.arange(-179.75,179.75+resolution,resolution)
y_range = np.arange(lat.min(), lat.max()+resolution, resolution)
y_range = np.arange(-89.75,89.75+resolution,resolution)
grid_x, grid_y = np.meshgrid(x_range, y_range)
# Interpolate data onto grid
points = list(zip(lon, lat))
grid_z = griddata(points, rgi, (grid_x, grid_y), method='nearest')
# Define transformation and CRS
transform = Affine.translation(grid_x[0][0]-resolution/2, grid_y[0][0]-resolution/2) * Affine.scale(resolution, resolution)
crs = CRS.from_epsg(4326)
# Write interpolated raster to file
interp_raster = rasterio.open('RGI.tif',
                              'w',
                              driver='GTiff',
                              height=grid_z.shape[0],
                              width=grid_z.shape[1],
                              count=1,
                              dtype=grid_z.dtype,
                              crs=crs,
                              transform=transform)
interp_raster.write(grid_z, 1)
interp_raster.close()
rgi_raster=np.squeeze(rioxarray.open_rasterio("RGI.tif").round())