1 Libraries & study IDs
library(tidyverse)
library(lubridate)
library(ctmm)
library(adehabitatHR)
library(sf)
library(sp)
library(terra)
library(MASS)
past_ids_vec <- c(
"AG082", "AG083", "AG088", "AG089", "AG096", "AG432",
"st2010-1048_2", "st2010-1052", "st2010-1053"
)
present_ids_vec <- c(
"171604 - adWBV05", "171605 - Matira", "171606 - Ilkeliani",
"171607 - Baldrick", "171608 - adWBV06", "171609 - Charlie",
"171612 - Olarro", "171613 - Tipilikwani"
)
all_ids <- c(past_ids_vec, present_ids_vec)
2 Load cleaned data
# Rename to t_/x_/y_ so all downstream chunks are consistent
gps <- readRDS("WBV_resampled_2h.rds") %>%
rename(t_ = timestamp, x_ = longitude, y_ = latitude)
WB_clean <- readRDS("WB_clean.rds")
wb_mig_past <- readRDS("WB_mig_past_proxy.rds")
3 Individual annual AKDEs
tel_list <- lapply(all_ids, function(nm) {
df <- gps %>% filter(individual.local.identifier == nm)
as.telemetry(data.frame(
individual.local.identifier = df$individual.local.identifier,
timestamp = as.POSIXct(df$t_, tz = "UTC"),
longitude = df$x_,
latitude = df$y_
))
})
names(tel_list) <- all_ids
cat("Telemetry objects created:", length(tel_list), "\n")
## Telemetry objects created: 17
# Run once only — takes 20–40 min
cat("Fitting movement models...\n")
fits <- lapply(all_ids, function(nm) {
cat(" Fitting:", nm, "\n")
guess <- ctmm.guess(tel_list[[nm]], interactive = FALSE)
ctmm.select(tel_list[[nm]], guess, verbose = FALSE)
})
names(fits) <- all_ids
saveRDS(fits, "WBV_fits.rds")
# Run once only
fits <- readRDS("WBV_fits.rds")
cat("Computing AKDEs...\n")
akdes <- lapply(all_ids, function(nm) {
cat(" AKDE:", nm, "\n")
akde(tel_list[[nm]], fits[[nm]])
})
names(akdes) <- all_ids
saveRDS(akdes, "WBV_akdes.rds")
cat("Saved: WBV_akdes.rds\n")
akdes <- readRDS("WBV_akdes.rds")
ud_areas <- do.call(rbind, lapply(all_ids, function(nm) {
do.call(rbind, lapply(c(0.25, 0.50, 0.75, 0.95), function(lev) {
est <- summary(akdes[[nm]], level.UD = lev)$CI["area (square kilometers)", "est"]
data.frame(
id = nm,
period = if_else(nm %in% past_ids_vec, "Past", "Present"),
level = paste0(lev * 100, "%"),
area = est
)
}))
}))
saveRDS(ud_areas, "WBV_ud_areas.rds")
ud_areas %>%
group_by(period, level) %>%
summarise(mean = round(mean(area)), median = round(median(area)),
n = n(), .groups = "drop") %>%
print()
## # A tibble: 8 × 5
## period level mean median n
## <chr> <chr> <dbl> <dbl> <int>
## 1 Past 25% 6114 3327 9
## 2 Past 50% 15944 9094 9
## 3 Past 75% 37076 19587 9
## 4 Past 95% 89211 43264 9
## 5 Present 25% 2045 1703 8
## 6 Present 50% 5721 5215 8
## 7 Present 75% 13793 13033 8
## 8 Present 95% 39534 30760 8
for (lev in c("25%", "50%", "75%", "95%")) {
p <- ud_areas %>% filter(level == lev, period == "Past") %>% pull(area)
r <- ud_areas %>% filter(level == lev, period == "Present") %>% pull(area)
wt <- wilcox.test(p, r)
cat(sprintf("%s UD — Past mean=%.0f vs Present mean=%.0f, W=%.0f, p=%.3f\n",
lev, mean(p), mean(r), wt$statistic, wt$p.value))
}
## 25% UD — Past mean=6114 vs Present mean=2045, W=43, p=0.541
## 50% UD — Past mean=15944 vs Present mean=5721, W=41, p=0.673
## 75% UD — Past mean=37076 vs Present mean=13793, W=42, p=0.606
## 95% UD — Past mean=89211 vs Present mean=39534, W=43, p=0.541
4 Population-level UD contours
template_ud <- terra::rast(
terra::ext(370000, 850000, 9550000, 9950000),
res = 5000,
crs = "+proj=utm +zone=36 +south +datum=WGS84 +units=m"
)
mean_ud_rast <- function(ids) {
rasters <- lapply(ids, function(nm) {
r <- terra::rast(raster::raster(akdes[[nm]], DF = "PDF"))
r <- terra::project(r, template_ud)
terra::resample(r, template_ud, method = "bilinear")
})
terra::mean(terra::rast(rasters), na.rm = TRUE)
}
mean_past_ud <- mean_ud_rast(past_ids_vec)
mean_present_ud <- mean_ud_rast(present_ids_vec)
ud_contours <- function(r) {
v <- terra::values(r); v[is.na(v) | v < 0] <- 0
sorted <- sort(v, decreasing = TRUE)
cumprob <- cumsum(sorted) / sum(sorted)
do.call(rbind, lapply(c(0.95, 0.75, 0.50, 0.25), function(lev) {
thresh <- sorted[which(cumprob >= lev)[1]]
r_bin <- r >= thresh; r_bin[r_bin == 0] <- NA
poly <- terra::as.polygons(r_bin) %>%
sf::st_as_sf() %>% sf::st_union() %>% sf::st_as_sf() %>%
sf::st_set_crs("+proj=utm +zone=36 +south +datum=WGS84 +units=m")
poly$level <- paste0(lev * 100, "%")
poly
}))
}
past_contours <- ud_contours(mean_past_ud) %>% mutate(period = "Past")
present_contours <- ud_contours(mean_present_ud) %>% mutate(period = "Present")
all_contours <- rbind(past_contours, present_contours) %>%
sf::st_transform(4326) %>%
mutate(level = factor(level, c("95%", "75%", "50%", "25%")))
saveRDS(all_contours, "WBV_pop_ud_contours.rds")
cat("Saved: WBV_pop_ud_contours.rds\n")
## Saved: WBV_pop_ud_contours.rds
5 Protected area overlap
pas <- sf::st_read(
"/Users/oli/University/Zoology/5th_Year_MSc/MSC_Thesis/GIS/FINAL_OUTPUTS.gpkg",
quiet = TRUE
)
pa_overlap <- do.call(rbind, lapply(all_ids, function(nm) {
do.call(rbind, lapply(c("50%", "95%"), function(lev) {
lev_num <- as.numeric(sub("%", "", lev)) / 100
ud_r <- terra::rast(raster::raster(akdes[[nm]], DF = "PDF"))
v <- terra::values(ud_r); v[is.na(v) | v < 0] <- 0
sorted <- sort(v, decreasing = TRUE)
cumprob <- cumsum(sorted) / sum(sorted)
thresh <- sorted[which(cumprob >= lev_num)[1]]
ud_bin <- ud_r >= thresh
pa_utm <- sf::st_transform(pas, terra::crs(ud_r))
inside <- terra::mask(ud_bin, terra::vect(pa_utm))
pct <- sum(terra::values(inside) == 1, na.rm = TRUE) /
sum(terra::values(ud_bin) == 1, na.rm = TRUE) * 100
data.frame(
id = nm,
period = if_else(nm %in% past_ids_vec, "Past", "Present"),
level = lev,
pct_in_pa = pct
)
}))
}))
saveRDS(pa_overlap, "WBV_pa_overlap.rds")
cat("Saved: WBV_pa_overlap.rds\n")
## Saved: WBV_pa_overlap.rds
pa_overlap %>%
group_by(period, level) %>%
summarise(mean = round(mean(pct_in_pa), 1), n = n(), .groups = "drop")
## # A tibble: 4 × 4
## period level mean n
## <chr> <chr> <dbl> <int>
## 1 Past 50% 63.2 9
## 2 Past 95% 41.9 9
## 3 Present 50% 63 8
## 4 Present 95% 57.3 8
for (lev in c("50%", "95%")) {
p <- pa_overlap %>% filter(level == lev, period == "Past") %>% pull(pct_in_pa)
r <- pa_overlap %>% filter(level == lev, period == "Present") %>% pull(pct_in_pa)
wt <- wilcox.test(p, r)
cat(sprintf("%s UD — Past=%.1f%% vs Present=%.1f%%, W=%.0f, p=%.3f\n",
lev, mean(p), mean(r), wt$statistic, wt$p.value))
}
## 50% UD — Past=63.2% vs Present=63.0%, W=38, p=0.912
## 95% UD — Past=41.9% vs Present=57.3%, W=22, p=0.200
6 Seasonal AKDEs (wet / dry)
avail <- gps %>%
mutate(month = month(as.POSIXct(t_, tz = "UTC"))) %>%
group_by(individual.local.identifier, period) %>%
summarise(
wet_months = n_distinct(month[month %in% c(11, 12, 1, 2, 3, 4)]),
dry_months = n_distinct(month[month %in% c(5, 6, 7, 8, 9, 10)]),
.groups = "drop"
)
wet_ids <- avail %>% filter(wet_months == 6) %>% pull(individual.local.identifier)
dry_ids <- avail %>% filter(dry_months >= 4) %>% pull(individual.local.identifier)
cat("Wet season qualifying:", length(wet_ids),
"— Past:", sum(wet_ids %in% past_ids_vec),
"Present:", sum(wet_ids %in% present_ids_vec), "\n")
## Wet season qualifying: 16 — Past: 8 Present: 8
cat("Dry season qualifying:", length(dry_ids),
"— Past:", sum(dry_ids %in% past_ids_vec),
"Present:", sum(dry_ids %in% present_ids_vec), "\n")
## Dry season qualifying: 17 — Past: 9 Present: 8
# Wet — run once
vul_wet <- gps %>%
filter(individual.local.identifier %in% wet_ids,
month %in% c(11, 12, 1, 2, 3, 4))
tel_wet <- lapply(wet_ids, function(nm) {
df <- vul_wet %>% filter(individual.local.identifier == nm)
as.telemetry(data.frame(
individual.local.identifier = df$individual.local.identifier,
timestamp = as.POSIXct(df$t_, tz = "UTC"),
longitude = df$x_,
latitude = df$y_
))
})
names(tel_wet) <- wet_ids
fits_wet <- lapply(wet_ids, function(nm) {
cat(" Wet:", nm, "\n")
ctmm.select(tel_wet[[nm]], ctmm.guess(tel_wet[[nm]], interactive = FALSE),
verbose = FALSE)
})
names(fits_wet) <- wet_ids
akdes_wet <- lapply(wet_ids, function(nm) akde(tel_wet[[nm]], fits_wet[[nm]]))
names(akdes_wet) <- wet_ids
saveRDS(akdes_wet, "WBV_std_akdes_wet.rds")
# Dry — run once
vul_dry <- gps %>%
filter(individual.local.identifier %in% dry_ids,
month %in% c(5, 6, 7, 8, 9, 10))
tel_dry <- lapply(dry_ids, function(nm) {
df <- vul_dry %>% filter(individual.local.identifier == nm)
as.telemetry(data.frame(
individual.local.identifier = df$individual.local.identifier,
timestamp = as.POSIXct(df$t_, tz = "UTC"),
longitude = df$x_,
latitude = df$y_
))
})
names(tel_dry) <- dry_ids
fits_dry <- lapply(dry_ids, function(nm) {
cat(" Dry:", nm, "\n")
ctmm.select(tel_dry[[nm]], ctmm.guess(tel_dry[[nm]], interactive = FALSE),
verbose = FALSE)
})
names(fits_dry) <- dry_ids
akdes_dry <- lapply(dry_ids, function(nm) akde(tel_dry[[nm]], fits_dry[[nm]]))
names(akdes_dry) <- dry_ids
saveRDS(akdes_dry, "WBV_std_akdes_dry.rds")
cat("Seasonal AKDEs saved.\n")
akdes_wet <- readRDS("WBV_std_akdes_wet.rds")
akdes_dry <- readRDS("WBV_std_akdes_dry.rds")
get_area <- function(akde_list) {
sapply(akde_list, function(a)
summary(a)$CI["area (square kilometers)", "est"])
}
areas_wet <- get_area(akdes_wet)
areas_dry <- get_area(akdes_dry)
cat("Wet 95% UD — Past mean:",
round(mean(areas_wet[names(areas_wet) %in% past_ids_vec])), "km²\n")
## Wet 95% UD — Past mean: 140316 km²
cat("Wet 95% UD — Present mean:",
round(mean(areas_wet[names(areas_wet) %in% present_ids_vec])), "km²\n")
## Wet 95% UD — Present mean: 44461 km²
cat("Dry 95% UD — Past mean:",
round(mean(areas_dry[names(areas_dry) %in% past_ids_vec])), "km²\n")
## Dry 95% UD — Past mean: 26271 km²
cat("Dry 95% UD — Present mean:",
round(mean(areas_dry[names(areas_dry) %in% present_ids_vec])), "km²\n")
## Dry 95% UD — Present mean: 29893 km²
wt_wet <- wilcox.test(areas_wet[names(areas_wet) %in% past_ids_vec],
areas_wet[names(areas_wet) %in% present_ids_vec])
wt_dry <- wilcox.test(areas_dry[names(areas_dry) %in% past_ids_vec],
areas_dry[names(areas_dry) %in% present_ids_vec])
cat(sprintf("Wet Past vs Present: W=%.0f, p=%.3f\n", wt_wet$statistic, wt_wet$p.value))
## Wet Past vs Present: W=40, p=0.442
cat(sprintf("Dry Past vs Present: W=%.0f, p=%.3f\n", wt_dry$statistic, wt_dry$p.value))
## Dry Past vs Present: W=24, p=0.277
7 Individual seasonal KDEs (fixed-bandwidth, h = 5 km)
make_kde_sf_utm <- function(lon, lat, level, h = 5000) {
pts <- sf::st_as_sf(data.frame(lon = lon, lat = lat),
coords = c("lon", "lat"), crs = 4326) %>%
sf::st_transform("+proj=utm +zone=36 +south +datum=WGS84")
coords_utm <- sf::st_coordinates(pts)
sp_pts <- SpatialPoints(coords_utm,
proj4string = CRS("+proj=utm +zone=36 +south +datum=WGS84"))
k <- kernelUD(sp_pts, h = h, grid = 500)
r <- raster::raster(k)
v <- raster::values(r); v[is.na(v)] <- 0
sorted <- sort(v, decreasing = TRUE)
thresh <- sorted[which(cumsum(sorted) / sum(sorted) >= level)[1]]
r_bin <- r >= thresh; r_bin[r_bin == 0] <- NA
terra::as.polygons(terra::rast(r_bin)) %>%
sf::st_as_sf() %>% sf::st_union() %>% sf::st_as_sf() %>%
sf::st_set_crs("+proj=utm +zone=36 +south +datum=WGS84") %>%
sf::st_transform(4326)
}
ind_smooth_contours <- do.call(rbind, lapply(all_ids, function(id) {
df <- gps %>% filter(individual.local.identifier == id)
period <- unique(df$period)
do.call(rbind, lapply(c("Wet", "Dry"), function(seas) {
d <- df %>% filter(season == seas)
if (nrow(d) < 50) return(NULL)
do.call(rbind, lapply(c(0.95, 0.75, 0.50, 0.25), function(lev) {
poly <- tryCatch(make_kde_sf_utm(d$x_, d$y_, lev), error = function(e) NULL)
if (is.null(poly)) return(NULL)
poly %>% mutate(id = id, period = period, season = seas,
level = paste0(lev * 100, "%"))
}))
}))
}))
ind_smooth_contours <- ind_smooth_contours %>%
mutate(level = factor(level, c("95%", "75%", "50%", "25%")),
season = factor(season, c("Wet", "Dry")),
period = factor(period, c("Past", "Present")))
saveRDS(ind_smooth_contours, "WBV_ind_seasonal_contours.rds")
cat("Saved: WBV_ind_seasonal_contours.rds —", nrow(ind_smooth_contours), "polygons\n")
## Saved: WBV_ind_seasonal_contours.rds — 136 polygons
ind_sc_utm <- ind_smooth_contours %>%
sf::st_transform("+proj=utm +zone=36 +south +datum=WGS84")
ind_sc_utm$area_km2 <- as.numeric(sf::st_area(ind_sc_utm)) / 1e6
area_df <- ind_sc_utm %>% sf::st_drop_geometry() %>% filter(level == "95%")
area_df %>%
group_by(period, season) %>%
summarise(median_km2 = round(median(area_km2), 1),
mean_km2 = round(mean(area_km2), 1),
n = n(), .groups = "drop") %>%
print()
## # A tibble: 4 × 5
## period season median_km2 mean_km2 n
## <fct> <fct> <dbl> <dbl> <int>
## 1 Past Wet 10105. 10396. 9
## 2 Past Dry 4020. 4580. 9
## 3 Present Wet 11317. 12778. 8
## 4 Present Dry 10599 9981. 8
for (per in c("Past", "Present")) {
w <- area_df %>% filter(period == per, season == "Wet") %>% pull(area_km2)
d <- area_df %>% filter(period == per, season == "Dry") %>% pull(area_km2)
wt <- wilcox.test(w, d)
cat(sprintf("%s — Wet mean=%.0f vs Dry mean=%.0f km², W=%.0f, p=%.3f\n",
per, mean(w), mean(d), wt$statistic, wt$p.value))
}
## Past — Wet mean=10396 vs Dry mean=4580 km², W=64, p=0.040
## Present — Wet mean=12778 vs Dry mean=9981 km², W=41, p=0.382
8 Seasonal centroid displacement
centroid_dist <- do.call(rbind, lapply(all_ids, function(id) {
df <- gps %>% filter(individual.local.identifier == id)
period <- unique(df$period)
wet <- df %>% filter(season == "Wet")
dry <- df %>% filter(season == "Dry")
if (nrow(wet) < 50 | nrow(dry) < 50) return(NULL)
wet_sf <- sf::st_as_sf(data.frame(lon = mean(wet$x_), lat = mean(wet$y_)),
coords = c("lon", "lat"), crs = 4326) %>%
sf::st_transform("+proj=utm +zone=36 +south +datum=WGS84")
dry_sf <- sf::st_as_sf(data.frame(lon = mean(dry$x_), lat = mean(dry$y_)),
coords = c("lon", "lat"), crs = 4326) %>%
sf::st_transform("+proj=utm +zone=36 +south +datum=WGS84")
data.frame(
id = id, period = period,
wet_lon = mean(wet$x_), wet_lat = mean(wet$y_),
dry_lon = mean(dry$x_), dry_lat = mean(dry$y_),
dist_km = as.numeric(sf::st_distance(wet_sf, dry_sf)) / 1000
)
}))
past_d <- centroid_dist %>% filter(period == "Past") %>% pull(dist_km)
present_d <- centroid_dist %>% filter(period == "Present") %>% pull(dist_km)
cat(sprintf("Past displacement: mean=%.1f km (n=%d)\n", mean(past_d), length(past_d)))
## Past displacement: mean=66.8 km (n=9)
cat(sprintf("Present displacement: mean=%.1f km (n=%d)\n", mean(present_d), length(present_d)))
## Present displacement: mean=48.0 km (n=8)
wt_cd <- wilcox.test(past_d, present_d)
cat(sprintf("Wilcoxon W=%.0f, p=%.3f\n", wt_cd$statistic, wt_cd$p.value))
## Wilcoxon W=42, p=0.606
saveRDS(centroid_dist, "WBV_centroid_dist.rds")
9 Net squared displacement
vul_nsd <- gps %>%
filter(individual.local.identifier %in% all_ids) %>%
group_by(individual.local.identifier, period) %>%
arrange(as.POSIXct(t_, tz = "UTC"), .by_group = TRUE) %>%
mutate(
timestamp = as.POSIXct(t_, tz = "UTC"),
first_lon = first(x_),
first_lat = first(y_),
nsd_km2 = ((x_ - first_lon) * 111.32 * cos(first_lat * pi / 180))^2 +
((y_ - first_lat) * 110.57)^2,
day_of_track = as.numeric(difftime(timestamp, min(timestamp), units = "days"))
) %>%
ungroup() %>%
group_by(individual.local.identifier) %>%
mutate(
season_change = season != lag(season, default = first(season)),
bout_id = cumsum(season_change)
) %>%
ungroup()
saveRDS(vul_nsd, "WBV_nsd.rds")
cat("Saved: WBV_nsd.rds\n")
## Saved: WBV_nsd.rds
10 Site fidelity
annual_centroids <- gps %>%
filter(individual.local.identifier %in% all_ids) %>%
mutate(year = year(as.POSIXct(t_, tz = "UTC"))) %>%
group_by(individual.local.identifier, period, year) %>%
summarise(cx = mean(x_, na.rm = TRUE), cy = mean(y_, na.rm = TRUE),
n = n(), .groups = "drop") %>%
filter(n >= 50)
site_fidelity <- annual_centroids %>%
arrange(individual.local.identifier, year) %>%
group_by(individual.local.identifier, period) %>%
filter(n() >= 2) %>%
mutate(
dist_km = sqrt(
((cx - lag(cx)) * 111.32 * cos(cy * pi / 180))^2 +
((cy - lag(cy)) * 110.57)^2
)
) %>%
filter(!is.na(dist_km)) %>%
ungroup()
sf_past <- site_fidelity %>% filter(period == "Past") %>% pull(dist_km)
sf_present <- site_fidelity %>% filter(period == "Present") %>% pull(dist_km)
cat(sprintf("Past site fidelity: mean=%.1f km (n=%d)\n", mean(sf_past), length(sf_past)))
## Past site fidelity: mean=80.6 km (n=7)
cat(sprintf("Present site fidelity: mean=%.1f km (n=%d)\n", mean(sf_present), length(sf_present)))
## Present site fidelity: mean=26.5 km (n=24)
wt_sf <- wilcox.test(sf_past, sf_present)
cat(sprintf("Wilcoxon W=%.0f, p=%.3f\n", wt_sf$statistic, wt_sf$p.value))
## Wilcoxon W=150, p=0.001
saveRDS(site_fidelity, "WBV_site_fidelity.rds")
11 Southern zone use
subregions <- sf::st_read(
"/Users/oli/University/Zoology/5th_Year_MSc/MSC_Thesis/Thesis_R/mahony.subregions.1.kml",
quiet = TRUE
)
zones <- subregions %>%
mutate(zone = case_when(
Name %in% c("north-west", "central") ~ "North",
Name %in% c("south-east", "south-west") ~ "South",
TRUE ~ NA_character_
)) %>%
filter(!is.na(zone)) %>%
group_by(zone) %>%
summarise(geometry = sf::st_union(geometry)) %>%
sf::st_set_crs(4326)
vul_sf <- gps %>%
filter(individual.local.identifier %in% all_ids) %>%
sf::st_as_sf(coords = c("x_", "y_"), crs = 4326, remove = FALSE)
vul_zones <- sf::st_join(vul_sf, dplyr::select(zones, zone), join = sf::st_within)
zone_pct <- vul_zones %>%
sf::st_drop_geometry() %>%
mutate(zone = if_else(is.na(zone), "North", zone)) %>%
group_by(individual.local.identifier, period, season) %>%
summarise(
total = n(),
n_south = sum(zone == "South", na.rm = TRUE),
pct_south = n_south / total * 100,
.groups = "drop"
)
zone_pct %>%
group_by(period, season) %>%
summarise(mean_south = round(mean(pct_south), 1), n = n(), .groups = "drop") %>%
print()
## # A tibble: 4 × 4
## period season mean_south n
## <chr> <chr> <dbl> <int>
## 1 Past Dry 0.2 9
## 2 Past Wet 14.5 9
## 3 Present Dry 1.6 8
## 4 Present Wet 16.2 8
for (per in c("Past", "Present")) {
w <- zone_pct %>% filter(period == per, season == "Wet") %>% pull(pct_south)
d <- zone_pct %>% filter(period == per, season == "Dry") %>% pull(pct_south)
wt <- wilcox.test(w, d)
cat(sprintf("%s — Wet=%.1f%% vs Dry=%.1f%%, W=%.0f, p=%.3f\n",
per, mean(w), mean(d), wt$statistic, wt$p.value))
}
## Past — Wet=14.5% vs Dry=0.2%, W=74, p=0.001
## Present — Wet=16.2% vs Dry=1.6%, W=40, p=0.442
for (seas in c("Wet", "Dry")) {
p <- zone_pct %>% filter(period == "Past", season == seas) %>% pull(pct_south)
r <- zone_pct %>% filter(period == "Present", season == seas) %>% pull(pct_south)
wt <- wilcox.test(p, r)
cat(sprintf("%s — Past=%.1f%% vs Present=%.1f%%, W=%.0f, p=%.3f\n",
seas, mean(p), mean(r), wt$statistic, wt$p.value))
}
## Wet — Past=14.5% vs Present=16.2%, W=34, p=0.912
## Dry — Past=0.2% vs Present=1.6%, W=8, p=0.004
saveRDS(zone_pct, "WBV_zone_pct.rds")
12 Vulture–wildebeest spatial overlap
make_wb_rast <- function(df, h = 5000, grid = 300) {
coords <- df %>% filter(!is.na(UTM_X), !is.na(UTM_Y))
sp_obj <- SpatialPoints(
cbind(coords$UTM_X, coords$UTM_Y),
proj4string = CRS("+proj=utm +zone=36 +south +datum=WGS84")
)
ud <- kernelUD(sp_obj, h = h, grid = grid)
spdf <- as(ud, "SpatialPixelsDataFrame")
r <- terra::rast(spdf)
terra::crs(r) <- "+proj=utm +zone=36 +south +datum=WGS84"
r_ll <- terra::project(r, "EPSG:4326")
r_ll / terra::global(r_ll, "max", na.rm = TRUE)$max
}
wb_mig_wet <- make_wb_rast(wb_mig_past %>% filter(season == "Wet"))
wb_mig_dry <- make_wb_rast(wb_mig_past %>% filter(season == "Dry"))
wb_mig_dry_r <- terra::resample(wb_mig_dry, wb_mig_wet, method = "bilinear")
make_mask <- function(r) {
v <- terra::values(r); v[is.na(v)] <- 0
s <- sort(v, decreasing = TRUE)
thr <- s[which(cumsum(s) / sum(s) >= 0.50)[1]]
m <- r >= thr; m[m == 0] <- NA; m
}
mask_mig_wet <- make_mask(wb_mig_wet)
mask_mig_dry <- make_mask(wb_mig_dry_r)
saveRDS(wb_mig_wet, "WB_mig_wet_rast.rds")
saveRDS(wb_mig_dry_r, "WB_mig_dry_rast.rds")
cat("Wildebeest rasters and masks built.\n")
## Wildebeest rasters and masks built.
overlap_pct <- function(akde_obj, wb_mask) {
r_vul <- terra::rast(raster::raster(akde_obj, DF = "PDF"))
r_vul <- terra::project(r_vul, terra::crs(wb_mask))
r_vul <- terra::resample(r_vul, wb_mask, method = "bilinear")
v <- terra::values(r_vul); v[is.na(v) | v < 0] <- 0; v <- v / sum(v)
sum(v[!is.na(terra::values(wb_mask))]) * 100
}
ov_wet <- sapply(names(akdes_wet),
function(nm) overlap_pct(akdes_wet[[nm]], mask_mig_wet))
ov_dry <- sapply(names(akdes_dry),
function(nm) overlap_pct(akdes_dry[[nm]], mask_mig_dry))
ov_df <- bind_rows(
data.frame(id = names(ov_wet), overlap = ov_wet, season = "Wet",
period = if_else(names(ov_wet) %in% past_ids_vec, "Past", "Present")),
data.frame(id = names(ov_dry), overlap = ov_dry, season = "Dry",
period = if_else(names(ov_dry) %in% past_ids_vec, "Past", "Present"))
) %>%
mutate(period = factor(period, c("Past", "Present")),
season = factor(season, c("Wet", "Dry")))
ov_df %>%
group_by(period, season) %>%
summarise(mean_pct = round(mean(overlap), 1), n = n(), .groups = "drop") %>%
print()
## # A tibble: 4 × 4
## period season mean_pct n
## <fct> <fct> <dbl> <int>
## 1 Past Wet 5.9 8
## 2 Past Dry 43.6 9
## 3 Present Wet 9.7 8
## 4 Present Dry 33.7 8
for (per in c("Past", "Present")) {
w <- ov_df %>% filter(period == per, season == "Wet") %>% pull(overlap)
d <- ov_df %>% filter(period == per, season == "Dry") %>% pull(overlap)
wt <- wilcox.test(w, d)
cat(sprintf("%s — Wet=%.1f%% vs Dry=%.1f%%, W=%.0f, p=%.3f\n",
per, mean(w), mean(d), wt$statistic, wt$p.value))
}
## Past — Wet=5.9% vs Dry=43.6%, W=1, p=0.000
## Present — Wet=9.7% vs Dry=33.7%, W=5, p=0.003
saveRDS(ov_df, "WBV_wb_overlap.rds")
cat("Saved: WBV_wb_overlap.rds\n")
## Saved: WBV_wb_overlap.rds
13 Bhattacharyya coefficient — pairwise overlap
ids_bc <- names(akdes)
n_ids <- length(ids_bc)
bc_mat <- matrix(NA, nrow = n_ids, ncol = n_ids,
dimnames = list(ids_bc, ids_bc))
for (i in seq_len(n_ids)) {
for (j in seq_len(n_ids)) {
if (i == j) { bc_mat[i, j] <- 1; next }
bc_mat[i, j] <- tryCatch(
ctmm::overlap(list(akdes[[ids_bc[i]]], akdes[[ids_bc[j]]]))$CI[2],
error = function(e) NA
)
}
}
bc_df <- as.data.frame(as.table(bc_mat)) %>%
rename(id1 = Var1, id2 = Var2, BC = Freq) %>%
filter(id1 != id2) %>%
mutate(
period1 = if_else(id1 %in% past_ids_vec, "Past", "Present"),
period2 = if_else(id2 %in% past_ids_vec, "Past", "Present"),
pair_type = case_when(
period1 == "Past" & period2 == "Past" ~ "Within-Past",
period1 == "Present" & period2 == "Present" ~ "Within-Present",
TRUE ~ "Between"
)
)
bc_df %>%
group_by(pair_type) %>%
summarise(mean_BC = round(mean(BC, na.rm = TRUE), 3), n = n(), .groups = "drop") %>%
print()
## # A tibble: 3 × 3
## pair_type mean_BC n
## <chr> <dbl> <int>
## 1 Between NaN 144
## 2 Within-Past NaN 72
## 3 Within-Present NaN 56
saveRDS(bc_df, "WBV_bc_overlap.rds")
cat("Saved: WBV_bc_overlap.rds\n")
## Saved: WBV_bc_overlap.rds