Jak mogę wykonać wykres punktowy pokolorowany według gęstości w matplotlib?

83

Chciałbym sporządzić wykres punktowy, na którym każdy punkt jest pokolorowany według gęstości przestrzennej pobliskich punktów.

Natknąłem się na bardzo podobne pytanie, które pokazuje przykład tego przy użyciu R:

Wykres punktowy R: kolor symbolu przedstawia liczbę nakładających się punktów

Jaki jest najlepszy sposób osiągnięcia czegoś podobnego w Pythonie przy użyciu matplotlib?

2964502
źródło
4
Cześć! Ludzie przegłosowali Cię prawdopodobnie dlatego, że nie przepisałeś pytania ani nie podałeś żadnego kontekstu, ani nie pokazałeś żadnej próby zrobienia tego samemu. Rozważ edytowanie pytania, aby było samowystarczalne (a nie tylko link), a w przypadku przyszłych pytań prosimy o podjęcie próby przed wysłaniem.
askewchan

Odpowiedzi:

157

Oprócz hist2dlub hexbinzgodnie z sugestią @askewchan możesz użyć tej samej metody, której używa zaakceptowana odpowiedź w pytaniu, do którego prowadzi link.

Jeśli chcesz to zrobić:

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde

# Generate fake data
x = np.random.normal(size=1000)
y = x * 3 + np.random.normal(size=1000)

# Calculate the point density
xy = np.vstack([x,y])
z = gaussian_kde(xy)(xy)

fig, ax = plt.subplots()
ax.scatter(x, y, c=z, s=100, edgecolor='')
plt.show()

wprowadź opis obrazu tutaj

Jeśli chcesz, aby punkty były wykreślane w kolejności gęstości, tak aby najgęstsze punkty zawsze znajdowały się na górze (podobnie jak w połączonym przykładzie), po prostu posortuj je według wartości z. Zamierzam również użyć mniejszego rozmiaru markera, ponieważ wygląda trochę lepiej:

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde

# Generate fake data
x = np.random.normal(size=1000)
y = x * 3 + np.random.normal(size=1000)

# Calculate the point density
xy = np.vstack([x,y])
z = gaussian_kde(xy)(xy)

# Sort the points by density, so that the densest points are plotted last
idx = z.argsort()
x, y, z = x[idx], y[idx], z[idx]

fig, ax = plt.subplots()
ax.scatter(x, y, c=z, s=50, edgecolor='')
plt.show()

wprowadź opis obrazu tutaj

Joe Kington
źródło
4
Sprytne, zwłaszcza te „najgęstsze” na górze :)
askewchan
5
@Leszek - Ether call plt.colorbar(), a jeśli wolisz być bardziej dosadny, zrób cax = ax.scatter(...)i wtedy fig.colorbar(cax). Należy pamiętać, że jednostki są różne. Ta metoda szacuje rozkład prawdopodobieństwa dla punktów, więc wartości będą wynosić od 0 do 1 (i zazwyczaj nie będą bardzo zbliżone do 1). Możesz przekonwertować z powrotem na coś bliższego zliczeniom histogramu, ale zajmuje to trochę pracy (musisz znać parametry, które gaussian_kdeoszacowano na podstawie danych).
Joe Kington
1
Bardzo dobrze! Sprawdzanie innych KDE w Pythonie również może być przydatne: jakevdp.github.io/blog/2013/12/01/kernel-density-estimation i scikit-learn.org/stable/modules/density.html W moim przypadku scipy.stats „KDE trwało zbyt długo
Rems
1
Dlaczego jądro Gaussa jest wywoływane dwukrotnie z (xy)?
Arjan Groen
@ArjanGroen Pierwsze wywołanie tworzy nowy obiekt gaussian_kde, a drugie wywołanie ocenia szacunkowy plik PDF na zbiorze punktów (skrót do wywołania metody oceny).
qRTPCR
35

Możesz zrobić histogram:

import numpy as np
import matplotlib.pyplot as plt

# fake data:
a = np.random.normal(size=1000)
b = a*3 + np.random.normal(size=1000)

plt.hist2d(a, b, (50, 50), cmap=plt.cm.jet)
plt.colorbar()

2dhist

askewchan
źródło
Aby lepiej dopasować skalę rozwiązania Joe Kingtona, możesz wykreślić w skali logu  : plt.hist2d(…, norm = LogNorm())(z from matplotlib.colors import LogNorm).
Skippy le Grand Gourou
29

Ponadto, jeśli liczba punktów powoduje, że obliczenia KDE są zbyt wolne, kolor można interpolować w np.histogram2d [Aktualizacja w odpowiedzi na komentarze: Jeśli chcesz pokazać pasek kolorów, użyj plt.scatter () zamiast ax.scatter () by plt.colorbar ()]:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.colors import Normalize 
from scipy.interpolate import interpn

def density_scatter( x , y, ax = None, sort = True, bins = 20, **kwargs )   :
    """
    Scatter plot colored by 2d histogram
    """
    if ax is None :
        fig , ax = plt.subplots()
    data , x_e, y_e = np.histogram2d( x, y, bins = bins, density = True )
    z = interpn( ( 0.5*(x_e[1:] + x_e[:-1]) , 0.5*(y_e[1:]+y_e[:-1]) ) , data , np.vstack([x,y]).T , method = "splinef2d", bounds_error = False)

    #To be sure to plot all data
    z[np.where(np.isnan(z))] = 0.0

    # Sort the points by density, so that the densest points are plotted last
    if sort :
        idx = z.argsort()
        x, y, z = x[idx], y[idx], z[idx]

    ax.scatter( x, y, c=z, **kwargs )

    norm = Normalize(vmin = np.min(z), vmax = np.max(z))
    cbar = fig.colorbar(cm.ScalarMappable(norm = norm), ax=ax)
    cbar.ax.set_ylabel('Density')

    return ax


if "__main__" == __name__ :

    x = np.random.normal(size=100000)
    y = x * 3 + np.random.normal(size=100000)
    density_scatter( x, y, bins = [30,30] )

Guillaume
źródło
1
To świetna wskazówka, dziękuję. Rysowałem 100k punktów, a gaussian_kde był zbyt wolny.
Emanuel
2
Ostrzeżenie, zauważyłem, że w niektórych przypadkach generuje to NaN, a ponieważ "bounds_error = False" jest ciche. Punkty z c ustawionym na NaNs nie są drukowane. To nie jest problem z gaussian_kde.
Emanuel
Wielkie dzięki za tę odpowiedź. Zwykle potrzebujemy takiej mapy cieplnej, gdy mamy dużą liczbę punktów danych, a KDE działa w tym przypadku bardzo wolno. Jednak nadal istnieje otwarty problem. Chcę dołączyć kolorowy pasek wskazujący częstotliwość! Powoduje to zgłoszenie błędu: obiekt „AxesSubplot” nie ma atrybutu „autoscale_None”. Zrobiłem „plt.colorbar (scat, ax = ax)”
Vinod Kumar
@VinodKumar, czy dowiedziałeś się, jak wykreślić pasek kolorów?
Daniel
1
@Daniel tak, jest to możliwe, zobacz zredagowaną odpowiedź. Następnie podczas budowania histogramu musisz ustawić gęstość = True, w przeciwnym razie pasek kolorów zależy od rozmiaru bin. @ Emanuel, rzeczywiście! Zamieniłem wartości NaN na zero, aby mieć pewność, że wykreślę wszystkie punkty (NaN powinno się zdarzyć, gdy nie ma zbyt wielu danych, więc 0,0 powinno wystarczyć)
Guillaume
6

Wykreślasz> 100 tys. Punktów danych?

Odpowiedź akceptowana , używając gaussian_kde () zajmie dużo czasu. Na mojej maszynie 100 tys. Rzędów zajęło około 11 minut . Tutaj dodam dwie alternatywne metody ( mpl-scatter-density i datashader ) i porównam podane odpowiedzi z tym samym zestawem danych.

Poniżej użyłem testowego zestawu danych 100 tys. Wierszy:

import matplotlib.pyplot as plt
import numpy as np

# Fake data for testing
x = np.random.normal(size=100000)
y = x * 3 + np.random.normal(size=100000)

Porównanie czasu wyjściowego i obliczeniowego

Poniżej znajduje się porównanie różnych metod.

1: mpl-scatter-density

Instalacja

pip install mpl-scatter-density

Przykładowy kod

import mpl_scatter_density # adds projection='scatter_density'
from matplotlib.colors import LinearSegmentedColormap

# "Viridis-like" colormap with white background
white_viridis = LinearSegmentedColormap.from_list('white_viridis', [
    (0, '#ffffff'),
    (1e-20, '#440053'),
    (0.2, '#404388'),
    (0.4, '#2a788e'),
    (0.6, '#21a784'),
    (0.8, '#78d151'),
    (1, '#fde624'),
], N=256)

def using_mpl_scatter_density(fig, x, y):
    ax = fig.add_subplot(1, 1, 1, projection='scatter_density')
    density = ax.scatter_density(x, y, cmap=white_viridis)
    fig.colorbar(density, label='Number of points per pixel')

fig = plt.figure()
using_mpl_scatter_density(fig, x, y)
plt.show()

Rysowanie trwało 0,05 sekundy: przy użyciu gęstości rozproszenia mpl

A powiększenie wygląda całkiem nieźle: powiększ gęstość rozproszenia mpl

2: datashader

pip install "git+https://github.com/nvictus/datashader.git@mpl"

Kod (źródło dsshow tutaj ):

from functools import partial

import datashader as ds
from datashader.mpl_ext import dsshow
import pandas as pd

dyn = partial(ds.tf.dynspread, max_px=40, threshold=0.5)

def using_datashader(ax, x, y):

    df = pd.DataFrame(dict(x=x, y=y))
    da1 = dsshow(df, ds.Point('x', 'y'), spread_fn=dyn, aspect='auto', ax=ax)
    plt.colorbar(da1)

fig, ax = plt.subplots()
using_datashader(ax, x, y)
plt.show()
  • Narysowanie tego zajęło 0,83 s:

wprowadź opis obrazu tutaj

a powiększony obraz wygląda świetnie!

wprowadź opis obrazu tutaj

3: scatter_with_gaussian_kde

def scatter_with_gaussian_kde(ax, x, y):
    # https://stackoverflow.com/a/20107592/3015186
    # Answer by Joel Kington

    xy = np.vstack([x, y])
    z = gaussian_kde(xy)(xy)

    ax.scatter(x, y, c=z, s=100, edgecolor='')
  • Narysowanie tego zajęło 11 minut: scatter_with_gaussian_kde

4: using_hist2d

import matplotlib.pyplot as plt
def using_hist2d(ax, x, y, bins=(50, 50)):
    # https://stackoverflow.com/a/20105673/3015186
    # Answer by askewchan
    ax.hist2d(x, y, bins, cmap=plt.cm.jet)

  • Narysowanie tego kosza zajęło 0,021 s = (50,50): using_hist2d_50
  • Narysowanie tego kosza zajęło 0,173 s = (1000,1000): using_hist2d_1000
  • Wady: Powiększone dane nie wyglądają tak dobrze, jak w przypadku gęstości rozproszonej mpl lub modułu danych. Musisz także samodzielnie określić liczbę pojemników.

powiększony hist2d 1000bins

5: density_scatter

  • Kod jest jak w odpowiedzi przez Guillaume .
  • Narysowanie tego z koszami = (50,50) zajęło 0,073 s: density_scatter_50bins
  • Narysowanie tego z koszami = (1000,1000) zajęło 0,368 s: density_scatter_1000bins
np8
źródło