Gere um mapa de calor no MatPlotLib usando um conjunto de dados de dispersão

Eu tenho um conjunto de pontos de dados X, Y (cerca de 10k) que são fáceis de plotar como um gráfico de dispersão, mas que eu gostaria de representar como um mapa de calor.

Examinei os exemplos no MatPlotLib e todos parecem já começar com os valores das células do mapa de calor para gerar a imagem.

Existe um método que converte um monte de x, y, todos diferentes, em um mapa de calor (onde zonas com maior frequência de x, y seriam "mais quentes")?

Se você não deseja hexágonos, pode usar a histogram2dfunção numpy :

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()

Isso cria um mapa de calor de 50x50. Se você quiser, digamos, 512x384, pode fazer bins=(512, 384)a chamada para histogram2d.

Exemplo:

No léxico Matplotlib , acho que você quer uma plotagem hexbin .

Se você não estiver familiarizado com esse tipo de gráfico, é apenas um histograma bivariado no qual o plano xy é mosaico por uma grade regular de hexágonos.

Assim, a partir de um histograma, você pode apenas contar o número de pontos que caem em cada hexágono, discretizar a região de plotagem como um conjunto de janelas , atribuir cada ponto a uma dessas janelas; finalmente, mapeie as janelas em uma matriz de cores e você terá um diagrama hexbin.

Embora menos comumente usado do que, por exemplo, círculos ou quadrados, os hexágonos sejam a melhor opção para a geometria do recipiente de armazenamento é intuitivo:

• hexágonos têm simetria do vizinho mais próximo (por exemplo, caixas quadradas não, por exemplo, a distância de um ponto na borda de um quadrado até um ponto dentro desse quadrado não é igual em todo lugar) e

• hexágono é o polígono n mais alto que fornece mosaico plano regular (ou seja, você pode remodelar com segurança o piso da cozinha com ladrilhos em formato hexagonal, porque você não terá espaço vazio entre os ladrilhos quando terminar - não é verdade para todos os outros polígonos com n superior, n> = 7).

( Matplotlib usa o termo plotagem de hexbin ; assim como (AFAIK) todas as bibliotecas de plotagem para R ; ainda não sei se esse é o termo geralmente aceito para plotagens desse tipo, embora eu suspeite que seja provável, dado que a hexbin é curta para bin hexagonal , que descreve a etapa essencial na preparação dos dados para exibição.)

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()   

Edit: Para uma melhor aproximação da resposta de Alejandro, veja abaixo.

Sei que essa é uma pergunta antiga, mas queria acrescentar algo ao analisador de Alejandro: Se você deseja uma boa imagem suavizada sem usar o py-sphviewer, pode usar np.histogram2de aplicar um filtro gaussiano (de scipy.ndimage.filters) no mapa de calor:

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()

Produz:

O gráfico de dispersão es = 16 plotados em cima um do outro para Agape Gal'lo (clique para ver melhor):

Uma diferença que notei com minha abordagem de filtro gaussiano e a abordagem de Alejandro foi que o método dele mostra estruturas locais muito melhores que as minhas. Portanto, implementei um método vizinho mais próximo simples no nível de pixel. Este método calcula para cada pixel a soma inversa das distâncias dos npontos mais próximos nos dados. Esse método está em uma alta resolução bastante caro em termos de computação e acho que há uma maneira mais rápida; então, deixe-me saber se você tem alguma melhoria.

Atualização: Como eu suspeitava, existe um método muito mais rápido usando o Scipy's scipy.cKDTree. Vejo a resposta de Gabriel para a implementação.

Enfim, aqui está o meu código:

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()

Resultado:

Em vez de usar o np.hist2d, que geralmente produz histogramas bastante feios, eu gostaria de reciclar o py-sphviewer , um pacote python para renderizar simulações de partículas usando um kernel de suavização adaptável e que pode ser facilmente instalado a partir do pip (consulte a documentação da página da web). Considere o seguinte código, que é baseado no exemplo:

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([3, len(x)])
    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()

que produz a seguinte imagem:

Como você vê, as imagens parecem muito boas e somos capazes de identificar diferentes subestruturas nela. Essas imagens são construídas espalhando um determinado peso para cada ponto dentro de um determinado domínio, definido pelo comprimento de suavização, que por sua vez é dado pela distância do vizinho nb mais próximo (eu escolhi 16, 32 e 64 para os exemplos). Portanto, regiões de alta densidade geralmente estão espalhadas por regiões menores em comparação com regiões de menor densidade.

A função myplot é apenas uma função muito simples que eu escrevi para fornecer os dados x, y ao py-sphviewer para fazer a mágica.

Se você estiver usando o 1.2.x

import numpy as np
import matplotlib.pyplot as plt

x = np.random.randn(100000)
y = np.random.randn(100000)
plt.hist2d(x,y,bins=100)
plt.show()

O Seaborn agora tem a função de plotagem conjunta, que deve funcionar bem aqui:

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

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

sns.jointplot(x=x, y=y, kind='hex')
plt.show()

e a pergunta inicial era ... como converter valores de dispersão em valores de grade, certo? histogram2dconta a frequência por célula; no entanto, se você tiver outros dados por célula além da frequência, precisará de mais trabalho.

x = data_x # between -10 and 4, log-gamma of an svc
y = data_y # between -4 and 11, log-C of an svc
z = data_z #between 0 and 0.78, f1-values from a difficult dataset

Então, eu tenho um conjunto de dados com resultados Z para as coordenadas X e Y. No entanto, eu estava calculando alguns pontos fora da área de interesse (grandes lacunas) e montes de pontos em uma pequena área de interesse.

Sim, aqui fica mais difícil, mas também mais divertido. Algumas bibliotecas (desculpe):

from matplotlib import pyplot as plt
from matplotlib import cm
import numpy as np
from scipy.interpolate import griddata

Hoje, o pyplot é meu mecanismo gráfico, cm é uma variedade de mapas de cores com algumas opções interessantes. numpy para os cálculos e dados de grade para anexar valores a uma grade fixa.

O último é importante, especialmente porque a frequência dos pontos xy não é igualmente distribuída nos meus dados. Primeiro, vamos começar com alguns limites adequados aos meus dados e um tamanho de grade arbitrário. Os dados originais possuem pontos de dados também fora desses limites x e y.

#determine grid boundaries
gridsize = 500
x_min = -8
x_max = 2.5
y_min = -2
y_max = 7

Portanto, definimos uma grade com 500 pixels entre os valores mínimo e máximo de x e y.

Nos meus dados, existem muito mais que os 500 valores disponíveis na área de alto interesse; considerando que na área de baixo interesse não existem nem 200 valores na grade total; entre os limites gráficos dex_min e x_maxhá ainda menos.

Portanto, para obter uma boa imagem, a tarefa é obter uma média dos altos valores de juros e preencher as lacunas em outros lugares.

Eu defino minha grade agora. Para cada par xx-yy, quero ter uma cor.

xx = np.linspace(x_min, x_max, gridsize) # array of x values
yy = np.linspace(y_min, y_max, gridsize) # array of y values
grid = np.array(np.meshgrid(xx, yy.T))
grid = grid.reshape(2, grid.shape[1]*grid.shape[2]).T

Por que a forma estranha? scipy.griddata quer uma forma de (n, D).

Griddata calcula um valor por ponto na grade, por um método predefinido. Eu escolhi "mais próximo" - pontos de grade vazios serão preenchidos com valores do vizinho mais próximo. Parece que as áreas com menos informações têm células maiores (mesmo que não seja o caso). Pode-se optar por interpolar "linear", então áreas com menos informações parecem menos nítidas. Questão de gosto, realmente.

points = np.array([x, y]).T # because griddata wants it that way
z_grid2 = griddata(points, z, grid, method='nearest')
# you get a 1D vector as result. Reshape to picture format!
z_grid2 = z_grid2.reshape(xx.shape[0], yy.shape[0])

E pulamos, entregamos ao matplotlib para exibir o gráfico

fig = plt.figure(1, figsize=(10, 10))
ax1 = fig.add_subplot(111)
ax1.imshow(z_grid2, extent=[x_min, x_max,y_min, y_max,  ],
            origin='lower', cmap=cm.magma)
ax1.set_title("SVC: empty spots filled by nearest neighbours")
ax1.set_xlabel('log gamma')
ax1.set_ylabel('log C')
plt.show()

Em torno da parte pontiaguda da V-Shape, você vê que fiz muitos cálculos durante minha busca pelo ponto ideal, enquanto as partes menos interessantes em quase todos os lugares têm uma resolução mais baixa.

Aqui está a melhor abordagem de vizinho mais próximo de Jurgy, mas implementada usando o scipy.cKDTree . Nos meus testes, é cerca de 100x mais rápido.

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')

Faça uma matriz bidimensional que corresponda às células da sua imagem final, chamada diga heatmap_cellse instancie-a como todos os zeros.

Escolha dois fatores de escala que definem a diferença entre cada elemento da matriz em unidades reais, para cada dimensão, digamos x_scaleey_scale . Escolha esses itens de forma que todos os seus pontos de dados caiam dentro dos limites da matriz do mapa de calor.

Para cada ponto de dados bruto com x_valuee y_value:

heatmap_cells[floor(x_value/x_scale),floor(y_value/y_scale)]+=1

Aqui está um que fiz em um milhão de pontos definido com 3 categorias (vermelho, verde e azul). Aqui está um link para o repositório, se você quiser experimentar a função. Repositório do 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)

Muito semelhante à resposta de @ Piti , mas usando 1 chamada em vez de 2 para gerar os pontos:

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()

Resultado:

Receio estar um pouco atrasado para a festa, mas tive uma pergunta semelhante há um tempo atrás. A resposta aceita (por @ptomato) me ajudou, mas eu também gostaria de postar isso, caso seja útil para alguém.

''' I wanted to create a heatmap resembling a football pitch which would show the different actions performed '''

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

#fixing random state for reproducibility
np.random.seed(1234324)

fig = plt.figure(12)
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

#Ratio of the pitch with respect to UEFA standards 
hmap= np.full((6, 10), 0)
#print(hmap)

xlist = np.random.uniform(low=0.0, high=100.0, size=(20))
ylist = np.random.uniform(low=0.0, high =100.0, size =(20))

#UEFA Pitch Standards are 105m x 68m
xlist = (xlist/100)*10.5
ylist = (ylist/100)*6.5

ax1.scatter(xlist,ylist)

#int of the co-ordinates to populate the array
xlist_int = xlist.astype (int)
ylist_int = ylist.astype (int)

#print(xlist_int, ylist_int)

for i, j in zip(xlist_int, ylist_int):
    #this populates the array according to the x,y co-ordinate values it encounters 
    hmap[j][i]= hmap[j][i] + 1   

#Reversing the rows is necessary 
hmap = hmap[::-1]

#print(hmap)
im = ax2.imshow(hmap)

Aqui está o resultado