使用 data.table 总结每月序列（计算特定事件）

Question

我希望这是一个可以接受的 R/data.table 问题。

我有一个包含以下内容的 3 列表：

• id 地理位置 ID（303,453 个位置）
• month 1990-2014 年超过 25 年的月份
• spei 一个在 -7 和 7 之间变化的气候指数。

我需要计算整个 1990-2014 年期间每个地点发生的干旱。干旱事件被定义为“SPEI 持续为负且 SPEI 达到 -1.0 或更低的时期。干旱在 SPEI 第一次低于零时开始，并以第一个正 SPEI 值结束。 -1.0 或更小的值”。

我知道使用 shift() 和滚动连接应该是可行的，但非常欢迎一些帮助！

# Sample table structure
dt <- data.table(
  id = rep(1:303453, each=25*12),
  month = rep(seq(as.Date("1990-01-01"), as.Date("2014-12-31"), "month"), 303453),
  spei = runif(303453*25*12, -7, 7))

# A minimal example with 1 location over 12 months
library(data.table)
library(xts)

dt <- data.table(
  id = rep("loc1", each=12),
  month = seq(as.Date("2014-01-01"), as.Date("2014-12-31"), "month"),
  spei = c(-2, -1.1, -0.5, 1.2, -1.2, 2.3, -1.7, -2.1, 0.9, 1.2, -0.9, -0.2))

spei.ts <- xts(dt$spei, order.by=dt$month, frequency="month")
plot(spei.ts, type="bars")

这显示了 1 年期间的 3 次干旱事件。这是我需要识别和计数的。

希望你们中的一些人更习惯于使用时间序列。 非常感谢，--Mel。

Answer 1

根据评论更新...

如果只需要计数，那么

# Let 'sp' = starting point of potential drought
# Let 'dv' = drought level validation
# The cumsum just gives unique ids to group by.
dt[, sp := (spei <= 0) & (shift(spei, fill = 1) > 0), by = id]
dt[, dv := min(spei) <= -1, by = .(id, cumsum(sp))]
dt[sp & dv, .N, by = id]

然而，正如 cmets 中所述，您已经去过那里，所以您已经了解了如何使用 shift。因为您也喜欢识别日期的想法。为什么不在那里也使用shift？

# Extending the previous columns...
dt[, ep := (shift(spei, type = "lead", fill = 1) > 0) & (spei <= 0), by = id]
cbind(dt[sp & dv, .(start = month), by = id],
      dt[ep & dv, .(end = month), by = id][,id := NULL])

如果您希望日期如图中红线所示，只需添加一个月，除非是最后一个月。我们也可以得到长度...

# Extending the previous columns again...
dt[, end.month := shift(month, type = "lead", fill = month[.N]), by = id]
dt[, orig.id := .I]
starts <- dt[sp & dv][, did := .I]
ends <- dt[ep & dv][, did := .I]
starts[ends, on = "did"][
  ,.(id = id, length = 1 + i.orig.id - orig.id, start = month, end = i.end.month)]

会成功

     id length      start        end
1: loc1      3 2014-01-01 2014-04-01
2: loc1      1 2014-05-01 2014-06-01
3: loc1      2 2014-07-01 2014-09-01

而且它仍然快！与n=300

> microbenchmark(max = max.full(copy(dt))[, .(nDroughts = .N), by = id],
+                thellcounts = thell.counts(copy(dt)),
+                thell .... [TRUNCATED] 
Unit: milliseconds
        expr       min        lq      mean    median        uq        max neval
         max 218.19152 220.30895 342.18605 222.75507 250.36644 1350.15847    10
 thellcounts  20.36785  22.27349  28.45167  23.39313  24.38610   78.25046    10
  thelldates  28.24378  28.64849  30.59897  30.57793  31.25352   34.51569    10
 thelldates2  36.19724  39.79588  42.34457  41.52455  42.41872   57.28073    10

与n=3000

> microbenchmark(max = max.full(copy(dt))[, .(nDroughts = .N), by = id],
+                thellcounts = thell.counts(copy(dt)),
+                thell .... [TRUNCATED] 
Unit: milliseconds
        expr       min        lq      mean    median        uq       max neval
         max 2126.1138 2148.3453 2207.7801 2205.3536 2241.2848 2340.1203    10
 thellcounts  197.7312  202.4817  234.2949  205.4828  304.1556  309.1028    10
  thelldates  261.9889  264.5597  283.9970  266.1244  267.8603  374.6406    10
 thelldates2  320.6352  331.7558  374.4110  340.2668  439.1490  441.8473    10

Answer 2

这是获得所需结果的起点。 或许专家可以建议提高速度。

编辑：通过删除 paste 将速度提高了约 8 倍。

library(data.table)
set.seed(42)
n <- 300  # 303453 will be ~1000 times slower
dt <- data.table(
    id = rep(1:n, each=25*12),
    month = rep(seq(as.Date("1990-01-01"), as.Date("2014-12-31"), "month"), n),
    spei = runif(n*25*12, -7, 7))

system.time({
  dt[, `:=`(neg = (spei < 0), neg1 = (spei <= -1))]
  dt[, runid := ifelse(neg, rleid(neg), NA)]
  res <- dt[!is.na(runid), 
            .(length = .N[any(neg1)], start = min(month), end = max(month)), 
            by = .(id, runid)][!is.na(length)]

})
#    user  system elapsed 
#   0.345   0.000   0.344 

# counts of droughts per id:
res[, .(nDroughts = .N), by = id]

# list of droughts per id: (NB: don't include 1st positive value after) 
res[, .(droughtN = seq_len(.N), start, end), by = id]