计算从日落开始的天数。

Question

我正在分析一个复杂数据集中的时间模式，其中包括多个环境变量以及来自各种动物物种的活动数据。这些数据由多个实验设置收集，并且每个设置的数据每分钟存储一次。该项目已经运行了几年，因此我的数据集相当大。

我的一个数据集的前几行如下所示：

> head(setup_01)
DateTime                Film_number unused PIR Wheel Temperature LightOld LightDay LightNight LightUV IDnumbers    error mouse shrew vole rat frog rest extra_info odour
1 2015-03-10 12:27:10                  x   0       0       13.40  1471.34    -0.97    1331.29  700.42           no error     0     0    0   0    0    0                1
2 2015-03-10 12:28:10                  x   0       0       13.43  1471.38    -1.07    1291.11  731.32           no error     0     0    0   0    0    0                1
3 2015-03-10 12:29:10                  x   0       0       13.31  1471.24    -1.08    1368.57 1016.02           no error     0     0    0   0    0    0                1

我想将这些变量与不同的自然循环如季节内的日出和日落联系起来，因此我使用了maptools包来计算日出和日落时间。

library(maptools)
gpclibPermit()

#set coordinates
crds=c(4.4900,52.1610)

# download the sunrise/sunset/etc data
setup_01$sunrise=sunriset(matrix(crds,nrow=1),dateTime=as.POSIXct(setup_01$DateTime),POSIXct.out=TRUE,direction="sunrise")
setup_01$sunset=sunriset(matrix(crds,nrow=1),dateTime=as.POSIXct(setup_01$DateTime),POSIXct.out=TRUE,direction="sunset")

#create a variable that's 0 except at sunrise, and one that's 0 except at sunset
setup_01$sunrise_act=0
setup_01$sunset_act=0
setup_01[abs(unclass(setup_01[,"DateTime"])-unclass(setup_01[,"sunrise"]$time))<30,]$sunrise_act=1
setup_01[abs(unclass(setup_01[,"DateTime"])-unclass(setup_01[,"sunset"]$time))<30,]$sunset_act=1

由于大多数动物的行为因白天和黑夜而异，我使用日落/日出时间来计算一个新变量，该变量在夜间为0，在白天为1：

#create a variable that's 0 at night and 1 at daytime
setup_01$daytime=0
setup_01[setup_01[,"DateTime"]>setup_01[,"sunrise"]$time & setup_01[,"DateTime"]<setup_01[,"sunset"]$time,]$daytime=1

到目前为止一切都很好......使用maptools甚至可以使用民用/航海/天文黄昏和黎明的开始来代替日出和日落。

然而，这就是我的问题所在。我想为我的实验中所有的日子编号。与通常在午夜增加日期计数器不同的是，我希望在日落时（或可能在未来的实验中，像日出、航海黄昏和黎明等其他可移动的时间）增加日期计数器。由于日落不会在每天的同一时间发生，对我来说这不是一个简单的问题。

我只想到了一个for循环的方式，这不是一个好的方法。另外，考虑到我已经收集了超过6年的数据点，在多个设置中每分钟收集一次数据，当R运行类似这样的循环时，我可以去看板块移动。

setup_01$day=0
day<-1
for(i in 1:nrow(setup_01)){
    setup_01[i,]$day<-day
    if(setup_01[i,]$sunset_act==1){
        day<-day+1
    }
}

除了丑陋和缓慢，这段代码还有一个大问题：它不能处理缺失值。有时由于设备故障，数据在几个小时或几天内根本没有记录。如果在日落期间没有记录任何数据，则上述代码不会增加天数计数器。这意味着我需要 - 在某种程度上 - 同时纳入日期/时间代码。很容易创建一个从实验开始以来的天数变量：

setup_01$daynumber<-as.integer(ceiling(difftime(setup_01$DateTime, setup_01$DateTime[1], units = "days")))

也许这些数字可以被使用，可能与Heroka's的好的rle算法相结合。

我已经使用dput从一个设置中制作了几个月的数据，包括一些大块丢失的数据，以及新创建的变量（如本帖子和Heroka's的答案所述），可在here找到。

我一直在寻找更好、更漂亮、特别是更快的东西，但一直无法想出一个好的技巧。我试图调整我的数据框的子集，但得出结论这可能是一个愚蠢的方法。我看过maptools、lubridate和GeoLight。我搜索了Google、Stack Overflow和各种书籍，比如Hadley Wickham的精彩Advanced R。但都没有找到。也许我错过了一些显而易见的东西。我希望这里有人能帮助我。

- Yuri Robbers

Answer 1

我提出了一种解决方案，可以在已经生成了0和1的情况下，并且可以与运行长度（runlengths）一起使用。

  #sunset/sunrise is series of 0's and 1's indicating night and daytime, so solution that works for random sequence
#will work for OP's dataset
set.seed(10)
sunset <- c(1,rbinom(20,1,0.5))

#counter needs to be x for sequence of 11111 (day) and 0000(night), and then increase when 0 reappears
#counter starts at 1

#intermediate step: number each half-day
rle_sunset <- rle(sunset)
period <- rep(1:length(rle_sunset$lengths),rle_sunset$lengths)
#calculate day so that each two subsequent periods are one day

day <- ceiling(period/2)

> cbind(sunset,period,day)
      sunset period day
 [1,]      1      1   1
 [2,]      1      1   1
 [3,]      0      2   1
 [4,]      0      2   1
 [5,]      1      3   2
 [6,]      0      4   2
 [7,]      0      4   2
 [8,]      0      4   2
 [9,]      0      4   2
[10,]      1      5   3
[11,]      0      6   3
[12,]      1      7   4
[13,]      1      7   4
[14,]      0      8   4
[15,]      1      9   5
[16,]      0     10   5
[17,]      0     10   5
[18,]      0     10   5
[19,]      0     10   5
[20,]      0     10   5
[21,]      1     11   6

- Heroka

---

谢谢你，Heroka！这正是我所需要的，并且速度几乎是即时的...再次感谢！ - Yuri Robbers

还有一个小问题，可能我太快接受这个答案为完美了：如果没有缺失数据，这个方法绝对没问题。但是一旦缺失的数据中至少有一个日落时间为“1”的实例，就会出现问题……不幸的是，我的数据中有这样的实例。我在原帖中没有提到这一点，抱歉。 - Yuri Robbers

你能否编辑你的帖子以适应缺失的数据？它们是如何生成的？（你是否有缺失的日期？）也许你可以使用dput提供更多（相关的！）数据。 - Heroka

谢谢，Heroka。我已相应地编辑了我的帖子，并提供了一个相关数据子集的链接，使用dput提供。 - Yuri Robbers

Answer 2

我更喜欢基于预先计算的表格的解决方案。虽然速度较慢，但我认为更容易理解。然后我使用 dplyr 来整理我需要的信息。

让我举个例子来说明。为了举例，我创建了一个日落时间列表。当然，您需要计算实际的日落时间。

library(dplyr)
n.obs=1000
set.seed(10)
t0 <- as.POSIXct('2015-03-08 18:00:00')
artificial.sunsets <- data.frame(num.day= seq(0,n.obs+35)) %>% mutate(sunset=cumsum(rlnorm(length(num.day))*30)+t0 + 24*3600*num.day)

artificial.sunsets 包含日落的日期和时间，但也可能包括有关该天的更多信息。

以下是一些人工数据：

t0 <- as.POSIXct('2015-03-10 12:27:10')
test.data <- data.frame(DateTime=t0+ seq(0, n.obs*24*3600, by=3600), observation=rnorm(24*n.obs+1))

然后可以使用以下代码找到之前的日落时间：

find.sunset.before <- function(x){
  cbind(x,artificial.sunsets %>% filter(sunset < x$DateTime) %>% tail(.,n=1))
}

data.with.sunset=test.data %>% rowwise() %>% do(find.sunset.before(.)) %>% ungroup()%>% mutate(rel.time = DateTime-sunset)
head(data.with.sunset)

生成的表格将会包含三列，分别为1) 对应日期的编号 2) 对应日落时间，和3) 日落后的时间。

由于日期编号在另一个表格中，所以这个方法可以很好地处理缺失的测量数据。您还可以轻松修改算法以使用不同的时间或者多种时间。

更新

使用data.table可以更快地完成所有操作:

library(data.table)
dt1 <- data.table(artificial.sunsets)
dt2 <- data.table(test.data)

dt1[,DateTime:=sunset]

setkey(dt1, DateTime)
setkey(dt2, DateTime)

r <- dt1[dt2,roll=TRUE]
r[,time.diff:=DateTime-sunset]

我尝试使用system.time对1000个观测值进行计时 - 前者需要约1分钟，而data.table的解决方案只需要0.011秒。

- bdecaf

---

谢谢，bdecaf。这是一个优雅的、易于理解的解决方案，而且它有效。唯一的缺点是速度：在我发布的原始帖子中可用的测试数据集上，它仍然需要半个多小时的时间，这意味着它可能需要运行大约8-10天才能处理完整个数据集。 - Yuri Robbers

我看到了 - 真的不应该花这么长时间。在这个解决方案中, find.sunset.before 是最慢的部分. 我找到了这个有趣的问题，可能与之相关: https://dev59.com/rGXWa4cB1Zd3GeqPP8Ha 并将其添加到答案中。在我的笔记本电脑上，这个解决方案近乎瞬间完成（我的笔记本非常慢） - bdecaf

谢谢，@bdecaf。看起来真的很有前途。你的旧版本的计时似乎与我的结果相差不大。请注意，我每天有1440个观测值，共有6年的8个设置。这将需要近两周的时间。你的 data.table 版本应该能把这个时间缩短到几分钟。我目前周末外出（正在手机上写这篇文章），但我会在周一之前尝试 data.table 解决方案。 - Yuri Robbers

结果令人惊讶。这让我在极短的时间内得到了想要的结果。非常感谢@bdecaf!! - Yuri Robbers