Setting up the workspace

# ==============================================================================
# * 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)))



Kernel Density Estimation

# ==============================================================================
# * 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 50% and 95% KDE countours

# ** 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)

Bats with relatively smaller 95% KDE areas (96.88 to 1437.45 ha)

Bats with relatively medium 95% KDE areas (4861.94 to 18985.86 ha)

Bat with relatively larger 95% KDE area (124413.78 ha)



Save KDE area information and KDE plots

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)



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] 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