Examples from the Advanced eBird Status Data Products in R workshop

aos-bc-2022_examples.R

library(raster)
library(sf)
library(exactextractr)
library(rnaturalearth)
library(tidyverse)
library(lubridate)
library(ebirdst)

# averaging and subsetting  ----

# load prairie warbler weekly relative abundance
prawar_abd <- get_species_path("prawar") %>%
  load_raster(product = "abundance", resolution = "mr")

# existing non-breeding season dates
ebirdst_runs %>%
  filter(species_code == "prawar") %>%
  select(nonbreeding_start, nonbreeding_end, nonbreeding_quality)
# define your own dates
start_date <- as.Date("2020-01-01")
end_date <- as.Date("2020-02-28")
# dates for each raster layer
week_dates <- parse_raster_dates(prawar_abd)
selected_weeks <- which(week_dates >= start_date & week_dates <= end_date)

# subset raster cube to these weeks and average
abd_breeding <- mean(prawar_abd[[selected_weeks]], na.rm = TRUE)

# subset to a focal region
pr_boundary <- ne_countries(scale = 50,
                            country = "Puerto Rico",
                            returnclass = "sf") %>%
  st_geometry() %>%
  st_transform(st_crs(abd_breeding)) %>%
  as_Spatial()
abd_breeding_pr <- crop(abd_breeding, pr_boundary)
plot(abd_breeding_pr)


# multi-species download and trajectories ----

# download just low resolution seasonal % of population for several species
species <- c("American Avocet", "Hudsonian Godwit", "Willet")
species_codes <- get_species(species)
# look at available files without downloading
ebirdst_download(species_codes[1], dry_run = TRUE)
# we want the ones ending in "percent-population_median_lr_2020.tif"
ebirdst_download(species_codes[1],
                 pattern = "percent-population_median_lr_2020.tif",
                 dry_run = TRUE)
# looks good, let's download for all species
dl_paths <- map_chr(species_codes, ebirdst_download,
                    pattern = "percent-population_median_lr_2020.tif")
names(dl_paths) <- species_codes
# boundary of kansas
kansas_boundary <- ne_states(iso_a2 = "US", returnclass = "sf") %>% 
  filter(postal == "KS") %>% 
  st_geometry()

# summarize % of pop over kansas boundary
pct_pop_all_species <- NULL
for (s in species_code) {
  pct_pop <- load_raster(dl_paths[s],
                         product = "percent-population",
                         resolution = "lr")
  pct_pop_region <- exact_extract(pct_pop, kansas_boundary, fun = "sum")
  pct_pop_all_species <- pct_pop_region %>% 
    pivot_longer(cols = everything(),
                 names_to = "week", 
                 names_prefix = "sum.w",
                 values_to = "pct_pop") %>% 
    mutate(species_code = s) %>% 
    bind_rows(pct_pop_all_species, .)
}
pct_pop_all_species$week <- ymd(pct_pop_all_species$week)
ggplot(pct_pop_all_species,
       aes(x = week, y = pct_pop, color = species_code)) +
  geom_line()