我有一组X,Y数据点(大约10k),很容易绘制成散点图,但我想用热图来表示。

我查看了Matplotlib中的示例,它们似乎都已经从热图单元值开始生成图像。

有没有一种方法可以将一堆不同的x, y转换为热图(其中x, y频率较高的区域会“更温暖”)?


当前回答

这些解决方案都不适用于我的应用程序,所以我想出了这个解决方案。本质上,我在每个点上都放置了一个二维高斯分布:

import cv2
import numpy as np
import matplotlib.pyplot as plt

def getGaussian2D(ksize, sigma, norm=True):
    oneD = cv2.getGaussianKernel(ksize=ksize, sigma=sigma)
    twoD = np.outer(oneD.T, oneD)
    return twoD / np.sum(twoD) if norm else twoD

def pt2heat(pts, shape, kernel=16, sigma=5):
    heat = np.zeros(shape)
    k = getGaussian2D(kernel, sigma)
    for y,x in pts:
        x, y = int(x), int(y)
        for i in range(-kernel//2, kernel//2):
            for j in range(-kernel//2, kernel//2):
                if 0 <= x+i < shape[0] and 0 <= y+j < shape[1]:
                    heat[x+i, y+j] = heat[x+i, y+j] + k[i+kernel//2, j+kernel//2]
    return heat


heat = pts2heat(pts, img.shape[:2])
plt.imshow(heat, cmap='heat')

以下是在相关图像上叠加的点,以及生成的热图:

其他回答

下面是我在100万个点集上做的一个,有3个类别(红色、绿色和蓝色)。如果您想尝试这个功能,这里有一个到存储库的链接。Github回购

histplot(
    X,
    Y,
    labels,
    bins=2000,
    range=((-3,3),(-3,3)),
    normalize_each_label=True,
    colors = [
        [1,0,0],
        [0,1,0],
        [0,0,1]],
    gain=50)

编辑:为了更好地近似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()

结果:

这些解决方案都不适用于我的应用程序,所以我想出了这个解决方案。本质上,我在每个点上都放置了一个二维高斯分布:

import cv2
import numpy as np
import matplotlib.pyplot as plt

def getGaussian2D(ksize, sigma, norm=True):
    oneD = cv2.getGaussianKernel(ksize=ksize, sigma=sigma)
    twoD = np.outer(oneD.T, oneD)
    return twoD / np.sum(twoD) if norm else twoD

def pt2heat(pts, shape, kernel=16, sigma=5):
    heat = np.zeros(shape)
    k = getGaussian2D(kernel, sigma)
    for y,x in pts:
        x, y = int(x), int(y)
        for i in range(-kernel//2, kernel//2):
            for j in range(-kernel//2, kernel//2):
                if 0 <= x+i < shape[0] and 0 <= y+j < shape[1]:
                    heat[x+i, y+j] = heat[x+i, y+j] + k[i+kernel//2, j+kernel//2]
    return heat


heat = pts2heat(pts, img.shape[:2])
plt.imshow(heat, cmap='heat')

以下是在相关图像上叠加的点,以及生成的热图:

下面是Jurgy使用scipy.cKDTree实现的最近邻方法。在我的测试中,它快了大约100倍。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from scipy.spatial import cKDTree


def data_coord2view_coord(p, resolution, pmin, pmax):
    dp = pmax - pmin
    dv = (p - pmin) / dp * resolution
    return dv


n = 1000
xs = np.random.randn(n)
ys = np.random.randn(n)

resolution = 250

extent = [np.min(xs), np.max(xs), np.min(ys), np.max(ys)]
xv = data_coord2view_coord(xs, resolution, extent[0], extent[1])
yv = data_coord2view_coord(ys, resolution, extent[2], extent[3])


def kNN2DDens(xv, yv, resolution, neighbours, dim=2):
    """
    """
    # Create the tree
    tree = cKDTree(np.array([xv, yv]).T)
    # Find the closest nnmax-1 neighbors (first entry is the point itself)
    grid = np.mgrid[0:resolution, 0:resolution].T.reshape(resolution**2, dim)
    dists = tree.query(grid, neighbours)
    # Inverse of the sum of distances to each grid point.
    inv_sum_dists = 1. / dists[0].sum(1)

    # Reshape
    im = inv_sum_dists.reshape(resolution, resolution)
    return im


fig, axes = plt.subplots(2, 2, figsize=(15, 15))
for ax, neighbours in zip(axes.flatten(), [0, 16, 32, 63]):

    if neighbours == 0:
        ax.plot(xs, ys, 'k.', markersize=5)
        ax.set_aspect('equal')
        ax.set_title("Scatter Plot")
    else:

        im = kNN2DDens(xv, yv, resolution, neighbours)

        ax.imshow(im, origin='lower', extent=extent, cmap=cm.Blues)
        ax.set_title("Smoothing over %d neighbours" % neighbours)
        ax.set_xlim(extent[0], extent[1])
        ax.set_ylim(extent[2], extent[3])

plt.savefig('new.png', dpi=150, bbox_inches='tight')

如果你不想要六边形,你可以使用numpy的histogram2d函数:

import numpy as np
import numpy.random
import matplotlib.pyplot as plt

# Generate some test data
x = np.random.randn(8873)
y = np.random.randn(8873)

heatmap, xedges, yedges = np.histogram2d(x, y, bins=50)
extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]

plt.clf()
plt.imshow(heatmap.T, extent=extent, origin='lower')
plt.show()

这是一个50x50的热图。如果你想要,比如说512x384,你可以在调用histogram2d时放入bins=(512,384)。

例子: