Density-2d plot on top of a map with ggmap

Question

I want to produce a 2d-density plot based on spatial point data. In the background I want to show an open map (e.g. stamen terrain). Besides I want to plot the borders of Austria. Both datasets (data points and border) are shapefiles in EPSG 4326.

I managed to produce such a plot (see screenshots and code V1 below), but the problem is that there is a shift between the map in the background on the one side and the plotted points and the borders of Austria on the other side, as you can see below.

2D-Density Plot V1 - full
2D-Density Plot V1 - detail

Here is the code (V1):

library(sf)
library(rgdal)
library(ggplot2)
library(ggmap)

# Read point data (EPSG: 4326)
sk <- st_read("points.shp") 

# Read country border polygon (EPSG: 4326)
blogr4326 <- readOGR(<path>, <layer name>)
bl4326_df <- fortify(blogr4326)

# Austria box (extent)
# Longitude: 9.6 to 16.94504
# Latitude: 46.52694 to 48.81667

map <- get_map(c(left = +9.6, bottom = 46.52694, right = +16.90, top = 48.99), color = "color", crop = FALSE)

hm_sk <- ggmap(map, extent = "panel", maprange=FALSE, darken=0.0) + 
  geom_point(data = sk, aes(x=X_WGS84, y=Y_WGS84)) +
  stat_density2d(data = sk, aes(x=X_WGS84, y=Y_WGS84, fill = ..density.., alpha=cut(..density..,breaks=c(-Inf,0.08,Inf))), contour = FALSE, bins=16, geom = 'raster', n=500) + 
  ggtitle("Schwarzkiefer 2016/2020") + xlab("X_WGS84") + ylab("X_WGS84") +
  scale_fill_distiller(palette= "Spectral", direction=-1, limits = c(0.08, 8.50)) + 
  scale_alpha_manual(values=c(0,0.7), guide="none") +
  geom_polygon(data=bl4326_df, aes(long, lat, group=group), color='black', fill='NA', inherit.aes = TRUE) +
  coord_fixed(1.5)

hm_sk

I found out that the shift is caused by the fact that the map in the background is in the projection EPSG:3857 and my shapefiles are in the projection EPSG:4326, as explained in this post. So I projected my shapefiles to EPSG 3857 and inserted the provided code into my code, as you can see here (V2):

library(sf)
library(rgdal)
library(ggplot2)
library(ggmap)

# Read point data (EPSG: 3857)
sk <- st_read("points.shp")

# Read country border polygon (EPSG: 3857)
blogr3857 <- readOGR(<path>, <layer name>)
bl3857_df <- fortify(blogr3857)

# Austria box (extent)
# Longitude: 9.6 to 16.94504
# Latitude: 46.52694 to 48.81667

map <- get_map(c(left = +9.6, bottom = 46.52694, right = +16.90, top = 48.99), color = "color", crop = FALSE)

#-------------------------------------------------------------------------------------------------
# Following code according to this link to avoid the shift between map and country border polygon:
# https://stackoverflow.com/questions/47749078/how-to-put-a-geom-sf-produced-map-on-top-of-a-ggmap-produced-raster

# Define a function to fix the bbox to be in EPSG:3857
ggmap_bbox <- function(map) {
  if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
  # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector, 
  # and set the names to what sf::st_bbox expects:
  map_bbox <- setNames(unlist(attr(map, "bb")), 
                       c("ymin", "xmin", "ymax", "xmax"))
  
  # Convert the bbox to an sf polygon, transform it to 3857, 
  # and convert back to a bbox (convoluted, but it works)
  bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), 3857))
  
  # Overwrite the bbox of the ggmap object with the transformed coordinates 
  attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
  attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
  attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
  attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
  map
}

# Use the function:
map <- ggmap_bbox(map)
#-------------------------------------------------------------------------------------------------


hm_sk <- ggmap(map,extent = "device", maprange=FALSE) +#, extent = "panel", maprange=FALSE, darken=0.0) + 
  coord_sf(crs = st_crs(3857)) + # f867orce the ggplot2 map to be in 3857
  geom_point(data = sk, aes(x=X_PM, y=Y_PM)) +
  stat_density2d(data = sk, aes(x=X_PM, y=X_PM, fill = ..density.., alpha=cut(..density..,breaks=c(-Inf,0.08,Inf))), contour = FALSE, bins=16, geom = 'raster', n=500) + 
  scale_fill_distiller(palette= "Spectral", direction=-1, limits = c(0.08, 8.50)) +
  scale_alpha_manual(values=c(0,0.7), guide="none") +
  geom_polygon(data=bl3857_df, aes(long, lat, group=group), color='black', fill='NA', inherit.aes = FALSE) +
  ggtitle("Schwarzkiefer 2016/2020") + xlab("X_3857") + ylab("X_3857")
  
hm_sk

Now, the problem with the shift is solved, but the density plot is not visible anymore (only map, points and borders are plotted), as you can see here:

2D-Density Plot V2 - full
2D-Density Plot V2 - detail

Any suggestions, how I can produce a plot that is proper aligned AND includes the density plot? Thanks a lot in advance!