Как построить 3D-карту плотности на Python с помощью matplotlib

У меня есть большой набор данных о положениях белков (x, y, z), и я хотел бы построить области с высокой загруженностью в виде тепловой карты. В идеале результат должен быть похож на объемную визуализацию ниже, но я не уверен, как этого добиться с помощью matplotlib.

http://i.stack.imgur.com/nsNEL.jpg

My initial idea was to display my positions as a 3D scatter plot and color their density via a KDE. I coded this up as follows with test data:

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 и графика рассеяния становится чрезвычайно медленным.

Небольшой образец моих реальных данных:

http://i.stack.imgur.com/BFT5V.png

Мои исследования показывают, что лучшим вариантом является оценка гауссовского 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() 

person nv_wu    schedule 13.08.2014    source источник
comment
Для этого приложения, я думаю, вам может быть лучше использовать mayavi, который более эффективен для приложений 3D-визуализации. Вот пример из документации, которая должна вам помочь начал.   -  person mwaskom    schedule 13.08.2014


Ответы (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.

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()
person nv_wu    schedule 17.08.2014
comment
Отлично! Вы знаете, как изобразить это, используя цветовую карту, отличную от струйной? - person crypdick; 30.04.2016
comment
@RovingRichard. Используя аргумент ключевого слова color = в вызове Volume, вы можете изменить его на один цвет. Более сложные цветовые карты могут быть выполнены путем создания ColorTransferFunction, как в примере на docs.enoughtt.com/mayavi/mayavi/auto/ - person coderforlife; 05.05.2016
comment
Я не могу заставить это работать с python 3.5.4, есть ли способ заставить это работать с matplotlib, seaborn или bokeh? - person O.rka; 13.01.2017
comment
Mayavi использует привязки VTK C ++ под капотом, а последняя в настоящее время не поддерживает Python 3. Я обнаружил, что matplotlib - плохое решение для изоповерхностей или чего-либо сложного в 3D. Я уверен, что это можно было бы переписать с помощью Brohh. Большинство графических библиотек на основе JavaScript / браузера обладают впечатляющими возможностями рендеринга. - person nv_wu; 15.01.2017
comment
У вас есть и пример такого рода сюжетов, когда вы рассматриваете только слайд сюжета? (Например, график с z = 0) - person Cruz; 17.10.2020