Plotting e-obs acc and gps burst data
Set workspace
# ==============================================================================
# * Environment
# ==============================================================================
library(dplyr)
library(lubridate)
library(data.table)
library(plotly)
library(tidyr)
library(hrbrthemes)
library(cowplot)
source("R/00_fxns.R")
# Set directory for saving plots
out_dir = file.path("results", "figs", "acc_gps_plots")
if(!dir.exists(out_dir)) dir.create(out_dir)Load acc + gps data and process for plotting
# ==============================================================================
# * Read in acceleration and GPS data
# ==============================================================================
# read in the behavior-classified acc and GPS data
acc_gps_file = file.path("data", "acc_behavL", "yz_acc_gps_behavL.RDS")
acc_gps_behavL = readRDS(acc_gps_file)
# ==============================================================================
# * Get acc and GPS data ready for plotting
# ==============================================================================
allIDs = names(acc_gps_behavL)
discard_tags = c("K5314", "K5316")
tagIDs = allIDs[ ! allIDs %in% discard_tags ]
# ID = tagIDs[ 1 ]
plot_info = lapply(tagIDs, function(ID, acc_behavL) {
cat("****", ID, "***\n")
ind = acc_gps_behavL[[ ID ]]
acc_behav = ind$acc_behav
gps_behav = ind$gps_behav
burst_info = ind$burst_info
# table(duplicated(acc_behav$timestamp))
# ** Make a long-format acc dataframe & add acc axes data
# ------------------------------------------------------------------------------
# lapply(burst_info, dim)
rep_acc = ncol(burst_info$xax)
acc_df = data.frame(timestamp = rep(acc_behav$timestamp, each = rep_acc),
activity = rep(acc_behav$activity, each = rep_acc),
stringsAsFactors = FALSE)
acc_df$xax = c(t(burst_info$xax))
acc_df$yax = c(t(burst_info$yax))
acc_df$zax = c(t(burst_info$zax))
# add an index column
acc_df$n = seq_len(nrow(acc_df))
# bk = acc_df
# acc_df = bk
# ** Gather acc data by axes & ready for plot
# ------------------------------------------------------------------------------
acc_plot = acc_df %>%
gather(key = "axis", value = "value", -c(timestamp, activity, n))
# ** Extract gps postions + bat activity + timestamps & map to acc data
# ------------------------------------------------------------------------------
gps = gps_behav %>%
select(activity, timestamp)
rownames(gps) = NULL
# Get nearest acc point for each gps (by timestamp)
acc_near_gps = lapply(gps$timestamp, get_near_acc, acc_behav) %>%
data.table::rbindlist() %>%
select(minN)
# Tabulate how often an acc reading is the closest for gps fixes
acc_gps_tb = acc_near_gps %>%
group_by(minN) %>%
tally
# Subset acc timestamps that were closest to gps fixes
near_acc_time = acc_behav$timestamp[ unique(acc_near_gps$minN) ]
# Subset rows in acc_df corresponding to selected timestamps
acc_df_gps = acc_df[ acc_df$timestamp %in% near_acc_time,
c("timestamp", "activity", "n")]
# Map the gps points to the long-format acc data
# This step helps plot the gps points above the corresponding acc burst
# E.g. if 10 gps points have "z" as the closest timestamp, they'll be
# plotted above that acc burst, corresponding to the first 10 points
# within that burst reading. They're mapped to acc's "n" column
acc_df_gpsL = split(acc_df_gps, acc_df_gps$timestamp)
gps_acc_N = lapply(seq_along(acc_df_gpsL), function(i, acc_df_gpsL, acc_gps_tb) {
df = acc_df_gpsL[[ i ]]
N = acc_gps_tb$n[ i ]
df[ 1:N, ]
}, acc_df_gpsL, acc_gps_tb) %>%
rbindlist()
# ** Get day info
# ------------------------------------------------------------------------------
daynight = dplyr::left_join(acc_df,
acc_behav[ , c("timestamp",
"sunrise",
"sunset",
"DayTime"), ],
by = "timestamp")
get_day_info = function(daynight) {
daynight$day = as.numeric(as.factor(daynight$DayTime))
start = daynight$day - lag(daynight$day)
end = daynight$day - lead(daynight$day)
start[1] = ifelse(daynight$day[1] %in% 1, -1, 1)
end[length(end)] = ifelse(daynight$day[length(end)] %in% 1, -1, 1)
start_n = start %in% -1
end_n = end %in% -1
data.frame(start_n = daynight$n[start_n],
start_ts = daynight$timestamp[ start_n ],
end_n = daynight$n[end_n],
end_ts = daynight$timestamp[ end_n])
}
day_info = get_day_info(daynight)
list(acc_df = acc_df,
acc_plot = acc_plot,
day_info = day_info,
gps_acc_N = gps_acc_N)
}, acc_behavL)## **** K5309 ***
## **** K5310 ***
## **** K5311 ***
## **** K5312 ***
## **** K5313 ***
## **** K5315 ***
## **** K5317 ***
## **** K5319 ***names(plot_info) = tagIDs
saveRDS(plot_info, "data/acc-gps_plot-info_yz.RDS")Plot acc & gps data
if(FALSE) {
ID = tagIDs[1]
}
sparse = FALSE
silent = TRUE
all_plots = lapply(tagIDs, function(ID, plot_info, sparse, silent, out_dir) {
if( !silent) cat("****", ID, "***\n")
acc_df = plot_info[[ ID ]]$acc_df
acc_plot = plot_info[[ ID ]]$acc_plot
day_info = plot_info[[ ID ]]$day_info
gps_acc_N = plot_info[[ ID ]]$gps_acc_N
# Max value of any axis reading
max_y = max(acc_plot$value + 200, na.rm = TRUE)
# Get x axis labels
x_lab = data.frame(n = c(day_info$start_n, day_info$end_n),
timestamp = c(day_info$start_ts, day_info$end_ts))
x_lab = x_lab[ order(x_lab$n), ]
x_lab$label = paste0(lubridate::month(x_lab$timestamp, label = TRUE),
lubridate::day(x_lab$timestamp),
"\n",
format(x_lab$timestamp, "%H:%M"))
tdiff = difftime(x_lab$timestamp, lag(x_lab$timestamp), units = "hours")
x_lab = x_lab[ tdiff >= 10 | is.na(tdiff), ]
# Subset acc_df for test plots
if(sparse) {
s = sample(1:nrow(acc_plot), size = 100000)
acc_plot = acc_plot[ s, ]
}
p = ggplot() +
geom_rect(data = day_info,
aes(xmin = start_n, xmax = end_n,
ymin = 0, ymax = max_y,
text = "Day"),
fill = "#F2F2F2") +
geom_line(data = acc_plot,
aes(x = n, y = value,
group = axis,
color = axis,
text = paste0(activity, "\n", timestamp)),
alpha = 0.6) +
scale_color_manual(values=c("#D55E00", "#009E73", "#0072B2"),
labels = c("X", "Y", "Z")) +
theme_ipsum(base_size = 10,
base_family = "Helvetica",
grid = FALSE,
axis_col = "#919191", axis = TRUE, ticks = TRUE) +
scale_x_continuous(breaks = x_lab$n,
labels = x_lab$label) +
ggtitle(ID) +
labs(color = "Acceleration\naxis") +
ylab("Acceleration\nreading") +
xlab("Time") +
theme(axis.title.y = element_text(angle = 0, hjust = 0, vjust = 0.5))
static_p = p + geom_segment(data = gps_acc_N[ gps_acc_N$activity %in% "NotFlying"],
aes(x = n,
y = (max_y + 100),
xend = n,
yend = (max_y + 200),
text = paste0(activity, "\n", timestamp)),
color = "#F2CB70") +
geom_segment(data = gps_acc_N[ gps_acc_N$activity %in% "Flying"],
aes(x = n,
y = (max_y + 100),
xend = n,
yend = (max_y + 200),
text = paste0(activity, "\n", timestamp)),
color = "#8C0C53")
gps_l = data.frame(activity = c("GPS: Flying", "GPS: Not Flying"), n = 2) %>%
ggplot(aes(activity, fill = activity)) +
geom_bar() + scale_fill_manual(values=c("#8C0C53", "#F2CB70")) +
theme(legend.position = "top",
legend.justification='right',
legend.key.size = unit(0.2, "cm"),
legend.text=element_text(size = 8, family = "Helvetica")) +
labs(fill = "")
gps_l = get_legend(gps_l)
day_l = data.frame(activity = "Day time", n = 2) %>%
ggplot(aes(x = activity, fill = activity))+
geom_bar() + scale_fill_manual(values="#F2F2F2") +
theme(legend.position = "top",
legend.justification='left',
legend.key.size = unit(0.2, "cm"),
legend.key.width = unit(0.6, "cm"),
legend.text=element_text(size = 8, family = "Helvetica")) +
labs(fill = "")
day_l = get_legend(day_l)
final_p = static_p +
annotation_custom(grob = day_l, ymin = max_y + 300) +
annotation_custom(grob = gps_l, ymin = max_y + 300)
print(final_p)
ggsave(file.path(out_dir, paste0(ID, ".png")), final_p, dpi = 600)
return(NULL)
}, plot_info, sparse = sparse, silent = silent, out_dir = out_dir)
Session info
sessionInfo()## R version 3.6.0 (2019-04-26)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] cowplot_1.0.0 hrbrthemes_0.6.0 tidyr_1.0.2 plotly_4.7.1.9900
## [5] ggplot2_3.3.0.9000 data.table_1.12.8 lubridate_1.7.4 dplyr_0.8.4
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.3 pillar_1.4.3 compiler_3.6.0 tools_3.6.0
## [5] extrafont_0.17 digest_0.6.23 viridisLite_0.3.0 jsonlite_1.6.1
## [9] evaluate_0.14 tibble_2.1.3 lifecycle_0.1.0 gtable_0.3.0
## [13] pkgconfig_2.0.3 rlang_0.4.4 yaml_2.2.1 xfun_0.12
## [17] Rttf2pt1_1.3.7 withr_2.1.2 stringr_1.4.0 httr_1.4.1
## [21] knitr_1.26 systemfonts_0.1.1 gdtools_0.2.1 vctrs_0.2.2
## [25] htmlwidgets_1.5.1 grid_3.6.0 tidyselect_1.0.0 glue_1.3.1
## [29] R6_2.4.1 rmarkdown_1.17.1 farver_2.0.3 extrafontdb_1.0
## [33] purrr_0.3.3 magrittr_1.5 ellipsis_0.3.0 scales_1.1.0
## [37] htmltools_0.4.0 assertthat_0.2.1 colorspace_1.4-1 labeling_0.3
## [41] stringi_1.4.5 lazyeval_0.2.2 munsell_0.5.0 crayon_1.3.4