【问题标题】:Python - 2D Histogram plot in log scale -- Error: `cannot convert float NaN to integer`Python - 对数刻度的 2D 直方图 - 错误:`无法将浮点 NaN 转换为整数`
【发布时间】:2016-07-28 16:56:09
【问题描述】:

这是我之前想要转换为二维直方图的图。

mass_bh = (subhalos['SubhaloBHMass'] * 1e10 / 0.704) # in units of M_sol h^-1
vdisp = subhalos['SubhaloVelDisp']

nbins = 200
H, xedges, yedges = np.histogram2d(mass_bh,vdisp,bins=nbins)

fig2 = plt.figure()
plt.pcolormesh(xedges,yedges,Hmasked)
cbar = plt.colorbar()
cbar.ax.set_ylabel('g-r')

plt.ylabel(' $\log(\sigma)\quad$ [km s$^{-1}$] ')
plt.xlabel('$\log(M_{BH})\quad$ [M$_{\odot}$]')
plt.title('$M_{BH}-\sigma$ relation')

这反而给了我这个

我之前的绘图将其xy 值都转换为对数缩放。但是对于这种直方图转换,效果并不理想。

我该如何解决这个问题?

谢谢!

【问题讨论】:

  • 如果不知道您希望您的绘图是什么样子或您拥有什么样的数据,就很难知道您所说的“不正确”是什么意思。
  • @Lanery - 你是对的。让我快速详细地添加一张图片。
  • 问题似乎出在数据上。注意轴上的限制是完全不同的。
  • @armitita - 是的,你是对的。上一个绘图对我的 xy 绘图使用了对数缩放,创建了您看到的第一个绘图。对直方图实施对数缩放会导致绘图不显示。第二个是没有对数缩放。

标签: python matplotlib plot histogram


【解决方案1】:

@armatita 关于数据的问题是正确的。我认为这一切都取决于你如何在histogram2d 内进行分箱。看看这个带有随机对数正态分布的例子是否有帮助。

import numpy as np
import matplotlib.pyplot as plt

n = 1000

x = np.logspace(2, 10, n)
y = x**1.5
y = y * np.random.lognormal(10, 3, n)

x_bins = np.logspace(np.log10(x.min()), np.log10(x.max()), np.sqrt(n))
y_bins = np.logspace(np.log10(y.min()), np.log10(y.max()), np.sqrt(n))
H, xedges, yedges = np.histogram2d(x, y, bins=[x_bins, y_bins])

fig = plt.figure()
ax1 = fig.add_subplot(211)
ax1.plot(x, y, 'o')
ax1.set_xscale('log')
ax1.set_yscale('log')

ax2 = fig.add_subplot(212)
ax2.pcolormesh(xedges, yedges, H.T)
ax2.set_xscale('log')
ax2.set_yscale('log')

我得到了下面的图片,我相信这就是您正在寻找的。还要注意H 上的转置。

【讨论】:

  • 酷!谢谢你。一件事,我从np.log10(vdisp.max()))]) 得到这个错误cannot convert float NaN to integer。另外,log10 中有些东西被零除
  • 这些都是你的数据的问题,所以你可能想调查为什么你的数据集中有nans,但同时你可以使用np.nanmin忽略nan s。见this post
  • 那么我会在 2dhistogram bin 部分应用 np.min 吗?
  • 让我详细说明一下,我将如何从您提供的链接中应用 np.min?我的两个数组的形状都约为 120000。
  • 您的意思是解决nans 的问题吗?我不明白你为什么不想通过像here 所示的某种过滤来完全摆脱你的nans。我会在调用 histogram2d 之前进行过滤。
【解决方案2】:

只是为了激发您的好奇心而提出的建议。尽管@lanery 清楚地回答了这个问题,但我想分享一种在 python 中获得漂亮的 2d 直方图的不同方法。我不想使用通常会产生非常丑陋的直方图的 np.histogram2d,而是想回收 py-sphviewer,这是一个使用自适应平滑内核渲染粒子模拟的 python 包。考虑以下代码,它基于 lanery 的示例:

将 numpy 导入为 np 将 matplotlib.pyplot 导入为 plt 将 sphviewer 导入为 sph

def myplot(x, y, extent=None, nb=8, xsize=500, ysize=500):   
    if(extent == None):
        xmin = np.min(x)
        xmax = np.max(x)
        ymin = np.min(y)
        ymax = np.max(y)
    else:
        xmin = extent[0]
        xmax = extent[1]
        ymin = extent[2]
        ymax = extent[3]

    k, = np.where( (x <= xmax) & (x >= xmin) & 
                   (y <= ymax) & (y >= ymin) )

    pos = np.zeros([3, len(k)])
    pos[0,:] = (x[k]-xmin)/(xmax-xmin)
    pos[1,:] = (y[k]-ymin)/(ymax-ymin)
    w = np.ones(len(k))

    P = sph.Particles(pos, w, nb=nb)
    S = sph.Scene(P)
    S.update_camera(r='infinity', x=0.5, y=0.5, z=0, 
                    extent=[-0.5,0.5,-0.5,0.5],
                    xsize=xsize, ysize=ysize)
    R = sph.Render(S)
    R.set_logscale()
    img = R.get_image()

    return img, [xmin,xmax,ymin,ymax]    


n = 1000

x = np.logspace(2, 10, n)
y = x**1.5
y = y * np.random.lognormal(10, 3, n)

H, xedges, yedges = np.histogram2d(x, y, bins=[np.logspace(np.log10(x.min()), np.log10(x.max())),
                                               np.logspace(np.log10(y.min()), np.log10(y.max()))])


img, extent = myplot(np.log10(x), np.log10(y))   #Call the function to make the 2d-histogram

fig = plt.figure()
ax1 = fig.add_subplot(311)
ax1.plot(x, y, 'o')
ax1.set_xscale('log')
ax1.set_yscale('log')

ax2 = fig.add_subplot(312)
ax2.pcolormesh(xedges, yedges, H.T)
ax2.set_xscale('log')
ax2.set_yscale('log')

ax3 = fig.add_subplot(313)
ax3.imshow(img, origin='lower', extent=extent, aspect='auto')

plt.show()

产生以下输出:

函数 myplot() 只是我编写的一个非常简单的函数,用于规范化数据并将其作为 py-sphviewer 的输入。平滑内核的长度通常由参数 nb 给出,它指定了执行平滑的邻居的数量。虽然乍一看似乎很复杂,但想法和实现都非常简单,并且与 np.histogram2d 相比,结果远远优于 np.histogram2d。但当然,这取决于您是否能够传播数据,以及这样做对您的研究有何意义和后果。

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2018-04-30
    • 2020-08-14
    • 1970-01-01
    • 2021-11-16
    相关资源
    最近更新 更多