如何基于双RasterBrick的年度栅格数据绘制小提琴图+箱线图

背景说明 #

你拥有两个RasterBrick对象，每个包含30年的月度栅格数据（共360个图层），目标是基于1960-2010年的数据，按年度维度绘制小提琴图+箱线图（匹配你提供的随机数据集可视化效果）。以下是基于随机模拟RasterBrick的完整操作流程：

步骤1：生成随机RasterBrick模拟数据 #

先创建两个对应不同变量的RasterBrick，模拟1960-1989年的月度数据（共360层），实际可根据需求扩展到1960-2010年：

library(raster)
library(tidyverse)
library(introdataviz)

set.seed(42)

# 创建基础栅格模板（10×10栅格，模拟实际空间数据）
r_template <- raster(nrow = 10, ncol = 10, vals = NA)

# 生成两个包含360层的RasterBrick
brick_ne <- brick(lapply(1:360, function(i) setValues(r_template, runif(ncell(r_template), 0.1, 1.9))))
brick_theo <- brick(lapply(1:360, function(i) setValues(r_template, runif(ncell(r_template), 0.1, 1.9))))

# 为图层命名（格式：年份_月份）
layer_dates <- seq.Date(as.Date("1960-01-01"), as.Date("1989-12-01"), by = "month")
names(brick_ne) <- str_replace_all(layer_dates, "-", "_")
names(brick_theo) <- str_replace_all(layer_dates, "-", "_")

步骤2：提取栅格数据并转换为长格式 #

由于实际数据量较大，不建议直接全量转换为data.frame，这里采用抽样处理（抽取20%栅格单元），同时整理为ggplot所需的长格式：

# 随机抽取部分栅格单元
sample_cells <- sample(ncell(r_template), size = ncell(r_template)*0.2)

# 提取抽样单元数据并整理
data_ne <- as.data.frame(brick_ne, cells = sample_cells) %>%
  mutate(cell_id = sample_cells) %>%
  pivot_longer(-cell_id, names_to = "date", values_to = "value") %>%
  mutate(measure = "chi_2_ne")

data_theo <- as.data.frame(brick_theo, cells = sample_cells) %>%
  mutate(cell_id = sample_cells) %>%
  pivot_longer(-cell_id, names_to = "date", values_to = "value") %>%
  mutate(measure = "chi_2_theo")

# 合并两个数据集
plot_data <- bind_rows(data_ne, data_theo)

步骤3：提取年份信息并筛选目标范围 #

从图层名称中提取年份，筛选1960-2010年的数据：

plot_data <- plot_data %>%
  mutate(year = str_sub(date, 1, 4) %>% as.integer()) %>%
  filter(year >= 1960 & year <= 1989) # 模拟数据到1989，实际可调整为2010

步骤4：绘制小提琴图+箱线图 #

沿用你提供的可视化逻辑，适配年度分组：

ggplot(plot_data, aes(x = factor(year), y = value, fill = measure)) +
  introdataviz::geom_split_violin(trim = FALSE, alpha = 0.4) +
  geom_boxplot(width = 0.2, alpha = 0.6) +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_fill_manual(values = c("#E69F00", "#56B4E9"), name = "变量类型") +
  labs(x = "年份", y = "数值") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # 旋转x轴标签避免重叠

大数据优化建议 #

如果实际数据量极大，可进一步优化：

• 调整抽样比例或分块处理数据，避免内存溢出
• 使用tidyterra包替代基础raster函数，提升栅格数据的tidy转换效率
• 若无需完整分布细节，可先按年度计算栅格统计量（如均值、分位数）再可视化

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