Plot 2D histogram of 3D+ multivariate data

The following seems to slice it, but please further QC its correctness:

for j in 1:nbins2
    x = getindex(hist.edges, 3)
    x = (x[1:end-1] + x[2:end])/2
    y = getindex(hist.edges, 1)
    y = (y[1:end-1] + y[2:end])/2
    z = hist.weights[:,j,:]
    display(heatmap(x,y,z, title="(i, j=$j, k)"))
    println("[ENTER] to continue...")
    readline()
end