cropping stat_density_2d output to a polygon

Question

I am trying to use the stat_density_2d() function to create a heat map based on GPS coordinates and would like to crop the output to a polygon in ggplot(). Here is an example:

x_coord <- c(-83, -82, -82, -83, -83)
y_coord <- c(41, 41, 42, 42, 41)

xy_coords <- cbind(x_coord, y_coord)

library(sp)

poly <- Polygon(xy_coords)
polys <- Polygons(list(poly), 1)
sp.polys <- SpatialPolygons(list(polys))

set.seed(2)
lon <- c(rnorm(20, -82.2, 0.1), rnorm(20, -82.1, 0.04))
lat <- c(rnorm(20, 41.2, 0.1), rnorm(20, 41.1, 0.02))

lon_lat <- data.frame(lon = lon,
                      lat = lat)

ggplot() +
  geom_polygon(data = sp.polys, aes(x_coord, y_coord), fill = 'transparent', color = 'black') +
  stat_density_2d(data = lon_lat, aes(x = lon, y = lat, fill = ..level.., alpha = 0.3), geom = 
  'polygon') +
  scale_fill_gradientn(colours=rev(brewer.pal(7, "Spectral"))) +
  scale_x_continuous(limits = c(-83.25, -81.75)) +
  scale_y_continuous(limits = c(40.75, 42.25))

which produces this figure

Is there a way to crop the stat_density output to remove the portion outside the square? Any thoughts or suggestions are greatly appreciated.

Answer 1

I was unable to find a way to crop the stat_density_2d output to a polygon but using the code below, can mask the area outside the polygon. The code below assumes that you have already used the code posted in the questions to create the polygon and stat_density_2d output.

library(sp)
library(ggplot2)
library(RColorBrewer)
library(rgeos)
library(raster)

x_coord2 <- c(-83.2, -81.8, -81.8, -83.2, -83.2)
y_coord2 <- c(40.8, 40.8, 42.2, 42.2, 40.8)

xy_coords2 <- cbind(x_coord2, y_coord2)

poly2 <- Polygon(xy_coords2)
polys2 <- Polygons(list(poly2), 1)
sp.polys2 <- SpatialPolygons(list(polys2))

r <- raster(x = extent(sp.polys2))
res(r) <- 0.003
r <- setValues(r, 1)
r <- mask(r, sp.polys, inverse = T)

rdf <- data.frame(rasterToPoints(r))


ggplot() +
  stat_density_2d(data = lon_lat, aes(x = lon, y = lat, fill = ..level..), alpha = 
                  0.7, geom = 'polygon') +
  scale_fill_gradientn(colours=rev(brewer.pal(7, "Spectral"))) +
  scale_x_continuous(limits = c(-83.25, -81.75)) +
  scale_y_continuous(limits = c(40.75, 42.25)) +
  geom_tile(data = rdf, aes(x, y), fill = 'white') +
  geom_polygon(data = sp.polys, aes(x_coord, y_coord), fill = 'transparent', color = 
               'black') +
  theme_bw() +
  theme(panel.grid.minor = element_blank(), 
        panel.grid.major = element_blank())