Overlaying a background with hexagonal points on ggplot r

I'm desperately trying to add inventory data (available online, as in this reproducible example below) in ggplot with a background of France, so that I can visualize the areas where there are no records of my species.

However, I want to do it with hexagons and not points. So, I used stat_summary_hex which prevents me from using geom_polygon to add the background to my ggplot. Then, I tried to overlay two ggplots, without success.

An exemple with Abies alba

## Extract data
temp <- tempfile()
download.file("https://inventaire-forestier.ign.fr/dataifn/data/export_dataifn_2005_2022.zip",temp)
couvert <- read.csv(unz(temp, "COUVERT.csv"), sep=";")
placette <- read.csv(unz(temp, "PLACETTE.csv"), sep=";")
espar <- read.csv(unz(temp,"espar-cdref13.csv"), sep=";")
unlink(temp)

## e.g. for Abies alba
ID <- espar$X...espar[espar$lib_cdref=="Abies alba"]
couvert_sp <- couvert[couvert$ESPAR_C == ID ,]
df <- merge(placette[ , c("IDP", "XL", "YL")], couvert_sp, by="IDP", all.y=T)

# coordinates conversion (required)
crs193 <- "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs"
wgs84 <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

xyl93 <- df[, c("XL", "YL", "IDP")]
require(sp)
coordinates(xyl93) <- ~ XL + YL
proj4string(xyl93) <- CRS(crs193)
xyl <- as.data.frame(spTransform(xyl93, wgs84))

require(dplyr)
df2 <- merge(xyl, df, by="IDP", all.x=T, all.y=T) %>% distinct(IDP, .keep_all = TRUE)

# graphic
require(ggplot2)
AB <- map_data("world") %>% filter(region %in% "France")

# background 
bck <- ggplot() +
  geom_polygon(data = AB, aes(x=long, y = lat, group = group), fill="#ECECEC") + 
  theme_void() + coord_map() 
bck

# my hex points
require(paletteer)
p <- ggplot(df2, aes(x=coords.x1, y=coords.x2, z=TCA)) +
  stat_summary_hex(fun = mean, bins = 80) +
  scale_fill_gradientn(colours = paletteer_c("ggthemes::Classic Green", 30)) +
  theme_void() + coord_map() 
p

p + bck

I want to overlay these two plots:

Any suggestions? Thanks in advance!

Solution

Instead of adding the two ggplots you could add the layers for the hex bin plot to your background plot like so:

library(ggplot2)
library(paletteer)

bck +
  stat_summary_hex(
    data = df2, aes(x = coords.x1, y = coords.x2, z = TCA),
    fun = mean, bins = 80
  ) +
  scale_fill_gradientn(colours = paletteer_c("ggthemes::Classic Green", 30))