#' ---
#' title: Occupancy Models for African Bird Atlas Data
#' format: clean-revealjs
#' logo: Seec.png
#' slide-number: true
#' embed-resources: true
#' html-math-method:
#'   method: mathjax
#'   url: https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js
#' author:
#'   name: Res Altwegg
#'   orcid: 0000-0002-4083-6561
#'   email: res.altwegg@uct.ac.za
#'   affiliations: Centre for Statistics in Ecology, Environment and Conservation, University
#'     of Cape Town
#' date: 3 September 2026
#' ---
#'

#| eval: false
#| echo: true
#
# install.packages("remotes")
# remotes::install_github("AfricaBirdData/ABAP")

#| eval: true
#| echo: true

library(ABAP)
library(sf)
library(dplyr, warn.conflicts = FALSE)
library(ggplot2)

#| eval: true
#| echo: true

# We can search for all starling species
starlings <- searchAbapSpecies("Starling")
starlings$Common_species[26:35]

# Then we can extract the code we are interested in
starlings[starlings$Common_species == "Pied", "Spp"]

#| eval: true
#| echo: true

spp_code   <- 746 # pied starling
region     <- "Western Cape"
region_type <- "province"
years <- c(2023,2024)

#| eval: true
#| echo: true
#| cache: true

my_Abap_data <- getAbapData(
  .spp_code     = spp_code,
  .region_type  = region_type,
  .region       = region,
  .years = years
)


# Force R to print massive widths without truncating or wrapping columns and avoid scientific notation
options(width = 999, scipen = 999, digits = 3)

#| eval: true
#| echo: true
#| class-output: .hscroll

print(my_Abap_data, width = Inf)

#| eval: true
#| echo: true

length(unique(my_Abap_data$Pentad)) # How many unique pentads?

#| eval: true
#| echo: true

length(unique(my_Abap_data$ObserverNo)) # How many unique observers

#| eval: true
#| echo: false
#| output: false



#width = 5; height = 3.5; pointsize = 13

#svg("observers.svg", width = width, height = height, pointsize=pointsize)
svg("exploring.svg", width = 9, height = 7, pointsize=13)

par(mar=c(5,5,3,1), oma=rep(1,4), mfrow=c(2,2))

hist(table(my_Abap_data$ObserverNo), breaks=30, main="Number of checklist per observer", xlab="# checklists", las=1)
hist(my_Abap_data$TotalHours, main="Number of hours spent per checklist", xlab="# hours", las=1)
barplot(table(my_Abap_data$TotalHours[my_Abap_data$TotalHours<10]), main="Only looking at checklist with < 10 hours", xlab="# hours", las=1)
plot(TotalSpp ~ TotalHours, data=my_Abap_data, axes=F, xlab="Hours observed", ylab = "# species detected", main = "Species vs hours", cex=0.5)
axis(1)
axis(2)
dev.off()

#| eval: true
#| echo: true

# Eliminate checklists with <2 hours effort
my_Abap_data <- my_Abap_data[my_Abap_data$TotalHours>1,]

# Eliminate checklists with 0 species detected
my_Abap_data <- my_Abap_data[my_Abap_data$TotalSpp>0,]


#| eval: true
#| echo: true

WC_pentads <- getRegionPentads(.region_type = region_type, .region = region)

#| eval: true
#| echo: true

visit_data <- my_Abap_data %>%
  left_join(WC_pentads, by = c("Pentad" = "pentad")) %>%
  st_as_sf()

#| eval: true
#| echo: true

head(visit_data)

#| eval: true
#| echo: true

# aggregate pdata er pentad
visit_data_ag <- visit_data %>%
  group_by(Pentad) %>%
  summarise(
    totalLists = n(),
    detected = sum(Spp==spp_code)
   )

# calculate reporting rate
visit_data_ag$ReportingRate <- visit_data_ag$detected/visit_data_ag$totalLists

#| eval: true
#| echo: false

library(patchwork)

p1 <- ggplot(data = visit_data_ag) +
  geom_sf(aes(fill = log(totalLists))) +
  theme_minimal() +
  scale_fill_viridis_c() +
  labs(
    title = "log(Number of checklists)",
    fill = "log(Count)"
  )

p2 <- ggplot(data = visit_data_ag) +
  geom_sf(aes(fill = ReportingRate)) +
  theme_minimal() +
  scale_fill_viridis_c() +
  labs(
    title = "Reporting rate",
    fill = "Proportion"
  )

p1+p2


#| eval: true
#| echo: false
#| output: false

par(mar=c(5,5,0,1), oma=rep(0,4))
svg("rr.svg", width = 12, height = 12, pointsize=30)
hist(visit_data_ag$ReportingRate[visit_data_ag$totalLists>1], main="Reporting rate in pentads with > 1 checklist", xlab="Reporting rate", las=1)
dev.off()

#| eval: true
#| echo: true

site_covs <- read.csv("site_covs.csv")

# this drops any pentads not in the abap data of interest, and ensures the order matches
pentad_order_abap <- sort(unique(visit_data$Pentad))
site_covs_WC <- site_covs[match(pentad_order_abap, site_covs$pentad),]

#| eval: true
#| echo: false

require(patchwork)

site_covs_sf <- st_as_sf(site_covs_WC,
                  coords = c("lon", "lat"),
                  crs = 4326, remove=FALSE)

pa <- ggplot(data = site_covs_sf) +
  geom_sf(aes(color = elev), size = 0.5) +
  theme_minimal() +
  labs(title = "Elevation", color = "")

pb <- ggplot(data = site_covs_sf) +
  geom_sf(aes(color = prcp), size = 0.5) +
  theme_minimal() +
  labs(title = "Precipitation", color = "")

pc <- ggplot(data = site_covs_sf) +
  geom_sf(aes(color = tmp_max), size = 0.5) +
  theme_minimal() +
  labs(title = "Max Temp", color = "")

pd <- ggplot(data = site_covs_sf) +
  geom_sf(aes(color = NDVI), size = 0.5) +
  theme_minimal() +
  labs(title = "NDVI", color = "")

pe <- ggplot(data = site_covs_sf) +
  geom_sf(aes(color = tmp_min), size = 0.5) +
  theme_minimal() +
  labs(title = "Min Temp", color = "")

pf <- ggplot(data = site_covs_sf) +
  geom_sf(aes(color = log(hum.dens)), size = 0.5) +
  theme_minimal() +
  labs(title = "log(Hum Dens)", color = "")

(pe + pa + pd) / (pc + pf + pb)

#| eval: true
#| echo: false
#| cache: true
#| warning: false

library(GGally)

ggpairs(site_covs_WC[,4:9])


#| eval: true
#| echo: true
#| cache: true

source("00_helper_functions_dataprep.R")

#check
stopifnot(identical(site_covs_WC$pentad, pentad_order_abap))

#need to drop the identifier as unmarked does not use an ID column
site_covs_WC <- site_covs_WC[, setdiff(names(site_covs_WC), "pentad"), drop = FALSE]

mySpecies_um <- abapToUnmarked_single2(
  abap_data = my_Abap_data,
  pentads   = WC_pentads,
  cap_visits = 50,       # caps visits at 50
  siteCovs = site_covs_WC,
  siteCovs_clean = TRUE  # selects complete cases
)

#| eval: true
#| echo: true

str(mySpecies_um)


#| eval: true
#| echo: true

library(unmarked)

m1 <- occu(~ log(hours) ~ NDVI + tmp_max + tmp_min + log(hum.dens), data = mySpecies_um)
summary(m1)


#| eval: true
#| echo: true

WC_pentad_predictions <- WC_pentads %>%
  left_join(
    site_covs, by = c("pentad" = "pentad")) %>%
  st_as_sf()

pr <- predict(m1, type='state', newdata=as.data.frame(WC_pentad_predictions[,c("NDVI","tmp_min","tmp_max","hum.dens")]), appendData=F)

WC_pentad_predictions <- cbind(WC_pentad_predictions, pr)


#| eval: true
#| echo: false

p3 <- ggplot(data = WC_pentad_predictions) +
  geom_sf(aes(fill = Predicted)) +
  theme_minimal() +
  scale_fill_viridis_c() +
  labs(
    title = "Predicted occupancy",
    fill = ""
  )

p3 + p2


#| eval: true
#| echo: false
#| output: false


df_clean <- na.omit(WC_pentad_predictions)

svg("covariate_relationships.svg", width = 10, height = 10, pointsize=30)

par(mfrow=c(2,2), mar=c(4,1,1,2),oma=c(1,5,1,0))

NDVI.p <- data.frame(NDVI = seq(min(df_clean$NDVI), max(df_clean$NDVI), length=100),
                     tmp_min = mean(df_clean$tmp_min), tmp_max = mean(df_clean$tmp_max), hum.dens = mean(df_clean$hum.dens))

plot(NDVI.p$NDVI,predict(m1,type='state',newdata=NDVI.p)$Predicted, type='l', axes=F, xlab='NDVI',ylim=c(0,1),ylab='')
axis(1,seq(round(min(df_clean$NDVI),1), round(max(df_clean$NDVI),1), length=3))
axis(2,c(0,0.5,1),las=1)

tmp_min.p <- data.frame(NDVI = mean(df_clean$NDVI), tmp_min = seq(min(df_clean$tmp_min), max(df_clean$tmp_min), length=100), tmp_max = mean(df_clean$tmp_max), hum.dens = mean(df_clean$hum.dens))

plot(tmp_min.p$tmp_min,predict(m1,type='state',newdata=tmp_min.p)$Predicted, type='l',axes=F, xlab='Min Temp',ylim=c(0,1),ylab='')
axis(1, seq(round(min(df_clean$tmp_min),0), round(max(df_clean$tmp_min),0), length=3))

tmp_max.p <- data.frame(NDVI = mean(df_clean$NDVI), tmp_min = mean(df_clean$tmp_min), tmp_max = seq(min(df_clean$tmp_max), max(df_clean$tmp_max), length=100), hum.dens = mean(df_clean$hum.dens))

plot(tmp_max.p$tmp_max,predict(m1,type='state',newdata=tmp_max.p)$Predicted, type='l',axes=F, xlab='Max temp',ylim=c(0,1),ylab='')
axis(1,seq(round(min(df_clean$tmp_max),0), round(max(df_clean$tmp_max),0), length=3))
axis(2,c(0,0.5,1),las=1)

hum.dens.p <- data.frame(NDVI = mean(df_clean$NDVI), tmp_min = mean(df_clean$tmp_min), tmp_max = mean(df_clean$tmp_max), hum.dens = seq(min(df_clean$hum.dens), max(df_clean$hum.dens), length=100))

plot(hum.dens.p$hum.dens,predict(m1,type='state',newdata=hum.dens.p)$Predicted, type='l',axes=F, xlab='Hum density',ylim=c(0,1),ylab='')
axis(1,seq(round(min(df_clean$hum.dens),0), round(max(df_clean$hum.dens),0), length=3))
mtext("Occupancy probability",side=2,line=3,outer=T)

dev.off()

svg("covariate_relationships_det.svg", width = 5, height = 5, pointsize=20)

par(mfrow=c(1,1), mar=c(4,4,1,2),oma=c(0,0,0,0))

hours.p <- data.frame(hours = seq(2,55, length=100))

plot(hours.p$hours,predict(m1,type='det',newdata=hours.p)$Predicted, type='l',axes=F, xlab='Hours', ylab="Detection probability", ylim=c(0,1))
axis(1, seq(2,55, length=3))
axis(2,c(0,0.5,1),las=1)

dev.off()