使用ggplot可视化时间序列风速遇stat_contour_fill报错求助

想要复刻Stack Overflow上《How to create contour with wind animation using gganimate?》的风场等高线动画效果，但使用时间序列风速数据集时遇到以下错误：

Warning message:
Computation failed in stat_contour_fill():
Supplied 2 columns to be assigned an empty list (which may be an empty data.table or data.frame since they are lists too). To delete multiple columns use NULL instead. To add multiple empty list columns, use list(list()).

原代码 #

pacman::p_load(ggplot2, tidyr, dplyr, raster,rgdal,sf,tidyverse,sp, gganimate, maptools, plyr, transformr, readxl, ggmap)


scale.breaks = scales::fullseq(range(test$speed), size = 10)

tim <- get_stamenmap(bbox = c(left = -105.2545,
                              bottom = 35.58840,
                              right = -105.6946,
                              top = 36.19384),
                     maptype = "terrain", 
                     crop = FALSE,
                     zoom = 11)


map.annotation <- list(
  annotation_raster(tim %>% unlist() %>%
                      alpha(0.4) %>% # change alpha setting for map here
                      matrix(nrow = dim(tim)[1], 
                             byrow = TRUE),
                    xmin = attr(tim, "bb")$ll.lon,
                    xmax = attr(tim, "bb")$ur.lon,
                    ymin = attr(tim, "bb")$ll.lat,
                    ymax = attr(tim, "bb")$ur.lat),
  coord_quickmap(xlim = c(attr(tim, "bb")$ll.lon, attr(tim, "bb")$ur.lon),
                 ylim = c(attr(tim, "bb")$ll.lat, attr(tim, "bb")$ur.lat),
                 expand = FALSE))

p.base <- ggplot(test, aes(x = longitude, y = latitude, z = speed))

p.base + 
  geom_contour_fill(breaks = scale.breaks) +
  facet_wrap(~DateTime) +
  map.annotation +
  scale_fill_gradient(low = "green", high = "red",
                      aesthetics = c("colour", "fill"),
                      limits = range(scale.breaks)) +
  theme_minimal()

数据集 #

tibble::tribble(
  ~speed, ~dir, ~longitude,     ~latitude,     ~DateTime,
  8.6,    269,  -105.265514515, 35.4631886028, as.POSIXct("2022-04-01 20:11:00"),
  8.3,    271,  -105.25959775,  35.4632014724, as.POSIXct("2022-04-01 20:11:00"),
  8.5,    273,  -105.253680981, 35.463214052,  as.POSIXct("2022-04-01 20:11:00"),
  8.9,    272,  -105.247764209, 35.4632263417, as.POSIXct("2022-04-01 20:11:00"),
  9.3,    266,  -105.241847433, 35.4632383413, as.POSIXct("2022-04-01 20:11:00"),
  9.1,    259,  -105.235930653, 35.4632500508, as.POSIXct("2022-04-01 20:11:00"),
  7.9,    260,  -105.230013871, 35.4632614704, as.POSIXct("2022-04-01 20:11:00"),
  7.7,    271,  -105.224097085, 35.4632725999, as.POSIXct("2022-04-01 20:11:00"),
  8.4,    278,  -105.218180295, 35.4632834394, as.POSIXct("2022-04-01 20:11:00"),
  9,      276,  -105.212263503, 35.4632939889, as.POSIXct("2022-04-01 20:11:00"),
)

错误原因 #

1. 等高线breaks设置无效：你的scale.breaks用scales::fullseq(range(test$speed), size = 10)生成，但数据集的风速范围仅为7.7-9.3，size=10会导致fullseq只返回一个值（7.7）。geom_contour_fill至少需要2个breaks才能划分出填充区间，这直接导致计算失败，触发错误提示。
2. 数据与地图范围不匹配：数据集的纬度集中在35.46左右，但你设置的地图bbox的bottom=35.58840，数据点完全在地图可视范围之外，即使修复breaks也看不到任何内容。
3. 底图处理方式复杂且易出错：手动将ggmap对象转矩阵再用annotation_raster的方式容易出问题，且plyr与dplyr同时加载会引发函数冲突。

修复后的代码 #

# 加载包，移除plyr避免冲突
pacman::p_load(ggplot2, tidyr, dplyr, raster, rgdal, sf, tidyverse, sp, gganimate, maptools, transformr, readxl, ggmap)

# 构造测试数据
test <- tibble::tribble(
  ~speed, ~dir, ~longitude,     ~latitude,     ~DateTime,
  8.6,    269,  -105.265514515, 35.4631886028, as.POSIXct("2022-04-01 20:11:00"),
  8.3,    271,  -105.25959775,  35.4632014724, as.POSIXct("2022-04-01 20:11:00"),
  8.5,    273,  -105.253680981, 35.463214052,  as.POSIXct("2022-04-01 20:11:00"),
  8.9,    272,  -105.247764209, 35.4632263417, as.POSIXct("2022-04-01 20:11:00"),
  9.3,    266,  -105.241847433, 35.4632383413, as.POSIXct("2022-04-01 20:11:00"),
  9.1,    259,  -105.235930653, 35.4632500508, as.POSIXct("2022-04-01 20:11:00"),
  7.9,    260,  -105.230013871, 35.4632614704, as.POSIXct("2022-04-01 20:11:00"),
  7.7,    271,  -105.224097085, 35.4632725999, as.POSIXct("2022-04-01 20:11:00"),
  8.4,    278,  -105.218180295, 35.4632834394, as.POSIXct("2022-04-01 20:11:00"),
  9,      276,  -105.212263503, 35.4632939889, as.POSIXct("2022-04-01 20:11:00"),
)

# 生成有效的等高线breaks：按0.5间隔划分风速范围
scale.breaks <- seq(min(test$speed), max(test$speed), by = 0.5)

# 调整地图bbox，包含数据点的经纬度范围
tim <- get_stamenmap(bbox = c(left = min(test$longitude) - 0.01,
                              bottom = min(test$latitude) - 0.01,
                              right = max(test$longitude) + 0.01,
                              top = max(test$latitude) + 0.01),
                     maptype = "terrain", 
                     crop = FALSE,
                     zoom = 14)

# 用ggmap快速添加底图，替代复杂的annotation_raster处理
p.base <- ggplot(test, aes(x = longitude, y = latitude, z = speed))

ggmap(tim) +
  geom_contour_fill(breaks = scale.breaks) +
  facet_wrap(~DateTime) +
  scale_fill_gradient(low = "green", high = "red",
                      aesthetics = c("fill"), # 只设置fill美学，colour对应等高线边界，这里不需要
                      limits = range(scale.breaks)) +
  theme_minimal()

动画效果扩展 #

如果要实现风场等高线的动画效果，只需在上述代码基础上添加gganimate的动画层：

ggmap(tim) +
  geom_contour_fill(breaks = scale.breaks) +
  scale_fill_gradient(low = "green", high = "red",
                      aesthetics = c("fill"),
                      limits = range(scale.breaks)) +
  theme_minimal() +
  transition_time(DateTime) + # 按时间序列过渡
  labs(title = "Wind Speed Contour: {frame_time}")

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