Misalignment between a raster overlay and DEM in rayshader 3D maps

I am using the rayshader package to create a 3D terrain visualization of a certain area. My R code is as follows:

# ─────────────────────────────────────────────────────────
# 1. Setup
# ─────────────────────────────────────────────────────────
pacman::p_load(
  osmdata, rstac, tidyverse,
  terra, sf, elevatr, ggtern,
  ggspatial, rayshader
)

# ─────────────────────────────────────────────────────────
# 2. Retrieve Al Hoceima Province Boundary
# ─────────────────────────────────────────────────────────
alhoceima <- osmdata::opq(bbox = "Al Hoceima") %>% 
  osmdata::add_osm_feature(
    key = "name:en", value = "Al Hoceima Province"
  )

osm_poly <- osmdata::osmdata_sf(alhoceima)

province_sf <- if (nrow(osm_poly$osm_multipolygons) > 0) {
  osm_poly$osm_multipolygons
} else {
  osm_poly$osm_polygons
}

province_sf <- sf::st_transform(province_sf, 4326)
province_bbox <- sf::st_bbox(province_sf)

plot(sf::st_geometry(province_sf))

# ─────────────────────────────────────────────────────────
# 3. Query ESA WorldCover Data
# ─────────────────────────────────────────────────────────
ms_query <- rstac::stac("https://planetarycomputer.microsoft.com/api/stac/v1")

collections <- "esa-worldcover"

ms_esa_query <- rstac::stac_search(
  q = ms_query,
  collections = collections,
  datetime = "2021-01-01T00:00:00Z/2021-12-31T23:59:59Z",
  bbox = province_bbox
) |> rstac::get_request()

# ─────────────────────────────────────────────────────────
# 4. Download Land Cover Data
# ─────────────────────────────────────────────────────────
ms_query_signin <- rstac::items_sign(
  ms_esa_query, rstac::sign_planetary_computer()
)

main_dir <- getwd()

rstac::assets_download(
  items = ms_query_signin,
  asset_names = "map",
  output_dir = main_dir,
  overwrite = TRUE
)

# ─────────────────────────────────────────────────────────
# 5. Load and Crop Land Cover Raster
# ─────────────────────────────────────────────────────────
version <- "v200"
year <- "2021"
asset_name <- "map"

data_dir <- file.path(main_dir, collections, version, year, asset_name)
raster_files <- list.files(data_dir, full.names = TRUE)

land_cover_raster <- terra::vrt(raster_files)

province_sf <- sf::st_transform(province_sf, crs = sf::st_crs(land_cover_raster))

province_land_cover <- terra::crop(
  land_cover_raster,
  province_sf,
  snap = "in",
  mask = TRUE
)

terra::plot(province_land_cover)

# ─────────────────────────────────────────────────────────
# 6. Extract Colors from Raster
# ─────────────────────────────────────────────────────────
raster_color_table <- do.call(
  data.frame, terra::coltab(province_land_cover)
)

hex_code <- ggtern::rgb2hex(
  r = raster_color_table[, 2],
  g = raster_color_table[, 3],
  b = raster_color_table[, 4]
)

cols <- hex_code[!hex_code == "#000000"]
province_land_cover <- na.omit(province_land_cover)

# ─────────────────────────────────────────────────────────
# 7. Get and Resample DEM for Province
# ─────────────────────────────────────────────────────────
dem <- elevatr::get_elev_raster(
  locations = province_sf,
  z = 9,
  clip = "locations"
) |> terra::rast()

proj <- "EPSG:4668"

province_land_cover_resampled <- terra::resample(
  x = province_land_cover,
  y = dem,
  method = "near"
) |> terra::project(
  proj,
  method = "near"
)

# ─────────────────────────────────────────────────────────
# 8. Theme for Plotting
# ─────────────────────────────────────────────────────────
theme_for_the_win <- function(){
  theme_minimal() +
    theme(
      axis.line = element_blank(),
      axis.title = element_blank(),
      axis.text = element_blank(),
      panel.grid = element_blank(),
      plot.background = element_rect(fill = "white", color = NA),
      plot.title = element_text(size = 20, hjust = .5, vjust = -2),
      plot.caption = element_text(size = 9, hjust = .5, vjust = 3),
      plot.margin = unit(rep(0,4), "lines")
    )
}

# ─────────────────────────────────────────────────────────
# 9. Prepare Data for Plotting
# ─────────────────────────────────────────────────────────
values <- seq(10, 100, by = 10)
values <- append(values, 95, after = 9)

labels <- c(
  "Tree cover", "Shrubland", "Grassland",
  "Cropland", "Built-up", "Bare/sparse vegetation",
  "Snow and Ice", "Permanent water bodies",
  "Herbaceous wetland", "Mangroves", "Moss and lichen"
)

codebook <- data.frame(
  values = values,
  labels = labels,
  cols = cols
)

province_land_cover_resampled_df <- as.data.frame(
  province_land_cover_resampled,
  xy = TRUE, na.rm = TRUE
)

names(province_land_cover_resampled_df)[3] <- "values"
codebook$values <- as.numeric(as.character(codebook$values))

province_land_cover_df <- dplyr::left_join(
  province_land_cover_resampled_df,
  codebook,
  by = "values"
)

province_land_cover_df$values <- as.factor(province_land_cover_df$values)

labs <- labels[labels %in% unique(province_land_cover_df$labels)]
pal <- cols[cols %in% unique(province_land_cover_df$cols)]

# ─────────────────────────────────────────────────────────
# 10. Create 2D Land Cover Map
# ─────────────────────────────────────────────────────────
lc_map <- ggplot() +
  geom_raster(
    data = province_land_cover_df,
    aes(x = x, y = y, fill = as.factor(values))
  ) +
  scale_fill_manual(
    name = "Land cover classes",
    values = pal,
    labels = labs
  ) +
  annotation_scale(location = "bl", width_hint = 0.15) +
  annotation_north_arrow(location = "tr", style = north_arrow_fancy_orienteering) +
  coord_sf(crs = proj) +
  labs(
    title = "Land Cover | Al Hoceima Province",
    caption = "©2024 Adapted from ESA WorldCover data"
  ) +
  theme_for_the_win() +
  theme(
    legend.background = element_blank(),
    # legend.position = "topleft",
    legend.position = c(.15, .72),
    legend.title = element_text(size = 9),
    legend.text = element_text(size = 8)
  )

ggsave(
  "alhoceima-landcover-2d.png", lc_map,
  width = 7, height = 7, bg = "white"
)

# ─────────────────────────────────────────────────────────
# 11. Create DEM Map
# ─────────────────────────────────────────────────────────
dem_df <- as.data.frame(
  terra::project(dem, proj),
  xy = TRUE, na.rm = TRUE
)

names(dem_df)[3] <- "dem"

dem_map <- ggplot() +
  geom_raster(
    data = dem_df,
    aes(x = x, y = y, fill = dem)
  ) +
  scale_fill_gradientn(colors = "white") +
  guides(fill = "none") +
  coord_sf(crs = proj) +
  theme_for_the_win()

# ─────────────────────────────────────────────────────────
# 12. Render 3D Rayshader Map
# ─────────────────────────────────────────────────────────
rayshader::plot_gg(
  ggobj = lc_map,
  ggobj_height = dem_map,
  width = 7, height = 7,
  windowsize = c(600, 600),
  scale = 100,
  shadow = TRUE,
  shadow_intensity = 1,
  phi = 87,
  theta = 0,
  zoom = .55,
  multicore = TRUE
)

u <- "https://dl.polyhaven.org/file/ph-assets/HDRIs/hdr/4k/brown_photostudio_02_4k.hdr"
hdri_file <- basename(u)
download.file(url = u, destfile = hdri_file, mode = "wb")

rayshader::render_highquality(
  filename = "alhoceima-esa-3d.png",
  preview = TRUE,
  light = FALSE,
  environment_light = hdri_file,
  intensity = 1,
  rotate_env = 90,
  parallel = TRUE,
  width = 2000,
  height = 2000,
  interactive = FALSE
)

Problem:

In the resulting 3D map, the landcover layer is visibly misaligned with the elevation.

Question:

What is the correct workflow to ensure precise alignment between the ggplot2-based raster overlay and the rayshader heightmap?