Setting up the workspace

# ==============================================================================
# * Libraries
# ==============================================================================

library(move)
library(lubridate)
library(dplyr)
library(data.table)
library(leaflet)
library(htmlwidgets)
library(sf)
library(ggplot2)
library(patchwork)
library(ggsn)
library(hrbrthemes)
library(DT)
library(viridis)
library(tidyr)
library(readxl)
library(cowplot)
source("R/00_fxns.R")

extrafont::loadfonts()

# Directory for saving plots and KDE information
plot_out_dir = file.path("results", "figs")
gps_info_out_dir = file.path("results", "tables")


# ==============================================================================
# * Read in GPS data with flying behaviour + tagging location
# ==============================================================================
# read in the behavior-classified acc and GPS data
gps = readRDS("data/gps_forage.RDS")
hr = read.csv("results/tables/KDE-area_eobs.csv", stringsAsFactors = FALSE)
hr$X = NULL

bats = read.csv("data/raw/bats/eobs_info.csv", stringsAsFactors = FALSE)
bats = bats[ , c("tagID", "location", "tag_type", "sex", "body_wt", "date")]
bats$tagID = paste0("K", bats$tagID)



GPS information

# * GPS time range
# ==============================================================================
gps_info = bats

tr = lapply(gps, function(df)
    data.frame(ts_min = min(df$timestamp),
               ts_max = max(df$timestamp))) %>%
    rbindlist()

gps_info = gps_info %>%
    mutate(
        date_deployed = format(tr$ts_min, "%b-%d"),
        hours = round(time_length(tr$ts_max - tr$ts_min, "hours"), 2),
        days = round(time_length(tr$ts_max - tr$ts_min, "days"), 2))

gps_info$foraging_nights = sapply(gps, function(df)
    max(df$time_cut, na.rm = TRUE))

gps_info$gps_fixes = sapply(gps, nrow)


# ==============================================================================
# * Cumulative distances
# ==============================================================================
dist_info = sapply(gps, function(bat) {
    nightly_dist = sapply(split(bat, bat$time_cut), get_dist_km)
    dist_mean = round(mean(nightly_dist), 2)
    dist_range = paste0(round(range(nightly_dist), 2), collapse = ", ")
    cum_dist = paste0(dist_mean, " (", dist_range, ")")
})

gps_info$cum_dist_km = dist_info

nightly_dist = lapply(gps, function(bat) {
    dists = sapply(split(bat, bat$time_cut), get_dist_km)
    data.frame(night = seq_along(dists), dists = dists, total_nights = length(dists))
})


dist_df = mapply(cbind, nightly_dist, "tagID" = names(nightly_dist),
                 SIMPLIFY = FALSE) %>%
    rbindlist()

summary(dist_df$dists)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3339  5.5329 12.7594 26.1401 41.4167 97.5681
night_dist = dist_df %>%
    mutate(dists = round(dists, 2)) %>%
    spread(night, dists) %>%
    arrange(tagID)


color_palette = c("#7F3C8D", "#3969AC", "#E68310",
                  "#E73F74", "#008695", "#CF1C90",
                  "#f97b72", "#4b4b8f")


dist_p = dist_df %>%
    mutate(tagID = paste0(tagID, " (", total_nights, " nights)")) %>%
    ggplot( aes(x = night, y = dists, color = tagID)) +
    geom_line(size = 0.5) +
    geom_point(size = 1) +
    ggtitle("Nightly cumulative distance flown by bats") +
    xlab("Foraging night") +
    ylab("Distance (km)") +
    scale_color_manual(values = color_palette) +
    theme_ipsum(base_size = 12, axis_title_size = 12, axis_text_size = 10) +
    theme(panel.grid.minor.x = element_blank()) +
    scale_x_continuous(limits = c(0, 7), breaks = 0:7) +
    labs(color = "Bat (nights tracked)")

night_dist
##   total_nights tagID     1     2     3     4     5     6
## 1            5 K5309 37.65 40.23 50.48 46.88 44.98    NA
## 2            5 K5310  2.07  5.39 96.46 97.57 77.66    NA
## 3            2 K5311 22.31 23.06    NA    NA    NA    NA
## 4            3 K5312  0.73  7.89 28.82    NA    NA    NA
## 5            6 K5313 14.64  8.87  7.93  8.21 10.88 22.91
## 6            1 K5315  2.92    NA    NA    NA    NA    NA
## 7            4 K5317  0.33  2.28  2.45 52.14    NA    NA
## 8            2 K5319  5.58 10.60    NA    NA    NA    NA
dist_p

write.csv(night_dist, file.path(gps_info_out_dir, "eobs_night-dist-info.csv"))
cowplot::save_plot(file.path(plot_out_dir, "line_nightly_dist.png"),
                   dist_p, base_height = 6)
                             

# ==============================================================================
# Add home range information
# ==============================================================================
gps_info = left_join(gps_info, hr, by = c("tagID" = "id"))
# gps_info


# ==============================================================================
# Add forage info
# ==============================================================================
gps_forage_n = gps %>%
    rbindlist() %>%
    dplyr::select(tagID, forage_pts) %>%
    filter(forage_pts) %>%
    group_by(tagID) %>%
    summarize(forage_n = n()) %>%
    as.data.frame()

gps_info = left_join(gps_info, gps_forage_n, by = "tagID")

kilombero_lonlat = data.frame(site = "kilombero",
                           lon = 277974.9,
                           lat = 9152539)

kilombero = st_as_sf(kilombero_lonlat, coords = c("lon", "lat"), crs = 32737)


forage_dist = lapply(gps, function(df, kilombero) {
    df_sf = df %>%
        st_as_sf(coords = c("utm.easting", "utm.northing"), crs = 32737)
    # first_loc = df_sf[1, ]
    forage_loc = df_sf[ df$forage_pts, ]
    dists = st_distance(kilombero, forage_loc)/1000
    dists = as.numeric(dists)
    summary(dists)
    data.frame(tagID = df$tagID[1],
               max_straightline_forage_dist_km = round(max(dists), 2),
               stringsAsFactors = FALSE)
}, kilombero) %>% rbindlist()

gps_info = left_join(gps_info, forage_dist, by = "tagID")


# ==============================================================================
# Add max dist between day roosts
# ==============================================================================

roost_dist = lapply(gps, function(df) {
    df_sf = df %>%
        st_as_sf(coords = c("utm.easting", "utm.northing"), crs = 32737)
    roost_loc = df_sf[ df$DayTime %in% "day" , ]
    dist_df = data.frame(tagID = df$tagID[1],
                         max_straightline_roost_dis_km = NA,
                         stringsAsFactors = FALSE)
    if(nrow(roost_loc) > 1){
        dists = st_distance(roost_loc, roost_loc)/1000
        dists = as.numeric(dists)
        dist_df$max_straightline_roost_dis_km = round(max(dists), 2)
    }
    dist_df
}) %>% rbindlist()

gps_info = left_join(gps_info, roost_dist, by = "tagID")
        
    
gps_info
##   tagID  location tag_type sex body_wt     date date_deployed  hours days
## 1 K5309 Kilombero    e-obs   M     284 11/14/16        Nov-14 106.64 4.44
## 2 K5310 Kilombero    e-obs   M     310 11/14/16        Nov-14 104.03 4.33
## 3 K5311 Kilombero    e-obs   M     272 11/14/16        Nov-14  34.28 1.43
## 4 K5312 Kilombero    e-obs   M     293 11/15/16        Nov-15  57.75 2.41
## 5 K5313 Kilombero    e-obs   M     271 11/14/16        Nov-14 129.11 5.38
## 6 K5315 Kilombero    e-obs   M     293 11/15/16        Nov-15  13.19 0.55
## 7 K5317 Kilombero    e-obs   F     271 11/15/16        Nov-15  84.96 3.54
## 8 K5319 Kilombero    e-obs   M     300 11/15/16        Nov-15  36.01 1.50
##   foraging_nights gps_fixes          cum_dist_km KDE_50_ha KDE_95_ha forage_n
## 1               5      9212 44.04 (37.65, 50.48)   3326.90  11743.12      150
## 2               5     10806  55.83 (2.07, 97.57)  33513.36 124413.78       84
## 3               2      1670 22.68 (22.31, 23.06)   1514.94   4861.94       25
## 4               3      1957  12.48 (0.73, 28.82)   1527.30   6497.41       45
## 5               6      3804  12.24 (7.93, 22.91)    142.25   1437.45      108
## 6               1       340    2.92 (2.92, 2.92)     23.62     96.88       12
## 7               4      2846   14.3 (0.33, 52.14)   4193.79  18985.86       48
## 8               2      1027    8.09 (5.58, 10.6)    118.16    469.22       33
##   max_straightline_forage_dist_km max_straightline_roost_dis_km
## 1                           19.83                         10.50
## 2                           62.63                         62.86
## 3                           11.03                          2.50
## 4                           11.17                          0.14
## 5                            6.37                          0.16
## 6                            0.97                          0.42
## 7                           16.68                          0.09
## 8                            2.84                          0.05
write.csv(gps_info, file.path(gps_info_out_dir, "eobs_gps-info.csv"))



Hourly distance travelled and heat maps

# ==============================================================================
# * Remove GPS points taken just at 7 am
# ==============================================================================
# These points would corespond to when the logger was getting a fix just
# before 7 am and the corresponding timestamp is milliseconds into 7 am
gps %>% rbindlist %>%
    mutate(gps_hour = hour(timestamp)) %>%
    filter(gps_hour %in% 7) %>%
    dplyr::select(tagID, timestamp)
##    tagID           timestamp
## 1  K5309 2016-11-16 07:00:00
## 2  K5309 2016-11-16 07:00:01
## 3  K5309 2016-11-16 07:00:02
## 4  K5309 2016-11-16 07:00:03
## 5  K5309 2016-11-16 07:00:04
## 6  K5309 2016-11-16 07:00:05
## 7  K5309 2016-11-16 07:00:41
## 8  K5309 2016-11-16 07:00:42
## 9  K5309 2016-11-16 07:00:43
## 10 K5309 2016-11-16 07:00:44
## 11 K5309 2016-11-16 07:00:45
## 12 K5309 2016-11-16 07:00:46
## 13 K5309 2016-11-16 07:00:47
## 14 K5309 2016-11-16 07:00:48
## 15 K5309 2016-11-16 07:00:49
## 16 K5309 2016-11-16 07:00:50
## 17 K5317 2016-11-18 07:00:47
## 18 K5317 2016-11-18 07:00:48
## 19 K5317 2016-11-18 07:00:49
## 20 K5317 2016-11-18 07:00:50
## 21 K5317 2016-11-18 07:00:51
## 22 K5317 2016-11-18 07:00:52
## 23 K5317 2016-11-18 07:00:53
## 24 K5317 2016-11-18 07:00:54
## 25 K5317 2016-11-18 07:00:55
## 26 K5317 2016-11-18 07:00:56
gps = lapply(gps, function(df) {
    gps_hour = hour(df$timestamp)
    df[ !gps_hour %in% 7, ]
})

# ==============================================================================
# Calculate hourly distances
# ==============================================================================
hourly_dist_info = lapply(gps, function(bat, lat, lon) {
    bat$night_hour = paste0(bat$time_cut, "_", hour(bat$timestamp))
    night_hour_dist = round(sapply(split(bat, bat$night_hour), get_dist_km, lon, lat), 2)
    hours = sapply(strsplit(names(night_hour_dist), "_"), `[[`, 2)
    hourly_dist = data.frame(hour = as.numeric(hours),
                             dist = as.numeric(night_hour_dist))
    hourly_dist %>%
        group_by(hour) %>%
        summarize(mean_dist = mean(dist),
                  min_dist = min(dist),
                  max_dist = max(dist)) %>%
        mutate_at(vars(mean_dist, min_dist, max_dist), round, 2) %>%
        as.data.frame()
}, lon = "utm.easting", lat = "utm.northing")

hourly_dist_info = Map(cbind, hourly_dist_info, tagID = names(gps)) %>%
    rbindlist

hour_df = data.frame(hour = c(17:23, 0:6))
hour_df$hour_order = factor(hour_df$hour, levels = rev(hour_df$hour))
hour_df = left_join(hour_df, hourly_dist_info, by = "hour")

hour_dist_hm = ggplot(data = hour_df,
                      mapping = aes(x = tagID, y = hour_order, fill = mean_dist)) +
    geom_tile() +
    xlab("Bat ID") +
    ylab("Hour") +
    labs(fill = "Mean distance \ntravelled (km)") +
    theme_ipsum(base_size = 14,
                plot_title_size = 14,
                axis_title_size = 14,
                grid = FALSE,
                axis = FALSE, ticks = TRUE) +
    theme(axis.title.y = element_text(angle = 0, vjust = 1)) +
    scale_fill_viridis(direction = -1)


hour_tbl = hour_df %>%
    mutate(dist_info = paste0(mean_dist, " (", min_dist, ", ", max_dist, ")")) %>%
    dplyr::select(hour_order, tagID, dist_info) %>%
    spread(key = tagID, value = dist_info) %>%
    arrange(desc(hour_order))

hour_tbl
##    hour_order               K5309              K5310              K5311
## 1          17         0 (0, 0.01)        0 (0, 0.01)           0 (0, 0)
## 2          18   2.34 (0.03, 4.28)  0.34 (0.01, 1.29)  0.19 (0.19, 0.19)
## 3          19 12.81 (0.05, 18.95)   15.62 (0, 21.51)  0.01 (0.01, 0.01)
## 4          20   0.19 (0.01, 0.55)     0.92 (0, 4.49)  0.09 (0.01, 0.17)
## 5          21  4.53 (0.03, 11.44)   0.1 (0.01, 0.29) 5.25 (0.04, 10.47)
## 6          22  3.49 (0.06, 11.55)     0.48 (0, 1.74)        0 (0, 0.01)
## 7          23  7.16 (0.05, 16.51)    3.13 (0, 10.25)       0.1 (0, 0.2)
## 8           0    2.5 (0.08, 8.73) 6.84 (0.06, 24.93) 5.66 (0.09, 11.22)
## 9           1   1.56 (0.01, 4.08) 8.29 (0.01, 38.08)        0 (0, 0.01)
## 10          2   0.87 (0.01, 2.44)  1.62 (0.01, 4.88)  0.01 (0.01, 0.01)
## 11          3   2.67 (0.02, 7.44)  6.1 (0.03, 30.37)  0.07 (0.06, 0.08)
## 12          4  5.51 (0.03, 15.31)     7.54 (0, 25.1)           0 (0, 0)
## 13          5      1.75 (0, 3.49)    7.12 (0, 25.52)  9.04 (7.1, 10.97)
## 14          6   0.04 (0.02, 0.09)  0.04 (0.01, 0.07)  1.55 (0.03, 3.08)
##                K5312             K5313             K5315              K5317
## 1  0.01 (0.01, 0.01)       0 (0, 0.01) 0.07 (0.07, 0.07)     0.01 (0, 0.03)
## 2   0.06 (0.02, 0.1) 1.85 (0.96, 3.35) 0.07 (0.07, 0.07)   0.45 (0.02, 0.9)
## 3   2.46 (0.01, 4.9) 2.54 (1.21, 4.82)          0 (0, 0)     0.02 (0, 0.05)
## 4  2.21 (1.74, 2.67)    0.16 (0, 0.86)          0 (0, 0)        0 (0, 0.01)
## 5  2.58 (0.03, 7.68) 0.68 (0.02, 3.72) 1.12 (1.12, 1.12)  0.05 (0.02, 0.12)
## 6        0 (0, 0.01)    0.72 (0, 4.21) 0.01 (0.01, 0.01)    4.24 (0, 16.94)
## 7   0.06 (0.01, 0.1)    0.78 (0, 3.94)          0 (0, 0)    3.67 (0, 14.64)
## 8  0.58 (0.09, 1.41) 0.21 (0.03, 0.81) 0.02 (0.02, 0.02)  0.08 (0.01, 0.23)
## 9     0.01 (0, 0.02) 0.16 (0.01, 0.82) 0.08 (0.08, 0.08)  0.02 (0.01, 0.03)
## 10 0.94 (0.01, 2.72)    0.32 (0, 0.98)          0 (0, 0)     0.42 (0, 1.68)
## 11  1.7 (0.01, 5.07) 1.19 (0.02, 3.57) 0.01 (0.01, 0.01) 4.14 (0.01, 16.45)
## 12    2.16 (0, 6.47)      0.7 (0, 3.4) 1.15 (1.15, 1.15)     0.41 (0, 1.04)
## 13    0.08 (0, 0.23)    2.26 (0, 8.01) 0.09 (0.09, 0.09)     0.22 (0, 0.84)
## 14 0.02 (0.01, 0.02) 0.06 (0.01, 0.25) 0.02 (0.02, 0.02)  0.02 (0.01, 0.03)
##                K5319
## 1  0.02 (0.02, 0.02)
## 2     0.02 (0, 0.03)
## 3  0.59 (0.01, 1.17)
## 4  1.74 (0.01, 3.47)
## 5  0.04 (0.02, 0.05)
## 6   0.11 (0.02, 0.2)
## 7  0.03 (0.01, 0.05)
## 8     1 (0.06, 1.95)
## 9  1.12 (0.03, 2.22)
## 10       0 (0, 0.01)
## 11 0.02 (0.02, 0.02)
## 12 0.76 (0.04, 1.47)
## 13 1.42 (0.96, 1.89)
## 14    0.02 (0, 0.04)
hour_dist_hm

# ==============================================================================
# Hourly distances travelled by bats by nights
# ==============================================================================
night_hourly_dist_info = lapply(gps, function(bat, lat, lon) {
    bat$night_hour = paste0(bat$time_cut, "_", hour(bat$timestamp))
    night_hour_dist = round(sapply(split(bat, bat$night_hour), get_dist_km, lon, lat), 2)
    nights = sapply(strsplit(names(night_hour_dist), "_"), `[[`, 1)
    hours = sapply(strsplit(names(night_hour_dist), "_"), `[[`, 2)
    hourly_dist = data.frame(night = nights,
                             hour = as.numeric(hours),
                             dist = as.numeric(night_hour_dist))
}, lon = "utm.easting", lat = "utm.northing")

night_hourly_dist_info = Map(cbind, night_hourly_dist_info, tagID = names(gps)) %>%
    rbindlist

night_hour_df = data.frame(hour = c(17:23, 0:6))
night_hour_df$hour_order = factor(night_hour_df$hour, levels = rev(night_hour_df$hour))
night_hour_df = left_join(night_hour_df, night_hourly_dist_info, by = "hour")


night_hour_p = ggplot(data = night_hour_df,
                      aes(x = night, y = hour_order, fill = dist)) +
    geom_tile() +
    xlab("Foraging night") +
    ylab("Hour") +
    ggtitle(paste0("Hourly distance travelled by bats for each tracked night")) +
    labs(fill = "Distance \ntravelled (km)") +
    theme_ipsum(base_size = 12,
                plot_title_size = 14,
                axis_title_size = 12,
                grid = FALSE,
                axis = FALSE, ticks = TRUE) +
    scale_fill_viridis(direction = -1) +
    facet_wrap( . ~ tagID, ncol = 4, scales = "free")


hour_summary = night_hour_df %>%
    group_by(tagID) %>%
    summarize(mean = mean(dist),
              min = min(dist),
              max = max(dist)) %>%
    mutate_at(vars(mean, min, max), round, 2)

hour_summary
## # A tibble: 8 x 4
##   tagID  mean   min   max
##   <fct> <dbl> <dbl> <dbl>
## 1 K5309  3.16     0 19.0 
## 2 K5310  4.12     0 38.1 
## 3 K5311  1.75     0 11.2 
## 4 K5312  0.89     0  7.68
## 5 K5313  0.81     0  8.01
## 6 K5315  0.19     0  1.15
## 7 K5317  0.98     0 16.9 
## 8 K5319  0.51     0  3.47
night_hour_p

# ==============================================================================
# Save heatmaps and nightly dist information
# ==============================================================================
write.csv(hour_tbl, file.path(gps_info_out_dir, "eobs_hour-dist-info.csv"))
write.csv(hour_summary, file.path(gps_info_out_dir, "eobs_hour-dist-summary.csv"))
          
hm_fn = file.path(plot_out_dir, "heatmap_hour-dist.png")
hm_fn_pdf = file.path(plot_out_dir, "heatmap_hour-dist.pdf")
cowplot::save_plot(hm_fn, hour_dist_hm, base_height = 7, dpi = 600)
cowplot::save_plot(hm_fn_pdf, hour_dist_hm, base_height = 7, dpi = 600)

hm_night_fn = file.path(plot_out_dir, "heatmap_hour-dist-night.png")
cowplot::save_plot(hm_night_fn, night_hour_p, base_height = 10, dpi = 600)