2_ConvertToOccupancy.Rmd

---
title: "ConvertToOccupancy"
output: html_document
date: "2025-05-31"
---

This file generates the species occupancy reports for occupancy modeling
At the end, it also runs the occupancy models

```{r}

library(dplyr)
library(tidyr)
library(lubridate)
library(purrr)
library(rmarkdown)
library(stringr)
library(ggplot2)
library(forcats)
library(googlesheets4)

```

```{r}

# Specifying input and output paths
input_path <- "../Input Data/"
report_output_path <- "../Species Occupancy/"

if (dir.exists(report_output_path) == F) {
  dir.create(report_output_path)
}

# Import file, change filepath/name as necessary
df <- read.csv(paste0(input_path, "multi_city_data_cleaned.csv"))

# Import file, change filepath/name as necessary
active_dates <- read.csv(paste0(input_path, "all_city_active_dates.csv"))

# Import mapping from original site IDs to clustered site IDs
site_map <- read.csv(paste0(input_path, "camera_sites_with_impervious.csv"))

site_map <- site_map %>%
  select(City, locAbbr, locAbbr2) %>%
  mutate(
    locAbbr = ifelse(City == "PACA", sub("^[^-]+-", "", locAbbr), locAbbr),
    locAbbr2 = ifelse(City == "PACA", sub("^[^-]+-", "", locAbbr2), locAbbr2)
  ) %>%
  distinct(City, locAbbr, locAbbr2)

# Align photo site IDs to clustered site IDs used in active_dates
df <- df %>%
  left_join(site_map, by = c("City", "locAbbr"))

df <- df %>%
  mutate(locAbbr = coalesce(locAbbr2, locAbbr)) %>%
  select(-locAbbr2)

```

```{r}

# Format before was FALSE for inactive, TRUE for active. Converts FALSE to NA and TRUE to 0 (will update to 1 if present)
active_dates <- active_dates %>%
  mutate(across(starts_with("Day_"), ~ ifelse(. == FALSE, NA, 0))) %>%
  arrange(Season)

# If PACA, remove the first 4 characters to match the df locAbbr names
active_dates <- active_dates %>%
  mutate(locAbbr = ifelse(
    City == "PACA",
    substr(locAbbr, 5, nchar(locAbbr)),
    locAbbr
  ))

# Dataframe with season information
season_dates <- data.frame(
  Season = 1:41,
  start_date = as.Date(c("2015-12-18",
                         "2016-03-18", "2016-06-17", "2016-09-17", "2016-12-18",
                         "2017-03-18", "2017-06-17", "2017-09-17", "2017-12-18",
                         "2018-03-18", "2018-06-17", "2018-09-17", "2018-12-18", 
                         "2019-03-18", "2019-06-17", "2019-09-17", "2019-12-18", 
                         "2020-03-18", "2020-06-17", "2020-09-17", "2020-12-18", 
                         "2021-03-18", "2021-06-17", "2021-09-17", "2021-12-18", 
                         "2022-03-18", "2022-06-17", "2022-09-17", "2022-12-18", 
                         "2023-03-18", "2023-06-17", "2023-09-17", "2023-12-18", 
                         "2024-03-18", "2024-06-17", "2024-09-17", "2024-12-18",
                         "2025-03-18", "2025-06-17", "2025-09-17", "2025-12-18")),
  end_date = as.Date(c("2016-02-14", "2016-05-14", "2016-08-14", "2016-11-14",
                       "2017-02-14", "2017-05-14", "2017-08-14", "2017-11-14",
                       "2018-02-14", "2018-05-14", "2018-08-14", "2018-11-14",
                       "2019-02-13", "2019-05-14", "2019-08-14", "2019-11-14", 
                       "2020-02-13", "2020-05-14", "2020-08-14", "2020-11-14", 
                       "2021-02-13", "2021-05-14", "2021-08-14", "2021-11-14", 
                       "2022-02-13", "2022-05-14", "2022-08-14", "2022-11-14", 
                       "2023-02-13", "2023-05-14", "2023-08-14", "2023-11-14", 
                       "2024-02-13", "2024-05-14", "2024-08-14", "2024-11-14",
                       "2025-02-13", "2025-05-14", "2025-08-14", "2025-11-14",
                       "2026-02-13"))
)

# Converts Timestamp to proper date, combines with season_dates to add relative_day column (what day of the season)
df <- df %>%
  mutate(Timestamp = ymd_hms(Timestamp)) %>%
  left_join(season_dates, by = "Season") %>%
  mutate(relative_day = as.integer(difftime(Timestamp, start_date, units = "days")) + 1) %>%
  select(-start_date)

```

```{r}

# Create header
header_df <- season_dates %>%
  rename("Season Legend:" = "Season", "Start Date" = "start_date", "End Date" = "end_date")

# Adds empty line of header to make margin before data
end_header_text <- data.frame(
  Season = "",
  Start = "",
  End = ""
)

# Combines headers
header <- bind_rows(header_df, end_header_text)

# Outputs file
write.csv(header, paste0(report_output_path, "header.csv"), row.names = FALSE)

```

```{r}

# TODO: update species group file?
# Get species group mapping from Google Sheets "Species Groups"
species_groups_file <- range_read("https://docs.google.com/spreadsheets/d/1czht5zz8dRII25loY9RI9efTi6vIsvsioSqXBgIwuvg/edit?gid=0#gid=0", trim_ws = FALSE)

# Change format such that Subspecies (commonName in df) maps to SpeciesGroup
species_map <- species_groups_file %>%
  pivot_longer(
    cols = -1,
    names_to = "SubCol",
    values_to = "Subspecies"
  ) %>%
  filter(!is.na(Subspecies), Subspecies != "") %>%
  rename(SpeciesGroup = 1) %>%
  select(SpeciesGroup, Subspecies)

# Map commonName to SpeciesGroup value
df <- df %>%
  left_join(species_map, by = c("commonName" = "Subspecies")) %>%
  filter(SpeciesGroup != "DELETE" | is.na(SpeciesGroup)) %>%
  mutate(commonName = coalesce(SpeciesGroup, commonName)) %>%
  filter(!is.na(commonName)) %>%
  select(-SpeciesGroup)

species_list <- unique(df$commonName)

print(species_list)

write.csv(species_list, paste0(report_output_path, "species_list.csv"), row.names = FALSE)

```

```{r}

# Dataframe for tracking site-seasons for human vs. md for each species
species_counts <- data.frame(
  city = character(),
  species = character(),
  source = character(),
  count = integer(),
  stringsAsFactors = FALSE
)

# Track amount of images after MD "removes" them below the threshold
total_images <- data.frame(
  threshold = numeric(),
  count = integer()
)

# City names
cities <- unique(df$City)

# Confidence threshold
confs <- c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.99)

# Exclude these groups from occupancy reports
exclude <- c("Human", "Empty", "Unknown")

```

```{r}

# Iterate through each confidence threshold
for (thresh in confs) {
  # Remove images below confidence threshold
  thresh_df <- df %>%
    filter(maxDetectionConf >= thresh)
  
  # Count number of images at current threshold
  cur_count <- data.frame(
    threshold = thresh,
    count = nrow(thresh_df)
  )
  # Combine into whole dataframe
  total_images <- rbind(total_images, cur_count)
  
  # For each city
  for (city in cities) {
    # Filter down to current city
    city_active_dates <- active_dates %>%
      filter(City == city)
    
    city_df <- thresh_df %>%
      filter(City == city) %>%
      drop_na(Season)
    
    # Define path to directory
    folder_path <- paste0(report_output_path, city, "/")
    
    # Create directory if not already existing
    if (dir.exists(folder_path) == F) {
      dir.create(folder_path)
    }
    
    folder_path <- paste0(folder_path, "Conf. ", thresh, "/")
    
    if (dir.exists(folder_path) == F) {
      dir.create(folder_path)
    }
    
    # Iterate through each species
    for (species in species_list) {
      if (species %in% exclude) {
        next
      }
      
      # Subsets to only current species
      species_df <- city_df %>%
        filter(commonName == species)
      
      # Reset cur_grid (where each section is hosted and then appended to total results)
      cur_grid <- city_active_dates %>%
        mutate(commonName = species)
        
      # Ignore if species wasn't observed
      if (nrow(species_df) > 0) {
        # For each entry in dataframe
        for (i in 1:nrow(species_df)) {
          # Get values of interest
          season <- species_df$Season[i]
          relative_day <- species_df$relative_day[i]
          site <- species_df$locAbbr[i]
          location <- species_df$City[i]
          
          # Get name of day column
          day_column_name <- paste0("Day_", relative_day)
          
          # Put 1 for an observation
          cur_grid[cur_grid$Season == season & cur_grid$locAbbr == site & cur_grid$City == location, day_column_name] <- 1
        }
      }
        
      # TODO: update target species?
      # Select which cities to observe across thresholds
      if (species %in% c("Bird", "Black bear", "Bobcat", "Coyote", "Domestic cat", 
                         "Lizard", "Mule deer", "Rabbit", "Raccoon", "Rodent", "Squirrel/Chipmunk", 
                         "Striped Skunk", "Virginia opossum")) {
        # Count number of site-seasons where species was observed
        site_seasons <- apply(cur_grid[, grep("^Day_", names(cur_grid))], 1, 
                              function(row) as.integer(any(row == 1, na.rm = TRUE)))
        
        # Put count in correct column format
        species_counts <- rbind(species_counts, 
                                data.frame(city = city,
                                           species = species, 
                                           source = paste0("Conf. ", thresh), 
                                           count = sum(site_seasons)))
      }
    
      # Order and name columns
      complete_grid <- cur_grid %>%
        rename("Species" = "commonName", "SeasonNumber" = "Season", "Site" = "locAbbr") %>%
        select(Species, SeasonNumber, City, Site, starts_with("Day_"))
      
      # Filepath specification
      output_filename <- paste0(folder_path, "OccupancyReport - ", species, " - Seasonal.csv")
      
      # Output file for each species
      write.csv(complete_grid, output_filename, row.names = FALSE)
    }
  }
}

```

```{r}

occupancy_output_path <- "../Output/"

# Plot total number of images at each threshold
ggplot(total_images, aes(x = factor(threshold), y = count, group=1)) + 
  geom_line(color = "steelblue", size = 1) + 
  geom_point(color = "steelblue", size = 3) +
  labs(x = "Threshold", y = "Total Images", title = "Images per Threshold") + 
  theme_minimal()

ggsave(paste0(occupancy_output_path, "images_per_threshold.jpg"), height=5, width=5, units="in", dpi=600)

for (city_name in cities) {
  city_counts <- species_counts %>%
    filter(city == city_name)
  
  # Order by highest to least count of species site-season observations
  species_order <- city_counts %>%
    filter(source == "Conf. 0") %>%
    arrange(desc(count)) %>%
    pull(species)

  # Apply order to species_counts, reverse for plotting
  city_counts <- city_counts %>%
    mutate(species = factor(species, levels = rev(species_order)))

  # Specify order for key
  desired_order <- c("Conf. 0", paste0("Conf. ", sprintf("%.1f", seq(0.1, 0.9, 0.1))), "Conf. 0.99")

  city_counts <- city_counts %>%
    mutate(source = factor(source, levels = rev(desired_order)))  # rev() makes df_human appear at the top
  
  # Plot site-seasons against species across confidence thresholds.
  count_plot <- ggplot(city_counts, aes(x = count, y = species, fill = source)) +
    geom_col(position = position_dodge(width = 0.7)) +
    labs(
      title = "Number of Site-Seasons with Detections by Source",
      x = "Number of Site-Seasons",
      y = "Species",
      fill = "Source"
    ) +
    guides(fill = guide_legend(reverse = TRUE)) +
    theme_bw(base_size = 12) +
    theme(
      plot.background = element_rect(fill = "white", color = "black"),
      axis.title = element_text(color = "black", size = 20),
      axis.text = element_text(color = "black", size = 12),
      legend.position = "top",
      legend.text = element_text(size = 12),
      legend.title = element_text(size = 16),
      plot.title = element_text(hjust = 0.5, size = 24),
      plot.margin = margin(t = 10, r = 20, b = 10, l = 10)
    )
  
  # Save plot
  ggsave(filename = paste0(occupancy_output_path, city_name, "/detection_by_treatment.png"), 
         plot = count_plot, 
         width = 16, 
         height = 16,
         dpi = 300
  )
}

```

```{r}

exclude <- c("Human", "Empty", "Unknown")

filepath = "../"

# TODO: remove day/week & all/middle logic?
# For each confidence threshold
for (threshold in confs) {
  # For each city
  for (location in cities) {
    # For considering by day and by week
    for (data_by in c("Day", "Week")) {
      if (data_by == "Day") {
        next
      }
      # For considering all dates or just middle dates
      for (data_range in c("All", "Middle")) {
        if (data_range == "Middle") {
          next
        }
        # For each species...
        for (name in species_list) {
          # Skip if part of exclude list
          if (name %in% exclude) {
            next
          }
          # Print current process
          print(paste(name, data_range, data_by, location, threshold))
          
          # Get desired data for city, threshold, species
          species_data <- read.csv(paste0(report_output_path, location, 
                                          "/Conf. ", threshold, 
                                          "/OccupancyReport - ", name, " - Seasonal.csv"))
          # Get day columns from data
          day_cols <- grep("^Day_", names(species_data), value = TRUE)
          # If no observations made for species, skip
          if (!any(species_data[, day_cols] == 1, na.rm = TRUE)) {
            message(paste("Skipping", name, "- no detections in any day column"))
            next
          }
          
          # Run "Occupancy Modeling.Rmd" with current species and specification
          tryCatch({
            render(paste0(filepath, "R Files/3_OccupancyModeling.Rmd"), 
                   params = list(path = filepath, sp_name = name, range = data_range, 
                                 by = data_by, city = location, conf = threshold))
          }, error = function(e) {
            message(paste("Error in rendering for", name, "in", location, "with threshold", threshold, ":", e$message))
          })
        }
      }
    }
  }
}

```