Genere un mapa de calor usando un conjunto de datos dispersos

14 minutos de lectura

avatar de usuario
gris

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

Revisé los ejemplos en Matplotlib y todos parecen comenzar con valores de celda de mapa de calor para generar la imagen.

¿Existe 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”)?

  • Igualmente relevante: método eficiente para calcular la densidad de puntos espaciados irregularmente

    – La Importancia De Ser Ernesto

    23/07/2019 a las 15:00

avatar de usuario
tomate

Si no quieres hexágonos, puedes usar numpy’s histogram2d función:

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 quieres, digamos, 512×384, puedes poner bins=(512, 384) en la llamada a histogram2d.

Ejemplo: Ejemplo de mapa de calor de Matplotlib

  • No pretendo ser un idiota, pero ¿cómo tiene realmente esta salida en un archivo PNG/PDF en lugar de mostrarla solo en una sesión interactiva de IPython? Estoy tratando de conseguir esto como algo normal axes instancia, donde puedo agregar un título, etiquetas de eje, etc. y luego hacer lo normal savefig() como lo haría con cualquier otra trama típica de matplotlib.

    – gotgenes

    15 de julio de 2011 a las 19:19


  • @gotgenes: no plt.savefig('filename.png') ¿trabajar? Si desea obtener una instancia de ejes, use la interfaz orientada a objetos de Matplotlib: fig = plt.figure() ax = fig.gca() ax.imshow(...) fig.savefig(...)

    – tomate

    16 de julio de 2011 a las 17:05

  • De hecho, ¡gracias! Supongo que no entiendo completamente eso imshow() está en la misma categoría de funciones que scatter(). sinceramente no entiendo porque imshow() convierte una matriz 2d de flotadores en bloques de color apropiado, mientras que entiendo lo que scatter() se supone que tiene que ver con tal matriz.

    – gotgenes

    21 de julio de 2011 a las 19:10


  • Una advertencia sobre el uso de imshow para trazar un histograma 2d de valores x/y como este: por defecto, imshow traza el origen en la esquina superior izquierda y transpone la imagen. Lo que haría para obtener la misma orientación que un diagrama de dispersión es plt.imshow(heatmap.T, extent=extent, origin = 'lower')

    – Jaime

    18 de noviembre de 2013 a las 13:29

  • Para aquellos que quieran hacer una barra de colores logarítmica, consulte esta pregunta stackoverflow.com/questions/17201172/… y simplemente haga from matplotlib.colors import LogNorm plt.imshow(heatmap, norm=LogNorm()) plt.colorbar()

    – tommy.carstensen

    16/03/2015 a las 20:25


avatar de usuario
doug

En matplotlib léxico, creo que quieres un hexágono gráfico.

Si no está familiarizado con este tipo de trama, es solo una histograma bivariado en el que el plano xy está dividido por una cuadrícula regular de hexágonos.

Entonces, a partir de un histograma, puede simplemente contar la cantidad de puntos que caen en cada hexágono, discretizar la región de trazado como un conjunto de ventanas, asigne cada punto a una de estas ventanas; finalmente, asigne las ventanas a un gama de coloresy tienes un diagrama hexadecimal.

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

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

  • hexágono es el n-polígono más alto que da teselado plano regular (es decir, puede remodelar de manera segura el piso de su cocina con mosaicos de forma hexagonal porque no tendrá ningún espacio vacío entre los mosaicos cuando haya terminado, lo que no es cierto para todos los demás polígonos de n superior, n >= 7 ).

(matplotlib utiliza el término hexágono gráfico; también (AFAIK) todos los bibliotecas de trazado por R; Todavía no sé si este es el término generalmente aceptado para tramas de este tipo, aunque sospecho que es probable dado que hexágono es corto para agrupamiento hexagonalque 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 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()   

ingrese la descripción de la imagen aquí

  • ¿Qué significa que “los hexágonos tienen simetría de vecino más cercano”? Usted dice que “la distancia desde un punto en el borde de un cuadrado y un punto dentro de ese cuadrado no es igual en todas partes”, pero ¿distancia a qué?

    – Jan

    11 de abril de 2014 a las 16:04

  • Para un hexágono, la distancia desde el centro hasta un vértice que une dos lados también es más larga que desde el centro hasta la mitad de un lado, solo que la relación es menor (2/sqrt(3) ≈ 1,15 para hexágono vs. sqrt(2) ≈ 1,41 por cuadrado). La única forma donde la distancia desde el centro hasta cada punto del borde es igual es el círculo.

    – Jan

    25 mayo 2014 a las 18:46

  • @Jaan Para un hexágono, todos los vecinos están a la misma distancia. No hay problema con 8 vecindarios o 4 vecindarios. No hay vecinos diagonales, solo un tipo de vecino.

    – isarandi

    8 de marzo de 2015 a las 16:06

  • @doug ¿Cómo eliges el gridsize= parámetro. Me gustaría elegirlo así, de modo que los hexágonos se toquen sin superponerse. Me di cuenta que gridsize=100 produciría hexágonos más pequeños, pero ¿cómo elegir el valor adecuado?

    – Alejandro Cska

    19 de abril de 2016 a las 9:05


  • El problema con estas gráficas (como con las gráficas de algunas otras respuestas) es que no queda claro dónde apuntan los datos y dónde está el fondo vacío.

    – Igor apoya a Ucrania

    18 de febrero a las 12:29

avatar de usuario
jurado

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

Sé que esta es una vieja pregunta, pero quería agregar algo a la respuesta de Alejandro: si desea una buena imagen suavizada sin usar py-sphviewer, puede usar en su lugar np.histogram2d y aplicar un filtro gaussiano (desde scipy.ndimage.filters) al 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 gráfico 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 enfoque de Alejandro fue que su método muestra las estructuras locales mucho mejor que el mío. Por lo tanto, implementé un método simple de vecino más cercano a nivel de píxel. Este método calcula para cada píxel la suma inversa de las distancias de los n puntos más cercanos en los datos. Este método tiene una alta resolución y es bastante costoso desde el punto de vista computacional y creo que hay una forma más rápida, así que avíseme si tiene alguna mejora.

Actualización: como sospechaba, hay un método mucho más rápido usando Scipy’s scipy.cKDTree. Vea la respuesta de Gabriel para la implementación.

De todos modos, aquí está mi 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:

Suavizado del vecino más cercano

  • Me encanta esto. Graph es tan bueno como la respuesta de Alejandro, pero no se requieren nuevos paquetes.

    – Nathan Clemente

    30 de noviembre de 2017 a las 21:16

  • Muy agradable ! Pero generas una compensación con este método. Puedes ver esto comparando un gráfico de dispersión normal con el coloreado. ¿Podrías agregar algo para corregirlo? ¿O simplemente para mover el gráfico en valores de x e y?

    – Ágape Gal’lo

    6 sep 2018 a las 20:09

  • Ágape Gal’lo, ¿a qué te refieres con offset? Si los trazas uno encima del otro, coinciden (ver la edición de mi publicación). Tal vez te desanime porque el ancho de la dispersión no coincide exactamente con los otros tres.

    – Jurgia

    7 de septiembre de 2018 a las 7:29

  • ¡Muchas gracias por trazar el gráfico solo para mí! Entendí mi error: había modificado la “extensión” para definir los límites x e y. Ahora entiendo que modificó el origen del gráfico. Entonces, tengo una última pregunta: ¿cómo puedo expandir los límites del gráfico, incluso para el área donde no hay datos existentes? Por ejemplo, entre -5 y +5 para x e y.

    – Ágape Gal’lo

    8 de septiembre de 2018 a las 13:59

  • Digamos que quiere que el eje x vaya de -5 a 5 y el eje y de -3 a 4; en el myplot función, agregue el range parámetro a np.histogram2d: np.histogram2d(x, y, bins=bins, range=[[-5, 5], [-3, 4]]) y en el bucle for establezca los límites x e y del eje: ax.set_xlim([-5, 5]) ax.set_ylim([-3, 4]). Además, por defecto, imshow mantiene la relación de aspecto idéntica a la relación de sus ejes (en mi ejemplo, una relación de 10: 7), pero si desea que coincida con su ventana de trazado, agregue el parámetro aspect='auto' a imshow.

    – Jurgia

    10 de septiembre de 2018 a las 7:38


avatar de usuario
Alejandro

En lugar de usar np.hist2d, que en general produce histogramas bastante feos, me gustaría reciclar py-sphviewer, un paquete de python para renderizar simulaciones de partículas utilizando un núcleo de suavizado adaptativo 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([len(x),3])
    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:

ingrese la descripción de la imagen aquí

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

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

avatar de usuario
Piti Ongmongkolkul

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

avatar de usuario
palabras para el sabio

Seaborn ahora tiene la función de diagrama conjunto 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 de demostración

avatar de usuario
gabriel

Aquí está el enfoque del gran vecino más cercano de Jurgy pero implementado usando scipy.cKDTree. En mis pruebas, es aproximadamente 100 veces más rápido.

ingrese la descripción de la imagen aquí

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

  • Sabía que mi implementación era muy ineficiente pero no sabía sobre cKDTree. ¡Bien hecho! Te mencionaré en mi respuesta.

    – Jurgia

    18 de marzo de 2020 a las 10:00

¿Ha sido útil esta solución?