底图和密度图答案

【问题标题】：Basemap and density plots底图和密度图
【发布时间】：2012-07-15 11:50:40
【问题描述】：

我想使用底图创建密度图。我拥有的数据是无网格的并且是重复的或彼此非常接近。我尝试对数据进行网格化，然后使用 pcolor 绘制 bin 的数量，但即使所有数据集的长度相同，我仍然会收到缓冲区大小错误。我最初的想法是使用下面的底图脚本，但我只能让 scatter 工作，尽管这并没有给我一个密度图。

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])

关于如何获取此数据的密度图的任何想法？感谢帮助。 /M

【问题讨论】：

希望我的回答有帮助。如果您愿意并且能够分享您的数据，那么知道它总是很有趣:-)

标签： python numpy matplotlib python-2.7 matplotlib-basemap

【解决方案1】：

如果您有兴趣对网格框中的每个纬度进行频率计数，可以使用 numpy 函数 histogram2d。

以下代码可能是您正在寻找的：

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()

编辑：我添加了更多内联 cmets 以帮助理解/可读性。

【讨论】：

看起来很棒！谢谢。这比我想象的要好。我对您的解决方案有几个问题。为什么要转置密度？为什么使用 numpy meshgrid xs, ys = m(*numpy.meshgrid(xs, ys))？你能解释一下这条线的作用吗？一个很好的解决方案，我什至从未想过。想象一下，一个二维直方图。酷！
您的权利。我将在回答中提供更多解释。
我已经添加了一些进一步的内联 cmets。转置问题是不幸的（而不是转置，我改变了纬度/经度的顺序）并且是一个奇怪的设计决定（尽管我确信它必须有充分的理由）。
感谢此次更新。它看起来更准确，尤其是在菲律宾。我昨晚很晚才玩这个，并注意到了不同之处。由于我将在同行评审文章中使用完整数据，因此我无法在此处发布。但它看起来很像你的例子。我不太确定如何发布图像。我真的很喜欢你的解决方案。我想知道你是否尝试过 pcolor 而不是 histogram_2d？
请注意，底图对象有它自己的 pcolormesh 和 colorbar 方法，它们也是从 matplotlib 派生的。它们可能更适合底图布局或更好地支持它。