You could interpolate your data onto a 2D grid. There are lots of ways to do this - probably the closest analogy to kernel density estimation would be to interpolate using a radial basis function:
import numpy as np
from scipy.interpolate import Rbf
from matplotlib import pyplot as plt
def f(x, y):
return np.sin(x) + np.cos(2 * y)
# 1D arrays of points
x = np.random.rand(100) * 2 * np.pi
y = np.random.rand(100) * 2 * np.pi
z = f(x, y)
# initialize radial basis function
rb = Rbf(x, y, z)
# interpolate onto a 100x100 regular grid
X, Y = np.mgrid[:2*np.pi:100j, :2*np.pi:100j]
Z = rb(X.ravel(), Y.ravel()).reshape(X.shape)
# plotting
fig, ax = plt.subplots(1, 1)
ax.set_aspect('equal')
ax.hold(True)
m = ax.contourf(X, Y, Z, 20, cmap=plt.cm.Greens)
ax.scatter(x, y, c=z, s=60, cmap=m.cmap, vmin=m.vmin, vmax=m.vmax)
cb = fig.colorbar(m)
cb.set_label('$f(x, y)$', fontsize='xx-large')
ax.set_xlabel('$x$', fontsize='xx-large')
ax.set_ylabel('$y$', fontsize='xx-large')
ax.margins(0.05)
fig.tight_layout()
plt.show()
