How to set alpha by bivariate density and fill by z in geom_raster?

I make a heatplot with geom_raster, something like this:

mround <- function(x,base){base*round(x/base)} # round to nearest 5
df <- data.frame(x= mround(rnorm(10000, 100, 10), 5),
                 y= mround(rnorm(10000, 100, 10), 5))
df$z <- df$x + df$y
library(ggplot2)
ggplot(df) +
  geom_raster(aes(x= x, y= y, fill= z)) +
  theme_classic()

I would like to set the alpha argument for the raster and it should depent on the bivariate density. We can visualize the bivariate density for raster like this:

ggplot() +
  stat_density_2d(data= df, aes(x= x, y= y, alpha = after_stat(density)), geom = "raster", contour= FALSE)

Basically, I want the first plot to have the alpha of the second. I tried all kind of thing, like geom_raster(aes(x= x, y= y, fill= z, alpha = after_stat(density))), but nothing works so far. How can I set the fill by z and alpha by 2d density at the same time?

This is difficult, because in order to calculate density, ggplot will have to aggregate your data, which loses the value of variables other than x and y at each unique point.

One option is to just precalculate density, which will work if z is indeed just a simple function of x and y:

ranges <- c(min(df$x, df$y), max(df$x, df$y))

n <- 1 + diff(ranges)/5

dens <- MASS::kde2d(df$x, df$y, n = n, lims = c(ranges, ranges))

cbind(expand.grid(x = dens$x, y = dens$y), 
      dens = c(dens$z), z = dens$x + dens$y) |>
  ggplot(aes(x, y)) +
  geom_tile(aes(fill = z, alpha = dens)) +
  theme_classic() +
  scale_alpha_continuous(range = c(0, 1)) +
  coord_equal()

If, in your actual data, z is an independent variable, then you may need to interpolate its value using something like interp::interp:

cbind(expand.grid(x = dens$x, y = dens$y), 
      dens = c(dens$z), 
      z = c(interp::interp(x = df$x, y = df$y, z = df$z, 
                           duplicate = "mean", xo = dens$x, yo = dens$y)$z)) |>
  ggplot(aes(x, y)) +
  geom_tile(aes(fill = z, alpha = dens)) +
  theme_classic() +
  scale_alpha_continuous(range = c(0, 1)) +
  coord_equal()

Which will give the same output in this particular instance.