Not a ggplot2 user, but basically you want to estimate a weighted 2d-density and make an image out of it. Your linked answer indicates that ggplot2::geom_density2d internally uses MASS::kde2d, but which only computes unweighted 2d-densities.
Inflate observations
Similar to @AllanCameron's suggestion (but without the need to use tidyr) we could inflate the data frame simply by copying each row by the number of milliseconds duration,
dfa <- df[rep(seq_len(nrow(df)), times=df$dur), -3]
and calculate the kde2d by hand.
n <- 1e3
system.time(
dens1 <- MASS::kde2d(dfa$x, dfa$y, n=n) ## this runs a while!
)
# user system elapsed
# 2253.285 2325.819 661.632
The n= argument denotes the number of grid points in each direction, the greater we choose it, the smoother the granularity will look in the heatmap image.
system.time(
dens1 <- MASS::kde2d(dfa$x, dfa$y, n=n) ## this runs a while
)
# user system elapsed
# 2253.285 2325.819 661.632
image(dens1, col=heat.colors(n, rev=TRUE))

This almost runs forever, though with n=1000...
Weighted 2d density estimation
In a comment on the answer mentioned above, @IRTFM links an ancient r-help post that provides a kde2d.weighted function which is lightning fast and that we could try (see code at the bottom).
dens2 <- kde2d.weighted(x=df$x, y=df$y, w=proportions(df$dur), n=n)
image(dens2, col=heat.colors(n, rev=TRUE))

However, the two versions look quite different, and I can't tell which is right, since I am not really an expert with this method. But at least there is a noticeable difference to the unweighted image:
Unweighted image
dens0 <- MASS::kde2d(df$x, df$y, n=n)
image(dens0, col=heat.colors(n, rev=TRUE))

Points
Still adding the points might be pointless, but you could run this line after image:
points(y ~ x, df, cex=proportions(dur)*2e3, col='green')
Taken from r-help (Ort 2006):
kde2d.weighted <- function(x, y, w, h, n=n, lims=c(range(x), range(y))) {
nx <- length(x)
if (length(y) != nx)
stop("data vectors must be the same length")
gx <- seq(lims[1], lims[2], length=n) ## gridpoints x
gy <- seq(lims[3], lims[4], length=n) ## gridpoints y
if (missing(h))
h <- c(MASS::bandwidth.nrd(x), MASS::bandwidth.nrd(y))
if (missing(w))
w <- numeric(nx) + 1
h <- h/4
ax <- outer(gx, x, "-")/h[1] ## distance of each point to each grid point in x-direction
ay <- outer(gy, y, "-")/h[2] ## distance of each point to each grid point in y-direction
z <- (matrix(rep(w,n), nrow=n, ncol=nx, byrow=TRUE)*
matrix(dnorm(ax), n, nx)) %*%
t(matrix(dnorm(ay), n, nx))/(sum(w)*h[1]*h[2]) ## z is the density
return(list(x=gx, y=gy, z=z))
}