R语言中基于SPI数据的干旱事件统计函数修正求助

我仔细看了你的需求和原函数的问题，咱们先梳理下核心问题出在哪，再给出修正后的方案：

原函数的问题分析 #

你的原函数主要有两个关键问题：

1. 时长计算逻辑错误：你把rle返回的段索引当成了实际的月份索引，用段索引相减得到的结果根本不是实际的干旱月数，这就是为什么你得到的总时长是1而不是7。
2. 干旱范围定义偏差：你只统计了连续SPI≤-1的月份，但根据你的规则，一旦干旱触发（连续2个月≤-1），之后直到SPI转正前的所有月份（包括SPI在(-1,0)之间的月份）都应该算入干旱时长，原函数漏掉了这部分。

修正后的函数 #

我重新梳理了逻辑，按照你的规则（连续≥2个月SPI≤-1触发干旱，SPI转正时结束）写了修正后的函数：

d.frequency <- function(x, na.rm = TRUE) {
  # 先处理缺失值（如果需要）
  if(na.rm) x <- x[!is.na(x)]
  
  # 标记每个月份的状态：是否为干旱候选（SPI≤-1）、是否为正值（SPI≥0）
  is_drought_candidate <- x <= -1
  is_positive <- x >= 0
  
  # 用rle找出连续的干旱候选段
  rle_candidate <- rle(is_drought_candidate)
  
  # 计算每个候选段在原数据中的起始和结束索引
  candidate_starts <- c(1, cumsum(rle_candidate$lengths)[-length(rle_candidate$lengths)] + 1)
  candidate_ends <- cumsum(rle_candidate$lengths)
  
  # 筛选出长度≥2的有效干旱候选段（满足触发条件）
  valid_candidates <- which(rle_candidate$values & rle_candidate$lengths >= 2)
  
  # 如果没有符合条件的干旱事件，直接返回结果
  if(length(valid_candidates) == 0) {
    return(list(num_droughts = 0, total_duration = 0, average_duration = 0))
  }
  
  drought_durations <- c()
  used_indices <- c() # 记录已经被统计过的月份，避免重复计算
  
  # 遍历每个有效候选段，确定实际的干旱时长
  for(i in valid_candidates) {
    start_idx <- candidate_starts[i]
    # 如果这个起始点已经被之前的干旱事件包含，跳过
    if(start_idx %in% used_indices) next
    
    # 找到从start_idx开始，第一个出现正值的位置
    first_positive <- which(is_positive[start_idx:length(is_positive)])[1]
    
    if(is.na(first_positive)) {
      # 如果后续没有正值，干旱持续到数据最后一个月
      end_idx <- length(x)
    } else {
      # 干旱结束在第一个正值的前一个月
      end_idx <- start_idx + first_positive - 2
    }
    
    # 计算当前干旱的时长
    duration <- end_idx - start_idx + 1
    drought_durations <- c(drought_durations, duration)
    
    # 标记这段索引为已使用
    used_indices <- c(used_indices, start_idx:end_idx)
  }
  
  # 计算最终统计值
  num_droughts <- length(drought_durations)
  total_duration <- sum(drought_durations)
  average_duration <- if(num_droughts > 0) total_duration / num_droughts else 0
  
  # 返回结果列表
  return(list(
    num_droughts = num_droughts,
    total_duration = total_duration,
    average_duration = average_duration
  ))
}

测试验证 #

用你提供的示例数据测试：

data <- c(-1,-1,-1,-0.5,0,-0.5,-1,0,-1,-1,-0.8,0,0,0)
d.frequency(data)

返回结果完全符合你的预期：

$num_droughts
[1] 2

$total_duration
[1] 7

$average_duration
[1] 3.5

逻辑说明 #

这个修正后的函数核心逻辑是：

1. 先标记每个月的状态，区分干旱候选和正值月份
2. 找到所有满足触发条件（连续≥2个月SPI≤-1）的起始点
3. 对每个起始点，追踪到第一个正值月份，确定干旱的结束点（正值的前一个月）
4. 避免重复统计重叠的干旱事件
5. 最终计算事件数量、总时长和平均时长

备注：内容来源于stack exchange，提问作者Rahmat Hidayat