如何在R中用ggplot2和sf在阿拉斯加地图上叠加密度图？

问题分析 #

你的核心问题在于坐标系不一致和密度图未按阿拉斯加边界裁剪：

• 美国底图usa仍使用WGS84（EPSG:4326），但密度图使用了阿拉斯加阿尔伯斯投影（EPSG:3467），两者坐标系不匹配导致叠加失败
• 直接用手动bbox生成的密度图会覆盖整个矩形范围，而非仅阿拉斯加区域

针对你的问题的解决方案 #

1. 如何将密度图叠加到阿拉斯加区域？ #

需要确保底图和密度图坐标系统一，同时用阿拉斯加的边界裁剪密度图，步骤如下：

• 将美国底图转换为与密度图相同的投影（EPSG:3467）
• 从美国地图中提取阿拉斯加的单独边界
• 用阿拉斯加边界裁剪密度图，只保留阿拉斯加范围内的密度值
• 按「底图→密度图→随机点」的顺序绘制图层

2. 是否遗漏步骤或投影错误？ #

你选择的EPSG:3467是阿拉斯加的专用投影，这个选择是正确的。但遗漏了两个关键步骤：

• 未将美国底图转换到目标投影
• 未用阿拉斯加边界约束密度图的范围

修正后的完整代码 #

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

# Load world data and filter for USA, then transform to Alaska's projection (EPSG:3467)
world <- ne_countries(scale='medium', returnclass = 'sf')
usa <- subset(world, admin == "United States of America") %>%
  st_transform(3467)  # 统一坐标系为阿拉斯加投影

# 提取阿拉斯加单独的边界
alaska_boundary <- usa %>%
  filter(subregion == "Alaska")

# 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)  # 转换到阿拉斯加投影

# Convert sf object to ppp object, 使用阿拉斯加边界作为窗口而非手动bbox
ppp_points <- as.ppp(sf_points)
Window(ppp_points) <- as.owin(alaska_boundary)  # 用阿拉斯加边界定义分析窗口

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

# Convert density to stars, 裁剪到阿拉斯加边界，再转为sf object
density_stars <- st_as_stars(density_spatstat) %>%
  st_crop(alaska_boundary)  # 只保留阿拉斯加范围内的密度
density_sf <- st_as_sf(density_stars) %>%
  st_set_crs(3467)

# Plot Alaska map with density overlay
(alaska_plot <- ggplot() +
    geom_sf(data = alaska_boundary, fill = "lightgray", color = "black") +  # 先画阿拉斯加底图
    geom_sf(data = density_sf, aes(fill = v), color = NA, alpha = 0.7) +  # 叠加密度图，调整透明度
    geom_sf(data = sf_points, size = 1, color = "red") +  # 叠加随机点（可选）
    scale_fill_viridis_c(name = "Point Density") +  # 用更美观的颜色标尺
    coord_sf(xlim = st_bbox(alaska_boundary)[c(1,3)], 
             ylim = st_bbox(alaska_boundary)[c(2,4)], 
             expand = FALSE, datum = NA) +
    theme_minimal() +
    labs(title = "Alaska Random Point Density"))

代码关键修改说明 #

• 统一坐标系：将美国底图和随机点都转换为EPSG:3467，确保所有图层坐标系一致
• 用阿拉斯加边界定义分析窗口：替代手动bbox，让密度计算仅针对阿拉斯加区域
• 裁剪密度图：用st_crop保留阿拉斯加范围内的密度值，避免密度图超出边界
• 调整图层顺序：底图→密度图→点，确保视觉层级正确
• 优化可视化：添加颜色标尺、调整透明度，提升图表可读性

内容的提问来源于stack exchange，提问作者Ali Roghani