基于R语言的单波段多栅格季节性探索性数据分析（EDA）方法问询

一、前期准备：栅格读取与数据整理 #

1. 加载核心工具包 #

优先选高效的terra处理栅格，配合数据处理和可视化包：

# 首次运行安装依赖包
install.packages(c("terra", "dplyr", "lubridate", "ggplot2", "tidyr"))

# 加载包
library(terra)
library(dplyr)
library(lubridate)
library(ggplot2)
library(tidyr)

2. 批量读取栅格并提取时间信息 #

假设栅格文件名包含可解析的日期（如20090105.tif对应2009年第1周），批量读取后拆分年份、周数：

# 获取所有栅格文件路径
raster_paths <- list.files(path = "你的栅格文件夹路径", pattern = "\\.tif$", full.names = TRUE)

# 从文件名提取日期并拆解年份、周数（根据实际命名规则调整）
file_dates <- basename(raster_paths) %>% 
  gsub("\\.tif$", "", .) %>% 
  ymd()

year_info <- year(file_dates)
week_info <- week(file_dates)

3. 按年份分组栅格 #

用terra的栅格集合工具按年份拆分数据：

# 创建栅格集合
raster_stack <- sprc(raster_paths)

# 按年份分组存储栅格索引
year_groups <- split(1:length(raster_stack), year_info)

二、季节性EDA核心操作 #

1. 计算年度周度均值 #

针对每个年份，计算每周栅格的全局均值（如需分区均值可调整代码）：

# 定义函数：计算单年份的周度均值
calc_year_weekly_mean <- function(year_idx) {
  year_rasters <- raster_stack[year_idx]
  weekly_means <- sapply(year_rasters, function(r) global(r, mean, na.rm = TRUE)[[1]])
  data.frame(
    Year = unique(year_info[year_idx]),
    Week = week_info[year_idx],
    Mean_Value = weekly_means
  )
}

# 批量计算所有年份的周度均值
year_weekly_data <- lapply(year_groups, calc_year_weekly_mean) %>% 
  bind_rows()

2. 季节性可视化 #

（1）年度周度趋势对比折线图

直观展示各年份的季节性波动规律：

ggplot(year_weekly_data, aes(x = Week, y = Mean_Value, color = factor(Year))) +
  geom_line(alpha = 0.7) +
  labs(x = "周数", y = "栅格均值", color = "年份") +
  theme_minimal() +
  scale_x_continuous(breaks = seq(0, 52, 4))

（2）周度数值分布箱线图

识别不同周数的数值分布特征，定位稳定季节性区间：

ggplot(year_weekly_data, aes(x = factor(Week), y = Mean_Value)) +
  geom_boxplot(fill = "lightblue", alpha = 0.6) +
  labs(x = "周数", y = "栅格均值") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

（3）年度季节性偏差热力图

对比各年份与整体均值的偏差，快速定位异常年份：

# 计算所有年份每周的整体均值
overall_weekly_mean <- year_weekly_data %>% 
  group_by(Week) %>% 
  summarise(Overall_Mean = mean(Mean_Value, na.rm = TRUE))

# 计算偏差并绘图
year_deviation <- year_weekly_data %>% 
  left_join(overall_weekly_mean, by = "Week") %>% 
  mutate(Deviation = Mean_Value - Overall_Mean)

ggplot(year_deviation, aes(x = Week, y = factor(Year), fill = Deviation)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
  labs(x = "周数", y = "年份", fill = "与整体均值的偏差") +
  theme_minimal()

三、进阶分析（可选） #

• 空间季节性分析：按周数计算多年空间均值，对比季节性空间分布：

# 按周数分组计算像元均值
seasonal_stack <- tapp(raster_stack, week_info, mean, na.rm = TRUE)
# 绘制空间对比图
plot(seasonal_stack)

• 时间序列季节性分解：拆分单年份的趋势、季节性、残差成分：

# 以2023年为例
year_2023 <- year_weekly_data %>% filter(Year == 2023)
ts_data <- ts(year_2023$Mean_Value, frequency = 52)
stl_result <- stl(ts_data, s.window = "periodic")
plot(stl_result)

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