Dimensions of C (1801, 3600) are incompatible with X (3600) and/or Y (1801); see help(pcolormesh)

Question

I'm looking to show surface wind data from "https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land?tab=form"

I'm using the code below to unpack the data in the form of a diagram "https://confluence.ecmwf.int/display/CKB/How+to+plot+GRIB+files+with+Python+and+matplotlib"

However I get the error 'TypeError: Dimensions of C (1801, 3600) are incompatible with X (3600) and/or Y (1801); see help(pcolormesh)', which should work as C is (rows, columns) and X represents columns and Y represents rows, but this error suggests the data does not fit?

The code is below, please any suggestions would be massively appreciated, thanks!

import pygrib
import matplotlib.pyplot as plt
import matplotlib.colors as colors
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import shiftgrid
import numpy as np
 
plt.figure(figsize=(12,8))
 
grib = 'adaptor.mars.internal-1669570066.0148444-9941-17-702abb2d-e37e-4ef7-a19e-a04fb24e5a20.grib' # Set the file name of your input GRIB file
grbs = pygrib.open(grib)
 
grb = grbs.select()[0]
data = grb.values
 
# need to shift data grid longitudes from (0..360) to (-180..180)
lons = np.linspace(float(grb['longitudeOfFirstGridPointInDegrees']), \
float(grb['longitudeOfLastGridPointInDegrees']), int(grb['Ni']) )
lats = np.linspace(float(grb['latitudeOfFirstGridPointInDegrees']), \
float(grb['latitudeOfLastGridPointInDegrees']), int(grb['Nj']) )
data, lons = shiftgrid(180., data, lons, start=False)
grid_lon, grid_lat = np.meshgrid(lons, lats) #regularly spaced 2D grid
 
m = Basemap(projection='cyl', llcrnrlon=-180, \
    urcrnrlon=180.,llcrnrlat=lats.min(),urcrnrlat=lats.max(), \
    resolution='c')
    
m.drawcoastlines()
m.drawmapboundary()
m.drawparallels(np.arange(-90.,120.,30.),labels=[1,0,0,0])
m.drawmeridians(np.arange(-180.,180.,60.),labels=[0,0,0,1])
    
x, y = m(grid_lon, grid_lat)

#x = np.transpose(x)
#y = np.transpose(y)

cs = m.pcolormesh(x,y,data,shading='flat',cmap=plt.cm.gist_stern_r)
 
plt.colorbar(cs,orientation='vertical', shrink=0.5)
plt.title('CAMS AOD forecast') # Set the name of the variable to plot
plt.savefig(grib+'.png') # Set the output file name

However, even if I transpose the x and y data I get the same error type but: X (1801) and/or Y (3600). What ideally should happen is a map loads with North-South and East-West wind surface speeds.

Numpy pcolormesh: TypeError: Dimensions of C are incompatible with X and/or Y

Above is a similar question. However, the suggested transposing data solution which worked in that case did not work here.