Convert lat,lon,data points to matrix (2D grid) at 0.5 degree resolution in Python

Question

I have a geodataframe which I load in as follows:

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

It looks as:

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   

I want to create a 2D matrix (world map) of the column "01REGION" at a 0.5 degree resolution (720x360 world map) with the mean as the aggregation method. How can I do this (preferably with cartopy?)

Answer 1

For anyone looking the same:

## 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())