## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(warning = FALSE, message = FALSE, collapse = TRUE, comment = "#>", out.width = "100%", fig.height = 4, fig.width = 7, fig.align = "center") # only build vignettes locally and not for R CMD check knitr::opts_chunk$set(eval = nzchar(Sys.getenv("BUILD_VIGNETTES"))) ## ----packages----------------------------------------------------------------- # library(dplyr) # library(ebirdst) # library(fields) # library(ggplot2) # library(lubridate) # library(rnaturalearth) # library(scico) # library(sf) # library(terra) # library(tidyr) # extract <- terra::extract ## ----map-load----------------------------------------------------------------- # # download seasonal relative abundance data # ebirdst_download_status("wesmea", # pattern = "abundance_seasonal_mean") # # # load seasonal mean relative abundance at 3km resolution # abd_seasonal <- load_raster("wesmea", # product = "abundance", # period = "seasonal", # metric = "mean", # resolution = "3km") # # # extract just the breeding season relative abundance # abd_breeding <- abd_seasonal[["breeding"]] ## ----map-simple, echo=-1------------------------------------------------------ # par(mar = c(0.25, 0.25, 0.25, 2)) # plot(abd_breeding, axes = FALSE) ## ----map-extent, echo=-1------------------------------------------------------ # par(mar = c(0.25, 0.25, 0.25, 0.25)) # # region boundary # region_boundary <- ne_states(iso_a2 = "US") |> # filter(name == "Montana") # # # project boundary to match raster data # region_boundary_proj <- st_transform(region_boundary, st_crs(abd_breeding)) # # # crop and mask to boundary of montana # abd_breeding_mask <- crop(abd_breeding, region_boundary_proj) |> # mask(region_boundary_proj) # # # map the cropped data # plot(abd_breeding_mask, axes = FALSE) ## ----map-projection, echo=-1-------------------------------------------------- # par(mar = c(0.25, 0.25, 0.25, 2)) # # find the centroid of the region # region_centroid <- region_boundary |> # st_geometry() |> # st_transform(crs = 4326) |> # st_centroid() |> # st_coordinates() |> # round(1) # # # define projection # crs_laea <- paste0("+proj=laea +lat_0=", region_centroid[2], # " +lon_0=", region_centroid[1]) # # # transform to the custom projection using nearest neighbor resampling # abd_breeding_laea <- project(abd_breeding_mask, crs_laea, method = "near") |> # # remove areas of the raster containing no data # trim() # # # map the cropped and projected data # plot(abd_breeding_laea, axes = FALSE, breakby = "cases") ## ----map-bins, echo=-1-------------------------------------------------------- # par(mar = c(0.25, 0.25, 0.25, 2)) # # quantiles of non-zero values # v <- values(abd_breeding_laea, na.rm = TRUE, mat = FALSE) # v <- v[v > 0] # breaks <- quantile(v, seq(0, 1, by = 0.1)) # # add a bin for 0 # breaks <- c(0, breaks) # # # status and trends palette # pal <- ebirdst_palettes(length(breaks) - 2) # # add a color for zero # pal <- c("#e6e6e6", pal) # # # map using the quantile bins # plot(abd_breeding_laea, breaks = breaks, col = pal, axes = FALSE) ## ----map-basemap, echo=-1----------------------------------------------------- # par(mar = c(0.25, 0.25, 0.25, 0.25)) # # natural earth boundaries # countries <- ne_countries(returnclass = "sf") |> # st_geometry() |> # st_transform(crs_laea) # states <- ne_states(iso_a2 = "US") |> # st_geometry() |> # st_transform(crs_laea) # # # define the map plotting extent with the region boundary polygon # region_boundary_laea <- region_boundary |> # st_geometry() |> # st_transform(crs_laea) # plot(region_boundary_laea) # # add basemap # plot(countries, col = "#cfcfcf", border = "#888888", add = TRUE) # # add relative abundance # plot(abd_breeding_laea, # breaks = breaks, col = pal, # maxcell = ncell(abd_breeding_laea), # legend = FALSE, add = TRUE) # # add boundaries # lines(vect(countries), col = "#ffffff", lwd = 3) # lines(vect(states), col = "#ffffff", lwd = 1.5, xpd = TRUE) # lines(vect(region_boundary_laea), col = "#ffffff", lwd = 3, xpd = TRUE) # # # add legend using the fields package # # label the bottom, middle, and top # labels <- quantile(breaks, c(0, 0.5, 1)) # label_breaks <- seq(0, 1, length.out = length(breaks)) # image.plot(zlim = c(0, 1), breaks = label_breaks, col = pal, # smallplot = c(0.90, 0.93, 0.15, 0.85), # legend.only = TRUE, # axis.args = list(at = c(0, 0.5, 1), # labels = round(labels, 2), # col.axis = "black", fg = NA, # cex.axis = 0.9, lwd.ticks = 0, # line = -0.5)) ## ----chron-------------------------------------------------------------------- # region_boundary <- ne_states(iso_a2 = "US") |> # filter(name == "Montana") ## ----chron-single-dl---------------------------------------------------------- # # download data if they haven't already been downloaded # # only weekly 3km relative abundance, median and confidence limits # ebirdst_download_status("Western Meadowlark", # pattern = "abundance_(median|upper|lower)_3km") # # # load the median weekly relative abundance and lower/upper confidence limits # abd_median <- load_raster("wesmea", product = "abundance", metric = "median") # abd_lower <- load_raster("wesmea", product = "abundance", metric = "lower") # abd_upper <- load_raster("wesmea", product = "abundance", metric = "upper") # # # project region boundary to match raster data # region_boundary_proj <- st_transform(region_boundary, st_crs(abd_median)) ## ----chron-single-region------------------------------------------------------ # # extract values within region and calculate the mean # abd_median_region <- extract(abd_median, region_boundary_proj, # fun = "mean", na.rm = TRUE, ID = FALSE) # abd_lower_region <- extract(abd_lower, region_boundary_proj, # fun = "mean", na.rm = TRUE, ID = FALSE) # abd_upper_region <- extract(abd_upper, region_boundary_proj, # fun = "mean", na.rm = TRUE, ID = FALSE) # # # transform to data frame format with rows corresponding to weeks # chronology <- data.frame(week = as.Date(names(abd_median)), # median = as.numeric(abd_median_region), # lower = as.numeric(abd_lower_region), # upper = as.numeric(abd_upper_region)) ## ----chron-single-chart------------------------------------------------------- # ggplot(chronology) + # aes(x = week, y = median) + # geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) + # geom_line() + # scale_x_date(date_labels = "%b", date_breaks = "1 month") + # labs(x = "Week", # y = "Mean relative abundance in Montana", # title = "Migration chronology for Western Meadowlark in Montana") ## ----chron-multi-chron-------------------------------------------------------- # grassland_species <- c("Baird's Sparrow", # "Bobolink", # "Chestnut-collared Longspur", # "Sprague's Pipit", # "Upland Sandpiper", # "Western Meadowlark") # # chronologies <- NULL # for (species in grassland_species) { # # download weekly 27km relative abundance, median and confidence limits # ebirdst_download_status(species, # pattern = "abundance_(median|upper|lower)_3km") # # # load the median weekly relative abundance and lower/upper confidence limits # abd_median <- load_raster(species) # abd_lower <- load_raster(species, metric = "lower") # abd_upper <- load_raster(species, metric = "upper") # # # total relative abundance across the entire modeled range of the species # abd_total <- global(abd_median, fun = sum, na.rm = TRUE)$sum # # # total abundance within the region of interest # abd_median_region <- extract(abd_median, region_boundary_proj, # fun = "sum", na.rm = TRUE, ID = FALSE) # abd_lower_region <- extract(abd_lower, region_boundary_proj, # fun = "sum", na.rm = TRUE, ID = FALSE) # abd_upper_region <- extract(abd_upper, region_boundary_proj, # fun = "sum", na.rm = TRUE, ID = FALSE) # # # proportion of population within the region of interest # prop_pop_median <- as.numeric(abd_median_region) / abd_total # prop_pop_lower <- as.numeric(abd_lower_region) / abd_total # prop_pop_upper <- as.numeric(abd_upper_region) / abd_total # # # transform to data frame format with rows corresponding to weeks # chronology <- data.frame(species = species, # week = as.Date(names(abd_median)), # median = prop_pop_median, # lower = prop_pop_lower, # upper = pmin(prop_pop_upper, 1)) # # # combine with other species # chronologies <- bind_rows(chronologies, chronology) # } ## ----chron-multi-chart-------------------------------------------------------- # ggplot(chronologies) + # aes(x = week, y = median, color = species, fill = species) + # geom_ribbon(aes(ymin = lower, ymax = upper), color = NA, alpha = 0.2) + # geom_line(linewidth = 1) + # scale_x_date(date_labels = "%b", date_breaks = "1 month") + # scale_y_continuous(labels = scales::label_percent()) + # scale_color_brewer(palette = "Set1") + # scale_fill_brewer(palette = "Set1") + # labs(x = NULL, # y = "Percent of population in Montana", # title = "Migration chronologies for grassland birds in Montana", # color = NULL, fill = NULL) + # theme(legend.position = "bottom") ## ----stats-dl----------------------------------------------------------------- # ebirdst_download_status("Golden Eagle") ## ----stats-seasonal-load------------------------------------------------------ # # seasonal proportion of population # prop_pop_seasonal <- load_raster("goleag", # product = "proportion-population", # period = "seasonal") # # # state boundaries, excluding hawaii # states <- ne_states(iso_a2 = "US") |> # filter(name != "Hawaii") |> # select(state = name) |> # # transform to match projection of raster data # st_transform(crs = st_crs(prop_pop_seasonal)) ## ----stats-seasonal-extract--------------------------------------------------- # state_prop_pop <- extract(prop_pop_seasonal, states, # fun = "sum", na.rm = TRUE, weights = TRUE, # bind = TRUE) |> # as.data.frame() |> # # sort in descending order or breeding proportion of population # arrange(desc(breeding)) # head(state_prop_pop) ## ----stats-relative-prep------------------------------------------------------ # # seasonal relative abundance # abd_seasonal <- load_raster("goleag", # product = "abundance", # period = "seasonal") # # # load country polygon, union into a single polygon, and project # noram <- ne_countries(country = c("United States of America", # "Canada", "Mexico")) |> # st_union() |> # st_transform(crs = st_crs(abd_seasonal)) |> # # vect converts an sf object to terra format for mask() # vect() # # # mask seasonal abundance # abd_seasonal_noram <- mask(abd_seasonal, noram) # # # total north american relative abundance for each season # abd_noram_total <- global(abd_seasonal_noram, fun = "sum", na.rm = TRUE) # # # proportion of north american population # prop_pop_noram <- abd_seasonal_noram / abd_noram_total$sum ## ----stats-relative-calc------------------------------------------------------ # state_prop_noram_pop <- extract(prop_pop_noram, states, # fun = "sum", na.rm = TRUE, weights = TRUE, # bind = TRUE) |> # as.data.frame() |> # # sort in descending order or breeding proportion of population # arrange(desc(breeding)) # head(state_prop_noram_pop) ## ----stats-custom-calc-------------------------------------------------------- # # weekly relative abundance, masked to north america # abd_weekly_noram <- load_raster("goleag", # product = "abundance", # resolution = "27km") |> # mask(noram) # # # total north american relative abundance for each week # abd_weekly_total <- global(abd_weekly_noram, fun = "sum", na.rm = TRUE) # # # proportion of north american population # prop_pop_weekly_noram <- abd_weekly_noram / abd_weekly_total$sum # # # proportion of weekly population in california # california <- filter(states, state == "California") # cali_prop_noram_pop <- extract(prop_pop_weekly_noram, california, # fun = "sum", na.rm = TRUE, # weights = TRUE, ID = FALSE) # prop_pop_weekly_noram <- data.frame( # week = as.Date(names(cali_prop_noram_pop)), # prop_pop = as.numeric(cali_prop_noram_pop[1, ])) # head(prop_pop_weekly_noram) ## ----stats-custom-month------------------------------------------------------- # prop_pop_weekly_noram |> # filter(month(week) == 1) |> # summarize(prop_pop = mean(prop_pop)) ## ----stats-coastal-wrong------------------------------------------------------ # # download only the season proportion of population layer # ebirdst_download_status("Surf Scoter", # pattern = "proportion-population_seasonal_mean_3km") # # # breeding season proportion of population # abd_nonbreeding <- load_raster("Surf Scoter", # product = "proportion-population", # period = "seasonal") |> # subset("nonbreeding") # # # load a polygon for the boundary of Mexico # mexico <- ne_countries(country = "Mexico") |> # st_transform(crs = st_crs(abd_nonbreeding)) # # # proportion in mexico # extract(abd_nonbreeding, mexico, fun = "sum", na.rm = TRUE, ID = FALSE) ## ----stats-coastal-buffer----------------------------------------------------- # # buffer by 5000m = 5km # mexico_buffer <- st_buffer(mexico, dist = 5000) # # # proportion in mexico # extract(abd_nonbreeding, mexico_buffer, fun = "sum", na.rm = TRUE, # touches = TRUE, ID = FALSE) ## ----aoi-species-------------------------------------------------------------- # # species list # grassland_species <- c("Baird's Sparrow", # "Bobolink", # "Chestnut-collared Longspur", # "Sprague's Pipit", # "Upland Sandpiper", # "Western Meadowlark") # # # region boundary # region_boundary <- ne_states(iso_a2 = "US") |> # filter(name == "Montana") |> # st_transform(st_crs(abd_breeding)) |> # vect() ## ----aoi-richness------------------------------------------------------------- # range_rasters <- list() # for (species in grassland_species) { # # download seasonal abundance at 3km # ebirdst_download_status(species, pattern = "abundance_seasonal_mean_3km") # # # load breeding season relative abundance # abd <- load_raster(species, period = "seasonal") |> # subset("breeding") # # crop and mask to region # abd_masked <- mask(crop(abd, region_boundary), region_boundary) # # convert to binary, presence-absence # range_rasters[[species]] <- abd_masked > 0 # } # # sum across species to calculate richness # richness <- sum(rast(range_rasters), na.rm = TRUE) ## ----aoi-richness-map--------------------------------------------------------- # # make a simple map # plot(richness, axes = FALSE) ## ----aoi-importance-pop------------------------------------------------------- # prop_pop <- list() # for (species in grassland_species) { # # download seasonal abundance at 3km # ebirdst_download_status(species, # pattern = "proportion-population_seasonal_mean_3km") # # # load breeding season proportion of population # pp <- load_raster(species, # product = "proportion-population", # period = "seasonal") |> # subset("breeding") # # crop and mask to region # prop_pop[[species]] <- mask(crop(pp, region_boundary), region_boundary) # } # # take mean across species # importance <- mean(rast(prop_pop), na.rm = TRUE) ## ----aoi-importance-simple---------------------------------------------------- # plot(importance, axes = FALSE) ## ----aoi-importance-clean----------------------------------------------------- # # drop zeros # importance <- ifel(importance == 0, NA, importance) # # drop anything below the median # cutoff <- global(importance, quantile, probs = 0.5, na.rm = TRUE) |> # as.numeric() # importance <- ifel(importance > cutoff, importance, NA) # # make a simple map # plot(importance, axes = FALSE) # plot(region_boundary, col = "grey", axes = FALSE, add = TRUE) # plot(importance, axes = FALSE, legend = FALSE, add = TRUE) ## ----aoi-importance-nice------------------------------------------------------ # # reproject # importance_proj <- trim(project(importance, crs_laea)) # region_boundary_proj <- project(region_boundary, crs_laea) # # basemap # par(mar = c(0, 0, 0, 0)) # plot(region_boundary_proj, col = "grey", axes = FALSE, # main = "Areas of importance for grassland birds in Montana") # # add importance raster # plot(importance_proj, legend = FALSE, add = TRUE) # # add legend # fields::image.plot(zlim = c(0, 1), legend.only = TRUE, # col = viridisLite::viridis(100), # breaks = seq(0, 1, length.out = 101), # smallplot = c(0.15, 0.85, 0.12, 0.15), # horizontal = TRUE, # axis.args = list(at = c(0, 0.5, 1), # labels = c("Low", "Medium", "High"), # fg = "black", col.axis = "black", # cex.axis = 0.75, lwd.ticks = 0.5, # padj = -1.5), # legend.args = list(text = "Relative Importance", # side = 3, col = "black", # cex = 1, line = 0)) ## ----ppms-quality------------------------------------------------------------- # horlar_review <- filter(ebirdst_runs, species_code == "horlar") |> # select(breeding_quality, breeding_start, breeding_end) # print(horlar_review) ## ----ppms-dl------------------------------------------------------------------ # # download and load ppm # ebirdst_download_status("horlar", download_ppms = TRUE) # bernoulli_dev <- load_ppm("horlar", ppm = "occ_bernoulli_dev") # print(bernoulli_dev) ## ----ppms-subset, echo=-1----------------------------------------------------- # par(mar = c(0, 0, 0, 0)) # # subset to weeks in breeding season and average # breeding_dates <- c(horlar_review$breeding_start, horlar_review$breeding_end) |> # format("%m-%d") # in_breeding <- names(bernoulli_dev) >= breeding_dates[1] & # names(bernoulli_dev) <= breeding_dates[2] # bernoulli_dev_breeding <- mean(bernoulli_dev[[in_breeding]], na.rm = TRUE) # # # mask to just canada and the united states # us_ca <- ne_countries(country = c("United States of America", "Canada")) |> # st_transform(st_crs(bernoulli_dev_breeding)) # bernoulli_dev_breeding_us_ca <- bernoulli_dev_breeding |> # crop(us_ca) |> # mask(us_ca) |> # trim() # # # make a map # ppm_cols <- rev(scico(100, palette = "vik")) # max_val <- global(abs(bernoulli_dev_breeding_us_ca), fun = max, na.rm = TRUE) |> # as.numeric() # plot(bernoulli_dev_breeding_us_ca, # range = c(-max_val, max_val), # col = ppm_cols, # axes = FALSE, box = TRUE) # plot(st_geometry(us_ca), add = TRUE)