R:计算在日落时开始的日子
我正在分析由多个环境变量组成的复杂数据集中的时间模式以及来自各种动物物种的活动数据。 这些数据已经通过多个实验设置收集,每个设置的数据每分钟存储一次。 这个项目已经运行好几年了,所以我的数据集相当大。
其中一个数据集的前几行如下所示:
> 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
-loop,这不是一个很好的做事方式。 另外,考虑到我在一些设置中每分钟收集一次超过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的很好结合rle
-algorithm。
我已经使用dput
从一个设置中创建了几个月的数据,其中包括一些缺少数据的大块数据,以及新创建的变量(如本文中和Heroka的答案中所述)。
我一直在寻找更好,更好,更快的东西,但一直未能提出一个好方法。 我已经调整了我的数据框的子集,但得出的结论是这可能是一个愚蠢的做法。 我看过maptools
, lubridate
和GeoLight
。 我搜索了Google,Stack Overflow和各种书籍,如Hadley Wickham的梦幻般的Advanced R.都无济于事。 也许我错过了很明显的事情。 我希望这里有人能帮助我。
我更喜欢基于预先计算表的解决方案。 这很慢,但我发现它更清楚地理解。 然后我使用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个观察值进行计时 - 以前需要约1m,data.table解决方案为0.011s。
我想出了一个生成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
链接地址: http://www.djcxy.com/p/87339.html
上一篇: R: count days that start at sunset
下一篇: Set initial height of parallax image in CollapsingToolbarLayout