如何按年份与气候分区（CLIMDIV）统计不同时长连续高温天数？

我需要从数据框中按气候分区（CLIMDIV）和年份统计热浪（即连续高温日）的发生次数。该数据集已从完整的气候分区与温度年份列表中过滤，仅保留最高热指数超过32℃的日期。我希望统计每个气候分区每年出现2天、3天、4天连续高温的次数，即不同时长的热浪发生频次。

我尝试过使用rle函数，但无法实现按气候分区和年份分组统计连续天数，也不知道如何生成不同连续天数的统计列表。我还尝试了以下代码，输出格式符合需求，但手动核对时发现结果数值过小且不准确。

library(dplyr) 

df18 <- data.frame(
  CLIMDIV = c(101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 102,
              102, 102, 102, 102, 102, 102),
  Year = c(1970, 1970, 1970, 1970, 1970, 1970, 1971, 1971, 1971, 1971, 1971, 1970, 1970, 1970,
           1972, 1972, 1972, 1972),
  Date = as.Date(c("1970-01-01", "1970-01-02", "1970-01-03", "1970-01-06", "1970-01-07", "1970-01-09",
                   "1971-01-01", "1971-01-03", "1971-01-04", "1971-01-05", "1971-01-06",
                   "1970-01-01", "1970-01-02", "1970-01-03",
                   "1972-01-01", "1972-01-02", "1972-01-03", "1972-01-04")),
  Max_Heat_Index = c(32.2, 32.3, 32.4, 32.5, 32.6, 32.7, 32.8, 32.9, 33.0, 33.1, 33.2, 33.3,
                     33.4, 33.5, 33.6, 33.7, 33.8, 33.9)
)

# 2 day heatwaves
grouped_data2day <- df18 %>% 
  group_by(CLIMDIV, Year) %>%
  mutate(Consecutive = cumsum(c(1, diff(Date) != 1))) %>%
  group_by(CLIMDIV, Year, Consecutive) %>%
  summarize(Count = n()) %>%
  filter(Count >= 2) %>%
  summarise(days2over32C = n()) 

# 3 day heatwaves
grouped_data3day <- df18 %>% 
  group_by(CLIMDIV, Year) %>%
  mutate(Consecutive = cumsum(c(1, diff(Date) != 1))) %>%
  group_by(CLIMDIV, Year, Consecutive) %>%
  summarize(Count = n()) %>%
  filter(Count >= 3) %>%
  summarise(days3over32C = n()) 

你之前的代码问题在于：

• 若需求是统计恰好n天的热浪事件数，filter(Count >= n)会把更长的连续段也纳入统计，导致结果不符合预期；
• 重复编写代码统计不同天数，效率较低。

以下是更简洁准确的实现，同时支持统计2、3、4天的热浪频次：

library(dplyr)
library(tidyr)

df18 <- data.frame(
  CLIMDIV = c(101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 102,
              102, 102, 102, 102, 102, 102),
  Year = c(1970, 1970, 1970, 1970, 1970, 1970, 1971, 1971, 1971, 1971, 1971, 1970, 1970, 1970,
           1972, 1972, 1972, 1972),
  Date = as.Date(c("1970-01-01", "1970-01-02", "1970-01-03", "1970-01-06", "1970-01-07", "1970-01-09",
                   "1971-01-01", "1971-01-03", "1971-01-04", "1971-01-05", "1971-01-06",
                   "1970-01-01", "1970-01-02", "1970-01-03",
                   "1972-01-01", "1972-01-02", "1972-01-03", "1972-01-04")),
  Max_Heat_Index = c(32.2, 32.3, 32.4, 32.5, 32.6, 32.7, 32.8, 32.9, 33.0, 33.1, 33.2, 33.3,
                     33.4, 33.5, 33.6, 33.7, 33.8, 33.9)
)

# 统计不同时长的热浪频次
heatwave_stats <- df18 %>%
  arrange(CLIMDIV, Year, Date) %>%  # 确保日期按顺序排列，避免识别错误
  group_by(CLIMDIV, Year) %>%
  mutate(
    # 标记每个连续日期组：日期差不为1时，组号递增
    consecutive_group = cumsum(c(TRUE, diff(Date) != 1))
  ) %>%
  group_by(CLIMDIV, Year, consecutive_group) %>%
  summarise(heatwave_length = n(), .groups = "drop_last") %>%  # 计算每个连续段的长度
  count(heatwave_length, name = "frequency") %>%  # 统计各长度热浪的出现次数
  pivot_wider(
    names_from = heatwave_length,
    names_prefix = "heatwave_",
    values_from = frequency,
    values_fill = 0  # 缺失的长度频次填0
  ) %>%
  # 确保2、3、4天的列存在，无数据则填0
  mutate(
    heatwave_2 = coalesce(heatwave_2, 0),
    heatwave_3 = coalesce(heatwave_3, 0),
    heatwave_4 = coalesce(heatwave_4, 0)
  ) %>%
  select(CLIMDIV, Year, heatwave_2, heatwave_3, heatwave_4)  # 调整列顺序

print(heatwave_stats)

代码说明 #

1. 排序：先按CLIMDIV、Year、Date排序，保证日期连续性识别准确；
2. 标记连续组：用cumsum(diff(Date) != 1)生成连续日期的分组ID；
3. 计算段长度：统计每个连续组的天数，即热浪持续时长；
4. 统计频次：按分区和年份统计不同时长热浪的发生次数；
5. 宽表转换：将结果转为宽表格式，方便查看2、3、4天的频次，缺失值补0。

结果示例 #

运行后输出如下（符合手动统计的准确结果）：

# A tibble: 4 × 5
# Groups:   CLIMDIV, Year [4]
  CLIMDIV  Year heatwave_2 heatwave_3 heatwave_4
    <dbl> <dbl>      <int>      <int>      <int>
1     101  1970          1          1          0
2     101  1971          0          0          1
3     102  1970          0          1          0
4     102  1972          0          0          1

用rle实现的版本 #

如果你偏好使用rle函数，也可以用以下代码实现相同效果：

library(dplyr)
library(tidyr)

heatwave_stats_rle <- df18 %>%
  arrange(CLIMDIV, Year, Date) %>%
  group_by(CLIMDIV, Year) %>%
  summarise(
    # 用rle提取连续日期段的长度
    heatwave_lengths = rle(c(TRUE, diff(Date) == 1))$lengths[rle(c(TRUE, diff(Date) == 1))$values]
  ) %>%
  unnest(heatwave_lengths) %>%
  count(heatwave_lengths, name = "frequency") %>%
  pivot_wider(
    names_from = heatwave_lengths,
    names_prefix = "heatwave_",
    values_from = frequency,
    values_fill = 0
  ) %>%
  mutate(
    heatwave_2 = coalesce(heatwave_2, 0),
    heatwave_3 = coalesce(heatwave_3, 0),
    heatwave_4 = coalesce(heatwave_4, 0)
  ) %>%
  select(CLIMDIV, Year, heatwave_2, heatwave_3, heatwave_4)

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