Plot BNG grid data on Matplotlib Basemap

Question

I would like to plot UK population density data in Matplotlib Basemap.

I got data in ESRI ASCII using BNG OSGB63 from https://catalogue.ceh.ac.uk/documents/61f10c74-8c2c-4637-a274-5fa9b2e5ce44

and then used the following code from https://stackoverflow.com/questions/23264404/gridded-data-basemap-mismatch-what-have-i-done-wrong#comment75355592_23264404

which in my case looks like this: https://pastebin.com/52GH9FnH

Here is an excerpt:

# Define a BNG coordinate grid for the data, based on information in the header:
   # LLcorner eastings = -200000
   # LLcorner northings = -200000
   # grid box size = 5000
   # All units in metres
lle = -10000 + 2500 # because we want to plot the box point at the centre of each box, not its lower left corner.
lln = 0 + 2500
delta = 1000
nn,ne = grid.shape

ns = np.arange(0.,nn,1)*delta + lln
es = np.arange(0.,ne,1)*delta + lle

# convert to mesh of eastings/northings.
eastings,northings = np.meshgrid(es,ns)

# convert mesh to lon/lat (wgs84) using pyproj:
bng = pyproj.Proj(init='epsg:27700')
wgs84 = pyproj.Proj(init='epsg:4326')
lons,lats = pyproj.transform(bng,wgs84,eastings,northings)

lats = lats[::-1] # because otherwise the map is upside down (differences between imshow and pcolormesh)
fig = plt.figure(figsize=(15, 15))
# Get a basemap of the UK
m = Basemap(projection='merc',llcrnrlon=-15 ,llcrnrlat=47.8,urcrnrlon=7.2 ,urcrnrlat=62.5,resolution='h')

# project the lon/lat grids onto map coordinates
x,y = m(lons,lats)

# Draw the map!
m.drawcoastlines()
m.pcolormesh(x,y,grid,cmap=plt.cm.YlOrRd)
plt.show()

This is the result:

Any idea, what I'm doing wrong and how to make the basemap outline and my grid data align?