###########
# Set up #
###########

# Load libraries and define focal year (2022). Most outputs below are filtered and specific to the focal year.

# load packages
library(sf)
library(here)
library(tidyverse)
library(auk)
library(ggthemes)
library(DT)

# custom functions
`%!in%` = Negate(`%in%`)

# resolve namesake conflicts
select <- dplyr::select

# set focal year
focal_year <- 2022

########################
# Load park boundaries #
########################

# Load park boundaries from TRPD spatial database. eBird data will be filtered to only retain records from within these boundaries.

# where to save clean parks sf file
f_parks <- file.path(here("Data", "Wildlife", "eBird", "parks_outlines_ManualEdits_30mBufferDissolve.RDS")) 

if (!file.exists(f_parks)) {

spatial.boundaries.parks <- st_read(here("Data", "Geospatial", "General", "Park_Bndry_Public_ManualEdits_30mBufferDissolve.shp")) %>%
  st_transform(crs = 4326)

parks <- spatial.boundaries.parks
saveRDS(parks, f_parks)} else {

parks <- readRDS(f_parks)
}

st_crs(parks) <- 4326

parks <- parks %>%
    select(LabelName) %>%
    rename(Park_Name = LabelName)

#mapview::mapview(parks, col.regions = "darkgreen", legend = FALSE)

###################
# Load eBird data #
###################

# Raw data were obtained directly from eBird. An account and simple pre-approval are needed to [download the data](https://science.ebird.org/en/use-ebird-data). 

# read in eBird data
ebd <- auk_ebd(file = file.path("S:\\Shared drives\\Blair Lab\\Data_External_General\\eBird\\EBD\\ebd_US-MN_relDec-2022\\ebd_US-MN_relDec-2022.txt"),
               file_sampling = file.path("S:\\Shared drives\\Blair Lab\\Data_External_General\\eBird\\EBD\\ebd_sampling_relDec-2022\\ebd_sampling_relDec-2022.txt"))

# define some attributes of starting EDB version and planned filtering to populate output file names
ebd_version <- "ebd_relDec-2022" # [UPDATE FOR RUN]
filter_region <- "TRPD_parks_manual_30mBuff" # [UPDATE FOR RUN]
filter_sp <- "all_sp" # [UPDATE FOR RUN]
filter_season <- "all_years" # [UPDATE FOR RUN]

# define filters to apply to the EBD
ebd_filters <- ebd %>%
  # select only records from study area
  auk_bbox(parks) #%>%
  # # data from any year but, could also filter
  # auk_date(date = c("2018-01-01", "2022-12-31"))

# define where filtered outputs will live
f_ebd <- file.path(here("Data", "Wildlife", "eBird", paste0(ebd_version,"_",filter_region,"_",filter_sp,"_",filter_season,".txt"))) 

f_sampling <- file.path(here("Data", "Wildlife", "eBird", paste0(ebd_version,"_","checklists","_",filter_region,"_",filter_sp,"_",filter_season,".txt")))

# implement filter, but only run if the files don't already exist
if (!file.exists(f_ebd)) {
  auk_filter(ebd_filters, file = f_ebd, file_sampling = f_sampling)
}

# define where park-point intersection will live (yielding only records within park boundaries)
f_int <- file.path(here("Data", "Wildlife", "eBird", paste0(ebd_version,"_",filter_region,"_",filter_sp,"_",filter_season,"_within_park_boundaries_int.RDS"))) 

# define where park-point filter will live (interim step without the park info)
f_filter <- file.path(here("Data", "Wildlife", "eBird", paste0(ebd_version,"_",filter_region,"_",filter_sp,"_",filter_season,"_within_park_boundaries_filter.RDS"))) 

# read filtered ebd and and clip to park boundaries, but only run if the files don't already exist
if (!file.exists(f_filter)) {
  dat <- read_ebd(f_ebd, unique = FALSE,  rollup = TRUE) %>%
    st_as_sf(coords = c("longitude", "latitude"), crs = 4326) 
  
  dat <- dat %>%
    select(global_unique_identifier, common_name, observation_count, locality, locality_id, locality_type, observation_date, time_observations_started, observer_id, sampling_event_identifier, protocol_type, duration_minutes, effort_distance_km, number_observers, all_species_reported, group_identifier)
  
  dat_parks <- st_filter(dat, parks) %>%
    mutate(Year = format(as.Date(observation_date, format="%Y-%m-%d"),"%Y")) # add column for year

 saveRDS(dat_parks, f_filter)
 
} else {
  dat_parks <- readRDS(f_filter)
}

# dat_parks_int <- st_intersection(dat_parks, parks)

### st_intersection taking too long. this loop is faster....

if (!file.exists(f_int)) {
# convert polygons to multipart polygons (one per park)
parks_mp <- parks %>%
  group_by(Park_Name) %>%
  summarise(geometry = st_combine(geometry)) 
  
dat_parks_int <- data.frame()
for (i in 1:nrow(parks_mp)) {
  park_name <- parks_mp$Park_Name[[i]]
  park<-parks_mp[i,]
  park_dat <- st_filter(dat_parks,park)
  park_dat$Park_Name <- park_name
  dat_parks_int <- rbind(dat_parks_int, park_dat)
}

saveRDS(dat_parks_int, f_int)
} else {
  dat_parks_int <- readRDS(f_int)
}

dat_parks <- dat_parks_int
rm(dat_parks_int)

# save checklist locations with total number as sf for later before dropping geometry
checklist_locations <-dat_parks %>%
  select(locality_id, locality_type, locality, Park_Name) %>%
  unique()

checklist_locations_visit_count <-dat_parks %>%
  select(locality_id, locality_type, locality, sampling_event_identifier, Park_Name) %>%
  unique() %>%
  group_by(locality_id) %>%
  summarise(total_visits = n()) %>%
  st_drop_geometry()

checklist_locations <- left_join(checklist_locations, checklist_locations_visit_count)

# Expand locality type codes
checklist_locations$locality_type[which(checklist_locations$locality_type == "H")] <- "Hotspot"
checklist_locations$locality_type[which(checklist_locations$locality_type == "P")] <- "Personal point"
checklist_locations$locality_type[which(checklist_locations$locality_type == "PC")] <- "Postal code"

dat_parks <- dat_parks %>%
  st_drop_geometry()

###########################
# Generate data summaries #
###########################

# This code summarizes data from the focal year (`r focal_year`) for each observer, park, and species. Visualizations and other displays of these summaries are below.

# Convert count of X to 1
dat_parks$observation_count[which(dat_parks$observation_count == "X")] <- 1

dat_parks$observation_count <- as.numeric(dat_parks$observation_count)

# Focal year summary by observer
observer_summaries <- dat_parks %>%
  st_drop_geometry() %>%
  filter(Year == focal_year) %>%
  group_by(observer_id) %>%
  summarize(total_species = length(unique(common_name)),
            total_birds = sum(as.numeric(observation_count)), # count of individual birds
            total_parks_visited = length(unique(Park_Name)),
            total_visits = length(unique(paste0(Park_Name,"_",observation_date))),
            total_checklists = length(unique(sampling_event_identifier))) %>%
  arrange(desc(total_species)) %>%
  ungroup()

observer_summaries_effort <- dat_parks %>%
  filter(Year == focal_year) %>%
  select(observer_id, sampling_event_identifier, Year, duration_minutes, effort_distance_km) %>%
  unique() %>%
  group_by(observer_id) %>%
  summarize(total_hours_birded = round(sum(duration_minutes, na.rm = TRUE)/60,1),
  total_miles_birded = round(sum(effort_distance_km, na.rm = TRUE)*0.621371,1))

observer_summaries <- left_join(observer_summaries, observer_summaries_effort)

# top observers, focal year
top_observers <- slice_head(observer_summaries, n = 10) %>%
  pull(observer_id)

## Get observer names from eBird webiste
library(rvest)
get_name_from_checklist <- function(sampling_event){
  read_html(paste0("https://ebird.org/checklist/",sampling_event)) %>%
  html_elements("meta") %>%
  html_attr("content") %>%
  .[10]
}

get_name_from_checklist <- Vectorize(get_name_from_checklist)

library(ratelimitr)
get_name_from_checklist_lim <- limit_rate(get_name_from_checklist,  rate(n = 1, period = 6))

f_obs_names <- here("Data", "Wildlife", "eBird", paste0("top_observer_names_",focal_year,".RDS"))
  
observer_names <- dat_parks %>%
  filter(Year == focal_year) %>%
  select(observer_id, sampling_event_identifier) %>%
  group_by(observer_id) %>%
  filter(sampling_event_identifier == first(sampling_event_identifier)) %>%
  slice(1) %>%
  filter(observer_id %in% top_observers)

if (file.exists(f_obs_names)) {
  observer_names_existing <- readRDS(f_obs_names)
  observer_names <- left_join(observer_names, observer_names_existing, by = c("observer_id"))}  

if (!file.exists(f_obs_names)) {
  observer_names <- observer_names %>%
    mutate(observer_name = NA)
  observer_names_existing <- data.frame(matrix(nrow=0, ncol=3))
  colnames(observer_names_existing) <- colnames(observer_names)
  observer_names_existing <- observer_names_existing %>%
    mutate(observer_id = as.character(observer_id),
           observer_name = as.character(observer_name))
  }

observer_names_needed <- filter(observer_names, is.na(observer_name)) %>%
  mutate(observer_name = get_name_from_checklist_lim(sampling_event_identifier)) %>%
  select(observer_id, observer_name)

if (nrow(observer_names_needed) != 0) {
observer_names <- bind_rows(observer_names_existing, observer_names_needed) %>%
 select(observer_id, observer_name)} else {
  observer_names <- observer_names_existing
}

saveRDS(observer_names, f_obs_names)

observer_summaries <- observer_summaries %>%
  left_join(observer_names) %>%
  relocate(observer_name) %>%
  mutate(observer_name = ifelse(observer_name == "Anonymous eBirder", observer_id, observer_name), # replace anonymous people with their ID
         observer = ifelse(!is.na(observer_name), observer_name, observer_id))

observer_summaries$observer <- reorder(observer_summaries$observer, -observer_summaries$total_species)

# Focal year summary by park
park_summaries <- dat_parks %>%
  filter(Year == focal_year) %>%
  group_by(Park_Name) %>%
  summarize(total_species = length(unique(common_name)),
            total_birds = sum(as.numeric(observation_count)), # count of individual birds
            total_visitors = length(unique(observer_id)),
            total_visits = length(unique(paste0(observer_id,"_",observation_date))),
            total_checklists = length(unique(sampling_event_identifier))) %>%
  arrange(desc(total_species))

park_summaries_effort <- dat_parks %>%
  filter(Year == focal_year) %>%
  select(observer_id, sampling_event_identifier, Year, duration_minutes, effort_distance_km, Park_Name) %>%
  unique() %>%
  group_by(Park_Name) %>%
  summarize(total_hours_birded = round(sum(duration_minutes, na.rm = TRUE)/60,1),
  total_miles_birded = round(sum(effort_distance_km, na.rm = TRUE)*0.621371,1))

park_summaries <- left_join(park_summaries, park_summaries_effort)

# Focal year summary by species (max count by park)
species_list <- dat_parks %>%
  filter(Year == focal_year) %>%
  pivot_wider(id_cols = common_name, names_from = Park_Name, values_from = observation_count, values_fn = max) %>%
  left_join(rebird:::tax %>% select(comName, taxonOrder), by = c(common_name = "comName")) %>%
  arrange(taxonOrder) %>%
  select(-taxonOrder)

# Focal year summary by species (summary)
species_summary <- dat_parks %>%
  filter(Year == focal_year) %>%
  group_by(common_name) %>%
  summarize(total_counted = sum(as.numeric(observation_count)),
            high_count = as.numeric(max(observation_count)),
            total_parks = length(unique(Park_Name)),
            total_observers = length(unique(observer_id)),
            total_days_observed = length(unique(observation_date))) %>%
  left_join(rebird:::tax %>% select(comName, taxonOrder), by = c(common_name = "comName")) %>%
arrange(taxonOrder) %>%
  select(-taxonOrder)

# Summaries by year across all parks & observers
year_summaries <- dat_parks %>%
  group_by(Year) %>%
  summarize(total_species = length(unique(common_name)),
            total_birds = sum(as.numeric(observation_count)), # count of individual birds
            total_visitors = length(unique(observer_id)),
            total_visits = length(unique(paste0(observer_id,"_",Park_Name,"_",observation_date))),
            total_checklists = length(unique(sampling_event_identifier)))

year_summaries_effort <- dat_parks %>%
  select(observer_id, sampling_event_identifier, Year, duration_minutes, effort_distance_km, Park_Name) %>%
  unique() %>%
  group_by(Year) %>%
  summarize(total_hours_birded = round(sum(duration_minutes, na.rm = TRUE)/60,1),
  total_miles_birded = round(sum(effort_distance_km, na.rm = TRUE)*0.621371,1))

year_summaries <- left_join(year_summaries, year_summaries_effort) %>%
  arrange(desc(Year))

# Species list for "winning" observer
winning_observer <- observer_summaries$observer_id[which(observer_summaries$total_species == max(observer_summaries$total_species))]

focal_observer <- winning_observer
observer_species_list <- dat_parks %>%
  filter(observer_id == winning_observer) %>%
  filter(Year == focal_year) %>%
  pivot_wider(id_cols = common_name, names_from = Park_Name, values_from = observation_count, values_fn = max) %>%
  left_join(rebird:::tax %>% select(comName, taxonOrder), by = c(common_name = "comName")) %>%
  arrange(taxonOrder) %>%
  select(-taxonOrder)
#observer_species_list[is.na(observer_species_list)] <- ""

dat_parks_focal_year <- dat_parks %>% filter(Year == focal_year)

# biggest day
biggest_days <- dat_parks_focal_year %>%
  group_by(observer_id, observation_date) %>%
  summarize(n_sp = length(unique(common_name)),
            length(unique(Park_Name))) %>%
  arrange(desc(n_sp))

# largest counts ever
high_counts <- dat_parks %>%
  slice_max(observation_count, n = 20)

# new high counts
new_high_counts <- dat_parks %>%
  group_by(common_name, Year) %>%
  summarize(high_count = max(observation_count)) %>%
  filter(high_count == max(high_count)) %>%
  filter(Year == min(Year)) %>%
  filter(Year == focal_year) %>%
  left_join(rebird:::tax %>% select(comName, taxonOrder), by = c(common_name = "comName")) %>%
arrange(taxonOrder) %>%
  select(-taxonOrder)

# SGCN
# determine SGCN sp in system
SGCN <- read_csv(here("Data", "Wildlife", "sgcn_2015-25_list_excel.csv")) %>%
  filter(Taxa == "Birds")
SGCN$SGCN <- 1 # make new column indicating species is on list
# align with eBird common names
SGCN$ComName[which(SGCN$ComName == "Le Conte's Sparrow")] <- "LeConte's Sparrow"
SGCN$ComName[which(SGCN$ComName == "Nelson's  Sparrow")] <- "Nelson's Sparrow"
SGCN$ComName[which(SGCN$ComName == "Rufa Red Knot")] <- "Red Knot"
SGCN$ComName[which(SGCN$ComName == "Belted kingfisher")] <- "Belted Kingfisher"
SGCN$ComName[which(SGCN$ComName == "Black-crowned Night-heron")] <- "Black-crowned Night-Heron"
SGCN$ComName[which(SGCN$ComName == "Greater Prairie-chicken")] <- "Greater Prairie-Chicken"
SGCN <- SGCN %>% 
  left_join(rebird:::tax %>% select(comName, taxonOrder), by = c(ComName = "comName")) %>%
  rename(common_name = "ComName")

SGCN_records <- dat_parks_focal_year %>%
  filter(common_name %in% unique(SGCN$common_name)) %>%
  group_by(common_name) %>%
  summarize(total_counted = sum(as.numeric(observation_count)),
            high_count = as.numeric(max(observation_count)),
            total_parks = length(unique(Park_Name)),
            total_observers = length(unique(observer_id)),
            total_days_observed = length(unique(observation_date))) %>%
  left_join(SGCN)

# how many new birders in focal year?
obs_focal_year <- dat_parks %>%
  filter(Year == focal_year)
obs_focal_year <- unique(obs_focal_year$observer_id)
  
obs_pre_focal_year <- dat_parks %>%
  filter(Year != focal_year)
obs_pre_focal_year <- unique(obs_pre_focal_year$observer_id)
  
obs_new_focal_year <- setdiff(obs_focal_year, obs_pre_focal_year)