Wygeneruj mapę cieplną w MatPlotLib przy użyciu zestawu danych rozproszonych

Mam zestaw punktów danych X, Y (około 10k), które są łatwe do wykreślenia jako wykres rozproszenia, ale które chciałbym przedstawić jako mapę termiczną.

Przejrzałem przykłady w MatPlotLib i wszystkie wydają się już zaczynać od wartości komórek heatmap w celu wygenerowania obrazu.

Czy istnieje metoda, która przekształca wiązkę x, y, wszystko inne, w mapę cieplną (gdzie strefy o wyższej częstotliwości x, y byłyby „cieplejsze”)?

Jeśli nie chcesz sześciokątów, możesz użyć histogram2dfunkcji 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()

To tworzy mapę cieplną 50x50. Jeśli chcesz, powiedzmy, 512x384, możesz bins=(512, 384)zadzwonić do histogram2d.

Przykład:

W leksykonie Matplotlib myślę, że chcesz wykresu szesnastkowego .

Jeśli nie znasz tego typu wykresu, jest to tylko dwuwymiarowy histogram, na którym płaszczyzna xy jest mozaikowana regularną siatką sześciokątów.

Tak więc z histogramu można po prostu policzyć liczbę punktów przypadających na każdy sześciokąt, dyskretyzować obszar kreślenia jako zestaw okien , przypisać każdy punkt do jednego z tych okien; na koniec zamapuj okna na tablicę kolorów , a otrzymasz schemat szesnastkowy.

Chociaż rzadziej używane niż np. Koła lub kwadraty, sześciokąty są lepszym wyborem dla geometrii pojemnika binningowego, jest intuicyjne:

• sześciokąty mają symetrię najbliższego sąsiada (np. kwadratowe kosze nie, np. odległość od punktu na granicy kwadratu do punktu wewnątrz tego kwadratu nie jest wszędzie równa) i

• sześciokąt jest najwyższym n-wielokątem, który zapewnia regularną teselację w płaszczyźnie (tzn. możesz bezpiecznie ponownie modelować podłogę w kuchni za pomocą płytek w kształcie sześciokąta, ponieważ po zakończeniu nie będziesz mieć pustej przestrzeni między płytkami - nie dotyczy to wszystkie inne wyższe-n, n> = 7, wielokąty).

( Matplotlib używa terminu wykres heksbinowy ; podobnie jak (AFAIK) wszystkie biblioteki kreślące dla R ; nadal nie wiem, czy jest to ogólnie przyjęty termin dla wykresów tego typu, chociaż podejrzewam, że jest prawdopodobne, biorąc pod uwagę, że hexbin jest krótki dla binoksu heksagonalnego , który opisuje niezbędny krok w przygotowaniu danych do wyświetlenia.)

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

Edycja: Dla lepszego przybliżenia odpowiedzi Alejandro, patrz poniżej.

Wiem, że to stare pytanie, ale chciałem dodać coś do odpowiedzi Alejandro: Jeśli chcesz uzyskać ładny wygładzony obraz bez użycia py-sphviewer, możesz zamiast tego użyć np.histogram2di zastosować filtr gaussowski (od scipy.ndimage.filters) do mapy cieplnej :

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

Produkuje:

Wykres rozproszenia is = 16 wykreślony jeden na drugim dla Agape Gal'lo (kliknij, aby uzyskać lepszy widok):

Jedną różnicę, którą zauważyłem przy moim podejściu z filtrem gaussowskim i podejściu Alejandro, było to, że jego metoda pokazuje struktury lokalne znacznie lepiej niż moja. Dlatego zaimplementowałem prostą metodę najbliższego sąsiada na poziomie pikseli. Ta metoda oblicza dla każdego piksela odwrotną sumę odległościn najbliższych punktów w danych. Ta metoda jest dość kosztowna obliczeniowo w wysokiej rozdzielczości i myślę, że jest szybszy sposób, więc daj mi znać, jeśli masz jakieś ulepszenia.

Aktualizacja: jak podejrzewałam, metoda Scipy jest znacznie szybsza scipy.cKDTree. Zobacz odpowiedź Gabriela na wdrożenie.

Tak czy inaczej, oto mój kod:

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

Wynik:

Zamiast używać np.hist2d, który generalnie generuje dość brzydkie histogramy, chciałbym poddać recyklingowi py-sphviewer , pakiet Pythona do renderowania symulacji cząstek przy użyciu adaptacyjnego wygładzania jądra, który można łatwo zainstalować z pip (patrz dokumentacja strony). Rozważ następujący kod oparty na przykładzie:

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

co daje następujący obraz:

Jak widać, obrazy wyglądają całkiem ładnie, a my jesteśmy w stanie zidentyfikować na nim różne podstruktury. Te obrazy są konstruowane rozkładając podaną wagę dla każdego punktu w określonej dziedzinie, zdefiniowanej przez długość wygładzania, która z kolei wynika z odległości do bliższego sąsiada nb (na przykład wybrałem 16, 32 i 64). Zatem regiony o większej gęstości są zazwyczaj rozłożone na mniejsze regiony w porównaniu do regionów o mniejszej gęstości.

Funkcja myplot jest po prostu bardzo prostą funkcją, którą napisałem, aby przekazać dane x, y py-sphviewerowi, aby wykonał magię.

Jeśli korzystasz z 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()

Seaborn ma teraz funkcję jointplot, która powinna tu dobrze działać:

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

a początkowe pytanie brzmiało ... jak przekonwertować wartości rozproszenia na wartości siatki, prawda? histogram2dzlicza częstotliwość na komórkę, jednak jeśli masz inne dane na komórkę niż tylko częstotliwość, potrzebujesz dodatkowej pracy.

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

Mam więc zestaw danych z wynikami Z dla współrzędnych X i Y. Jednak obliczałem kilka punktów poza obszarem zainteresowania (duże luki) i stosy punktów w małym obszarze zainteresowania.

Tak, tutaj staje się trudniejsze, ale także przyjemniejsze. Niektóre biblioteki (przepraszam):

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

pyplot jest dziś moim silnikiem graficznym, cm to szereg map kolorów z pewnym nieciekawym wyborem. numpy do obliczeń i griddata do dołączania wartości do stałej siatki.

Ten ostatni jest ważny, szczególnie, że częstotliwość punktów xy nie jest równomiernie rozłożona w moich danych. Najpierw zacznijmy od pewnych granic pasujących do moich danych i dowolnego rozmiaru siatki. Oryginalne dane mają punkty danych również poza tymi granicami xiy.

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

Zdefiniowaliśmy więc siatkę z 500 pikselami między wartościami min i max x i y.

W moich danych jest o wiele więcej niż 500 wartości dostępnych w obszarze o dużym zainteresowaniu; mając na uwadze, że w obszarze niskiego zainteresowania nie ma nawet 200 wartości w całej sieci; między granicami graficznymi x_minix_max jest jeszcze mniej.

Aby uzyskać ładny obraz, zadaniem jest uzyskanie średniej dla wysokich wartości zainteresowania i wypełnienie luk gdzie indziej.

Teraz definiuję swoją siatkę. Dla każdej pary xx-yy chcę mieć kolor.

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

Skąd ten dziwny kształt? scipy.griddata chce mieć kształt (n, D).

Griddata oblicza jedną wartość na punkt w siatce za pomocą predefiniowanej metody. Wybieram „najbliższy” - puste punkty siatki zostaną wypełnione wartościami od najbliższego sąsiada. Wygląda na to, że obszary z mniejszą ilością informacji mają większe komórki (nawet jeśli tak nie jest). Można wybrać interpolację „liniową”, a obszary z mniejszą ilością informacji wyglądają mniej ostro. Naprawdę sprawa gustu.

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

I hop, przekazujemy matplotlib, aby wyświetlić fabułę

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

Wokół spiczastej części Kształtu V widać, że wykonałem wiele obliczeń podczas poszukiwania słodkiego miejsca, podczas gdy mniej interesujące części prawie wszędzie indziej mają niższą rozdzielczość.

Oto podejście najbliższego sąsiada Jurgy'ego, ale zaimplementowane przy użyciu scipy.cKDTree . W moich testach jest około 100 razy szybszy.

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

Utwórz dwuwymiarową tablicę, która odpowiada komórkom na ostatecznym obrazie, zwaną say heatmap_cellsi utwórz ją jako wszystkie zera.

Wybierz dwa współczynniki skalowania, które określają różnicę między każdym elementem tablicy w jednostkach rzeczywistych, dla każdego wymiaru, powiedz x_scalei y_scale. Wybierz je w taki sposób, aby wszystkie twoje punkty danych mieściły się w granicach tablicy Heatmap.

Dla każdego nieprzetworzonego punktu danych za pomocą x_valuei y_value:

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

Oto jeden, który wykonałem na zestawie 1 milion punktów z 3 kategoriami (w kolorze czerwonym, zielonym i niebieskim). Oto link do repozytorium, jeśli chcesz wypróbować tę funkcję. 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)

Bardzo podobny do odpowiedzi @ Piti , ale do wygenerowania punktów używa 1 połączenia zamiast 2:

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

Wynik:

Obawiam się, że jestem trochę spóźniony na imprezę, ale już kiedyś miałem podobne pytanie. Zaakceptowana odpowiedź (autor @ptomato) pomogła mi, ale chciałbym również zamieścić to na wypadek, gdyby było dla kogoś przydatne.

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

Oto wynik