r ggplot2 stat_density2d - how to change the bandwidth calculation to Sheather and Jones method

Question

I'm creating two maps, one with contour lines and the other is a heatmap. I'm using the geom_density2d function to create the one with contour lines, and stat_density2d to create the heatmap.

Note that geom_density2d does not seem to have the same parameters as geom_density, which now includes a bw argument to specify the bandwidth.

I understand (after a lot of searching) that both of these functions are defaulting to "Silverman's rule of thumb" bw.nrd0 to calculate the bandwidth, as described in this blog post. However I would like to change the bandwidth calculation to use the Sheather and Jones method instead bw = "SJ" but have no idea how to do this as this parameter is not accepted by ggplot2 (neither the geom nor the stat functions).

Example data is provided here.

The code I'm using to create the contour map is below:

# Create the contour map with geom_density2d:
contourmap <- ggplot2::ggplot() + 
  geom_sf(data = phecll, 
          colour = "black", 
          fill = "#004529",
          alpha = 0.5,
          size = 0.75) +
  coord_sf() + 
  geom_density2d(data = rdmdata, 
                 mapping = aes(x = home_long, 
                               y = home_lat, 
                               alpha = 0.5),
                 inherit.aes = FALSE, 
                 contour_var = "density")

And for the heatmap:

# Create the heatmap with stat_density2d:
heatmap <- ggplot2::ggplot() + 
  geom_sf(data = phecll, 
          colour = "black", 
          fill = "#004529",
          alpha = 0.5,
          size = 0.75) +
  coord_sf() + 
  stat_density2d(data = rdmdata, 
                 mapping = aes(x = home_long, 
                               y = home_lat, 
                               fill = ..level.., 
                               alpha = ..level..), 
                 size = 0.01,  
                 bins = 50,
                 geom = "polygon", 
                 inherit.aes = FALSE) + 
  scale_fill_gradient(low = "blue", 
                      high = "red", 
                      guide = "none") + 
  scale_alpha(range = c(0, 0.5),
              guide = "none")

How can I change the bandwidth calculation from nrd0 to sj?