R：计算从日落开始的天数
R: count days that start at sunset
我正在分析复杂数据集中的时间模式，该数据集包含多个环境变量以及 activity 来自不同动物物种的数据。这些数据由多个实验设置收集，每个设置的数据每分钟存储一次。该项目已经 运行 好几年了，所以我的数据集相当大。
我的一个数据集的前几行如下所示：
> 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
由于大多数动物的行为不同，取决于是白天还是晚上，我使用了 sunset/sunrise 次数据来计算一个新变量，该变量在夜间为 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 甚至可以使用 civil/nautical/astronomical dusk 的开始和黎明而不是日出和日落。
然而，这就是我的问题开始的地方。我想对实验中的所有日子进行编号。而不是像往常一样容易地在 午夜 增加日计数器，我想在 sunset 增加日计数器（或者可能在未来的实验中，一天中的另一个可移动时间，如日出、航海 dusk 和黎明，...）。由于日落不会每天都在同一时间发生，所以对我来说，这不是一个可以直接解决的问题。
我只提出了一个 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
}
}
除了丑陋和缓慢之外，这段代码还有一个大问题：它不处理缺失值。有时，由于设备故障，数小时或数天根本没有记录数据。如果在日落期间没有记录数据，则上述代码不会增加日计数器。这意味着我需要 - 以某种方式 - 合并 date/time 代码。很容易创建自实验开始以来天数的变量：
setup_01$daynumber<-as.integer(ceiling(difftime(setup_01$DateTime, setup_01$DateTime[1], units = "days")))
也许可以使用这些数字，可能与 Heroka's 不错的 rle-算法结合使用。
我已经使用 dput 从一个设置中获取了几个月的数据，包括一些大块丢失的数据，以及新创建的变量（如本 [=62= 中所述） ] 并在 Heroka's 答案中）可用 此处 .
我一直在寻找更好、更好、尤其是更快的东西，但一直无法想出一个好的技巧。我一直在摆弄我的数据框的子集，但得出的结论是这可能是一种愚蠢的方法。我查看了 maptools、lubridate 和 GeoLight。我搜索了 Google、Stack Overflow 和各种书籍，例如 Hadley Wickham 的精彩 Advanced R。都无济于事。也许我错过了一些明显的东西。我希望这里有人可以帮助我。
我想出了一个关于生成的 0 和 1 的解决方案（因为你已经生成了它们），它适用于运行长度。
#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
我更喜欢基于预先计算的 tables 的解决方案。这比较慢，但我发现它更容易理解。然后我用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)
结果 table 将包含另外三列 1) 相应的日期 2) 相应的日落时间，以及 3) 日落后的时间。
这应该能够防止丢失测量值，因为日期编号出现在另一个 table 中。您还可以轻松修改算法以使用不同的时间，甚至应用多个。
更新
所有这些都可以使用 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 次观察 - 前一个大约需要 1m，data.table 解决方案是 0.011s。
我正在分析复杂数据集中的时间模式，该数据集包含多个环境变量以及 activity 来自不同动物物种的数据。这些数据由多个实验设置收集，每个设置的数据每分钟存储一次。该项目已经 运行 好几年了，所以我的数据集相当大。
我的一个数据集的前几行如下所示：
> 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
由于大多数动物的行为不同，取决于是白天还是晚上，我使用了 sunset/sunrise 次数据来计算一个新变量，该变量在夜间为 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 甚至可以使用 civil/nautical/astronomical dusk 的开始和黎明而不是日出和日落。
然而，这就是我的问题开始的地方。我想对实验中的所有日子进行编号。而不是像往常一样容易地在 午夜 增加日计数器，我想在 sunset 增加日计数器（或者可能在未来的实验中，一天中的另一个可移动时间，如日出、航海 dusk 和黎明，...）。由于日落不会每天都在同一时间发生，所以对我来说，这不是一个可以直接解决的问题。
我只提出了一个 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
}
}
除了丑陋和缓慢之外，这段代码还有一个大问题：它不处理缺失值。有时，由于设备故障，数小时或数天根本没有记录数据。如果在日落期间没有记录数据，则上述代码不会增加日计数器。这意味着我需要 - 以某种方式 - 合并 date/time 代码。很容易创建自实验开始以来天数的变量：
setup_01$daynumber<-as.integer(ceiling(difftime(setup_01$DateTime, setup_01$DateTime[1], units = "days")))
也许可以使用这些数字，可能与 Heroka's 不错的 rle-算法结合使用。
我已经使用 dput 从一个设置中获取了几个月的数据，包括一些大块丢失的数据，以及新创建的变量（如本 [=62= 中所述） ] 并在 Heroka's 答案中）可用 此处 .
我一直在寻找更好、更好、尤其是更快的东西，但一直无法想出一个好的技巧。我一直在摆弄我的数据框的子集，但得出的结论是这可能是一种愚蠢的方法。我查看了 maptools、lubridate 和 GeoLight。我搜索了 Google、Stack Overflow 和各种书籍，例如 Hadley Wickham 的精彩 Advanced R。都无济于事。也许我错过了一些明显的东西。我希望这里有人可以帮助我。
我想出了一个关于生成的 0 和 1 的解决方案（因为你已经生成了它们），它适用于运行长度。
#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
我更喜欢基于预先计算的 tables 的解决方案。这比较慢，但我发现它更容易理解。然后我用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)
结果 table 将包含另外三列 1) 相应的日期 2) 相应的日落时间，以及 3) 日落后的时间。
这应该能够防止丢失测量值，因为日期编号出现在另一个 table 中。您还可以轻松修改算法以使用不同的时间，甚至应用多个。
更新
所有这些都可以使用 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 次观察 - 前一个大约需要 1m，data.table 解决方案是 0.011s。