Creating Scatter Plots Colored by Density: A Comprehensive Guide with Python and Matplotlib

Introduction and Problem Context

In data visualization, scatter plots are a common tool for displaying relationships between two variables. However, when the number of data points is large or densely distributed, traditional scatter plots may fail to effectively convey spatial density information, obscuring important patterns. Therefore, scatter plots colored by density have emerged, using color mapping to reflect the spatial density around each point, thereby enhancing the informational content of visualizations. This article aims to explore various methods for implementing such plots in the Python environment using the Matplotlib library, analyzing their advantages and disadvantages.

Basic Method: Using Gaussian Kernel Density Estimation

A straightforward approach involves using scipy.stats.gaussian_kde to compute density values for each point and map them to colors. The following complete example code demonstrates how to generate and visualize a density-colored scatter plot.

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

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

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

# Create figure
fig, ax = plt.subplots()
ax.scatter(x, y, c=z, s=100, edgecolor='none')
plt.colorbar(ax.collections[0], label='Density Value')
plt.show()

In this example, the gaussian_kde function estimates the probability density of the data points, with the result z being an array representing the density estimate for each point. By passing z to the c parameter of the scatter function, Matplotlib automatically colors the points based on density. To optimize visual effects, data points can be sorted by density to ensure that the densest points are displayed on top, avoiding occlusion.

# Sort points by density
idx = z.argsort()
x, y, z = x[idx], y[idx], z[idx]

fig, ax = plt.subplots()
ax.scatter(x, y, c=z, s=50, cmap='viridis', edgecolor='none')
plt.colorbar(ax.collections[0], label='Density Value')
plt.show()

This method is simple and easy to use, but for large-scale datasets (e.g., over 100,000 points), computing kernel density estimates can be very time-consuming, taking approximately 11 minutes for 100,000 points on standard hardware. Therefore, more efficient alternatives need to be considered.

Efficient Alternatives: For Large-Scale Datasets

When dealing with large-scale data, the following methods offer better performance and visual quality.

Using the mpl-scatter-density Library

mpl-scatter-density is a library specifically designed for high-density scatter plots, achieving fast rendering through pixel-level density calculations. After installation, its custom projection can be used to create density plots.

import mpl_scatter_density
from matplotlib.colors import LinearSegmentedColormap

# Define custom colormap
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)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection='scatter_density')
density = ax.scatter_density(x, y, cmap=white_viridis)
plt.colorbar(density, label='Points per Pixel')
plt.show()

This method requires only about 0.05 seconds to render 100,000 points in tests and maintains high quality during zooming, making it suitable for interactive visualizations.

Using the Datashader Library

Datashader is another powerful tool designed for big data visualization, supporting integration with Matplotlib. It efficiently computes density through rasterization techniques.

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

# Convert data to DataFrame
df = pd.DataFrame(dict(x=x, y=y))
fig, ax = plt.subplots()
dsartist = dsshow(
    df,
    ds.Point("x", "y"),
    ds.count(),
    vmin=0,
    vmax=35,
    norm="linear",
    aspect="auto",
    ax=ax,
)
plt.colorbar(dsartist, label='Count')
plt.show()

Datashader processes 100,000 points in about 0.83 seconds and offers flexible coloring options, such as coloring based on a third variable, making it suitable for complex data analysis.

Using the hist2d Method

Matplotlib's built-in hist2d function computes density through two-dimensional histograms, providing a fast and simple approach.

fig, ax = plt.subplots()
ax.hist2d(x, y, bins=(50, 50), cmap=plt.cm.jet)
plt.colorbar(ax.collections[0], label='Density')
plt.show()

With 50x50 bins, rendering takes only about 0.021 seconds, but zoomed-in effects may not be as detailed as the previous methods. Increasing the number of bins (e.g., 1000x1000) improves resolution but increases time to about 0.173 seconds and requires manual parameter adjustment.

Density Interpolation Based on np.histogram2d

Another method combines np.histogram2d with interpolation techniques to compute density values for each point. The following function implements this process.

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):
    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)
    z[np.isnan(z)] = 0.0
    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 = plt.colorbar(cm.ScalarMappable(norm=norm), ax=ax)
    cbar.ax.set_ylabel('Density')
    return ax

# Usage example
density_scatter(x, y, bins=[30, 30], cmap='plasma')
plt.show()

This method processes 100,000 points in about 0.073 seconds (50x50 bins) or 0.368 seconds (1000x1000 bins) in tests, balancing speed and quality but relying on interpolation accuracy.

Performance and Quality Comparison

Comparing the various methods comprehensively, gaussian_kde is simple and effective for small datasets but computationally expensive; mpl-scatter-density and Datashader excel in large-scale scenarios, offering fast rendering and high zoom quality; hist2d and np.histogram2d-based methods provide a good balance for medium-sized data. Selection should consider data scale, computational resources, and visualization needs. For example, mpl-scatter-density may be the best choice for real-time analysis, while gaussian_kde or interpolation methods are more suitable for research requiring precise density estimates.

Conclusion and Best Practices

Creating scatter plots colored by density is a crucial means of enhancing data visualization effectiveness. In Python, Matplotlib combined with extension libraries offers multiple implementation pathways. For beginners or small datasets, starting with the gaussian_kde method is recommended; for large-scale data, consider mpl-scatter-density or Datashader to improve performance. In practical applications, it is advisable to first assess data characteristics, then select an appropriate method, and optimize output by adjusting parameters (e.g., colormap, bin size). In the future, with advances in computing technology, these methods are expected to further integrate, providing more efficient and flexible visualization solutions.