-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path03_wavelet-analysis.R
156 lines (132 loc) · 8.97 KB
/
03_wavelet-analysis.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
## -----------------------------------------------------
##
## Script name: 03_wavelet-analysis.R
##
## Purpose of script: Perform a wavelet analysis on hourly spurdog mean depth and add variable DVM
##
## Author: Antonia Kloecker
##
## Institution: Institute of Marine Research (IMR), Norway
##
## E-Mail: antonia.kloecker@hi.no
##
## Date Created: 2023-09-19
#
## Input
# - `archives-hourly-cleaned.csv`
## Output
# - `archive_hourly_nested_wavelet.rds`
# - `archive_hourly_24hwavelet.rds`
## -----------------------------------------------------
## Load packages ########################################
library(tidyverse)
library(WaveletComp)
## Read Data ############################################
hdata <- readr::read_rds("data/spurdog/processed_data/archives-hourly-cleaned.rds")
## Wavelet Analysis #####################################
#depth 1/24
###################
# nest wavelet analysis in tibble
wt24 <- hdata %>% rename(date = DateTime) %>% nest(.by = tagID) %>%
mutate(wavelet = map(.x = data, ~WaveletComp::analyze.wavelet(my.data = .x,
my.series = "Depth_median",
loess.span = 0, # no detrending
dt = 1/24, # every hour
# here we choose 1/12 so we have an hour around the 24h and respective signf. levels
dj = 1/12, # "resolution" (no of suboctaves) - should correspond with col n.levels
lowerPeriod = 1/4, # 6 h,
upperPeriod = 2^7, # 128 d
# null model is AR1 with p = 0.5
method = "AR",
params = list(p = 0.5),
#make.pval = F,
n.sim = 1000))) # choose bigger value min. 100 better 1000
readr::write_rds(x = wt24, file = "data/spurdog/processed_data/archive_hourly_nested_wavelet_depth_1_24.rds", compress = "gz")
wt24_red <- wt24 %>% select(-data)
readr::write_rds(x = wt24_red , file = "data/spurdog/processed_data/archive_hourly_nested_wavelet_depth_1_24_wodata.rds", compress = "gz")
# wt24 <- readr::read_rds(file = "data/spurdog/processed_data/archive_hourly_nested_wavelet_depth_1_24.rds")
# for comparative analysis on same time period
###################
hdata_comp <- hdata %>% mutate(comp = case_when(sharkID %in% 1:4 & Date >= "2019-12-14" & Date <= "2020-05-23" ~ T,
sharkID %in% c(6:8,10) & Date >= "2020-12-14" & Date <= "2021-05-24" ~ T,
sharkID %in% c(11,12,14,15) & Date >= "2021-12-14" & Date <= "2022-05-24" ~ T,
sharkID %in% 16:19 & Date >= "2022-12-14" & Date <= "2023-05-24" ~ T,
.default = F))
#hdata_comp %>% filter(comp ==T) %>% group_by(sharkID) %>% summarise(n()) # check that they are same length
wt_comp <- hdata_comp %>% filter(comp ==T) %>% rename(date = DateTime) %>% nest(.by = tagID) %>%
mutate(wavelet = map(.x = data, ~WaveletComp::analyze.wavelet(my.data = .x,
my.series = "Depth_median",
loess.span = 0, # no detrending
dt = 1/24, #
dj = 1/12, # "resolution" (no of suboctaves) - should correspond with col n.levels
lowerPeriod = 1/4, # 6 h,
upperPeriod = 2^7, # 128 d
# null model is AR1 with p = 0.5
method = "AR",
params = list(p = 0.5),
#make.pval = F,
n.sim = 100))) # choose bigger value min. 100 better 1000
readr::write_rds(x = wt_comp, file = "data/spurdog/processed_data/archive_hourly_nested_wavelet_comp_1_24.rds", compress = "gz")
#wt_comp <- readr::read_rds(file = "data/spurdog/processed_data/archive_hourly_nested_wavelet_comp_1_24.rds")
#wavelet for fast starts
###########################
wtFS <- hdata %>% rename(date = DateTime) %>% nest(.by = tagID) %>%
mutate(wavelet = map(.x = data, ~WaveletComp::analyze.wavelet(my.data = .x,
my.series = "FS95_cum",
loess.span = 0, # no detrending
dt = 1/24, # every hour
# here we choose 1/12 so we have an hour around the 24h and respective signf. levels
dj = 1/12, # "resolution" (no of suboctaves) - should correspond with col n.levels
lowerPeriod = 1/4, # 6 h,
upperPeriod = 2^7, # 128 d
# null model is AR1 with p = 0.5
method = "AR",
params = list(p = 0.5),
#make.pval = F,
n.sim = 1000))) # choose bigger value min. 100 better 1000
readr::write_rds(x = wtFS, file = "data/spurdog/processed_data/archive_hourly_nested_wavelet_FS95cum_1_24.rds", compress = "gz")
#wtFS <- readr::read_rds(file = "data/spurdog/processed_data/archive_hourly_nested_wavelet_FS95cum_1_24.rds")
wtFS_red <- wtFS %>% select(-data)
readr::write_rds(x = wtFS_red , file = "data/spurdog/processed_data/archive_hourly_nested_wavelet_FS95cum_1_24_wodata.rds", compress = "gz")
## Generate scalogramms ####################################
#wt <- readr::read_rds(file = "data/spurdog/processed_data/archive_hourly_nested_wavelet.rds", compress = "gz")
## Generate avg power df to prepare plots ##############################
get_wtwave <- function(wavelet){return(wavelet$Wave)}
get_wtperiod <- function(wavelet){return(wavelet$Period)}
get_wtpw <- function(wavelet){return(wavelet$Power.avg)}
get_wtpw_sig <- function(wavelet){return(wavelet$Power.avg.pval)}
wtavg <- wt24 %>%
#extract average power across periods and respective significance
mutate(period = map(.x = wavelet, get_wtperiod),
pw_avg = map(.x = wavelet, get_wtpw),
pw_avg_sig = map(.x = wavelet, get_wtpw_sig)) %>%
# remove wavelet list to avoid replicating them when unnesting
select(-wavelet) %>% select(-data) %>%
unnest(c(period, pw_avg, pw_avg_sig))
readr::write_rds(x = wtavg, file = "data/spurdog/processed_data/archive_hourly_depth_avgwavelet.rds", compress = "gz")
#readr::write_rds(x = wtavg, file = "data/spurdog/processed_data/archive_hourly_depth_FS95cum_avgwavelet.rds", compress = "gz")
### Extract comparive wave data for dendrogramme###########
wtc <- wt_comp %>% mutate(wave = map(.x = wavelet, get_wtwave)) %>% select(-wavelet) %>% select(-data)
readr::write_rds(x = wtc, file = "data/spurdog/processed_data/archive_hourly_depth_compWave.rds", compress = "gz")
## Assess if DVM is present ##############################
#define functions to extract 23-25 h power and significance level
get_24h <- function(wavelet) {
poi <- which(wavelet$Period > 0.958 & wavelet$Period < 1.042)
return(wavelet$Power[poi,])
}
get_24h_sig <- function(wavelet) {
poi <- which(wavelet$Period > 0.958 & wavelet$Period < 1.042)
return(wavelet$Power.pval[poi,])
}
hdata_wtDVM <- wt24 %>%
# extract power & resp. significance level at period 23-25 h
mutate(power24 = map(wavelet, get_24h),
power24_sig = map(wavelet, get_24h_sig)) %>%
# remove wavelet list to avoid replicating them when unnesting
select(-wavelet) %>%
unnest(c(data, power24, power24_sig)) %>%
# add logical variable if considered as DVM or not (sig.level <= 0.05)
mutate(DVM = if_else(power24_sig > 0.05, "nonDVM", "DVM")) %>%
# reverse the renaming for WaveletComp
rename(DateTime = date)
readr::write_rds(x = hdata_wtDVM, file = "data/spurdog/processed_data/archive_hourly_DVM_wavelet.rds", compress = "gz")