-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path2-define_decision_points_add_SPL.R
399 lines (309 loc) · 14.4 KB
/
2-define_decision_points_add_SPL.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
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
## ---------------------------
##
## Script name: 2-define_decision_points_add_SPL
##
## Purpose of script: This script defines decision points, splits possible trajectories
## into 'segments', and matches these to modelled microbarom infrasound
##
## Author: Dr. Lucia Martina Martin Lopez (with support by Dr Natasha Gillies)
##
## Email: gilliesne@gmail.com
##
## ---------------------------
### 0.0 Load the packages ------------------------------------------------------
# 0.0.0 Define the packages
packages <- c("dplyr", "ncdf4", "lubridate", "birk", "pracma")
# Install packages not yet installed - change lib to library path
#installed_packages <- packages %in% rownames(installed.packages())
#if (any(installed_packages == FALSE)) {
# install.packages(packages[!installed_packages], lib = "C:/Users/libraryPath")
#}
# 0.0.1 Load packages
invisible(lapply(packages, library, character.only = TRUE))
# 0.0.2 Load boutFinder custom function
source("FUNCTION_boutFinder.R")
## 0.1 Load the data -----------------------------------------------------------
gps_2013 <- read.csv("Data_inputs/WAAL_2013_gps_labelled.csv", stringsAsFactors = F)
### 1.0 Create travel bouts ----------------------------------------------------
# 1.0.0 Index all rows in each trip for each bird
gps_2013Trav <- gps_2013 %>%
group_by(TripID) %>%
mutate(idx = seq(1:n())) %>%
# 1.0.1 Filter fixes labelled as 'travel'
filter(State == "Travel") %>%
# 1.0.2 Label individual bouts
mutate(travbout = boutFinder(idx)) %>%
data.frame()
gps_2013Trav$travbout <- as.factor(as.character(gps_2013Trav$travbout))
## 1.2 Estimate distance travelled in each bout & filter for travel bouts > 20km ----
gps_2013Trav20 <- gps_2013Trav %>%
group_by(TripID, travbout) %>%
mutate(TotdisttravBout = sum(DistTrav)) %>%
filter(TotdisttravBout > 20) %>%
data.frame()
# 1.2.0 Output dataframe - used for sensitivity analysis
#write.csv(gps_2013Trav20, "Data_inputs/WAAL_2013_gps_Trav20.csv", row.names = F)
## 1.3 Identify the first point of each travelling period (decision point) -----
gps_2013Trav20_1stpoint <- gps_2013Trav20 %>%
group_by(TripID, travbout) %>%
slice_head(n = 1) %>%
ungroup()
### 2.0 Match each travel bout to SPL map --------------------------------------
#### NOTE: This script may take several hours to run depending on computing power.
## 2.1 Point to soundscape data ------------------------------------------------
path_to_IS_maps <- "E:/Soundscapes/"
IS_folder_maps <- dir(path_to_IS_maps, pattern = "2013")
GPS_ID_segments <- list()
## 2.2 Define segment parameters -----------------------------------------------
aperture <- 60
transectlength <- 2025
segmentno <- 360/aperture
angDiffs <- c(seq(0, 180, by = aperture),
seq((-180 + aperture),-aperture, aperture))
for (x in 1:length(IS_folder_maps)) {
## 2.3 Point to maps and isolate bird data -----------------------------------
# 2.3.1 Point to IS map relevant to bird
birdmapid <- IS_folder_maps[x]
path_to_IS_files <- paste0(path_to_IS_maps, birdmapid, "/Interp/")
# 2.3.2 Load infrasound map files
IS_files <- list.files(path = path_to_IS_files)
# 2.3.3 Extract bird ID from the folder name
IDbirdmap <- sapply(strsplit(IS_folder_maps[x], "_", fixed = TRUE),
function(i) paste(head(i, -1), collapse = "_"))
if (grepl("_", IDbirdmap, fixed = TRUE) == FALSE) { IDbirdmap <- paste0(IDbirdmap, "_1")}
# Change underscore to '.'
IDbirdmap <- gsub("_", ".", IDbirdmap)
# 2.3.4 Select the bird trip for which the soundscape is being loaded
gps_2013_ID1 <- gps_2013Trav20_1stpoint %>%
filter(TripID == IDbirdmap)
# If there are no GPS data for that map, skip
if (nrow(gps_2013_ID1) == 0) next
# Check data formatting
gps_2013_ID1$TripID <- as.factor(as.character(gps_2013_ID1$TripID))
gps_2013_ID1$DateTime <- as.POSIXct(gps_2013_ID1$DateTime, tz = "UTC")
# 2.3.5 Get closest hour for each GPS fix to match with soundscapes
gps_2013_ID1$maptobe <-
format(round(gps_2013_ID1$DateTime, units = "hours"), format = "%Y-%m-%d %H:%M")
## 2.4 Match each GPS fix to the closest SPL map in time and compute segments ----
for (i in 1:nrow(gps_2013_ID1)) {
# 2.4.1 Find the closest hour for each GPS point in the IS_files
TG <- gps_2013_ID1$maptobe[i]
Maptobe <- paste0(sub("\\_Soundscape_.*",
"",
IS_files[i]),
"_Soundscape_",
substr(TG, 1, 4),
substr(TG, 6, 7),
substr(TG, 9, 10),
"_",
substr(TG, 12, 13),
".nc")
IDX <- which(IS_files == Maptobe)
# If no map is found then fill it with NAs.
if (length(IDX) == 0) {
segments = data.frame(
segment_ID = as.numeric(NA),
segment_n = as.numeric(NA),
segment_vert_lef = as.numeric(NA),
segment_vert_rig = as.numeric(NA),
abs_SPL_2000dB = as.numeric(NA),
abs_SPL_2000 = as.numeric(NA),
abs_SPL_2000_std = as.numeric(NA),
abs_SPL_2000dB_std = as.numeric(NA),
birdID = as.character(birdmapid),
TripID = as.character(NA),
mapID = as.character(TG),
counter = as.numeric(NA),
Sex = as.character(gps_2013_ID1$Sex),
x_lon = as.numeric(NA),
y_lat = as.numeric(NA),
State = as.character(NA),
Dist_cro_shelf = as.numeric(NA),
WindSp = as.numeric(NA), # wind speed
Dev.wind2 = as.numeric(NA), # wind direction relative to track direction with directionality removed i.e. 0 to 180 rather than -180 to 180)
relDir_adj.bearing = as.numeric(NA),
stringsAsFactors = FALSE
)
} else {
ID = IS_files[IDX]
# 2.4.2 Extract the dateTime from the SPL map
a <- paste0(substr(ID, 19, 22),
":",
substr(ID, 23, 24),
":",
substr(ID, 25, 26),
":",
substr(ID, 28, 29),
":00:00")
mapDateTime <- parse_date_time(a, 'ymd HMS')
mapDateTime <- as.POSIXct(mapDateTime, tz = "UTC")
# 2.4.3 Load the nc file - extract SPL for each latitude/longitude
nc_file <- nc_open(paste0(path_to_IS_files, ID))
latitude <- ncvar_get(nc_file, "lat") #deg
longitude <- ncvar_get(nc_file, "lon") # deg
OSWIPa <- ncvar_get(nc_file, "OSWI_interp_spatial_time") #Interpolated OSWI + TLloss + NormGrid + Time in Pa
Pref <- 20e-6
OSWIdB <- 10 * log10(OSWIPa/(Pref^2))
latitude <- as.vector(latitude)
longitude <- as.vector(longitude)
SPL_dB <- as.vector(OSWIdB)
SPL_Pa <- as.vector(OSWIPa)
baz <- ncvar_get(nc_file,"baz") # back azimuth angle, clockwise to north
baz <- as.vector(baz)
Gdist <- ncvar_get(nc_file,"dist") # Geodesic distance in degrees
X <- as.data.frame(cbind(latitude,longitude,SPL_dB,SPL_Pa, Gdist,baz))
X2 <- X[complete.cases(X * 0), , drop = FALSE]
nc_close(nc_file)
# 2.4.4 Set Back zenith angles (BAZ) in accordance with the GPS Bearing
X2$baz_converted <- ifelse(
X2$baz <= 90 & X2$baz >= 0,
X2$baz - 180,
ifelse(
X2$baz > 90 & X2$baz <= 180,
X2$baz - 180,
ifelse(
X2$baz >= -180 & X2$baz <= -90,
(X2$baz + 180),
ifelse(X2$baz > -90 &
X2$baz < 0, (X2$baz + 180), NA)
)
)
)
## 2.5 Within each soundscape map divide the area into 6 segments of 60 deg each ----
## Get the baz angles for the 12 segments starting from the left side of
## the focal segment, and then clockwise.
## Always the same relative angles from the ontrack one!
# 2.5.1 Set the focal segment as +/- 30 degrees from bird's bearing
focal_segment_deg <-
c( round(gps_2013_ID1$Bearing[i]) - aperture / 2,
round(gps_2013_ID1$Bearing[i]) + aperture / 2 )
if (any(abs(focal_segment_deg) > 180)) {
a <- which(abs(focal_segment_deg) > 180)
if (focal_segment_deg[a] < 0) {
focal_segment_deg[a] <- 360 + focal_segment_deg[a]
} else{
focal_segment_deg[a] <- 360 - focal_segment_deg[a]
}
}
# 2.5.2 Create the other segments from the focal segment
if (focal_segment_deg[1] <= (-180 + aperture)) {
segment_vert_lef <-
(seq(from = focal_segment_deg[1], to = 180, by = aperture))
segment_vert_rig <-
c(segment_vert_lef[c(2:length(segment_vert_lef), 1)])
} else {
segment_vert_lef <-
c(seq(from = focal_segment_deg[1], to = 180, by = aperture),
rev(
seq(
from = focal_segment_deg[1] - aperture,
to = -180,
by = -aperture
)
))
segment_vert_rig <-
c(segment_vert_lef[c(2:length(segment_vert_lef), 1)])
}
# 2.5.3 Add the exception to remove the 0 degree segment that will be created
# when segment_vert_lef = 180 and segment_vert_rig = -180
if (any(abs(diff(segment_vert_lef)) == 360)) {
segment_vert_lef <-
segment_vert_lef[-which(diff(segment_vert_lef) != aperture)]
segment_vert_rig <-
c(segment_vert_lef[c(2:length(segment_vert_lef), 1)])
}
if (length(segment_vert_lef) > segmentno) {
segment_vert_lef <- c(segment_vert_lef[c(1:segmentno)])
segment_vert_rig <- c(segment_vert_rig[c(1:segmentno)])
}
# 2.5.4 Label the segments as focal vs non-focal and numeric ID (in
# clockwise direction)
segment_ID <- c(0,rep(1, segmentno - 1))
segment_n <- seq(1,segmentno)
# 2.5.5 Make a dataframe containing all the segments
segments <- as.data.frame(cbind(segment_ID, # segment identity (focal vs non)
segment_n, # segment number
segment_vert_lef, # left bound
segment_vert_rig)) # right bound
## 2.6 For each segment estimate the abs & standarized SPL and gdist to 45dB ----
# 2.6.0 Add new empty column to dataframe to be filled
segments$abs_SPL_2000dB <- NA
segments$abs_SPL_2000 <- NA
# 2.6.1 Loop through segments and calculate integrated SPL in each
for (c in 1:nrow(segments)) {
if (segment_vert_lef[c] > 0 & segment_vert_rig[c] < 0) {
X2_1 <- X2 %>%
filter(X2$baz_converted >= segment_vert_lef[c] &
X2$baz_converted <= 180)
X2_2 <- X2 %>%
filter(X2$baz_converted <= segment_vert_rig[c] &
X2$baz_converted >= -180)
newX <- rbind(X2_1, X2_2)
abs_SPL_2000 <- newX %>%
filter(Gdist <= transectlength) %>%
summarise(x = sum(SPL_Pa))
segments$abs_SPL_2000[c] <- as.numeric(abs_SPL_2000)
abs_SPL_2000dB <- newX %>%
filter(Gdist <= transectlength) %>%
summarise(x = 10 * log10(sum(SPL_Pa) / (Pref ^ 2)))
segments$abs_SPL_2000dB[c] <- as.numeric(abs_SPL_2000dB)
} else {
abs_SPL_2000<- X2 %>%
filter(X2$baz_converted >= segment_vert_lef[c] & X2$baz_converted <= segment_vert_rig[c] & Gdist <= transectlength) %>%
summarise(x = sum(SPL_Pa))
segments$abs_SPL_2000[c]<-as.numeric(abs_SPL_2000)
abs_SPL_2000dB <- X2 %>%
filter( X2$baz_converted >= segment_vert_lef[c] &
X2$baz_converted <= segment_vert_rig[c] &
Gdist <= transectlength ) %>%
summarise(x = 10 * log10(sum(SPL_Pa) / (Pref ^ 2)))
segments$abs_SPL_2000dB[c] <- as.numeric(abs_SPL_2000dB)
}
}
segments$abs_SPL_2000_std <- scale(segments$abs_SPL_2000)
segments$abs_SPL_2000dB_std <- scale(segments$abs_SPL_2000dB)
segments$birdID <- birdmapid
segments$TripID <- gps_2013_ID1$TripID[i]
segments$mapID <- TG
segments$counter = gps_2013_ID1$counter[i]
segments$Sex = gps_2013_ID1$Sex[i]
segments$x_lon = gps_2013_ID1$x_lon[i]
segments$y_lat = gps_2013_ID1$y_lat[i]
segments$State = gps_2013_ID1$State[i]
segments$Dist_cro_shelf = gps_2013_ID1$Dist_cro_shelf[i]
segments$WindSp = gps_2013_ID1$WindSp[i] # wind speed
segments$Dev.wind2 = gps_2013_ID1$Dev.wind2[i] # wind direction relative to track direction with directionality removed i.e. 0 to 180 rather than -180 to 180)
## 2.7 Find wind direction for each segment -------------------------------
# 2.7.1 Find angular differences between each segment
segments$segment_vert_lef.DIFF <- rep(angDiffs, nrow(segments) / segmentno)
# 2.7.2 Calculate relative wind direction adjusted for segment and convert to c(-180, 180)
segments$relDir_adj <- segments$Dev.wind2 + segments$segment_vert_lef.DIFF
segments$relDir_adj.bearing <- abs(ifelse(segments$relDir_adj > 180, -360 + segments$relDir_adj,
segments$relDir_adj))
# 2.7.3 Remove columns used for calculation
segments$segment_vert_lef.DIFF <- NULL
segments$relDir_adj <- NULL
segments$relDir <- NULL
}
if (x == 1) {
GPS_ID_segments = segments
} else{
GPS_ID_segments = rbind(GPS_ID_segments, segments)
}
}
print(x)
}
### 3.0 Process and output segment dataframe ----------------------------------
# 3.0.0 Remove NA variables
GPS_ID_segments <- subset(GPS_ID_segments, !is.na(segment_ID))
# 3.0.1 Rename relDir columns
GPS_ID_segments <- rename(GPS_ID_segments, relDir = relDir_adj.bearing)
# 3.0.2 Make (decision) Point ID variable
GPS_ID_segments$pointID <- paste(GPS_ID_segments$TripID, GPS_ID_segments$counter, sep = ".")
# 3.0.3 Make focal segment '1' and non-focal '0'
GPS_ID_segments$segment_ID <- abs(GPS_ID_segments$segment_ID - 1)
# 3.0.4 Remove columns not relevant for analysis
cols_rm <- c("segment_vert_lef", "segment_vert_rig", "mapID", "counter", "Dev.wind2")
GPS_ID_segments[,cols_rm] <- NULL
## 3.1 Write to csv ------------------------------------------------------------
write.csv(GPS_ID_segments, paste0("./Data_inputs/WAAL_2013_gps_processed_aperture", aperture, "deg.csv"), row.names = F)