如果你想要的是装箱，那么binned_statistic_2d可能适合你。下面是一个例子：

from scipy.stats import binned_statistic_2d
import numpy as np

x = np.random.uniform(0, 10, 1000)
y = np.random.uniform(10, 20, 1000)
z = np.exp(-(x-3)**2/5 - (y-18)**2/5) + np.random.random(1000)

x_bins = np.linspace(0, 10, 10)
y_bins = np.linspace(10, 20, 10)

ret = binned_statistic_2d(x, y, z, statistic=np.mean, bins=[x_bins, y_bins])

fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(12, 4))
ax0.scatter(x, y, c=z)
ax1.imshow(ret.statistic.T, origin='bottom', extent=(0, 10, 10, 20))

@Andrea的回答非常清楚和有帮助，但我想提一个更快的替代方案，不使用scipy库。
我们的想法是做一个x和y的二维直方图，用z变量加权（它有每个bin中z变量的总和），然后根据没有权重的直方图进行归一化（它有每个bin中的计数）。这样，你将计算每个bin中z变量的平均值。
代码：

import numpy as np
import matplotlib.pyplot as plt

x = np.random.uniform(0, 10, 10**7)
y = np.random.uniform(10, 20, 10**7)
z = np.exp(-(x-3)**2/5 - (y-18)**2/5) + np.random.random(10**7)

x_bins = np.linspace(0, 10, 50)
y_bins = np.linspace(10, 20, 50)

H, xedges, yedges = np.histogram2d(x, y, bins = [x_bins, y_bins], weights = z)
H_counts, xedges, yedges = np.histogram2d(x, y, bins = [x_bins, y_bins]) 
H = H/H_counts

plt.imshow(H.T, origin='lower',  cmap='RdBu',
            extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])
plt.colorbar()

在我的计算机中，这种方法比使用scipy的binned_statistic_2d快**5倍。