如何实现单元格大小与密度成反比的热力图（每格至少10个数据源）

核心思路 #

要复刻Donley等人2012年论文的图8效果，不能用固定分箱的np.histogram2d，得采用自适应空间划分的方式：先把数据按密度拆分成每个区域至少包含10个样本的单元格，再基于每个单元格的样本密度调整大小，用z值的平均值映射颜色。

具体实现步骤 #

1. 自适应划分单元格 #

推荐用四分树（Quadtree）递归拆分，它能根据数据密度自动调整区域大小，直到每个子区域的样本数≥10：

• 从整个数据的x、y范围作为初始区域开始，不断将区域拆分为四个象限；
• 当某个子区域的样本数不足10时停止拆分，反之继续递归拆分。
• 也可以用sklearn.cluster.MeanShift聚类，调整带宽参数控制聚类大小，确保每个聚类样本数达标。

2. 计算单元格参数 #

对每个划分好的单元格：

• 统计样本数count，计算密度=count/单元格面积，设置单元格大小与密度成反比（比如用1/sqrt(density)作为边长，保证面积与密度成反比）；
• 计算该区域内z值的平均值z_mean，作为颜色映射的依据。

3. 绘制可变大小热力图 #

用matplotlib的Rectangle逐个绘制单元格，将z_mean映射到指定色卡，完成可视化。

代码示例 #

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

# 替换为你的真实数据
np.random.seed(42)
x = np.random.normal(0, 1, 1000)
y = np.random.normal(0, 1, 1000)
z = np.sin(x) + np.cos(y) + np.random.normal(0, 0.1, 1000)

# 自定义四分树拆分函数，确保每个区域样本数≥10
def quadtree_split(x, y, z, min_samples=10):
    cells = []
    x_min, x_max = x.min(), x.max()
    y_min, y_max = y.min(), y.max()
    
    def split_region(x_range, y_range, x_data, y_data, z_data):
        if len(x_data) <= min_samples:
            z_mean = np.mean(z_data)
            count = len(x_data)
            area = (x_range[1]-x_range[0])*(y_range[1]-y_range[0])
            density = count / area
            # 边长与密度平方根成反比，0.1为缩放系数，可按需调整
            side_length = np.sqrt(1/density) * 0.1
            cells.append({
                'x_min': x_range[0], 'y_min': y_range[0],
                'side': side_length, 'z_mean': z_mean
            })
            return
        # 拆分区域为四个象限
        x_mid = (x_range[0]+x_range[1])/2
        y_mid = (y_range[0]+y_range[1])/2
        idx1 = (x_data <= x_mid) & (y_data <= y_mid)
        idx2 = (x_data > x_mid) & (y_data <= y_mid)
        idx3 = (x_data <= x_mid) & (y_data > y_mid)
        idx4 = (x_data > x_mid) & (y_data > y_mid)
        # 递归拆分各子区域
        split_region([x_range[0],x_mid], [y_range[0],y_mid], x_data[idx1], y_data[idx1], z_data[idx1])
        split_region([x_mid,x_range[1]], [y_range[0],y_mid], x_data[idx2], y_data[idx2], z_data[idx2])
        split_region([x_range[0],x_mid], [y_mid,y_range[1]], x_data[idx3], y_data[idx3], z_data[idx3])
        split_region([x_mid,x_range[1]], [y_mid,y_range[1]], x_data[idx4], y_data[idx4], z_data[idx4])
    
    split_region([x_min,x_max], [y_min,y_max], x, y, z)
    return cells

# 获取划分后的单元格
cells = quadtree_split(x, y, z, min_samples=10)

# 绘制热力图
fig, ax = plt.subplots(figsize=(8,8))
norm = plt.Normalize(min(c['z_mean'] for c in cells), max(c['z_mean'] for c in cells))
cmap = plt.get_cmap('viridis')

for cell in cells:
    rect = Rectangle((cell['x_min'], cell['y_min']), cell['side'], cell['side'],
                     color=cmap(norm(cell['z_mean'])), alpha=0.8)
    ax.add_patch(rect)

# 配置坐标轴与颜色条
ax.set_xlim(x.min(), x.max())
ax.set_ylim(y.min(), y.max())
sm = plt.cm.ScalarMappable(norm=norm, cmap=cmap)
sm.set_array([])
fig.colorbar(sm, ax=ax, label='Average z value')

plt.xlabel('X')
plt.ylabel('Y')
plt.title('Adaptive Density Heatmap (≥10 samples per cell)')
plt.show()

注意事项 #

• 四分树的缩放系数需要根据你的数据范围手动调整，保证图面布局美观；
• 如果数据量极大，自定义四分树效率较低，可以用pysal的空间划分模块优化；
• 若需要非正方形单元格，可修改拆分逻辑，允许x、y方向的拆分步长独立调整。

内容的提问来源于stack exchange，提问作者Kat0mi