I’m trying to plot a density map of some randomly generated points in Alaska using ggplot2 and sf. I want to overlay the density on a map of the USA including Alaska. However, I am unable to make geom_density work on my map. Here is what I have done so far:

I generated 100 random points within the latitude and longitude ranges for Alaska.
I converted these points to an sf object and then to a ppp object for density calculation.
I used spatstat to calculate the density and then converted the result to an sf object.
Finally, I plotted the map using ggplot2.
My main problem is integrating the density plot with the map of the USA, particularly getting geom_density or an equivalent to work correctly. I am not sure how to overlay the density on the Alaska map properly.

Here are my specific questions:

How can I overlay the density plot on the Alaska region within the USA map?
Am I missing any steps or using incorrect projections?
Any guidance or suggestions would be greatly appreciated!

Here is the code:

# Load necessary libraries
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)
library(dplyr)
library(ggplot2)
library(spatstat)
library(stars)

# Load world data and filter for USA
world <- ne_countries(scale='medium', returnclass = 'sf')
usa <- subset(world, admin == "United States of America")

# Generate random latitude and longitude coordinates within specified ranges for Alaska
set.seed(123)  # For reproducibility
Latitude_Alaska <- runif(100, min = 54, max = 71)
Longitude_Alaska <- runif(100, min = -170, max = -130)

# Combine latitude and longitude into a data frame and convert to sf object
alaska_points <- data.frame(Latitude = Latitude_Alaska, Longitude = Longitude_Alaska)
sf_points <- alaska_points %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(3467)  # Adjust the projection accordingly

# Define the bounding box for Alaska in the projected CRS
alaska_bbox <- st_bbox(c(xmin = -2500000, xmax = 2000000, ymin = 200000, ymax = 2500000), crs = st_crs(3467))

# Convert sf object to ppp object
ppp_points <- as.ppp(sf_points)
Window(ppp_points) <- as.owin(alaska_bbox)

# Plot the points
par(mar = rep(0, 4))
plot(ppp_points, main = "")

# Calculate density
density_spatstat <- density(ppp_points, dimyx = 256)

# Convert density to stars and then to sf object
density_stars <- st_as_stars(density_spatstat)
density_sf <- st_as_sf(density_stars) %>%
  st_set_crs(3467)


# Plot USA map with Alaska
(alaska <- ggplot(data = usa) +
     geom_sf(data = density_sf, aes(fill = v), col= NA) +
     coord_sf(crs = st_crs(3467), xlim = c(-2400000, 1600000), ylim = c(200000, 
         2500000), expand = FALSE, datum = NA))