consider using a gaussian KDE to estimate spatial rate

[lepmik]

scipy.stats.gaussian_kde

[lepmik]

Working example, but has issues with edges, consider using some these boundary conditions

from scipy.interpolate import interp1d
from scipy.stats import gaussian_kde

sx = interp1d(t, x)(spike_times)
sy = interp1d(t, y)(spike_times)
    
X, Y = np.meshgrid(xbins, ybins)
positions = np.vstack([X.ravel(), Y.ravel()])

xy = np.vstack([x, y])
s_xy = np.vstack([sx, sy])

bw = 0.2

kernel = gaussian_kde(s_xy, bw_method=bw)
smap = np.reshape(kernel(positions), X.shape)

kernel = gaussian_kde(xy, bw_method=bw)
omap = np.reshape(kernel(positions), X.shape)

omap = np.flipud(np.rot90(omap))
smap = np.flipud(np.rot90(smap))
rmap = smap / omap

fig, (ax1, ax2, ax3) = plt.subplots(1,3)
ax1.imshow(omap)
ax2.imshow(smap)
ax3.imshow(rmap)