1

I am trying to understand why a hexbin plot in a north or south polar stereo projection shows squashed hexagons, even though the area of the grid is square and the projection is approximately equal area.

I've tried both north and south polar stereo projections using basemap.

import numpy as np
from numpy.random import uniform
import matplotlib.pyplot as plt 
from mpl_toolkits.basemap import Basemap

fig = plt.figure(figsize=(12,10)) # width, height in inches
ax =fig.add_axes([-0.02,0.1,0.74,0.74]) 

m = Basemap(epsg='3413',lon_0=0.,resolution='l',width=6000000,height=6000000)

m.drawcoastlines()

m.drawmapscale(0.,90.,0.,90.,1000)

npts=2000
lats = uniform(60.,80.,size=npts)
lons = uniform(0.,360.,size=npts)
data = uniform(0.,4800.,size=npts)

x,y=m(lons, lats)

thiscmap=plt.cm.get_cmap('viridis')

p=m.hexbin(x,y,C=data,gridsize=[10,10],cmap=thiscmap)

plt.show()

plot output

finefoot
  • 9,914
  • 7
  • 59
  • 102
Alan
  • 13
  • 2
  • Not quite sure why but if I do p=m.hexbin(x,y,C=data,gridsize=10,cmap=thiscmap) instead of p=m.hexbin(x,y,C=data,gridsize=[10,10],cmap=thiscmap) then I get more regular sized hexagons – Alan Jan 09 '19 at 23:58

1 Answers1

2

I don't know why you get squashed hexagons. But you can adjust the hexagon shape by setting appropriate values of gridsize. Here I modify your code and get better plot.

import numpy as np
from numpy.random import uniform
import matplotlib.pyplot as plt 
from mpl_toolkits.basemap import Basemap

fig = plt.figure(figsize=(12,10)) # width, height in inches
ax =fig.add_axes([-0.02, 0.1, 0.74, 0.74]) 

# North polar stereographic projection epsg='3413'; ***large areal distortion***
#m = Basemap(epsg='3413', lon_0=0., resolution='c', width=6000000, height=6000000)

# 'laea':  Lambert Azimuthal Equal Area
# Thematic mapping with ground surface data should be plotted on 'equal-area' projection
m = Basemap(projection='laea', lon_0=0., lat_0=90, resolution='l', width=6000000, height=6000000)

m.drawcoastlines(linewidth=0.5)

m.drawmapscale(0.,90.,0.,90.,1000)  # 1000 km?

npts = 2000
lats = uniform(60.,80.,size=npts)  # not cover N pole
lons = uniform(0.,360.,size=npts)  # around W to E
data = uniform(0.,4800.,size=npts)

x,y = m(lons, lats)

thiscmap = plt.cm.get_cmap('viridis')

# To get 'rounded' hexagons, gridsize should be specified appropriately
# need some trial and error to get them right
#p=m.hexbin(x, y, C=data, gridsize=[10,10], cmap=thiscmap)  # original code
m.hexbin(x, y, C=data, gridsize=[16,11], cmap=thiscmap)     # better

plt.colorbar()  # useful on thematic map

plt.show()

The projection you use (epsg:3413) is stereographic projection which has large areal distortion. More appropriate projection for thematic mapping in this case is Lambert Azimuthal Equal Area.

The resulting plot:

enter image description here

swatchai
  • 17,400
  • 3
  • 39
  • 58
  • Thanks for your help. I also found that I could manually adjust the grid size parameter to get more regular hexagons. I guess that the hexagon density will be greater in the 'y' direction because of the way they are packed, and hence a 16:11 ratio is about correct. If you only give a single grid size parameter then hex bin works it out itself. This is more obvious when looking at your plot than my original one. I guess it would be interesting to know what the exact ratio should be (maybe 16:11 is it?) – Alan Jan 10 '19 at 15:34