# ==============================================================================
# * Libraries
# ==============================================================================
library(move)
library(adehabitatHR)
library(dplyr)
library(data.table)
library(sf)
library(ggplot2)
library(ggsn)
library(hrbrthemes)
library(patchwork)
source("R/00_fxns.R")
# Directory for saving plots and KDE information
plot_out_dir = file.path("results", "maps", "kde")
kde_info_out_dir = file.path("results", "tables")
invisible(lapply(c(plot_out_dir, kde_info_out_dir), function(out_dir)
if( !dir.exists(out_dir)) dir.create(out_dir, recursive = TRUE)))
# ==============================================================================
# * Read in the behavior-classified gps data
# ==============================================================================
acc_gps_behavL = readRDS("data/acc_behavL/yz_acc_gps_behavL.RDS")
kilombero_lonlat = data.frame(site = "kilombero",
lon = 277974.9,
lat = 9152539)
kilombero = st_as_sf(kilombero_lonlat, coords = c("lon", "lat"), crs = 32737)
gpsALL = lapply(acc_gps_behavL, `[[`, "gps_behav") %>%
rbindlist() %>%
as.data.frame()
# Subset "not flying" gps fixes
# ------------------------------------------------------------------------------
gpsALl = gpsALL[ gpsALL$activity %in% "NotFlying", ]
# ==============================================================================
# * Home range Estimation (using KDE method)
# ==============================================================================
bats_df = gpsALL[, c("tagID", "utm.easting", "utm.northing")]
bats_sf = st_as_sf(bats_df, coords = c("utm.easting", "utm.northing"), crs = 32737)
bats_sp = as(bats_sf, 'Spatial')
bats_ud <- kernelUD(bats_sp, h = "href") # href = the reference bandwidth
if(FALSE) {
mcp95 <- mcp(bats_sp[ ,"tagID"], percent = 95)
mcp50 <- mcp(bats_sp[ ,"tagID"], percent = 50)
}
ud50 = getverticeshr(bats_ud, percent = 50)
ud95 = getverticeshr(bats_ud, percent = 95)
# Make data frames for ggplot
ud50_df = fortify(ud50)
ud95_df = fortify(ud95)
# Split to make individual plots
ud50L = split(ud50_df, ud50_df$id)
ud95L = split(ud95_df, ud95_df$id)
# KDE areas (in ha)
ud50_area = ud50 %>%
as.data.frame() %>%
mutate(area = round(area, 2)) %>%
rename(KDE_50_ha = area)
ud95_area = ud95 %>%
as.data.frame() %>%
mutate(area = round(area, 2)) %>%
rename(KDE_95_ha = area)
ud_area = dplyr::left_join(ud50_area, ud95_area, by = "id")
DT::datatable(ud_area, options = list(dom = 't'))
# ** Plot KDE
# ==============================================================================
tagIDs = unique(bats_sp$tagID)
color_palette = c("#E58606", "#5D69B1", "#52BCA3", "#99C945",
"#CC61B0", "#24796C", "#DAA51B", "#2F8AC4",
"#764E9F", "#ED645A", "#CC3A8E", "#A5AA99")
tag_colors = color_palette[ 1:length(tagIDs) ]
names(tag_colors) = tagIDs
small_area_IDs = c("K5313", "K5315", "K5319")
large_area_IDs = "K5310"
medium_area_IDs = tagIDs[ ! tagIDs %in% c(small_area_IDs, large_area_IDs) ]
small_extent = extent(bats_sp[ bats_sp$tagID %in% small_area_IDs, ]) %>%
as('SpatialPolygons') %>%
st_as_sf %>%
st_set_crs(st_crs(bats_sp)) %>%
st_buffer(1500)
medium_extent = extent(bats_sp[ bats_sp$tagID %in% medium_area_IDs, ]) %>%
as('SpatialPolygons') %>%
st_as_sf %>%
st_set_crs(st_crs(bats_sp)) %>%
st_buffer(1000)
large_extent = extent(bats_sp[ bats_sp$tagID %in% large_area_IDs, ]) %>%
as('SpatialPolygons') %>%
st_as_sf %>%
st_set_crs(st_crs(bats_sp)) %>%
st_buffer(100)
plot_kde = function(ID, bats_df, ud50L, ud95L, tag_colors, map_extent, kilombero) {
ind_bat = bats_df[ bats_df$tagID %in% ID, ]
ind_ud50 = ud50L[[ ID ]]
ind_ud95 = ud95L[[ ID ]]
col = tag_colors[ ID ]
g = ggplot() +
geom_sf(data = map_extent, fill = NA, color = NA) +
geom_polygon(data = ind_ud50, aes(x = long, y = lat),
alpha = 0.6, fill = col) +
geom_polygon(data = ind_ud95, aes(x = long, y = lat),
alpha = 0.4, fill = col) +
geom_path(data = ind_bat, aes(x = utm.easting, y = utm.northing),
color = "black", size = .05) +
geom_point(data = ind_bat, aes(x = utm.easting, y = utm.northing),
color = "black", size = 0.1, alpha = 0.3) +
geom_sf(data = kilombero, color = "white", size = 2, aes(shape = site)) +
scale_shape_manual(values = 8) +
ggthemes::scale_fill_gdocs() +
theme_ipsum(base_size = 12,
plot_title_size = 14,
axis_title_size = 12,
axis_text_size = 10,
grid = FALSE,
axis = FALSE, ticks = TRUE) +
theme(legend.position="none") +
coord_sf() +
scalebar(map_extent, dist = 4, dist_unit = "km",
location = "topleft",
height = 0.01,
st.dist = 0.05,
st.size = 3,
box.fill = c("grey", "black"),
box.color = NA,
transform = FALSE) +
scale_x_continuous(
labels=label_fill) +
scale_y_continuous(
labels=label_fill) +
xlab("Longitude") +
ylab("Latitude") +
ggtitle(ID)
g
}
p_sm = lapply(small_area_IDs, plot_kde,
bats_df, ud50L, ud95L, tag_colors, map_extent = small_extent, kilombero)
names(p_sm) = small_area_IDs
p_sm_all = (p_sm$K5313 | p_sm$K5315 | p_sm$K5319)
p_med = lapply(medium_area_IDs, plot_kde,
bats_df, ud50L, ud95L, tag_colors, map_extent = medium_extent, kilombero)
names(p_med) = medium_area_IDs
p_med_all = (p_med$K5309 | p_med$K5311) / (p_med$K5312 | p_med$K5317)
p_large = plot_kde(large_area_IDs, bats_df, ud50L, ud95L,
tag_colors, map_extent = large_extent, kilombero)
write.csv(ud_area, file.path(kde_info_out_dir, "KDE-area_eobs.csv"))
cowplot::save_plot(file.path(plot_out_dir, "KDE_small_area.png"),
p_sm_all, base_width = 12, base_height = 8, dpi = 600)
cowplot::save_plot(file.path(plot_out_dir, "KDE_med_area.png"),
p_med_all, base_width = 12, base_height = 8, dpi = 600)
cowplot::save_plot(file.path(plot_out_dir, "KDE_large_area.png"),
p_large, base_width = 12, base_height = 8, dpi = 600)
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] adehabitatHR_0.4.16 adehabitatLT_0.3.24 CircStats_0.2-6
## [4] boot_1.3-22 MASS_7.3-51.4 adehabitatMA_0.3.13
## [7] ade4_1.7-13 deldir_0.1-21 finalfit_1.0.0
## [10] kableExtra_1.1.0 knitr_1.26 maptools_0.9-5
## [13] cowplot_1.0.0 readxl_1.3.1 tidyr_1.0.2
## [16] viridis_0.5.1 viridisLite_0.3.0 DT_0.11.5
## [19] hrbrthemes_0.6.0 ggsn_0.5.0 patchwork_1.0.0
## [22] ggplot2_3.3.0.9000 sf_0.8-0 htmlwidgets_1.5.1
## [25] leaflet_2.0.2 data.table_1.12.8 dplyr_0.8.4
## [28] lubridate_1.7.4 move_3.2.0 rgdal_1.4-4
## [31] raster_2.9-5 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] ggmap_3.0.0 gtable_0.3.0 glue_1.3.1 ggthemes_4.2.0
## [33] Rcpp_1.0.3 cellranger_1.1.0 vctrs_0.2.2 nlme_3.1-139
## [37] extrafontdb_1.0 crosstalk_1.0.0 xfun_0.12 stringr_1.4.0
## [41] lme4_1.1-21 rvest_0.3.4 mime_0.9 lifecycle_0.1.0
## [45] pan_1.6 scales_1.1.0 hms_0.5.2 promises_1.1.0
## [49] parallel_3.6.0 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