Base maps and density graphs

Question

I would like to create a density graph using a basemap. The data that I have remains unchanged and is repeated or very close to each other. I tried gridding the data and then plotting the number of boxes using pcolor, but I keep getting a buffer size error even if all the data sets are of equal length. My initial idea was to use the base map script below, but I can only get the scatter for work, although this does not give me a density graph.

m = Basemap(resolution='f',projection='merc', lon_0=160, llcrnrlat=-30.0, urcrnrlat=30.0, llcrnrlon=100., urcrnrlon=270.0, lat_ts=0.0) m.drawmapboundary(fill_color='white') m.fillcontinents(color='#F5DEB3',lake_color='#85A6D9') m.drawcoastlines(color='black', linewidth=.4) m.drawcountries(color='#6D5F47', linewidth=.4) m.drawmeridians(np.arange(0, 360, 20), labels=[0,0,0,1], color='black', dashes=[1,0], labelstyle='+/-', linewidth=0.2) m.drawparallels(np.arange(-90, 90, 10), labels=[1,0,0,0], color='black', dashes=[1,0], labelstyle='+/-', linewidth=0.2) lats = array([ 14.375, 14.125, 14.125, 9.375, 13.625, 14.375, 5.625, 13.875, 14.625, 5.875, 8.875, 5.625, 8.875, 13.375, 11.125, 8.375, 12.375, 6.125, 5.375, 8.375, 7.375, 7.875, 14.375, 14.875, 9.875, 11.125, 14.875, 7.875, 9.125, 11.625, 5.125, 10.875, 5.125, 12.125, 12.625, 8.625, 5.125, 8.375, 11.625, 11.375, 12.875, 14.375, 8.875, 8.375, 6.375, 8.625, 5.875, 8.125, 9.375, 5.875, 8.125, 8.875, 5.375, 8.875, 5.625, 11.875, 9.875, 9.875, 10.875, 11.375, 9.875, 9.375, 13.125, 14.125, 8.125, 14.875, 9.875, 9.625, 10.625, 12.125, 9.375, 5.625, 13.625, 6.375, 10.125, 14.875, 8.875, 5.125, 6.125, 8.625, 6.875, 9.375, 9.125, 9.375, 13.625, 6.125, 6.125, 6.875, 11.375, 13.375, 10.375, 6.875, 7.625, 7.625, 13.875, 5.125, 6.125, 14.125, 7.375, 5.375]) lons = array([ 122.125, 125.875, 122.375, 132.125, 124.375, 122.625, 206.125, 122.625, 120.375, 187.625, 243.625, 161.625, 220.375, 121.875, 125.625, 130.125, 219.625, 223.875, 223.375, 190.125, 132.875, 185.875, 122.625, 120.875, 259.375, 125.875, 204.875, 131.875, 130.875, 125.375, 234.375, 241.375, 188.125, 124.375, 124.375, 140.375, 205.625, 130.875, 256.875, 239.375, 124.125, 123.125, 131.625, 126.375, 241.125, 130.875, 233.625, 184.375, 205.125, 169.375, 150.375, 183.375, 168.875, 239.625, 180.125, 241.375, 131.375, 244.875, 244.375, 125.625, 131.375, 150.125, 203.625, 122.125, 243.125, 159.125, 125.625, 243.375, 125.875, 127.375, 130.125, 146.625, 203.125, 185.125, 204.625, 120.625, 130.125, 233.875, 131.875, 258.625, 130.125, 173.375, 258.125, 129.125, 137.875, 215.625, 214.125, 234.625, 125.375, 136.625, 231.125, 145.625, 128.625, 232.875, 125.125, 145.125, 239.875, 138.625, 241.375, 169.625]) data = array([ 84.8125 , 81.6875 , 79.75 , 78.5625 , 70.8125 , 70.75 , 69.0625 , 69. , 67.3125 , 63.65625, 63.375 , 62.5625 , 62.46875, 61.09375, 60.4375 , 59.21875, 58.84375, 58.8125 , 58.75 , 58.4375 , 57.8125 , 57.375 , 57.1875 , 57.15625, 57.15625, 57.09375, 56.96875, 56.875 , 56.3125 , 56.28125, 56.1875 , 55.59375, 55.5 , 55.46875, 55.4375 , 55.375 , 55.3125 , 55.21875, 55.21875, 55.15625, 54.90625, 54.75 , 54.6875 , 54.65625, 54.625 , 54.59375, 54.4375 , 54.34375, 54.21875, 54.0625 , 53.9375 , 53.75 , 53.65625, 53.21875, 53.15625, 53.0625 , 52.78125, 52.75 , 52.4375 , 52.25 , 52.0625 , 51.71875, 51.71875, 51.59375, 51.4375 , 51.3125 , 51.28125, 51.28125, 51.25 , 51.03125, 51. , 50.96875, 50.9375 , 50.875 , 50.8125 , 50.75 , 50.71875, 50.6875 , 50.5 , 50.4375 , 50.4375 , 50.28125, 50.21875, 50.125 , 50. , 49.5 , 49.5 , 49.40625, 49.21875, 49.21875, 49.0625 , 48.96875, 48.53125, 48.46875, 48.375 , 48.28125, 48.21875, 48.0625 , 47.875 , 47.84375]) 

Any ideas on how to get a graph of the density of this data? Appreciate the help. / M

Answer 1

If you are interested in performing a frequency for each lat-long in a gridbox, you can use the numpy2 function.

The following code might be what you are looking for:

 nx, ny = 10, 3 # compute appropriate bins to histogram the data into lon_bins = numpy.linspace(lons.min(), lons.max(), nx+1) lat_bins = numpy.linspace(lats.min(), lats.max(), ny+1) # Histogram the lats and lons to produce an array of frequencies in each box. # Because histogram2d does not follow the cartesian convention # (as documented in the numpy.histogram2d docs) # we need to provide lats and lons rather than lons and lats density, _, _ = numpy.histogram2d(lats, lons, [lat_bins, lon_bins]) # Turn the lon/lat bins into 2 dimensional arrays ready # for conversion into projected coordinates lon_bins_2d, lat_bins_2d = numpy.meshgrid(lon_bins, lat_bins) # convert the xs and ys to map coordinates xs, ys = m(lon_bins_2d, lat_bins_2d) plt.pcolormesh(xs, ys, density) plt.colorbar(orientation='horizontal') # overlay the scatter points to see that the density # is working as expected plt.scatter(*m(lons, lats)) plt.show() 

EDIT . I added additional inline comments to help with understanding / readability.