netCDF polygons conversion to fixed-sized raster image

Question

I'm trying to convert data from netCDF files; the data contains polygons and values for methane concentrations.

from matplotlib.patches import Polygon
import cartopy.crs as ccrs
import cartopy.feature as cfeature
# %% Plotting the polygons
fig, ax = plt.subplots(figsize=(10, 10), subplot_kw={'projection': ccrs.PlateCarree()})
ax.set_extent([-9, -5, 32, 34], crs=ccrs.PlateCarree())

# Add a base map
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=':')

# Normalize the ch4 values for the colormap
norm = mcolors.Normalize(vmin=1850, vmax=1880)
cmap = matplotlib.colormaps['rainbow']

# Plot the grid cells and color them based on ch4 values
for lat_corners, lon_corners, ch4 in zip(lat_corners, lon_corners, ch4_values):
    color = cmap(norm(ch4))
    poly = Polygon(list(zip(lon_corners, lat_corners)), facecolor=color,
                   edgecolor='grey', linewidth=0.5, alpha=0.7)
    ax.add_patch(poly)

# Add a colorbar
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([])
cbar = plt.colorbar(sm, ax=ax, orientation='horizontal', fraction=0.046, pad=0.04)
cbar.set_label('Methane Concentration (ppb)')
plt.show()

This image is very specific, and I want to apply this to generate a dataset with the shape of [N,M,32x32] pixel images, where N is the number of images and M contains different forms of data (methane values, wind, etc.).

One of the big issues I've been having is treating spaces between polygons (where one is not there) as a pixel.

I apologize for any lack of information; this is something very new to me.

def project_to_image(lat_corners, lon_corners, ch4_values):
    # Create a 32x32 empty image
    #image = np.zeros((32, 32, 3), dtype=np.uint8)
    #image = np.full((32, 32, 3), 255, dtype=np.uint8)
    norm = mcolors.Normalize(vmin=1850, vmax=1880)
    cmap = matplotlib.colormaps['rainbow']
    #min_value_color = cmap(norm(1850))  # Color for the minimum value in the range
    #min_value_color = (np.array(min_value_color[:3]) * 255).astype(np.uint8)  # Convert to RGB
    min_value_color = 255
    # Create a 32x32 empty image with the background color set to the minimum value color
    image = np.full((32, 32, 3), min_value_color, dtype=np.uint8)

    # Create a Basemap instance
    m = Basemap(projection='cyl', llcrnrlat=31.5, urcrnrlat=34.0, llcrnrlon=-9.0, urcrnrlon=-5.0, resolution='h')

    for lat_corner, lon_corner, ch4_value in zip(lat_corners, lon_corners, ch4_values):
        if np.any(~lon_corner.mask) and np.any(~lat_corner.mask):
            color = cmap(norm(ch4_value))
            x, y = m(lon_corner[~lon_corner.mask], lat_corner[~lat_corner.mask])

            # Limit coordinates to Basemap range
            x = np.clip(x, m.llcrnrx, m.urcrnrx)
            y = np.clip(y, m.llcrnry, m.urcrnry)


            # Adjust coordinates to fit 32x32 image
            # Map coordinates to 0-31 range
            x = (31*(x - m.llcrnrx) / (m.urcrnrx - m.llcrnrx) ).astype(int)
            y = (31*(y - m.llcrnry) / (m.urcrnry - m.llcrnry) ).astype(int)
            # Filter out of range coordinates
            #valid_idx = (x >= 0) & (x < 32) & (y >= 0) & (y < 32)
            #x = x[valid_idx]
            #y = y[valid_idx]

            # Draw polygon
            for i in range(len(x)):
                image[int(y[i]), int(x[i])] = (np.array(color[:3]) * 255).astype(np.uint8)

    return image

what i am trying to do:

The research article that I am trying to replicate the data from can be found here:

Research article