library(tidyverse)
## ── Attaching packages ────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.1.0       ✔ purrr   0.3.1  
## ✔ tibble  2.0.1       ✔ dplyr   0.8.0.1
## ✔ tidyr   0.8.3       ✔ stringr 1.4.0  
## ✔ readr   1.3.1       ✔ forcats 0.4.0
## ── Conflicts ───────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(viridis)
## Loading required package: viridisLite
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(tictoc)

generate data

a <- tibble( x=rnorm(20000, 8.2, 1.3), y=rnorm(20000, 9, 1.6) )
b <- tibble( x=rnorm(20000, 20.5, 2), y=rnorm(20000, 16.5, 2.2) )
c <- tibble( x=rnorm(20000, 9.4, 1.3), y=rnorm(20000, 14.4, 2.6) )
xy <- rbind(a,b,c)

dot plot

ggplot(xy) + 
  geom_point(aes(x, y), size=0.2, pch=16, alpha=0.5) + 
  #geom_density_2d(aes(x, y), col="cyan") +
  coord_fixed() +
  theme_classic()

viridis plot using base graphics

z <- kde2d(xy$x, xy$y, n=200)
par(pty="s")
image(z, col=viridis_pal()(200))

#points(xy$x, xy$y, col="white", pch=".")

slower version

# make_xyz <- function(xy, n=200){
#   z <- kde2d(xy[[1]], xy[[2]], n=n)
#   d <- tibble()
#   for (i in 1:length(z$x)){
#     for (j in 1:length(z$y)){
#       d <- rbind(d, tibble(x=z$x[i], y=z$y[j], z=z$z[i,j]))
#     }
#   }
#   d
# }
# tic()
# d <- make_xyz(xy, 200)
# toc()
# 
# # 67.2 sec elapsed

faster version

make_xyz <- function(x, y, n=200){
  z <- MASS::kde2d(x, y, n=n)
  d <- matrix(nrow=n*n, ncol=3)
  k = 1
  for (i in 1:n){
    for (j in 1:n){
      d[k,] <- c(z$x[i], z$y[j], z$z[i,j])
      k <- k + 1
    }
  }
  colnames(d) <- c("x", "y", "z")
  d <- as_tibble(d)
  d
}
# tic()
d <- make_xyz(xy$x, xy$y, 200)
# toc()
# 3.735 sec elapsed
plt <- ggplot(d) +
  geom_raster(aes(x=x, y=y, fill=z)) +
  coord_fixed() +  
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(legend.position='none') +
  scale_fill_viridis()

plt

xy_sub <- xy[seq(1, dim(xy)[1], 100),]
plt + geom_point(aes(x, y), alpha=0.5, pch=1, data=xy_sub, col="white")

plt + 
  geom_density_2d(aes(x,y), data=xy, alpha=0.2, col="white")

plt +
  geom_point(aes(x, y), alpha=0.5, pch=1, data=xy_sub, col="white", size=0.5) +
  geom_density_2d(aes(x,y), data=xy, alpha=0.2, col="white")