Identify roosts and further subset not-flying points for foraging/feeding roost locations
Setting up the workspace
# ==============================================================================
# * Libraries
# ==============================================================================
library(move)
library(maptools)
library(sp)
library(rgdal)
library(rgeos)
library(lubridate)
library(dplyr)
library(data.table)
library(leaflet)
library(htmlwidgets)
library(adehabitatHR)
library(sf)
library(ggplot2)
library(ggsn)
library(hrbrthemes)
library(DT)
library(tidyr)
library(ggmap)
library(patchwork)
library(cowplot)
source("R/00_fxns.R")
gps_info_out_dir = file.path("results", "tables")
plot_out_dir = "results/maps/gps/eobs"
plot_facet = file.path(plot_out_dir, "facet")
if( !dir.exists(plot_facet)) dir.create(plot_facet, recursive = TRUE)
# ==============================================================================
# * Read in GPS data with flying behaviour + tagging location
# ==============================================================================
# read in the behavior-classified acc and GPS data
acc_gps_behavL = readRDS("data/acc_behavL/yz_acc_gps_behavL.RDS")
gps = lapply(acc_gps_behavL, `[[`, "gps_behav")
gps = lapply(gps, break_nights_fxn)Find day roosts
# ==============================================================================
# * Find day roosts
# ==============================================================================
# Kilombero roost where tagging took place
kilombero_colony = data.frame(xlim = c(36.98629, 36.99085),
ylim = c(-7.66273, -7.66052))
# Determine day roosts as those gps fixes collected during the day time
get_roosts = function(df, kilombero_colony) {
is_roost = ifelse(df$DayTime %in% "day" & df$activity %in% "NotFlying",
TRUE, FALSE)
xlim = kilombero_colony$xlim
ylim = kilombero_colony$ylim
is_kilombero = df$location.long > xlim[1] &
df$location.long < xlim[2] &
df$location.lat > ylim[1] &
df$location.lat < ylim[2]
new_roost = is_roost & !is_kilombero
kilombero_roost = is_roost & is_kilombero
df$roost[ new_roost ] = "new"
df$roost[ kilombero_roost ] = "kilombero"
df$roost_date = date(df$timestamp)
df
}
gps = lapply(gps, get_roosts, kilombero_colony)
# lapply(gps, function(df) table(df$roost))
gps_lines = lapply(gps, make_lines)
if(FALSE) {
# call and save maps
htmltools::tagList(
lapply(names(gps), map_gps, gps, gps_lines, map_roosts = TRUE)
)
}
roost_info = lapply(gps, function(df) {
cat("*** ", df$tagID[1], " ***\n")
info = df %>%
filter(roost == "new") %>%
dplyr::select(tagID, roost, roost_date) %>%
group_by(tagID, roost_date) %>%
summarize(n = n()) %>%
as.data.frame()
}) %>% rbindlist()## *** K5309 ***
## *** K5310 ***
## *** K5311 ***
## *** K5312 ***
## *** K5313 ***
## *** K5315 ***
## *** K5317 ***
## *** K5319 ***roost_info## tagID roost_date n
## 1: K5309 2016-11-15 30
## 2: K5310 2016-11-17 30
## 3: K5311 2016-11-15 55Further subset of “not flying” points
# ** Get time differences
gps_tdiff = lapply(gps, function(df) {
bat_nights = split(df, df$time_cut)
time_diff = lapply(bat_nights, function(df_night) {
df_night$time_diff = difftime(lead(df_night$timestamp),
df_night$timestamp, units = "sec")
df_night$time_diff = as.numeric(df_night$time_diff)
df_night
})
time_diff = rbindlist(time_diff)
as.data.frame(time_diff)
})
# ** Forage points
# ==============================================================================
find_forage_pts = function(df, min_diff = 30, kilombero_colony) {
avg_speed = mean(df$ground.speed)
# avg_speed = mean(df$ground.speed[ df$activity %in% "Flying"])
sep_m = (avg_speed * min_diff)
dist = spDists(as.matrix(df[ , c("utm.easting", "utm.northing")]),
longlat = FALSE)
colnames(dist) = 1:ncol(dist)
rownames(dist) = 1:nrow(dist)
coloc = data.table(which(dist <= sep_m, arr.ind=T))
colnames(coloc) = c("pt1", "pt2")
coloc_diff = coloc$pt2 - coloc$pt1
coloc_consec = which(coloc_diff == 1)
coloc_raw = coloc[ coloc_consec, ]
coloc_pts = unique(unlist(coloc_raw[ , c(1:2)], use.names = FALSE))
coloc_time = df$time_diff[ coloc_pts ]
# coloc_time[ is.na(coloc_time) ] = 9000 # spuriously large number
coloc_keep = coloc_pts[ coloc_time >= min_diff ]
df$forage_pts = FALSE
df$forage_pts[ coloc_keep ] = TRUE
df$forage_pts[ df$activity %in% "Flying" ] = FALSE
df$forage_pts[ df$DayTime %in% "day" ] = FALSE
xlim = kilombero_colony$xlim
ylim = kilombero_colony$ylim
is_kilombero = df$location.long > xlim[1] &
df$location.long < xlim[2] &
df$location.lat > ylim[1] &
df$location.lat < ylim[2]
df$forage_pts[ is_kilombero ] = FALSE
# Include distances
df$dist = NA
for(i in 1:(nrow(dist) - 1)) {
df$dist[i + 1] = dist[i, i + 1]
}
# Aggregating subsequent foraging points to one group
x = df$forage_pts
y = split(x, df$time_cut)
z = lapply(y, function(v) {
c(head(v, 1), tail(v, -1) - head(v, -1) == 1)
})
z = unlist(z)
names(z) = NULL
z = cumsum(c(head(z, 1), tail(z, -1) - head(z, -1) == 1))
df$forage_cluster = ifelse(x, z, NA)
forage_time = aggregate(df$time_diff,
by = list(df$forage_cluster),
FUN = sum, na.rm = TRUE)
names(forage_time) = c("forage_cluster", "forage_time")
df = dplyr::left_join(df, forage_time, by = "forage_cluster")
}
gps_forage = lapply(gps_tdiff, find_forage_pts, min_diff = 30, kilombero_colony)
# lapply(gps_forage, function(df) table(df$forage_pts))
# lapply(gps_forage, function(df) table(df$forage_pts, df$activity))
forage_info = rbindlist(gps_forage) %>%
dplyr::select(tagID, forage_pts) %>%
filter(forage_pts) %>%
group_by(tagID) %>%
summarize(n = n()) %>%
as.data.frame()
forage_info_night = rbindlist(gps_forage) %>%
dplyr::select(tagID, forage_pts, time_cut) %>%
filter(forage_pts) %>%
group_by(tagID, time_cut) %>%
summarize(n = n()) %>%
as.data.frame() %>%
spread(time_cut, n)
if(FALSE){
htmltools::tagList(
lapply(names(gps_forage), map_gps, gps_forage, gps_lines,
map_roosts = TRUE, map_forage = TRUE, map_notflying = TRUE)
)
}
saveRDS(gps_forage, "data/gps_forage.RDS")
gps_forage_df = gps_forage %>%
rbindlist()
write.csv(gps_forage_df, "data/gps_forage.csv")Static maps of gps points
# Map gps points
# ------------------------------------------------------------------------------
kilombero_lonlat = data.frame(site = "kilombero",
lon = 277974.9,
lat = 9152539)
kilombero_wgs = data.frame(site = "kilombero",
lon = 36.98726,
lat = -7.662083)
kilombero = st_as_sf(kilombero_lonlat, coords = c("lon", "lat"), crs = 32737)
tagIDs = names(gps_forage)
color_palette = c("#7F3C8D", "#3969AC", "#E68310",
"#E73F74", "#008695", "#CF1C90",
"#f97b72", "#4b4b8f")
tag_colors = color_palette[ 1:length(tagIDs) ]
names(tag_colors) = tagIDs
# ID = tagIDs[1]
small_area_IDs = c("K5313", "K5315", "K5319")
large_area_IDs = "K5310"
medium_area_IDs = tagIDs[ ! tagIDs %in% c(small_area_IDs, large_area_IDs) ]
area_IDs = list(small = small_area_IDs,
medium = medium_area_IDs,
large = large_area_IDs)
bats_sf = gps_forage %>%
rbindlist() %>%
st_as_sf(coords = c("utm.easting", "utm.northing"), crs = 32737)
map_extents = lapply(area_IDs, function(IDs, bats_sf) {
bats = bats_sf[ bats_sf$tagID %in% IDs, ]
bats_sp = as(bats, 'Spatial')
map_extent = extent(bats_sp) %>%
as('SpatialPolygons') %>%
st_as_sf %>%
st_set_crs(st_crs(bats_sp)) %>%
st_buffer(3000)
}, bats_sf)
names(map_extents) = c("small", "medium", "large")
basemaps = lapply(map_extents, function(map_extent) {
extent = map_extent %>%
st_transform(crs = 4326) %>%
st_bbox()
get_map(as.vector(extent),
zoom=13,
maptype = 'satellite', source = 'google')
})
p_sm = lapply(small_area_IDs, plot_gps,
"small", gps_forage, tag_colors, basemaps,
kilombero_wgs,
area_IDs,
plot_out_dir)
p_med = lapply(medium_area_IDs, plot_gps,
"medium", gps_forage, tag_colors, basemaps,
kilombero_wgs,
area_IDs,
plot_out_dir)
p_large = lapply(large_area_IDs, plot_gps,
"large", gps_forage, tag_colors, basemaps,
kilombero_wgs,
area_IDs,
plot_out_dir)
# ** All plots
ind_bat = gps_forage %>% rbindlist
bats_sp = as(bats_sf, 'Spatial')
all_map_extent = extent(bats_sp) %>%
as('SpatialPolygons') %>%
st_as_sf %>%
st_set_crs(st_crs(bats_sp)) %>%
st_buffer(5000)
all_extent = all_map_extent %>%
st_transform(crs = 4326) %>%
st_bbox()
all_basemap = get_map(as.vector(all_extent),
zoom=13,
maptype = 'satellite', source = 'google')
sc_info = extent(st_transform(bats_sf, crs = 4326))
g = ggmap(all_basemap) +
geom_point(data = kilombero_wgs, aes(x = lon, y = lat),
color = "#08B8FF", size = 4, alpha = 0.65) +
geom_sf(data = st_transform(map_extents$small, crs = 4326),
fill = "#6B8278", color = "#6B8278", size = 0.25, alpha = 0.5,
inherit.aes = FALSE) +
geom_sf(data = st_transform(map_extents$medium, crs = 4326),
fill = "#6B8278", color = "#6B8278", size = 0.25, alpha = 0.5,
inherit.aes = FALSE) +
geom_path(data = ind_bat, aes(x = location.long, y = location.lat, color = tagID),
size = 0.75) +
scale_color_manual(values = tag_colors) +
theme_ipsum(base_size = 10,
base_family = "Helvetica",
plot_title_size = 12,
grid = FALSE,
axis = FALSE, ticks = TRUE) +
theme(
plot.background = element_rect(fill = '#D1D6AB', colour = '#D1D6AB'),
legend.position="none",
plot.margin=grid::unit(c(1,1,1,1), "mm")
)
g_scale =
g +
ggsn::scalebar(
x.min = sc_info[1] - 0.001,
x.max = sc_info[2] + 0.001,
y.min = sc_info[3] - 0.001,
y.max = sc_info[4] + 0.001,
dist = 4,
dist_unit = "km",
transform = TRUE,
model = "WGS84",
location = "bottomright",
height = 0.01,
st.dist = 0.025,
st.size = 4,
box.fill = c("white", "black"),
box.color = "black",
border.size = 0.25) +
xlab("Longitude") +
ylab("Latitude") +
ggtitle("")
fn = file.path(plot_out_dir, paste0("all_tags", ".png"))
cowplot::save_plot(fn, g_scale, base_height = 5)
g_scale
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] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ggmap_3.0.0 rgeos_0.4-3 adehabitatHR_0.4.16
## [4] adehabitatLT_0.3.24 CircStats_0.2-6 boot_1.3-22
## [7] MASS_7.3-51.4 adehabitatMA_0.3.13 ade4_1.7-13
## [10] deldir_0.1-21 finalfit_1.0.0 kableExtra_1.1.0
## [13] knitr_1.26 maptools_0.9-5 cowplot_1.0.0
## [16] readxl_1.3.1 tidyr_1.0.2 viridis_0.5.1
## [19] viridisLite_0.3.0 DT_0.11.5 hrbrthemes_0.6.0
## [22] ggsn_0.5.0 patchwork_1.0.0 ggplot2_3.3.0.9000
## [25] sf_0.8-0 htmlwidgets_1.5.1 leaflet_2.0.2
## [28] data.table_1.12.8 dplyr_0.8.4 lubridate_1.7.4
## [31] move_3.2.0 rgdal_1.4-4 raster_2.9-5
## [34] sp_1.3-1 geosphere_1.5-10
##
## loaded via a namespace (and not attached):
## [1] minqa_1.2.4 colorspace_1.4-1 rjson_0.2.20 ellipsis_0.3.0
## [5] class_7.3-15 rstudioapi_0.11 mice_3.5.0 farver_2.0.3
## [9] xml2_1.2.2 codetools_0.2-16 splines_3.6.0 extrafont_0.17
## [13] jsonlite_1.6.1 nloptr_1.2.1 broom_0.5.2 Rttf2pt1_1.3.7
## [17] png_0.1-7 shiny_1.4.0 readr_1.3.1 compiler_3.6.0
## [21] httr_1.4.1 backports_1.1.5 assertthat_0.2.1 Matrix_1.2-17
## [25] fastmap_1.0.1 later_1.0.0 htmltools_0.4.0 tools_3.6.0
## [29] gtable_0.3.0 glue_1.3.1 ggthemes_4.2.0 Rcpp_1.0.3
## [33] cellranger_1.1.0 vctrs_0.2.2 nlme_3.1-139 extrafontdb_1.0
## [37] crosstalk_1.0.0 xfun_0.12 stringr_1.4.0 lme4_1.1-21
## [41] rvest_0.3.4 mime_0.9 lifecycle_0.1.0 pan_1.6
## [45] scales_1.1.0 hms_0.5.2 promises_1.1.0 parallel_3.6.0
## [49] curl_4.3 yaml_2.2.1 memoise_1.1.0 gridExtra_2.3
## [53] gdtools_0.2.1 rpart_4.1-15 stringi_1.4.5 highr_0.8
## [57] e1071_1.7-2 RgoogleMaps_1.4.3 rlang_0.4.4 pkgconfig_2.0.3
## [61] systemfonts_0.1.1 matrixStats_0.54.0 bitops_1.0-6 evaluate_0.14
## [65] lattice_0.20-38 purrr_0.3.3 labeling_0.3 tidyselect_1.0.0
## [69] plyr_1.8.5 magrittr_1.5 R6_2.4.1 generics_0.0.2
## [73] mitml_0.3-7 DBI_1.0.0 pillar_1.4.3 foreign_0.8-71
## [77] withr_2.1.2 units_0.6-3 survival_2.44-1.1 nnet_7.3-12
## [81] tibble_2.1.3 crayon_1.3.4 jomo_2.6-8 KernSmooth_2.23-15
## [85] rmarkdown_1.17.1 jpeg_0.1-8.1 forcats_0.4.0 digest_0.6.23
## [89] classInt_0.4-2 webshot_0.5.2 xtable_1.8-4 httpuv_1.5.2
## [93] munsell_0.5.0