How to add contours to a density scatter plot in R?

Question

I have a code that generates a density scatter plot like this (found here How to add an ellipse to a scatter plot colored by density in R? ) :

## Data in a data.frame
x1 <- rnorm(n=1000, sd=2)
x2 <- x1*1.2 + rnorm(n=1000, sd=2)
df <- data.frame(x1,x2)

## Use densCols() output to get density at each point
x <- densCols(x1,x2, colramp=colorRampPalette(c("black", "white")))
df$dens <- col2rgb(x)[1,] + 1L

## Map densities to colors
cols <-  colorRampPalette(c("#FF3100", "#FF9400", "#FCFF00", 
                            "#45FE4F", "#00FEFF", "#000099"))(6)
df$col <- ifelse(df$dens >= 250, cols[1], ifelse(df$dens >= 200, cols[2], ifelse(df$dens >= 150, cols[3], ifelse(df$dens >= 100, cols[4], ifelse(df$dens >= 50, cols[5], cols[6])))))

## Plot it, reordering rows so that densest points are plotted on top
plot(x2~x1, data=df[order(df$dens),], pch=20, col=col, cex=1)

I'd like to add contours exactly like in the image below (found here https://fr.mathworks.com/matlabcentral/fileexchange/8577-scatplot) :

I have to do this without ggplot2.

Does anyone have an idea please?

Answer 1

These contour lines show the estimated density of the data at given levels. You may compute the density estimate using MASS::kde2d() and then plot using contour():

biv_kde <- MASS::kde2d(x1, x2)
contour(biv_kde, add = T)

Edit:

Here's another approach to get colored density contours using KernSmooth:

(Note that I used set.seed(123) before generating the data.)

# estimate bandwidths
h <- c(dpik(df$x1), dpik(df$x2))

# obtain density estimatte
f1 <- bkde2D(df[, 1:2], bandwidth = h, gridsize = c(200, 200))

# setup vector of density levels to obtain contours for
contour_levels <- pretty(f1$fhat, 9)

# setup color palette
crp <- colorRampPalette(rev(c("#FF3100", "#FF9400", "#FCFF00", 
                              "#45FE4F", "#00FEFF", "#000099")))

# density based colors 
df$col <- densCols(x1, x2, bandwidth = h, colramp = crp)

# plot
plot(x2 ~ x1, data = df, pch = 20, col = col, cex = 0.5)
contour(f1$x1, f1$x2, f1$fhat, levels = contour_levels, col = crp(9),
        lwd = 2, add = T)