如何修复ggplot2 stat_density_2d在数据集边界框边缘的多边形异常？

问题描述 #

我正在使用ggplot2::stat_density_2d分析sfc_POINT类型的数据集，但生成的多边形在数据集边界框边缘出现变形——ggplot2似乎无法绘制超出源数据点构成的边界框之外的多边形。请问如何修复这些密度多边形？是否有办法扩展stat_density_2d的边界框？

可复现代码：

library(dplyr)
library(sf)
library(geofi)
library(ggplot2)

# 芬兰市政数据
muns <- geofi::get_municipalities(year = 2022)

# 创建样本点
points <- sf::st_sample(muns, 50) %>% as.data.frame()
points[c("x", "y")] <- sf::st_coordinates(points$geometry)

p4 <- ggplot() +
  geom_sf(data = muns) +
  coord_sf(default_crs = sf::st_crs(3067)) +
  stat_density_2d(geom = "polygon", 
                  data = points, 
                  aes(x = x, y = y, fill = ..level..), 
                  alpha = 0.75) +
  scale_fill_viridis_c() +
  geom_point(data = points, aes(x = x, y = y), alpha = 0.5)
p4

问题表现：密度多边形在数据点最外侧构成的边界框边缘会出现扭曲变形。

解决方案 #

方法1：添加虚拟点扩展数据范围 #

在现有数据的边界外添加几个不影响密度计算的虚拟点，强制stat_density_2d将计算范围扩展到目标区域（比如芬兰地图的边界）。

library(dplyr)
library(sf)
library(geofi)
library(ggplot2)

# 加载芬兰市政数据
muns <- geofi::get_municipalities(year = 2022)
# 获取地图边界框
mun_bbox <- st_bbox(muns)

# 创建样本点
points <- sf::st_sample(muns, 50) %>% as.data.frame()
points[c("x", "y")] <- sf::st_coordinates(points$geometry)

# 在地图四个角外侧添加虚拟点（偏移量根据坐标系单位调整，这里是米）
dummy_points <- tibble(
  x = c(mun_bbox[["xmin"]] - 10000, mun_bbox[["xmax"]] + 10000, 
        mun_bbox[["xmin"]] - 10000, mun_bbox[["xmax"]] + 10000),
  y = c(mun_bbox[["ymin"]] - 10000, mun_bbox[["ymin"]] - 10000, 
        mun_bbox[["ymax"]] + 10000, mun_bbox[["ymax"]] + 10000)
)

# 合并原始点与虚拟点
extended_points <- bind_rows(points, dummy_points)

# 绘图
p4 <- ggplot() +
  geom_sf(data = muns) +
  coord_sf(default_crs = sf::st_crs(3067)) +
  stat_density_2d(geom = "polygon", 
                  data = extended_points, 
                  aes(x = x, y = y, fill = ..level..), 
                  alpha = 0.75) +
  scale_fill_viridis_c() +
  geom_point(data = points, aes(x = x, y = y), alpha = 0.5)
p4

方法2：手动调用密度函数控制计算范围 #

使用MASS::kde2d直接计算密度，指定计算范围为地图边界框，再转换成多边形数据，完全避免边界限制问题。

library(dplyr)
library(sf)
library(geofi)
library(ggplot2)
library(MASS)
library(tibble)

# 加载数据
muns <- geofi::get_municipalities(year = 2022)
mun_bbox <- st_bbox(muns)

points <- sf::st_sample(muns, 50) %>% as.data.frame()
points[c("x", "y")] <- sf::st_coordinates(points$geometry)

# 手动计算2D密度，指定范围为地图边界框外侧
dens <- kde2d(points$x, points$y, 
              lims = c(mun_bbox[["xmin"]]-10000, mun_bbox[["xmax"]]+10000,
                       mun_bbox[["ymin"]]-10000, mun_bbox[["ymax"]]+10000),
              n = 200) # n值越大，密度图分辨率越高

# 转换为ggplot可用的数据格式
dens_df <- expand.grid(x = dens$x, y = dens$y) %>%
  mutate(z = as.vector(dens$z))

# 绘图
p4 <- ggplot() +
  geom_sf(data = muns) +
  coord_sf(default_crs = sf::st_crs(3067)) +
  geom_contour_filled(data = dens_df, aes(x = x, y = y, z = z), alpha = 0.75) +
  scale_fill_viridis_c() +
  geom_point(data = points, aes(x = x, y = y), alpha = 0.5)
p4

注意：单纯扩展绘图范围无效 #

如果仅通过coord_sf的xlim/ylim或expand参数扩大显示范围，只会让画布变大，但密度多边形的计算范围还是受限于原始数据点的边界，依然会出现变形，因此不推荐这种方式。

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