I am trying to produce a heat map based on point locations using stat_density2d_filled. Using stat_density2d results in ugly artifacts. stat_density2d works well with scale_fill_gradient, but when trying the same things with stat_density2d_filled I get an error: “Error: Discrete value supplied to continuous scale”.

Please see my example code below. The first two plots work, but the third one does not.

rm(list = ls())

library(tidyverse)
library(tidygeocoder)

no <- ne_states(geounit = "Norway", returnclass = "sf") 

no_south <- no %>%
  filter(name == "Aust-Agder" | 
           name == "Vest-Agder" |
           name == "Telemark") 

towns <- tribble(
  ~town, 
  "Kristiansand, NO",   "Vest-Agder",       
  "Søgne, NO",     "Vest-Agder",                    
  "Vennesla, NO",
  "Arendal, NO",
  "Lillesand Town, NO",
  "Porsgrunn, NO",
  "Skien, NO") %>%
  tidygeocoder::geocode(town, method = "osm", 
                        lat = latitude , long = longitude) 

# Map with stat_density_2d and scale_fill_gradient works:
ggplot(data = no_south) +
  geom_sf() +
  stat_density_2d(data = towns,
                        aes(x = longitude, y = latitude,
                            fill = after_stat(level)),
                  geom = "polygon",
                  alpha = 0.5) +
  scale_fill_gradient() +
  geom_point(data = towns,
             aes(x = longitude, y = latitude),
             size = rel(3),
             shape = 20) 

# Map with stat_density_2d_filled and withouth scale_fill_gradient works:
ggplot(data = no_south) +
  geom_sf() +
  stat_density2d_filled(data = towns,
                        aes(x = longitude, y = latitude,
                            fill = after_stat(level)),
                        alpha = 0.5) +
  # scale_fill_gradient() +
  geom_point(data = towns,
             aes(x = longitude, y = latitude),
             size = rel(3),
             shape = 20) 

# Map with stat_density_2d_filled and scale_fill_gradient does not work:
ggplot(data = no_south) +
  geom_sf() +
  stat_density2d_filled(data = towns,
                        aes(x = longitude, y = latitude,
                            fill = after_stat(level)),
                        alpha = 0.5) +
  scale_fill_gradient() +
  geom_point(data = towns,
             aes(x = longitude, y = latitude),
             size = rel(3),
             shape = 20)