使用ggplot2与sf包可视化空间密度时的难题

问题描述 #

空间密度可视化中，采样点稀疏区域（如示例中的犹他州）无法在密度图中清晰呈现，即便使用scale_fill_viridis_c调整配色，视觉对比度仍不足以区分密度差异。

现有代码回顾 #

数据代码 #

sf_points <- data.frame(
  lat = c(39.52, 39.62, 40.03, 40.00, 39.93, 39.94, 40.12, 40.54, 35.78, 35.77),
  lon = c(-116.32, -116.00, -116.42, -116.40, -116.41, -116.41, -116.59, -116.56, -111.89, -111.51)
)

实现与绘图代码（修正原代码笔误） #

library(sf)
library(dplyr)
library(ggplot2)
library(spatstat)
library(stars)
library(ggthemes)

sf_points <- sf_points %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  st_transform(32650)

ppp_points <- as.ppp(sf_points)
Window(ppp_points) <- as.owin(usa_map)

par(mar = rep(0, 4))
plot(ppp_points, main = "")

density_spatstat <- density(ppp_points, dimyx = 256)
density_stars <- stars::st_as_stars(density_spatstat)
density_sf <- st_as_sf(density_stars) %>%
  st_set_crs(32650)

# 绘图
ggplot(usa_map) +
  geom_sf(data = density_sf, aes(fill = v), col = NA) +
  theme_bw() +
  ggthemes::theme_map() + 
  geom_sf(data = st_boundary(usa_map)) +
  coord_sf(crs = st_crs("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"), datum = NA) +
  scale_fill_viridis_c(option = "c")

优化方案 #

1. 调整配色方案：增强稀疏区域对比度 #

默认配色对低值区分度不足，可通过以下方式优化：

• 对密度值做对数转换，压缩高值区间，放大低值差异
• 手动设置配色断点，给稀疏区域分配更易区分的颜色
• 使用专为稀疏数据设计的渐变配色

示例代码：

# 对数转换+自定义断点
ggplot(usa_map) +
  geom_sf(data = density_sf, aes(fill = v), col = NA) +
  theme_bw() +
  ggthemes::theme_map() + 
  geom_sf(data = st_boundary(usa_map)) +
  coord_sf(crs = st_crs("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"), datum = NA) +
  scale_fill_viridis_c(
    option = "c",
    trans = "log1p",  # 避免0值报错的对数转换
    breaks = c(0, 1e-6, 1e-5, 1e-4),
    labels = c("0", "1e-6", "1e-5", "1e-4")
  )

# 稀疏数据专用渐变配色
ggplot(usa_map) +
  geom_sf(data = density_sf, aes(fill = v), col = NA) +
  theme_bw() +
  ggthemes::theme_map() + 
  geom_sf(data = st_boundary(usa_map)) +
  coord_sf(crs = st_crs("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"), datum = NA) +
  scale_fill_gradientn(
    colors = c("#f7fbff", "#deebf7", "#c6dbef", "#9ecae1", "#6baed6", "#4292c6", "#2171b5", "#08519c"),
    trans = "log1p",
    name = "密度"
  )

2. 优化密度计算参数：适配稀疏区域 #

spatstat::density()默认带宽可能过度平滑稀疏区域，可调整参数增强信号：

• 使用bw.diggle/bw.ppl等针对稀疏数据的带宽选择方法
• 手动缩小带宽，保留稀疏区域的密度细节

示例代码：

# 稀疏数据专用带宽计算
bw_opt <- bw.diggle(ppp_points)
density_spatstat <- density(ppp_points, dimyx = 256, bw = bw_opt)

# 手动设置带宽（单位：米，对应CRS 32650）
density_spatstat <- density(ppp_points, dimyx = 256, bw = 5000)

# 后续转换与绘图代码不变
density_stars <- stars::st_as_stars(density_spatstat)
density_sf <- st_as_sf(density_stars) %>% st_set_crs(32650)

3. 替代可视化技术：直接突出稀疏区域 #

若密度图仍无法满足需求，可尝试以下方法：

• 核密度+点叠加：在密度图上叠加原始采样点，直观展示稀疏区域的点分布
• 分级密度图：将密度值分箱，用离散颜色强化稀疏区域辨识度
• 热点分析（Getis-Ord Gi）*：识别统计意义上的冷点区域，直接标记稀疏区

示例代码：

核密度+点叠加

ggplot(usa_map) +
  geom_sf(data = density_sf, aes(fill = v), col = NA) +
  geom_sf(data = sf_points, color = "red", size = 1.5) +
  theme_bw() +
  ggthemes::theme_map() + 
  geom_sf(data = st_boundary(usa_map)) +
  coord_sf(crs = st_crs("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"), datum = NA) +
  scale_fill_viridis_c(trans = "log1p")

分级密度图

# 密度值分箱
density_sf <- density_sf %>%
  mutate(density_bin = cut(v, breaks = c(0, 1e-6, 1e-5, 1e-4, Inf),
                           labels = c("极低密度", "低密度", "中密度", "高密度")))

ggplot(usa_map) +
  geom_sf(data = density_sf, aes(fill = density_bin), col = NA) +
  theme_bw() +
  ggthemes::theme_map() + 
  geom_sf(data = st_boundary(usa_map)) +
  coord_sf(crs = st_crs("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"), datum = NA) +
  scale_fill_brewer(palette = "YlOrRd", name = "密度等级")

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