R 使用 data.table 中的条件查找波的频率和持续时间高于给定值
R Find the frequency and duration a wave is above a given value using conditional in data.table
下面贴了一个MRE
MRE
date<-c('2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02')
time<-c('07:00:00 GMT','08:00:00 GMT','09:00:00 GMT','10:00:00 GMT','11:00:00 GMT','12:00:00 GMT','13:00:00 GMT','14:00:00 GMT','15:00:00 GMT','16:00:00 GMT','17:00:00 GMT', '18:00:00 GMT','19:00:00 GMT','20:00:00 GMT','21:00:00 GMT','22:00:00 GMT','23:00:00 GMT','00:00:00 GMT', '01:00:00 GMT','02:00:00 GMT','03:00:00 GMT','04:00:00 GMT','05:00:00 GMT','06:00:00 GMT','07:00:00 GMT','08:00:00 GMT','09:00:00 GMT','10:00:00 GMT','11:00:00 GMT','12:00:00 GMT','13:00:00 GMT','14:00:00 GMT','15:00:00 GMT','16:00:00 GMT','17:00:00 GMT','18:00:00 GMT','19:00:00 GMT','20:00:00 GMT','21:00:00 GMT')
el<-c(0.257,0.687,1.861,3.288, 4.821,6.172,7.048,7.258,6.799,5.654,4.463,3.443,2.704,2.708,3.328,4.23,5.244,5.985,6.317,6.074,5.234,3.981,2.662,1.615,0.88,0.746,1.405,2.527,3.928,5.283,6.517,7.179,7.252,6.625,5.454,4.214,3.144,2.491,2.357)
Time<-as.POSIXct(paste(date, time),tz="GMT")
wave<-data.table(Time, el)
ggplot(wave, aes(wave$Time, wave$el)) + geom_point() + labs(x="time", y="elevation") + geom_hline(aes(yintercept=4))
我有一个波浪时间序列,我希望能够有一个函数能够告诉我波浪高于给定海拔的频率和 mean/median 持续时间。在我的示例中,我选择了 4。
我想对波在上升沿和下降沿到达4的时间进行插值,求出每个波的两点之间的时间差。
我可以使用 for 循环来完成此操作,但我认为我应该能够在 data.table 中更快地完成此操作。我有几个位置的 100 万+ 点,我认为 for 循环不会有效。
对于上升波我想做这样的事情:
wave[,timeIs4:=ifelse(elev<3 & elev[+1]>4,TRUE,FALSE )]
但在我的插值计算中输入的不是 TRUE。我不知道如何访问数据 table 中的前面和后面的值,例如在 for 循环 i+1 或 i-1 中。
期望的输出
上升腿
我想在第 4 点和第 5 点之间进行插值; 15 和 16; 29 和 30。
落腿
我想在第 11 点和第 12 点之间进行插值; 21 和 22; 36 和 37
大概结果
Rising Falling
10:28:00 17:27:00
21:45:00 3:59:00
11:03:00 18:12:00
然后我将能够使用 difftime() 从下降中减去上升以确定水位高于给定高度的时间量。
这会给我水高于给定海拔的频率和持续时间。
这是使用 devel
version from GH 的可能解决方案。 shift
函数(如@Jan 所述)和接受表达式的新 dcast
方法都需要它。此外,您的 MRE 中没有分钟数,因此不确定您从哪里获得预期输出中的分钟数。
无论如何,对于初学者,我们将创建一个索引(我们将其称为 Wave
,以便您知道#它来自哪个波浪),它会告诉我们波浪是上升还是下降使用 shift
。然后,我们将 dcast
匹配值,同时使用 na.omit
删除不匹配的值(如果您喜欢使用 setcolorder
函数,您可以稍后重新排序列名)
library(data.table) ## V 1.9.5+
dt[elev <= 4 & shift(elev, type = "lead") > 4, Wave := "Rising"]
dt[elev > 4 & shift(elev, type = "lead") <= 4, Wave := "Falling"]
dcast(na.omit(dt), cumsum(Wave == "Rising") ~ Wave, value.var = "time")
# Wave Falling Rising
# 1: 1 2001-01-01 17:00:00 2001-01-01 10:00:00
# 2: 2 2001-01-02 03:00:00 2001-01-01 21:00:00
# 3: 3 2001-01-02 18:00:00 2001-01-02 11:00:00
这是另一个可能的想法:
elev = 4
#a helper function to calculate elapsed time
ff = function(el1, el2, el, time1, time2)
time1 + ((el - el1) / (el2 - el1)) * (time2 - time1)
dif = diff(findInterval(wave$el, c(-Inf, elev, Inf)))
ris = which(dif == 1) #risings
fal = which(dif == -1) #fallings
ff(wave$el[ris], wave$el[ris + 1], elev, wave$Time[ris], wave$Time[ris + 1])
#[1] "2001-01-01 10:27:52 GMT" "2001-01-01 21:44:42 GMT" "2001-01-02 11:03:11 GMT"
ff(wave$el[fal], wave$el[fal + 1], elev, wave$Time[fal], wave$Time[fal + 1])
#[1] "2001-01-01 17:27:14 GMT" "2001-01-02 03:59:05 GMT" "2001-01-02 18:12:00 GMT"
下面贴了一个MRE
MRE
date<-c('2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-01','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02','2001-01-02')
time<-c('07:00:00 GMT','08:00:00 GMT','09:00:00 GMT','10:00:00 GMT','11:00:00 GMT','12:00:00 GMT','13:00:00 GMT','14:00:00 GMT','15:00:00 GMT','16:00:00 GMT','17:00:00 GMT', '18:00:00 GMT','19:00:00 GMT','20:00:00 GMT','21:00:00 GMT','22:00:00 GMT','23:00:00 GMT','00:00:00 GMT', '01:00:00 GMT','02:00:00 GMT','03:00:00 GMT','04:00:00 GMT','05:00:00 GMT','06:00:00 GMT','07:00:00 GMT','08:00:00 GMT','09:00:00 GMT','10:00:00 GMT','11:00:00 GMT','12:00:00 GMT','13:00:00 GMT','14:00:00 GMT','15:00:00 GMT','16:00:00 GMT','17:00:00 GMT','18:00:00 GMT','19:00:00 GMT','20:00:00 GMT','21:00:00 GMT')
el<-c(0.257,0.687,1.861,3.288, 4.821,6.172,7.048,7.258,6.799,5.654,4.463,3.443,2.704,2.708,3.328,4.23,5.244,5.985,6.317,6.074,5.234,3.981,2.662,1.615,0.88,0.746,1.405,2.527,3.928,5.283,6.517,7.179,7.252,6.625,5.454,4.214,3.144,2.491,2.357)
Time<-as.POSIXct(paste(date, time),tz="GMT")
wave<-data.table(Time, el)
ggplot(wave, aes(wave$Time, wave$el)) + geom_point() + labs(x="time", y="elevation") + geom_hline(aes(yintercept=4))
我有一个波浪时间序列,我希望能够有一个函数能够告诉我波浪高于给定海拔的频率和 mean/median 持续时间。在我的示例中,我选择了 4。
我想对波在上升沿和下降沿到达4的时间进行插值,求出每个波的两点之间的时间差。
我可以使用 for 循环来完成此操作,但我认为我应该能够在 data.table 中更快地完成此操作。我有几个位置的 100 万+ 点,我认为 for 循环不会有效。
对于上升波我想做这样的事情:
wave[,timeIs4:=ifelse(elev<3 & elev[+1]>4,TRUE,FALSE )]
但在我的插值计算中输入的不是 TRUE。我不知道如何访问数据 table 中的前面和后面的值,例如在 for 循环 i+1 或 i-1 中。
期望的输出
上升腿 我想在第 4 点和第 5 点之间进行插值; 15 和 16; 29 和 30。
落腿 我想在第 11 点和第 12 点之间进行插值; 21 和 22; 36 和 37
大概结果
Rising Falling
10:28:00 17:27:00
21:45:00 3:59:00
11:03:00 18:12:00
然后我将能够使用 difftime() 从下降中减去上升以确定水位高于给定高度的时间量。
这会给我水高于给定海拔的频率和持续时间。
这是使用 devel
version from GH 的可能解决方案。 shift
函数(如@Jan 所述)和接受表达式的新 dcast
方法都需要它。此外,您的 MRE 中没有分钟数,因此不确定您从哪里获得预期输出中的分钟数。
无论如何,对于初学者,我们将创建一个索引(我们将其称为 Wave
,以便您知道#它来自哪个波浪),它会告诉我们波浪是上升还是下降使用 shift
。然后,我们将 dcast
匹配值,同时使用 na.omit
删除不匹配的值(如果您喜欢使用 setcolorder
函数,您可以稍后重新排序列名)
library(data.table) ## V 1.9.5+
dt[elev <= 4 & shift(elev, type = "lead") > 4, Wave := "Rising"]
dt[elev > 4 & shift(elev, type = "lead") <= 4, Wave := "Falling"]
dcast(na.omit(dt), cumsum(Wave == "Rising") ~ Wave, value.var = "time")
# Wave Falling Rising
# 1: 1 2001-01-01 17:00:00 2001-01-01 10:00:00
# 2: 2 2001-01-02 03:00:00 2001-01-01 21:00:00
# 3: 3 2001-01-02 18:00:00 2001-01-02 11:00:00
这是另一个可能的想法:
elev = 4
#a helper function to calculate elapsed time
ff = function(el1, el2, el, time1, time2)
time1 + ((el - el1) / (el2 - el1)) * (time2 - time1)
dif = diff(findInterval(wave$el, c(-Inf, elev, Inf)))
ris = which(dif == 1) #risings
fal = which(dif == -1) #fallings
ff(wave$el[ris], wave$el[ris + 1], elev, wave$Time[ris], wave$Time[ris + 1])
#[1] "2001-01-01 10:27:52 GMT" "2001-01-01 21:44:42 GMT" "2001-01-02 11:03:11 GMT"
ff(wave$el[fal], wave$el[fal + 1], elev, wave$Time[fal], wave$Time[fal + 1])
#[1] "2001-01-01 17:27:14 GMT" "2001-01-02 03:59:05 GMT" "2001-01-02 18:12:00 GMT"