Generieren Sie eine Heatmap in MatPlotLib mithilfe eines Streudatensatzes

Ich habe eine Reihe von X- und Y-Datenpunkten (ca. 10.000), die sich leicht als Streudiagramm darstellen lassen, die ich aber gerne als Heatmap darstellen möchte.

Ich habe die Beispiele in MatPlotLib durchgesehen und sie scheinen alle bereits mit Heatmap-Zellenwerten zu beginnen, um das Bild zu generieren.

Gibt es eine Methode, die eine Reihe von x, y, die alle unterschiedlich sind, in eine Heatmap konvertiert (wobei Zonen mit einer höheren Frequenz von x, y "wärmer" wären)?

Wenn Sie keine Sechsecke möchten, können Sie die histogram2dFunktion von numpy verwenden:

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

Dies ergibt eine 50x50 Heatmap. Wenn Sie beispielsweise 512 x 384 möchten, können Sie bins=(512, 384)den Anruf bei tätigen histogram2d.

Beispiel:

Im Matplotlib- Lexikon möchten Sie einen Hexbin- Plot.

Wenn Sie mit dieser Art von Plot nicht vertraut sind, handelt es sich nur um ein bivariates Histogramm, bei dem die xy-Ebene durch ein regelmäßiges Sechseckgitter tesselliert wird.

Aus einem Histogramm können Sie also einfach die Anzahl der Punkte zählen, die in jedes Sechseck fallen, den Plotbereich als eine Reihe von Fenstern diskretisieren und jeden Punkt einem dieser Fenster zuweisen. Schließlich ordnen Sie die Fenster auf ein Farbarray , und Sie haben ein hexbin Diagramm bekam.

Obwohl Sechsecke weniger häufig verwendet werden als z. B. Kreise oder Quadrate, ist es eine intuitive Wahl für die Geometrie des Binning-Containers:

• Sechsecke haben Symmetrie zum nächsten Nachbarn (z. B. sind quadratische Bins nicht, z. B. ist der Abstand von einem Punkt an der Grenze eines Quadrats zu einem Punkt innerhalb dieses Quadrats nicht überall gleich) und

• Sechseck ist das höchste n-Polygon, das eine regelmäßige ebene Tessellation ergibt (dh Sie können Ihren Küchenboden mit sechseckigen Fliesen sicher neu modellieren, da Sie nach Fertigstellung keinen Hohlraum zwischen den Fliesen haben - nicht wahr alle anderen höheren n, n> = 7, Polygone).

( Matplotlib verwendet den Begriff hexbin Grundstück; so tun (AFAIK) alle der Plotten Bibliotheken für R , noch ich weiß nicht , ob dies die allgemein akzeptierte Bezeichnung für Grundstücke dieser Art ist, obwohl ich es wahrscheinlich gegeben vermuten , dass hexbin kurz ist für hexagonales Binning , das den wesentlichen Schritt bei der Vorbereitung der Daten für die Anzeige beschreibt.)

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

Bearbeiten: Für eine bessere Annäherung an Alejandros Antwort siehe unten.

Ich weiß, dass dies eine alte Frage ist, wollte aber Alejandros Antwort etwas hinzufügen: Wenn Sie ein schönes geglättetes Bild ohne Verwendung von py-sphviewer wünschen, können Sie stattdessen np.histogram2deinen Gauß-Filter (von scipy.ndimage.filters) verwenden und auf die Heatmap anwenden :

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

Produziert:

Das Streudiagramm und s = 16 sind für Agape Gal'lo übereinander aufgetragen (zur besseren Ansicht klicken):

Ein Unterschied, den ich bei meinem Gaußschen Filteransatz und Alejandros Ansatz bemerkte, war, dass seine Methode lokale Strukturen viel besser zeigt als meine. Daher habe ich eine einfache Methode für den nächsten Nachbarn auf Pixelebene implementiert. Diese Methode berechnet für jedes Pixel die inverse Summe der Abstände der nnächstgelegenen Punkte in den Daten. Diese Methode ist mit einer hohen Auflösung ziemlich rechenintensiv und ich denke, es gibt einen schnelleren Weg. Lassen Sie mich wissen, wenn Sie Verbesserungen haben.

Update: Wie ich vermutet habe, gibt es eine viel schnellere Methode mit Scipy's scipy.cKDTree. Siehe Gabriels Antwort für die Implementierung.

Wie auch immer, hier ist mein Code:

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

Ergebnis:

Anstatt np.hist2d zu verwenden, das im Allgemeinen ziemlich hässliche Histogramme erzeugt, möchte ich py-sphviewer recyceln , ein Python-Paket zum Rendern von Partikelsimulationen mit einem adaptiven Glättungskern, das einfach über pip installiert werden kann (siehe Webseiten-Dokumentation). Betrachten Sie den folgenden Code, der auf dem Beispiel basiert:

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

welches das folgende Bild erzeugt:

Wie Sie sehen, sehen die Bilder ziemlich gut aus, und wir können verschiedene Unterstrukturen darauf identifizieren. Diese Bilder werden so konstruiert, dass sie für jeden Punkt innerhalb eines bestimmten Bereichs ein bestimmtes Gewicht verteilen, das durch die Glättungslänge definiert wird, die wiederum durch den Abstand zum näheren nb- Nachbarn gegeben ist (ich habe für die Beispiele 16, 32 und 64 gewählt). Daher sind Regionen mit höherer Dichte im Vergleich zu Regionen mit niedrigerer Dichte typischerweise über kleinere Regionen verteilt.

Die Funktion myplot ist nur eine sehr einfache Funktion, die ich geschrieben habe, um py-sphviewer die x, y-Daten zu geben, um die Magie auszuführen.

Wenn Sie 1.2.x verwenden

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

Seaborn hat jetzt die Joint-Plot-Funktion, die hier gut funktionieren sollte:

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

und die erste Frage war ... wie man Streuwerte in Gitterwerte umwandelt, richtig? histogram2dzählt jedoch die Häufigkeit pro Zelle. Wenn Sie jedoch andere Daten pro Zelle als nur die Häufigkeit haben, müssen Sie einige zusätzliche Arbeiten ausführen.

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

Ich habe also einen Datensatz mit Z-Ergebnissen für X- und Y-Koordinaten. Ich berechnete jedoch nur wenige Punkte außerhalb des interessierenden Bereichs (große Lücken) und jede Menge Punkte in einem kleinen interessierenden Bereich.

Ja hier wird es schwieriger, aber auch lustiger. Einige Bibliotheken (sorry):

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

Pyplot ist heute meine Grafik-Engine. cm ist eine Reihe von Farbkarten mit einer interessanten Auswahl. numpy für die Berechnungen und griddata zum Anhängen von Werten an ein festes Gitter.

Letzteres ist besonders wichtig, weil die Häufigkeit von xy-Punkten in meinen Daten nicht gleichmäßig verteilt ist. Beginnen wir zunächst mit einigen Grenzen, die zu meinen Daten passen, und einer beliebigen Rastergröße. Die Originaldaten haben Datenpunkte auch außerhalb dieser x- und y-Grenzen.

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

Wir haben also ein Raster mit 500 Pixeln zwischen den Min- und Max-Werten von x und y definiert.

In meinen Daten sind viel mehr als die 500 Werte im Bereich von hohem Interesse verfügbar; in der Erwägung, dass es im Niedrigzinsbereich nicht einmal 200 Werte im Gesamtnetz gibt; zwischen den grafischen Grenzen von x_minund x_maxgibt es noch weniger.

Um ein schönes Bild zu erhalten, besteht die Aufgabe darin, einen Durchschnitt für die hohen Zinswerte zu erhalten und die Lücken an anderer Stelle zu schließen.

Ich definiere jetzt mein Raster. Für jedes xx-yy Paar möchte ich eine Farbe haben.

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

Warum die seltsame Form? scipy.griddata möchte eine Form von (n, D).

Griddata berechnet einen Wert pro Punkt im Raster nach einer vordefinierten Methode. Ich wähle "am nächsten" - leere Gitterpunkte werden mit Werten des nächsten Nachbarn gefüllt. Dies sieht so aus, als hätten die Bereiche mit weniger Informationen größere Zellen (auch wenn dies nicht der Fall ist). Man könnte wählen, "linear" zu interpolieren, dann sehen Bereiche mit weniger Informationen weniger scharf aus. Geschmackssache, wirklich.

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

Und hüpfen, wir übergeben an matplotlib, um die Handlung anzuzeigen

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

Um den spitzen Teil der V-Form herum haben Sie bei meiner Suche nach dem Sweet Spot viele Berechnungen durchgeführt, während die weniger interessanten Teile fast überall eine niedrigere Auflösung haben.

Hier ist Jurgys Ansatz für den nächsten Nachbarn, der jedoch mit scipy.cKDTree implementiert wurde . In meinen Tests ist es ungefähr 100x schneller.

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

Erstellen Sie ein zweidimensionales Array, das den Zellen in Ihrem endgültigen Bild entspricht heatmap_cells und instanziieren Sie es als alle Nullen.

Wählen Sie zwei Skalierungsfaktoren aus, die den Unterschied zwischen jedem Array-Element in realen Einheiten für jede Dimension definieren, z. B. x_scaleundy_scale . . Wählen Sie diese so aus, dass alle Ihre Datenpunkte innerhalb der Grenzen des Heatmap-Arrays liegen.

Für jeden Rohdatenpunkt mit x_valueund y_value:

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

Hier ist eine, die ich mit einem 1-Millionen-Punkte-Set mit 3 Kategorien (rot, grün und blau) erstellt habe. Hier ist ein Link zum Repository, wenn Sie die Funktion ausprobieren möchten. Github Repo

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)

Sehr ähnlich der Antwort von @ Piti , aber mit 1 Anruf anstelle von 2, um die Punkte zu generieren:

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

Ausgabe:

Ich fürchte, ich bin etwas spät zur Party, aber ich hatte vor einiger Zeit eine ähnliche Frage. Die akzeptierte Antwort (von @ptomato) hat mir geholfen, aber ich möchte sie auch posten, falls sie für jemanden von Nutzen ist.

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

Hier ist das Ergebnis