Generar un mapa de calor en MatPlotLib usando un conjunto de datos de dispersión

Tengo un conjunto de puntos de datos X, Y (alrededor de 10k) que son fáciles de trazar como un diagtwig de dispersión pero que me gustaría representar como un mapa de calor.

Miré los ejemplos en MatPlotLib y todos parecen comenzar con valores de celdas de mapa de calor para generar la imagen.

¿Hay algún método que convierta un grupo de x, y, todos diferentes, en un mapa de calor (donde las zonas con mayor frecuencia de x, y serían “más cálidas”)?

Si no quieres hexágonos, puedes usar la función histogram2d de 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() 

Esto hace un mapa de calor de 50×50. Si desea, digamos, 512×384, puede poner bins=(512, 384) en la llamada a histogram2d .

Ejemplo: Ejemplo de mapa de calor de Matplotlib

En el léxico de Matplotlib , creo que quieres una gráfica de hexbin .

Si no estás familiarizado con este tipo de gráfico, es solo un histogtwig bivariado en el que el plano xy está teselado por una cuadrícula regular de hexágonos.

Entonces, a partir de un histogtwig, puede simplemente contar el número de puntos que caen en cada hexágono, discretizar la región de trazado como un conjunto de ventanas , asignar cada punto a una de estas ventanas; finalmente, asigne las ventanas a una matriz de colores , y tendrá un diagtwig de hexbin.

Aunque se usa con menos frecuencia que, por ejemplo, los círculos o los cuadrados, los hexágonos son una mejor opción, ya que la geometría del contenedor de agrupación es intuitiva:

  • los hexágonos tienen simetría del vecino más cercano (p. ej., los cubos cuadrados no, p. ej., la distancia desde un punto en el borde de un cuadrado a un punto dentro de ese cuadrado no es igual en todas partes) y

  • el hexágono es el polígono n más alto que proporciona una teselación de plano regular (es decir, puede volver a modelar de manera segura el piso de su cocina con azulejos de forma hexagonal porque no tendrá ningún espacio vacío entre los azulejos cuando haya terminado, no es cierto para todos los demás polos superiores a n, n> = 7).

( Matplotlib usa el término gráfico de hexbin ; también lo hace (AFAIK) todas las bibliotecas de trazado para R ; aún no sé si este es el término generalmente aceptado para los gráficos de este tipo, aunque sospecho que es posible dado que hexbin es corto para el agrupamiento hexagonal , que describe el paso esencial en la preparación de los datos para su visualización.)


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

introduzca la descripción de la imagen aquí

En lugar de usar np.hist2d, que en general produce histogtwigs bastante feos, me gustaría reciclar py-sphviewer , un paquete de python para representar simulaciones de partículas utilizando un núcleo de suavizado adaptable y que se puede instalar fácilmente desde pip (consulte la documentación de la página web). Considere el siguiente código, que se basa en el ejemplo:

 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 produce la siguiente imagen:

introduzca la descripción de la imagen aquí

Como puede ver, las imágenes se ven muy bien y podemos identificar diferentes subestructuras en ellas. Estas imágenes se construyen extendiendo un peso determinado para cada punto dentro de un determinado dominio, definido por la longitud de suavizado, que a su vez viene dada por la distancia al vecino nb más cercano (he elegido 16, 32 y 64 para los ejemplos). Por lo tanto, las regiones de mayor densidad normalmente se extienden en regiones más pequeñas en comparación con las regiones de menor densidad.

La función myplot es solo una función muy simple que he escrito para dar los datos x, y a py-sphviewer para hacer la magia.

Si está utilizando 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() 

gaussian_2d_heat_map

Edición: para una mejor aproximación de la respuesta de Alejandro, vea a continuación.

Sé que esta es una pregunta antigua, pero quería agregar algo a la respuesta de Alejandro: si desea una buena imagen suavizada sin usar py-sphviewer, puede usar np.histogram2d y aplicar un filtro gaussiano (desde scipy.ndimage.filters ) a el 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() 

Produce:

Imágenes de salida

El diagtwig de dispersión y s = 16 trazados uno encima del otro para Agape Gal’lo (haga clic para ver mejor):

Encima del otro


Una diferencia que noté con mi enfoque de filtro gaussiano y el de Alejandro fue que su método muestra las estructuras locales mucho mejor que las mías. Por lo tanto, implementé un método vecino más cercano a nivel de píxeles. Este método calcula para cada píxel la sum inversa de las distancias de los n puntos más cercanos en los datos. Este método es de alta resolución y es computacionalmente costoso y creo que hay una manera más rápida, así que avíseme si tiene alguna mejora. De todos modos, aquí está el 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:

Vecino más cercano alisado

Seaborn ahora tiene la función jointplot que debería funcionar bien aquí:

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

imagen demo

Cree una matriz bidimensional que se corresponda con las celdas de su imagen final, llamada say heatmap_cells y heatmap_cells instancia como todos los ceros.

Elija dos factores de escala que definan la diferencia entre cada elemento de matriz en unidades reales, para cada dimensión, por ejemplo, x_scale y y_scale . Elija estos para que todos sus puntos de datos queden dentro de los límites de la matriz de mapa de calor.

Para cada punto de datos sin procesar con x_value y y_value :

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

y la pregunta inicial era … ¿cómo convertir valores de dispersión a valores de cuadrícula, verdad? histogram2d cuenta la frecuencia por celda, sin embargo, si tiene otros datos por celda que solo la frecuencia, necesitará algún trabajo adicional para realizar.

 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 

Entonces, tengo un conjunto de datos con los resultados Z para las coordenadas X e Y. Sin embargo, estaba calculando algunos puntos fuera del área de interés (grandes brechas), y montones de puntos en un área pequeña de interés.

Sí, aquí se vuelve más difícil pero también más divertido. Algunas bibliotecas (lo siento):

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

Pyplot es mi motor gráfico hoy en día, cm es una gama de mapas en color con algunas opciones interesantes. numpy para los cálculos y griddata para adjuntar valores a una cuadrícula fija.

El último es importante, especialmente porque la frecuencia de los puntos xy no se distribuye por igual en mis datos. Primero, comencemos con algunos límites que se ajusten a mis datos y un tamaño de cuadrícula arbitrario. Los datos originales tienen puntos de datos también fuera de esos límites x e y.

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

Así que hemos definido una cuadrícula con 500 píxeles entre los valores mínimo y máximo de x e y.

En mis datos, hay mucho más que los 500 valores disponibles en el área de alto interés; mientras que en el área de bajo interés, no hay ni siquiera 200 valores en la cuadrícula total; entre los límites gráficos de x_min y x_max hay incluso menos.

Entonces, para obtener una buena imagen, la tarea es obtener un promedio de los valores de interés alto y llenar los vacíos en otros lugares.

Defino mi cuadrícula ahora. Para cada par xx-yy, quiero tener un color.

 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 qué la forma extraña? scipy.griddata quiere una forma de (n, D).

Griddata calcula un valor por punto en la cuadrícula, por un método predefinido. Elijo “más cercano”: los puntos de la cuadrícula vacíos se llenarán con los valores del vecino más cercano. Esto parece que las áreas con menos información tienen celdas más grandes (incluso si no es el caso). Uno podría elegir interpolar “lineal”, luego las áreas con menos información se ven menos definidas. La cuestión del gusto, de verdad.

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

Y salto, entregamos a matplotlib para mostrar la ttwig.

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

Alrededor de la parte puntiaguda de la Forma en V, ves que hice muchos cálculos durante mi búsqueda del punto dulce, mientras que las partes menos interesantes en casi todas partes tienen una resolución más baja.

Heatmap de un SVC en alta resolución.

Muy similar a la respuesta de @ Piti , pero usando 1 llamada en lugar de 2 para generar los puntos:

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

Salida:

2d_gaussian_heatmap