|
6 | 6 |
|
7 | 7 | # create north polar stereographic basemap |
8 | 8 | m = Basemap(lon_0=270, boundinglat=20, projection='npstere',round=True) |
| 9 | +#m = Basemap(lon_0=-105,lat_0=40,projection='ortho') |
9 | 10 |
|
10 | | -# number of points to plot. |
| 11 | +# number of points, bins to plot. |
11 | 12 | npts = 10000 |
12 | 13 | bins = 40 |
13 | 14 | # generate random points on a sphere, |
|
16 | 17 | # http://mathworld.wolfram.com/SpherePointPicking.html |
17 | 18 | u = uniform(0.,1.,size=npts) |
18 | 19 | v = uniform(0.,1.,size=npts) |
19 | | -lons1 = 360.*u |
20 | | -lats1 = (180./np.pi)*np.arccos(2*v-1) - 90. |
| 20 | +lons = 360.*u |
| 21 | +lats = (180./np.pi)*np.arccos(2*v-1) - 90. |
21 | 22 | # toss points outside of map region. |
22 | | -lats = np.compress(lats1 > 20, lats1) |
23 | | -lons = np.compress(lats1 > 20, lons1) |
| 23 | +lats = np.compress(lats > 20, lats) |
| 24 | +lons = np.compress(lats > 20, lons) |
24 | 25 | # convert to map projection coordinates. |
25 | | -x, y = m(lons, lats) |
| 26 | +x1, y1 = m(lons, lats) |
| 27 | +# remove points outside projection limb. |
| 28 | +x = np.compress(np.logical_or(x1 < 1.e20,y1 < 1.e20), x1) |
| 29 | +y = np.compress(np.logical_or(x1 < 1.e20,y1 < 1.e20), y1) |
26 | 30 | # function to plot at those points. |
27 | 31 | xscaled = 4.*(x-0.5*(m.xmax-m.xmin))/m.xmax |
28 | 32 | yscaled = 4.*(y-0.5*(m.ymax-m.ymin))/m.ymax |
|
31 | 35 | # make plot using hexbin |
32 | 36 | fig = plt.figure(figsize=(12,5)) |
33 | 37 | ax = fig.add_subplot(121) |
34 | | -CS = plt.hexbin(x,y,C=z,gridsize=bins,cmap=plt.cm.jet) |
| 38 | +CS = m.hexbin(x,y,C=z,gridsize=bins,cmap=plt.cm.jet) |
35 | 39 | # draw coastlines, lat/lon lines. |
36 | 40 | m.drawcoastlines() |
37 | 41 | m.drawparallels(np.arange(0,81,20)) |
|
41 | 45 |
|
42 | 46 | # use histogram2d instead of hexbin. |
43 | 47 | ax = fig.add_subplot(122) |
44 | | -#m = Basemap(lon_0=270, boundinglat=20, projection='npstere',round=True) |
| 48 | +# remove points outside projection limb. |
45 | 49 | bincount, xedges, yedges = np.histogram2d(x, y, bins=bins) |
46 | 50 | mask = bincount == 0 |
47 | 51 | # reset zero values to one to avoid divide-by-zero |
48 | 52 | bincount = np.where(bincount == 0, 1, bincount) |
49 | 53 | H, xedges, yedges = np.histogram2d(x, y, bins=bins, weights=z) |
50 | 54 | H = np.ma.masked_where(mask, H/bincount) |
51 | | -# set color of masked values to white (hexbin does this by default) |
| 55 | +# set color of masked values to axes background (hexbin does this by default) |
52 | 56 | palette = plt.cm.jet |
53 | | -palette.set_bad('white', 1.0) |
| 57 | +palette.set_bad(ax.get_axis_bgcolor(), 1.0) |
54 | 58 | CS = m.pcolormesh(xedges,yedges,H.T,shading='flat',cmap=palette) |
55 | 59 | # draw coastlines, lat/lon lines. |
56 | 60 | m.drawcoastlines() |
|
0 commit comments