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- петля, которая не очень хороший способ делать вещи. Кроме того, учитывая, что у меня есть данные за 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-алгоритм.
я использовал dput сделать данные за несколько месяцев из одной установки, включая несколько больших кусков отсутствующих данных, а также недавно созданные переменные (как описано в этом посте и в Heroka), доступные здесь.
Я искал что-то лучше, лучше и особенно быстрее, но не смог придумать хороший трюк. Я возился с подмножеством своего информационного кадра, но пришел к выводу, что это, вероятно, глупый подход. Я смотрел на maptools, lubridate, а также GeoLight, Я искал Google, Stack Overflow и различные книги, такие как фантастический Advanced R Хэдли Уикхэма. Все безрезультатно. Возможно, я упускаю что-то очевидное, хотя. Я надеюсь, что кто-то здесь может мне помочь.
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 с.
Я придумал решение для сгенерированных 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