在R中使用gratia包draw函数绘制经纬度图时添加geom_sf报错求助

问题场景 #

用gratia包的draw()函数绘制含s(longitude, latitude)平滑项的多分类gam模型时，能生成经纬度效应的等高线图，但叠加通过giscoR获取的梵蒂冈矢量边界时出现报错：

Coordinate system already present. Adding a new coordinate system, which will replace the existing one.
Error in geom_sf():
! Problem while computing aesthetics.
ℹ Error occurred in the 5th layer.
Caused by error in .data[["longitude"]]:
! Column longitude not found in .data.

测试用模型代码如下：

sound <- sample(0:5, size=960, replace=T)
word <- as.factor(rep(c('alpha', 'bravo', 'charlie', 'delta', 'echo', 'foxtrot'), each=4))
age <- as.factor(rep(c('young', 'old'), times=480))
gender <- as.factor(rep(c('female', 'female', 'male', 'male'), times=240))
longitude <- rep(c(runif(40, min=41, max=42)), each=24)
latitude <- rep(c(runif(40, min=12, max=13)), each=24)
pronunciation <- data.frame(sound, word, age, gender, longitude, latitude)

library(mgcv)

model = gam(list(sound ~ word + s(longitude, latitude),
                       ~ word + s(longitude, latitude),
                       ~ word + s(longitude, latitude),
                       ~ word + s(longitude, latitude),
                       ~ word + s(longitude, latitude)),
                       data=pronunciation,
                       family=multinom(K=5))

报错原因 #

1. draw()生成的ggplot对象，其图层数据中没有longitude列，但geom_sf默认会继承主图的美学映射，导致它尝试从主数据中读取不存在的列，触发报错。
2. 主图的坐标系与sf边界数据的坐标系不一致，引发坐标系替换的警告。

解决方法 #

步骤1：获取并处理梵蒂冈边界数据 #

先通过giscoR获取边界，并转换到与模型经纬度匹配的WGS84坐标系（EPSG:4326）：

library(giscoR)
library(sf)

# 获取梵蒂冈边界并指定坐标系
vatican <- gisco_get_countries(country = "VAT", epsg = 4326)

步骤2：叠加边界时关闭美学继承并匹配坐标系 #

在调用geom_sf时，设置inherit.aes = FALSE避免继承主图的无效美学映射，同时明确指定坐标系消除冲突：

library(gratia)
library(ggplot2)

# 绘制模型平滑项图并叠加边界
draw(model, select = "s(longitude,latitude)") +
  geom_sf(data = vatican, inherit.aes = FALSE, color = "red", fill = NA) +
  coord_sf(crs = 4326)

补充适配方案 #

如果仍存在坐标系不匹配的问题，可以先提取draw生成图像的坐标系，再将边界数据转换为对应坐标系：

# 先生成基础图并提取其坐标系
p <- draw(model, select = "s(longitude,latitude)")
target_crs <- ggplot2::layer_scales(p)$x$crs

# 转换边界数据的坐标系
vatican_transformed <- st_transform(vatican, crs = target_crs)

# 叠加转换后的边界
p + geom_sf(data = vatican_transformed, inherit.aes = FALSE, color = "red", fill = NA)

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