use hexbin method.

Author: Jeff Whitaker

examples/hexbin_demo.py

@@ -6,8 +6,9 @@
 
 # create north polar stereographic basemap
 m = Basemap(lon_0=270, boundinglat=20, projection='npstere',round=True)
+#m = Basemap(lon_0=-105,lat_0=40,projection='ortho')
 
-# number of points to plot.
+# number of points, bins to plot.
 npts = 10000
 bins = 40
 # generate random points on a sphere,
@@ -16,13 +17,16 @@
 # http://mathworld.wolfram.com/SpherePointPicking.html
 u = uniform(0.,1.,size=npts)
 v = uniform(0.,1.,size=npts)
-lons1 = 360.*u
-lats1 = (180./np.pi)*np.arccos(2*v-1) - 90.
+lons = 360.*u
+lats = (180./np.pi)*np.arccos(2*v-1) - 90.
 # toss points outside of map region.
-lats = np.compress(lats1 > 20, lats1)
-lons = np.compress(lats1 > 20, lons1)
+lats = np.compress(lats > 20, lats)
+lons = np.compress(lats > 20, lons)
 # convert to map projection coordinates.
-x, y = m(lons, lats)
+x1, y1 = m(lons, lats)
+# remove points outside projection limb.
+x = np.compress(np.logical_or(x1 < 1.e20,y1 < 1.e20), x1)
+y = np.compress(np.logical_or(x1 < 1.e20,y1 < 1.e20), y1)
 # function to plot at those points.
 xscaled = 4.*(x-0.5*(m.xmax-m.xmin))/m.xmax
 yscaled = 4.*(y-0.5*(m.ymax-m.ymin))/m.ymax
@@ -31,7 +35,7 @@
 # make plot using hexbin
 fig = plt.figure(figsize=(12,5))
 ax = fig.add_subplot(121)
-CS = plt.hexbin(x,y,C=z,gridsize=bins,cmap=plt.cm.jet)
+CS = m.hexbin(x,y,C=z,gridsize=bins,cmap=plt.cm.jet)
 # draw coastlines, lat/lon lines.
 m.drawcoastlines()
 m.drawparallels(np.arange(0,81,20))
@@ -41,15 +45,16 @@
 
 # use histogram2d instead of hexbin.
 ax = fig.add_subplot(122)
+# remove points outside projection limb.
 bincount, xedges, yedges = np.histogram2d(x, y, bins=bins)
 mask = bincount == 0
 # reset zero values to one to avoid divide-by-zero
 bincount = np.where(bincount == 0, 1, bincount)
 H, xedges, yedges = np.histogram2d(x, y, bins=bins, weights=z)
 H = np.ma.masked_where(mask, H/bincount)
-# set color of masked values to white (hexbin does this by default)
+# set color of masked values to axes background (hexbin does this by default)
 palette = plt.cm.jet
-palette.set_bad('white', 1.0)
+palette.set_bad(ax.get_axis_bgcolor(), 1.0)
 CS = m.pcolormesh(xedges,yedges,H.T,shading='flat',cmap=palette)
 # draw coastlines, lat/lon lines.
 m.drawcoastlines()