我有一组X,Y数据点(大约10k)，很容易绘制成散点图，但我想用热图来表示。
我查看了Matplotlib中的示例，它们似乎都已经从热图单元值开始生成图像。
有没有一种方法可以将一堆不同的x, y转换为热图(其中x, y频率较高的区域会“更温暖”)?
当前回答
编辑:为了更好地近似Alejandro的答案，请看下面。
我知道这是一个老问题，但想在Alejandro的回答中添加一些东西:如果你想要一个漂亮的平滑图像而不使用py-sphviewer，你可以使用np。Histogram2d，并对热图应用高斯滤波器(from scipy. nmage .filters):
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from scipy.ndimage.filters import gaussian_filter
def myplot(x, y, s, bins=1000):
heatmap, xedges, yedges = np.histogram2d(x, y, bins=bins)
heatmap = gaussian_filter(heatmap, sigma=s)
extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
return heatmap.T, extent
fig, axs = plt.subplots(2, 2)
# Generate some test data
x = np.random.randn(1000)
y = np.random.randn(1000)
sigmas = [0, 16, 32, 64]
for ax, s in zip(axs.flatten(), sigmas):
if s == 0:
ax.plot(x, y, 'k.', markersize=5)
ax.set_title("Scatter plot")
else:
img, extent = myplot(x, y, s)
ax.imshow(img, extent=extent, origin='lower', cmap=cm.jet)
ax.set_title("Smoothing with $\sigma$ = %d" % s)
plt.show()
生产:
Agape Gal'lo的散点图和s=16相互叠加(点击查看更好的视图):
我注意到我的高斯滤波方法和亚历杭德罗的方法的一个区别是，他的方法显示局部结构比我的好得多。因此，我在像素级上实现了一个简单的最近邻方法。该方法为每个像素计算数据中n个最近点距离的逆和。这种方法的分辨率很高，计算成本很高，我认为有更快的方法，所以如果你有任何改进，请告诉我。
更新:正如我所怀疑的，有一个更快的方法使用Scipy的Scipy . ckdtree。关于实现，请参阅Gabriel的回答。
总之，这是我的代码:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
def data_coord2view_coord(p, vlen, pmin, pmax):
dp = pmax - pmin
dv = (p - pmin) / dp * vlen
return dv
def nearest_neighbours(xs, ys, reso, n_neighbours):
im = np.zeros([reso, reso])
extent = [np.min(xs), np.max(xs), np.min(ys), np.max(ys)]
xv = data_coord2view_coord(xs, reso, extent[0], extent[1])
yv = data_coord2view_coord(ys, reso, extent[2], extent[3])
for x in range(reso):
for y in range(reso):
xp = (xv - x)
yp = (yv - y)
d = np.sqrt(xp**2 + yp**2)
im[y][x] = 1 / np.sum(d[np.argpartition(d.ravel(), n_neighbours)[:n_neighbours]])
return im, extent
n = 1000
xs = np.random.randn(n)
ys = np.random.randn(n)
resolution = 250
fig, axes = plt.subplots(2, 2)
for ax, neighbours in zip(axes.flatten(), [0, 16, 32, 64]):
if neighbours == 0:
ax.plot(xs, ys, 'k.', markersize=2)
ax.set_aspect('equal')
ax.set_title("Scatter Plot")
else:
im, extent = nearest_neighbours(xs, ys, resolution, neighbours)
ax.imshow(im, origin='lower', extent=extent, cmap=cm.jet)
ax.set_title("Smoothing over %d neighbours" % neighbours)
ax.set_xlim(extent[0], extent[1])
ax.set_ylim(extent[2], extent[3])
plt.show()
结果:
其他回答
创建一个与最终图像中的单元格对应的二维数组，称为say heatmap_cells，并将其实例化为全零。
选择两个比例因子来定义每个数组元素在实际单位中的差异，对于每个维度，例如x_scale和y_scale。选择这些，使所有数据点都在热图数组的范围内。
对于每个带x_value和y_value的原始数据点:
heatmap_cells[地板(x_value / x_scale),地板(y_value / y_scale)] + = 1
在Matplotlib词典，我认为你需要一个hexbin plot。
如果你不熟悉这种类型的图，它只是一个二元直方图，其中xy平面由一个规则的六边形网格镶嵌。
在直方图中，你可以数出每个六边形中的点的数量，将绘图区域离散化为一组窗口，将每个点分配给这些窗口中的一个;最后，将窗口映射到一个颜色数组上，你就得到了一个hexbin图。
虽然不像圆形或正方形那样常用，但直觉上，六边形是装箱容器的几何形状的更好选择:
六边形具有最近邻对称性(例如，方形容器没有， 例如，从正方形边界上的一点到另一点的距离 正方形内部并非处处相等)和 六边形是给出正平面的最高n多边形 镶嵌(例如，你可以安全地用六边形瓷砖重新设计厨房地板，因为当你完成时，瓷砖之间不会有任何空隙——而不是所有其他高n, n >= 7的多边形)。
(Matplotlib使用术语hexbin plot;所以(AFAIK)所有的绘图库的R;我仍然不知道这是否是这种类型的图表的普遍接受术语，尽管我怀疑它很可能是六角形装箱的缩写，这描述了准备数据显示的基本步骤。)
from matplotlib import pyplot as PLT
from matplotlib import cm as CM
from matplotlib import mlab as ML
import numpy as NP
n = 1e5
x = y = NP.linspace(-5, 5, 100)
X, Y = NP.meshgrid(x, y)
Z1 = ML.bivariate_normal(X, Y, 2, 2, 0, 0)
Z2 = ML.bivariate_normal(X, Y, 4, 1, 1, 1)
ZD = Z2 - Z1
x = X.ravel()
y = Y.ravel()
z = ZD.ravel()
gridsize=30
PLT.subplot(111)
# if 'bins=None', then color of each hexagon corresponds directly to its count
# 'C' is optional--it maps values to x-y coordinates; if 'C' is None (default) then
# the result is a pure 2D histogram
PLT.hexbin(x, y, C=z, gridsize=gridsize, cmap=CM.jet, bins=None)
PLT.axis([x.min(), x.max(), y.min(), y.max()])
cb = PLT.colorbar()
cb.set_label('mean value')
PLT.show()
而不是用np。我想回收py-sphviewer，这是一个使用自适应平滑内核渲染粒子模拟的python包，可以很容易地从pip安装(见网页文档)。考虑以下基于示例的代码:
import numpy as np
import numpy.random
import matplotlib.pyplot as plt
import sphviewer as sph
def myplot(x, y, nb=32, xsize=500, ysize=500):
xmin = np.min(x)
xmax = np.max(x)
ymin = np.min(y)
ymax = np.max(y)
x0 = (xmin+xmax)/2.
y0 = (ymin+ymax)/2.
pos = np.zeros([len(x),3])
pos[:,0] = x
pos[:,1] = y
w = np.ones(len(x))
P = sph.Particles(pos, w, nb=nb)
S = sph.Scene(P)
S.update_camera(r='infinity', x=x0, y=y0, z=0,
xsize=xsize, ysize=ysize)
R = sph.Render(S)
R.set_logscale()
img = R.get_image()
extent = R.get_extent()
for i, j in zip(xrange(4), [x0,x0,y0,y0]):
extent[i] += j
print extent
return img, extent
fig = plt.figure(1, figsize=(10,10))
ax1 = fig.add_subplot(221)
ax2 = fig.add_subplot(222)
ax3 = fig.add_subplot(223)
ax4 = fig.add_subplot(224)
# Generate some test data
x = np.random.randn(1000)
y = np.random.randn(1000)
#Plotting a regular scatter plot
ax1.plot(x,y,'k.', markersize=5)
ax1.set_xlim(-3,3)
ax1.set_ylim(-3,3)
heatmap_16, extent_16 = myplot(x,y, nb=16)
heatmap_32, extent_32 = myplot(x,y, nb=32)
heatmap_64, extent_64 = myplot(x,y, nb=64)
ax2.imshow(heatmap_16, extent=extent_16, origin='lower', aspect='auto')
ax2.set_title("Smoothing over 16 neighbors")
ax3.imshow(heatmap_32, extent=extent_32, origin='lower', aspect='auto')
ax3.set_title("Smoothing over 32 neighbors")
#Make the heatmap using a smoothing over 64 neighbors
ax4.imshow(heatmap_64, extent=extent_64, origin='lower', aspect='auto')
ax4.set_title("Smoothing over 64 neighbors")
plt.show()
产生如下图像:
如你所见，这些图像看起来非常漂亮，我们能够识别出它上面不同的子结构。这些图像是在一个特定的域内为每个点扩展一个给定的权重，由平滑长度定义，而平滑长度又由到更近的nb邻居的距离给出(我选择了16,32和64作为示例)。因此，高密度区域通常分布在较小的区域，与低密度区域相比。
myplot函数是我写的一个非常简单的函数它是为了将x y数据交给py-sphviewer来完成这个魔术。
编辑:为了更好地近似Alejandro的答案，请看下面。
我知道这是一个老问题，但想在Alejandro的回答中添加一些东西:如果你想要一个漂亮的平滑图像而不使用py-sphviewer，你可以使用np。Histogram2d，并对热图应用高斯滤波器(from scipy. nmage .filters):
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from scipy.ndimage.filters import gaussian_filter
def myplot(x, y, s, bins=1000):
heatmap, xedges, yedges = np.histogram2d(x, y, bins=bins)
heatmap = gaussian_filter(heatmap, sigma=s)
extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
return heatmap.T, extent
fig, axs = plt.subplots(2, 2)
# Generate some test data
x = np.random.randn(1000)
y = np.random.randn(1000)
sigmas = [0, 16, 32, 64]
for ax, s in zip(axs.flatten(), sigmas):
if s == 0:
ax.plot(x, y, 'k.', markersize=5)
ax.set_title("Scatter plot")
else:
img, extent = myplot(x, y, s)
ax.imshow(img, extent=extent, origin='lower', cmap=cm.jet)
ax.set_title("Smoothing with $\sigma$ = %d" % s)
plt.show()
生产:
Agape Gal'lo的散点图和s=16相互叠加(点击查看更好的视图):
我注意到我的高斯滤波方法和亚历杭德罗的方法的一个区别是，他的方法显示局部结构比我的好得多。因此，我在像素级上实现了一个简单的最近邻方法。该方法为每个像素计算数据中n个最近点距离的逆和。这种方法的分辨率很高，计算成本很高，我认为有更快的方法，所以如果你有任何改进，请告诉我。
更新:正如我所怀疑的，有一个更快的方法使用Scipy的Scipy . ckdtree。关于实现，请参阅Gabriel的回答。
总之，这是我的代码:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
def data_coord2view_coord(p, vlen, pmin, pmax):
dp = pmax - pmin
dv = (p - pmin) / dp * vlen
return dv
def nearest_neighbours(xs, ys, reso, n_neighbours):
im = np.zeros([reso, reso])
extent = [np.min(xs), np.max(xs), np.min(ys), np.max(ys)]
xv = data_coord2view_coord(xs, reso, extent[0], extent[1])
yv = data_coord2view_coord(ys, reso, extent[2], extent[3])
for x in range(reso):
for y in range(reso):
xp = (xv - x)
yp = (yv - y)
d = np.sqrt(xp**2 + yp**2)
im[y][x] = 1 / np.sum(d[np.argpartition(d.ravel(), n_neighbours)[:n_neighbours]])
return im, extent
n = 1000
xs = np.random.randn(n)
ys = np.random.randn(n)
resolution = 250
fig, axes = plt.subplots(2, 2)
for ax, neighbours in zip(axes.flatten(), [0, 16, 32, 64]):
if neighbours == 0:
ax.plot(xs, ys, 'k.', markersize=2)
ax.set_aspect('equal')
ax.set_title("Scatter Plot")
else:
im, extent = nearest_neighbours(xs, ys, resolution, neighbours)
ax.imshow(im, origin='lower', extent=extent, cmap=cm.jet)
ax.set_title("Smoothing over %d neighbours" % neighbours)
ax.set_xlim(extent[0], extent[1])
ax.set_ylim(extent[2], extent[3])
plt.show()
结果:
非常类似于@Piti的答案，但使用1次调用而不是2次调用来生成点:
import numpy as np
import matplotlib.pyplot as plt
pts = 1000000
mean = [0.0, 0.0]
cov = [[1.0,0.0],[0.0,1.0]]
x,y = np.random.multivariate_normal(mean, cov, pts).T
plt.hist2d(x, y, bins=50, cmap=plt.cm.jet)
plt.show()
输出:
推荐文章
- 为什么我得到“'str'对象没有属性'读取'”当尝试使用' json。载入字符串?
- 不区分大小写的列表排序，没有降低结果?
- 排序后的语法(key=lambda:…)
- 在烧瓶中返回HTTP状态代码201
- 使用python创建一个简单的XML文件
- APT命令行界面式的yes/no输入?
- 如何打印出状态栏和百分比?
- 在Python中获取大文件的MD5哈希值
- 在Python格式字符串中%s是什么意思?
- 如何循环通过所有但最后一项的列表?
- python用什么方法避免默认参数为空列表?
- ValueError: numpy。Ndarray大小改变，可能表示二进制不兼容。期望从C头得到88，从PyObject得到80
- Anaconda /conda -安装特定的软件包版本
- 我在哪里调用Keras的BatchNormalization函数?
- 打印测试执行时间并使用py.test锁定缓慢的测试
