0

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: map

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

  • I think that m = Basemap(projection='merc',llcrnrlon=-15 ,llcrnrlat=47.8,urcrnrlon=7.2 ,urcrnrlat=62.5,resolution='h') is not going to be in lat/lon so may need reprojecting too? – Ian Turton May 26 '17 at 9:09
  • Hi @iant, thanks for the suggestion. I tried # use major and minor sphere radii from WGS84 ellipsoid. m = Basemap(llcrnrlon=-145.5,llcrnrlat=1.,urcrnrlon=-2.566,urcrnrlat=46.352,\ rsphere=(6378137.00,6356752.3142),\ resolution='l',area_thresh=1000.,projection='lcc',\ lat_1=50.,lon_0=-107.,ax=ax) as well as explicitly setting the ` ellps='WGS84',` - both without a difference. I have the suspicion that's not it. – Dan Schien May 26 '17 at 9:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.