# ==============================================================================
# * Reading in and sorting Lotek data
# ==============================================================================
# ** Importing lotek data
# ------------------------------------------------------------------------------
bats = read.csv("data/raw/bats/lotek_info.csv", stringsAsFactors = FALSE)
bats$argosID = as.character(bats$argosID)
bats$tagID = paste0(ifelse(bats$location %in% "Morogoro", "M", "K"), bats$argosID)
lotek = list.files("data/raw/lotek", full.names = TRUE)
lotek_list = lapply(lotek, function(x) {
fname = strsplit(x, "_")[[1]][2]
df = read.csv(x, stringsAsFactors = FALSE)
df$argosID = fname
df = df[ , c(ncol(df), 1:(ncol(df) - 1))]
df
})
gps = do.call(rbind, lotek_list)
# ** Subsetting valid rows, and sorting by lat long
# ------------------------------------------------------------------------------
gps = gps[ gps$CRC != "Fail" & gps$Fix %in% c("3D", "2D"), ]
# gps = gps[ order(-gps[ , "Latitude"], gps[ , "Longitude"]), ]
gps = gps[ gps$Latitude > -20, ] # there's one fix with a latitude of -63.025063
drop_cols = c("Fix", "Act1", "Act2", "Act3", "Act4")
gps = gps[ , ! names(gps) %in% drop_cols]
rownames(gps) = NULL
# ** Creating timestamp for Africa/Dar_es_Salaam
# ------------------------------------------------------------------------------
gps$lotek_timestamp = as.POSIXct(paste(gps$Date, gps$Time),
format="%d/%m/%Y %H:%M:%S", tz = "UTC")
gps$timestamp = with_tz(gps$lotek_timestamp, "Africa/Dar_es_Salaam")
attr(gps$timestamp, "tzone")
# ** Combining lotek data with bat-specific info
# ==============================================================================
gps = dplyr::left_join(gps, bats, "argosID")
gps$argosID = NULL
bats$argosID = NULL
# ** Create local (Tanzania timezone) timestamp
# ------------------------------------------------------------------------------
gps$Date = NULL
gps$date = date(gps$timestamp)
gps$day = day(gps$timestamp)
gps$month = month(gps$timestamp)
gps$year = year(gps$timestamp)
# * Break the GPS data by foraging nights
# ==============================================================================
gps_l = split(gps, gps$tagID)
break_nights_fxn = function(df, write = FALSE) {
dates = sort(unique(df$date))
date_seq = seq(min(dates) - 1, max(dates) + 1, by = "days")
if (length(unique(date_seq)) >= 3) {
time_cut = as.numeric(cut(df$timestamp,
breaks = as.POSIXct(paste0(date_seq, " 17:00:00"),
tz = "Africa/Dar_es_Salaam")))
if(min(time_cut) > 1) {time_cut = time_cut - 1}
df$time_cut = time_cut
} else {
df$time_cut = "<1"
}
df
}
gps_l = lapply(gps_l, break_nights_fxn)
# lapply(gps_l, function(df) table(df$time_cut))
# * Get GPS time info - min, max, length and units (day or hours)
# ==============================================================================
gps_times = lapply(gps_l, function(df) {
tmax = max(df$timestamp)
tmin = min(df$timestamp)
diff = difftime(tmax, tmin, units = "days")
tag_nights = df$time_cut
data.frame(tagID = df$tagID[1],
tmin = tmin,
tmax = tmax,
length = as.numeric(diff),
units = attributes(diff)$units,
nights = length(unique(tag_nights)),
stringsAsFactors = FALSE)
}) %>% rbindlist()
# gps_times = gps_times[ , c(ncol(gps_times), 1:(ncol(gps_times) - 1))]
# ** add GPS times info to "bat_info" df for article
bat_info = left_join(bats, gps_times, by = "tagID")
bat_info = bat_info[ !is.na(bat_info$nights), ]
# * Getting distances travelled
# ==============================================================================
get_dist_km = function(bat) {
df = as.data.frame(bat)
coordinates(df) = c("Longitude", "Latitude")
pp = list(Lines(list(Line(coordinates(df))), ID = "w"))
ps = SpatialLines(pp)
# plot(ps)
data.frame(tagID = bat$tagID[1],
dist = SpatialLinesLengths(ps, longlat = TRUE),
stringsAsFactors = FALSE)
}
# ** add distance travelled to "bat_info" df for article
dists = lapply(gps_l, get_dist_km) %>% rbindlist()
bat_info = left_join(bat_info, dists, by = "tagID")
# ==============================================================================
# * Cumulative distances
# ==============================================================================
get_dist_km = function(df, lon = "Longitude", lat = "Latitude") {
coordinates(df) = c(lon, lat)
pp = list(Lines(list(Line(coordinates(df))), ID = "w"))
ps = SpatialLines(pp)
# plot(ps)
SpatialLinesLengths(ps, longlat = TRUE)
}
dist_info = lapply(gps_l, 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 = ", ")
cumulative_dist = paste0(dist_mean, " (", dist_range, ")")
data.frame(tagID = bat$tagID[1], cumulative_dist_km = cumulative_dist,
stringsAsFactors = FALSE)
}) %>% rbindlist()
bat_info = left_join(bat_info, dist_info, by = "tagID")
cumulative_dist = sapply(gps_l, function(bat) {
nightly_dist = sapply(split(bat, bat$time_cut), get_dist_km)
}) %>% unlist()
# * Determining foraging and roosting points
# ==============================================================================
# We can't determine them from this data, so assigning all to NA
gps_l = lapply(gps_l, function(df) {
df$forage_pts = NA
xlim = c(37.66670, 37.66755)
ylim = c(-6.82482 ,-6.82405)
df$roost[ df$Longitude > xlim[1] &
df$Longitude < xlim[2] &
df$Latitude > ylim[1] &
df$Latitude < ylim[2] ] = TRUE
df$roost_loc = ifelse(df$roost, "nunge", NA)
forage_nunge = which(df$roost_loc %in% "nunge")
df$forage_pts[forage_nunge] = FALSE
df
})
# * Bat table for publication
# ==============================================================================
# lapply(gps_l,function(df) unique(df$tagID))
gps_fixes = lapply(gps_l, function(df)
data.frame(tagID = df$tagID[1], gps_n = nrow(df),
stringsAsFactors = FALSE)) %>% rbindlist()
bat_info = left_join(bat_info, gps_fixes, by = "tagID")
bat_info$foraging_pts = NA
bat_info$new_roosts = 0
# ==============================================================================
# * Home range Estimation (using KDE method)
# ==============================================================================
gpsALL = gps_l %>% rbindlist()
bats_df = gpsALL[, c("tagID", "Longitude", "Latitude")]
fixes_5 = bat_info$tagID[ bat_info$gps_n >= 5 ]
bats_df = bats_df[ bats_df$tagID %in% fixes_5, ]
bats_sf = st_as_sf(bats_df, coords = c("Longitude", "Latitude"), crs = 4326)
bats_sf = bats_sf %>%
st_transform(crs = 32737)
bats_sp = as(bats_sf, 'Spatial')
bats_ud <- kernelUD(bats_sp, h = "href", extent = 2.5)
ud50 = getverticeshr(bats_ud, percent = 50)
ud95 = getverticeshr(bats_ud, percent = 95)
# 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")
bat_info = dplyr::left_join(bat_info, ud_area, by = c("tagID" = "id"))
# ==============================================================================
# * Export results
# ==============================================================================
gps_export = do.call(rbind, gps_l)
row.names(gps_export) = NULL
saveRDS(gps_export, "data/gps_lotek.RDS")
write.csv(bat_info, file.path(gps_info_out_dir, "lotek_gps-info.csv"))
# ==============================================================================
# * Mapping
# ==============================================================================
dfm = gps_l %>%
rbindlist() %>%
filter(location %in% "Morogoro")
dfk = gps_l %>%
rbindlist() %>%
filter(location %in% "Kilombero")
# col_pal = c("#E58606", "#5D69B1", "#52BCA3", "#99C945", "#CC61B0", "#24796C", "#DAA51B", "#2F8AC4", "#764E9F", "#ED645A", "#CC3A8E", "#A5AA99")
col_pal = c("#7F3C8D", "#11A579", "#3969AC", "#F2B701", "#E73F74", "#80BA5A", "#E68310", "#008695", "#CF1C90", "#f97b72", "#4b4b8f", "#A5AA99")
plot_lotek = function(df, scale_location, location, col_pal = col_pal, plot_out_dir) {
df_sf = df %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326)
df_sp = as(st_transform(df_sf, crs = 32737), 'Spatial')
map_extent = extent(df_sp) %>%
as('SpatialPolygons') %>%
st_as_sf %>%
st_set_crs(st_crs(df_sp)) %>%
st_buffer(1000) %>%
st_transform(crs = 4326)
extent = map_extent %>%
st_bbox()
basemap = get_map(as.vector(extent),
zoom=13,
maptype = 'satellite', source = 'google')
g = ggmap(basemap, darken = 0) +
geom_path(data = df, aes(x = Longitude, y = Latitude, color = tagID),
# color = col,
size = .5) +
geom_point(data = df, aes(x = Longitude, y = Latitude, color = tagID),
# color = col,
size = 1.5, alpha = 1) +
# geom_sf(data = kilombero, color = "blue", size = 2) +
# geom_sf(data = kilombero, color = "white", size = 1.5, aes(shape = site)) +
scale_shape_manual(values = 8) +
scale_colour_manual(values = col_pal) +
ggthemes::scale_fill_gdocs() +
theme_ipsum(base_size = 12,
plot_title_size = 12,
axis_title_size = 12,
grid = FALSE,
axis = FALSE, ticks = TRUE) +
theme(legend.position="top",
# plot.margin=grid::unit(c(0,0,0,0), "mm")
) +
coord_sf(crs = st_crs(4326))
g = g + ggsn::scalebar(df_sf,
dist = 4,
dist_unit = "km",
transform = TRUE,
model = "WGS84",
location = scale_location,
height = 0.01,
st.dist = 0.025,
st.size = 4,
box.fill = c("grey", "black"),
box.color = NA) +
labs(color = paste0(location, " bats \n(Argos satellite tags)")) +
xlab("Longitude") +
ylab("Latitude")
fn = file.path(plot_out_dir, paste0("lotek_", location, ".png"))
ggsave(fn, g, dpi = 500)
fn_pdf = file.path(plot_out_dir, paste0("lotek_", location, ".pdf"))
ggsave(fn_pdf, g, dpi = 600)
g
}
plot_lotek(dfm, "topleft", "Morogoro", col_pal, plot_out_dir)
