如何用precintcon R包批量计算并保存多列年SPI值

问题分析 #

precintcon包的spi()函数仅支持单序列处理，直接用summarize(across(...))无法完成多站点的SPI计算——因为SPI需要基于时间序列的年降水聚合和分布拟合，不是简单的列运算。必须批量处理每个站点的降水序列，再将结果按年份合并。

解决方案代码 #

1. 加载依赖包 #

library(precintcon)
library(dplyr)
library(tidyr)
library(lubridate)

2. 自定义年SPI计算函数 #

这个函数负责将单站点的月降水数据转换为precintcon要求的格式，聚合年降水，再计算年SPI：

calc_annual_spi <- function(precip_vec, date_vec) {
  # 构造precintcon所需的月尺度数据框
  monthly_df <- tibble(
    year = year(date_vec),
    month = month(date_vec),
    precipitation = precip_vec
  )
  
  # 聚合年降水量
  annual_precip <- annual.precipitation(monthly_df)
  
  # 计算年SPI（将年数据转为"月=1"的序列，scale=1对应年尺度）
  spi_result <- spi(
    data = data.frame(
      year = annual_precip$year,
      month = rep(1, nrow(annual_precip)),
      precipitation = annual_precip$precipitation
    ),
    scale = 1
  ) %>%
    select(year, spi)
  
  return(spi_result)
}

3. 批量处理所有站点并合并结果 #

假设你的原始数据框df包含Date列（日期格式）和RainGauge1至RainGauge6列：

# 获取所有雨量站列名
gauge_cols <- colnames(df)[grepl("^RainGauge", colnames(df))]

# 对每个站点计算年SPI，得到结果列表
spi_results <- lapply(gauge_cols, function(col) {
  spi_df <- calc_annual_spi(df[[col]], df$Date)
  # 重命名SPI列，加上站点标识
  colnames(spi_df)[2] <- paste0("SPI_", col)
  spi_df
})

# 按年份合并所有站点的SPI结果
final_spi_table <- spi_results %>% reduce(left_join, by = "year")

替代方案（tidyverse长格式处理） #

如果习惯用tidyverse的长格式操作，也可以这样写：

final_spi_table <- df %>%
  # 宽转长，方便按站点分组处理
  pivot_longer(
    cols = starts_with("RainGauge"),
    names_to = "Station",
    values_to = "Precipitation"
  ) %>%
  mutate(
    year = year(Date),
    month = month(Date)
  ) %>%
  # 按站点和年份聚合年降水
  group_by(Station, year) %>%
  summarize(Annual_Precip = sum(Precipitation, na.rm = TRUE), .groups = "drop") %>%
  # 构造precintcon所需格式（每个年对应month=1）
  mutate(month = 1) %>%
  # 按站点分组计算SPI
  group_by(Station) %>%
  group_modify(function(data, station) {
    spi(data[, c("year", "month", "Annual_Precip")], scale = 1) %>%
      select(year, spi)
  }) %>%
  # 长转宽，得到最终表格
  pivot_wider(
    names_from = Station,
    values_from = spi,
    names_prefix = "SPI_"
  )

最终输出格式 #

运行后final_spi_table会生成如下结构的表格：

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