如何使用 matplotlib 在 python 中绘制 3D 密度图

Posted 2023-02-14

【中文标题】如何使用 matplotlib 在 python 中绘制 3D 密度图

【英文标题】：How to plot a 3D density map in python with matplotlib

【发布时间】：2014-10-06 20:25:44

【问题描述】：

我有一个 (x,y,z) 蛋白质位置的大型数据集，并希望将高占用率区域绘制为热图。理想情况下，输出应该类似于下面的体积可视化，但我不确定如何使用 matplotlib 实现这一点。

我最初的想法是将我的位置显示为 3D 散点图，并通过 KDE 为它们的密度着色。我用测试数据将其编码如下：

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D 

mu, sigma = 0, 0.1 
x = np.random.normal(mu, sigma, 1000)
y = np.random.normal(mu, sigma, 1000)
z = np.random.normal(mu, sigma, 1000)

xyz = np.vstack([x,y,z])
density = stats.gaussian_kde(xyz)(xyz) 

idx = density.argsort()
x, y, z, density = x[idx], y[idx], z[idx], density[idx]

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, c=density)
plt.show()

这很好用！但是，我的真实数据包含数千个数据点，计算 kde 和散点图变得非常缓慢。

我的真实数据的一个小样本：

我的研究表明，更好的选择是在网格上评估高斯 kde。我只是不确定如何在 3D 中做到这一点：

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D 

mu, sigma = 0, 0.1 
x = np.random.normal(mu, sigma, 1000)
y = np.random.normal(mu, sigma, 1000)

nbins = 50

xy = np.vstack([x,y])
density = stats.gaussian_kde(xy) 

xi, yi = np.mgrid[x.min():x.max():nbins*1j, y.min():y.max():nbins*1j]
di = density(np.vstack([xi.flatten(), yi.flatten()]))

fig = plt.figure()
ax = fig.add_subplot(111)
ax.pcolormesh(xi, yi, di.reshape(xi.shape))
plt.show() 

【问题讨论】：

对于这个应用程序，我认为您最好使用 mayavi，它对于 3D 可视化应用程序更强大。这是文档中的example，应该可以帮助您入门。

【参考方案1】：

感谢 mwaskon 推荐 mayavi 库。

我在 mayavi 中重新创建了密度散点图，如下所示：

import numpy as np
from scipy import stats
from mayavi import mlab

mu, sigma = 0, 0.1 
x = 10*np.random.normal(mu, sigma, 5000)
y = 10*np.random.normal(mu, sigma, 5000)
z = 10*np.random.normal(mu, sigma, 5000)

xyz = np.vstack([x,y,z])
kde = stats.gaussian_kde(xyz)
density = kde(xyz)

# Plot scatter with mayavi
figure = mlab.figure('DensityPlot')
pts = mlab.points3d(x, y, z, density, scale_mode='none', scale_factor=0.07)
mlab.axes()
mlab.show()

将 scale_mode 设置为“none”可防止字形与密度向量成比例缩放。此外，对于大型数据集，我禁用了场景渲染并使用遮罩来减少点数。

# Plot scatter with mayavi
figure = mlab.figure('DensityPlot')
figure.scene.disable_render = True

pts = mlab.points3d(x, y, z, density, scale_mode='none', scale_factor=0.07) 
mask = pts.glyph.mask_points
mask.maximum_number_of_points = x.size
mask.on_ratio = 1
pts.glyph.mask_input_points = True

figure.scene.disable_render = False 
mlab.axes()
mlab.show()

接下来，评估网格上的高斯 kde：

import numpy as np
from scipy import stats
from mayavi import mlab

mu, sigma = 0, 0.1 
x = 10*np.random.normal(mu, sigma, 5000)
y = 10*np.random.normal(mu, sigma, 5000)    
z = 10*np.random.normal(mu, sigma, 5000)

xyz = np.vstack([x,y,z])
kde = stats.gaussian_kde(xyz)

# Evaluate kde on a grid
xmin, ymin, zmin = x.min(), y.min(), z.min()
xmax, ymax, zmax = x.max(), y.max(), z.max()
xi, yi, zi = np.mgrid[xmin:xmax:30j, ymin:ymax:30j, zmin:zmax:30j]
coords = np.vstack([item.ravel() for item in [xi, yi, zi]]) 
density = kde(coords).reshape(xi.shape)

# Plot scatter with mayavi
figure = mlab.figure('DensityPlot')

grid = mlab.pipeline.scalar_field(xi, yi, zi, density)
min = density.min()
max=density.max()
mlab.pipeline.volume(grid, vmin=min, vmax=min + .5*(max-min))

mlab.axes()
mlab.show()

作为最后的改进，我通过并行调用 kde 函数加快了对 kensity 密度函数的评估。

import numpy as np
from scipy import stats
from mayavi import mlab
import multiprocessing

def calc_kde(data):
    return kde(data.T)

mu, sigma = 0, 0.1 
x = 10*np.random.normal(mu, sigma, 5000)
y = 10*np.random.normal(mu, sigma, 5000)
z = 10*np.random.normal(mu, sigma, 5000)

xyz = np.vstack([x,y,z])
kde = stats.gaussian_kde(xyz)

# Evaluate kde on a grid
xmin, ymin, zmin = x.min(), y.min(), z.min()
xmax, ymax, zmax = x.max(), y.max(), z.max()
xi, yi, zi = np.mgrid[xmin:xmax:30j, ymin:ymax:30j, zmin:zmax:30j]
coords = np.vstack([item.ravel() for item in [xi, yi, zi]]) 

# Multiprocessing
cores = multiprocessing.cpu_count()
pool = multiprocessing.Pool(processes=cores)
results = pool.map(calc_kde, np.array_split(coords.T, 2))
density = np.concatenate(results).reshape(xi.shape)

# Plot scatter with mayavi
figure = mlab.figure('DensityPlot')

grid = mlab.pipeline.scalar_field(xi, yi, zi, density)
min = density.min()
max=density.max()
mlab.pipeline.volume(grid, vmin=min, vmax=min + .5*(max-min))

mlab.axes()
mlab.show()

【讨论】：

这太棒了！你知道如何使用不是 jet 的颜色图来绘制吗？

@RovingRichard 通过在对音量的调用中使用关键字参数 color=，您可以将其更改为单一颜色。可以通过创建一个 ColorTransferFunction 来完成更高级的颜色映射，如docs.enthought.com/mayavi/mayavi/auto/…的示例中所示

我无法让它与 python 3.5.4 一起使用，是否有办法让它与 matplotlib、seaborn 或 bokeh 一起使用？

Mayavi 在后台使用 VTK C++ 绑定，而后者目前不支持 Python 3。我发现 matplotlib 对等值面或 3D 中的任何复杂事物都不是一个糟糕的解决方案。我敢肯定它可以用breakh重写。大多数基于 javascript/浏览器的图形库都具有令人印象深刻的渲染能力。

当您只考虑情节的幻灯片时，您有这种情节的例子吗？ （例如 z=0 的图）