This vignette describes how to recreate the seasonal map products found on the eBird Status and Trends pages. First, the vignette will cover how to average the abundance data seasonally, then it will proceed with examples of making abundance maps, aggregating and smoothing data for range maps, and, finally, provide examples of calculating regional summary statistic. Throughout this vignette we will use the simplified example dataset available through the ebirdst_download()
function, which consists of the Yellow-bellied Sapsucker data spatially subset to the state of Michigan. However, these methods can be applied to the full datasets for any of the eBird Status and Trends species.
Let’s begin by loading all the packages we’ll need for this vignette.
library(ebirdst)
library(raster)
library(sf)
library(smoothr)
library(rnaturalearth)
library(dplyr)
library(tidyr)
library(stringr)
library(ggplot2)
# resolve namespace conflicts
<- dplyr::select
select sf_use_s2(FALSE)
Before we do anything else, we’ll need to download and load the example abundance data for Yellow-bellied Sapsucker.
<- ebirdst_download(species = "example_data")
sp_path # load the abundance data
# this automatically labels layers with their dates
<- load_raster(sp_path, "abundance") abd
The abundance data now consist of a RasterStack
object with 52 layers, each corresponding to the relative abundance for a single week.
We’ll also need some additional spatial data (state and country borders, lakes, etc.) to provide context for the maps we’ll make. Natural Earth is the best source for free map data and there’s an associated R package (rnaturalearth
) for accessing the data. However, we’ll work with a prepared subset of the data stored in the ebirdst
GitHub repository. For all the maps in this vignette, we’ll use the Eckert IV projection, a good option for continental scale maps.
<- "+proj=eck4 +lon_0=-90 +x_0=0 +y_0=0 +ellps=WGS84"
proj
# download natural earth data
<- tempfile(fileext = ".gpkg")
temp_gpkg file.path("https://github.com/CornellLabofOrnithology/ebirdst/raw/",
"main/data-raw/ebirdst_gis-data.gpkg") %>%
download.file(temp_gpkg)
# land polygon
<- read_sf(temp_gpkg, layer = "ne_land") %>%
ne_land st_transform(crs = proj) %>%
st_geometry()
# country lines
<- read_sf(temp_gpkg, layer = "ne_country_lines") %>%
ne_country_lines st_transform(crs = proj) %>%
st_geometry()
# state lines
<- read_sf(temp_gpkg, layer = "ne_state_lines") %>%
ne_state_lines st_transform(crs = proj) %>%
st_geometry()
# rivers
<- read_sf(temp_gpkg, layer = "ne_rivers") %>%
ne_rivers st_transform(crs = proj) %>%
st_geometry()
# lakes
<- read_sf(temp_gpkg, layer = "ne_lakes") %>%
ne_lakes st_transform(crs = proj) %>%
st_geometry()
Generally, the seasons for eBird Status and Trends products are defined on a species-specific basis through expert review. For information on the details of defining seasons, please see the seasons section of the FAQ. While it is certainly possible to define your own seasons when making seasonal abundance and range maps, if you want to recreate the products with the same seasons as the website, you’ll need to use the definitions included in the ebirdst_runs
data frame included in this package.
Let’s start by getting the seasonal definitions for Yellow-bellied Sapsucker and transforming the data into a more usable format. For some species, expert review may have indicated that the models are poor in certain seasons. These problematic seasons are identified by missing dates in ebirdst_runs
for the season in question. Although all seasons passed review for Yellow-bellied Sapsucker, for generality I’ll add an additional column (passed
) indicating whether a seasonal model passed review.
# subset to the yellow-bellied sapsucker season definitions
<- filter(ebirdst_runs, species_code == "yebsap")
run_review <- run_review %>%
yebsap_dates # just keep the seasonal definition columns
select(setdiff(matches("(start)|(end)"), matches("year_round"))) %>%
# transpose
gather("label", "date") %>%
# spread data so start and end dates are in separate columns
separate(label, c("season", "start_end"), "_(?=s|e)") %>%
spread(start_end, date) %>%
select(season, start_dt = start, end_dt = end) %>%
filter(season != "resident")
# did the season pass review
<- run_review[paste0(yebsap_dates$season, "_quality")]
quality_rating $pass <- as.integer(quality_rating) > 1
yebsap_dates yebsap_dates
Now, we’ll use these season definitions to assign each of the weekly abundance layers to a season.
# dates for each abundance layer
<- parse_raster_dates(abd)
weeks # assign to seasons
<- rep(NA_character_, length(weeks))
weeks_season for (i in seq_len(nrow(yebsap_dates))) {
<- yebsap_dates[i, ]
s # skip seasona assignment if season failed
if (!s$pass) {
next()
}# handle seasons cross jan 1 separately
if (s$start_dt <= s$end_dt) {
<- weeks >= s$start_dt & weeks <= s$end_dt
in_season else {
} <- weeks >= s$start_dt | weeks <= s$end_dt
in_season
}<- s$season
weeks_season[in_season]
}table(weeks_season)
Finally, let’s average all the weeks within each season to produce seasonal relative abundance rasters.
# drop weeks not assigned to season
<- !is.na(weeks_season)
week_pass <- abd[[which(week_pass)]]
abd <- weeks[week_pass]
weeks <- weeks_season[week_pass]
weeks_season # average over weeks in season
<- function(s) {
mean_season calc(abd[[which(weeks_season == s)]], mean, na.rm = TRUE)
}<- unique(weeks_season)
seasons <- lapply(seasons, mean_season) %>%
abd_season stack() %>%
setNames(seasons)
abd_season
Before we create maps of seasonal relative abundance we need to account for two addition considerations that the online maps make: whether to split pre- and post-breeding migration and whether to explicitly show the year-round distribution as a separate color. For species that use different paths for their two migrations (less than 40% overlap) we show pre-breeding migration (green) and post-breeding migration (yellow) separately.
<- 0.4
migration_threshold <- c("prebreeding_migration", "postbreeding_migration")
mig_seasons if (all(mig_seasons %in% names(abd_season))) {
# identify areas with abundance in only one season
<- abd_season[[mig_seasons]] > 0
abd_nz <- mask(abd_nz[["prebreeding_migration"]],
just_pre "postbreeding_migration"]],
abd_nz[[maskvalue = 1)
<- mask(abd_nz[["postbreeding_migration"]],
just_post "prebreeding_migration"]],
abd_nz[[maskvalue = 1)
# count the number of cells with abundance in only one season
<- cellStats(stack(just_pre, just_post), sum)
n_just <- cellStats(abd_nz, sum)
n_all # is the proportion of one season cells above the 40% threshold
<- max(n_just / n_all, na.rm = TRUE) >= migration_threshold
split_migration else {
} <- FALSE
split_migration
}/ n_all
n_just split_migration
In this case, there is essentially complete overlap between pre- and post-breeding migration, so we won’t split them into separate colors on the map. Next, we’ll calculate the average annual abundance, which we’ll display as a separate color if at least 1% of the range is occupied in all four seasons.
<- 0.01
threshold_yearround # decide whether to show year-round layer
if (nlayers(abd_season) == 4) {
# annual abundance
<- calc(abd, fun = mean, na.rm = TRUE)
abd_yr # mask out cells that aren't occupied year-round
<- calc(abd_season > 0, fun = sum, na.rm = TRUE) == 4
year_round <- mask(abd_yr, year_round, maskvalue = 0)
abd_yr_mask # determine proportion of celss that are occupied year round
<- cellStats(abd_yr_mask > 0, sum)
n_yr <- cellStats(abd_yr > 0, sum)
n_an # only show year round abundance if it's above 1% of range threshold
<- ((n_yr / n_an) >= threshold_yearround)
show_yearround else {
} <- FALSE
show_yearround
} show_yearround
Mapping relative abundance across the full-annual cycle requires a specialized set of color bins in order to ensure the full range of abundance values is effectively displayed. For simplicity, we’ll use quantile bins based on the seasonal abundance, however, you could also use the function calc_bins()
to calculate the optimal break points of bins for Status and Trends abundance data.
<- abd_season %>%
vals getValues() %>%
na.omit() %>%
as.numeric()
<- quantile(vals[vals > 0], seq(0, 1, by = 0.05))
bin_breaks <- signif(bin_breaks[c("5%", "50%", "95%")], 3)
lbls rm(vals)
Now that everything is in place, we can actually make the seasonal relative abundance map! Note that the Status and Trends maps distinguish between regions of zero abundance and regions with no prediction, which are displayed in different shades of gray, and regions with non-zero abundance, shown in color. To account for this, we’ll need to generate a raster layer delineating the region within which predictions were made.
# project the abundance data to mollweide
# use nearest neighbour resampling to preserve true zeros
<- projectRaster(abd_season, crs = proj, method = "ngb")
abd_season_proj # determine spatial extent for plotting
<- calc_full_extent(abd_season_proj)
ext # set the plotting order of the seasons
<- c("postbreeding_migration", "prebreeding_migration",
season_order "nonbreeding", "breeding")
# prediction region, cells with predicted value in at least one week
<- calc(abd_season_proj, mean, na.rm = TRUE)
pred_region # mask to land area
<- st_buffer(ne_land, dist = max(res(pred_region)) / 2)
ne_land_buffer <- mask(pred_region, as_Spatial(ne_land_buffer))
pred_region
# remove zeros from abundance layers
<- subs(abd_season_proj, data.frame(from = 0, to = NA),
abd_no_zero subsWithNA = FALSE)
# set up plot area
par(mar = c(0 , 0, 0, 0))
plot(ne_land, col = "#cfcfcf", border = NA,
xlim = c(ext@xmin, ext@xmax),
ylim = c(ext@ymin, ext@ymax))
# prediction region and explicit zeros
plot(pred_region, col = "#e6e6e6", maxpixels = raster::ncell(pred_region),
legend = FALSE, add = TRUE)
# lakes
plot(ne_lakes, col = "#ffffff", border = "#444444", lwd = 0.5, add = TRUE)
# land border
plot(ne_land, col = NA, border = "#444444", lwd = 0.5, add = TRUE)
# seasonal layer
<- intersect(season_order, names(abd_no_zero))
plot_seasons for (s in plot_seasons) {
# handle splitting of migration seasons into different colors
if (!split_migration && s %in% c("prebreeding_migration",
"postbreeding_migration")) {
<- "migration"
pal_season
else {
} <- s
pal_season
}<- abundance_palette(length(bin_breaks) - 1, pal_season)
pal plot(abd_no_zero[[s]], col = pal, breaks = bin_breaks,
maxpixels = ncell(abd_no_zero[[s]]),
legend = FALSE, add = TRUE)
}# year round
if (show_yearround) {
<- projectRaster(year_round, crs = mollweide, method = "ngb")
year_round_proj plot(year_round_proj,
col = abundance_palette(length(bin_breaks$bins) - 1, "year_round"),
breaks = bin_breaks$bins,
maxpixels = ncell(year_round_proj),
legend = FALSE, add = TRUE)
}# linework
plot(ne_rivers, col = "#ffffff", lwd = 0.75, add = TRUE)
plot(ne_state_lines, col = "#ffffff", lwd = 1.5, add = TRUE)
plot(ne_country_lines, col = "#ffffff", lwd = 2, add = TRUE)
# legends
<- plot_seasons
legend_seasons if (!split_migration) {
%in% c("prebreeding_migration",
legend_seasons[legend_seasons "postbreeding_migration")] <- "migration"
<- unique(legend_seasons)
legend_seasons
}if (show_yearround) {
<- c(legend_seasons, "year_round")
legend_seasons
}# thin out labels
# plot legends
<- seq(0, 1, length.out = length(lbls))
lbl_at for (i in seq_along(legend_seasons)) {
<- abundance_palette(length(bin_breaks) - 1, legend_seasons[i])
pal if (i == 1) {
<- list(at = lbl_at, labels = lbls, line = -1,
axis_args cex.axis = 0.75, lwd = 0)
else {
} <- list(at = lbl_at, labels = rep("", length(lbls)),
axis_args cex.axis = 0.75, lwd = 0)
}<- legend_seasons[i] %>%
legend_title str_replace_all("_", " ") %>%
str_to_title()
::image.plot(zlim = c(0, 1),
fieldslegend.only = TRUE,
breaks = seq(0, 1, length.out = length(bin_breaks)),
col = pal,
smallplot = c(0.05, 0.35, 0.01 + 0.06 * i, 0.03 + 0.06 * i),
horizontal = TRUE,
axis.args = axis_args,
legend.args = list(text = legend_title, side = 3,
cex = 0.9, col = "black", line = 0.1))
}title("Yellow-bellied Sapsucker Relative Abundance",
line = -1, cex.main = 1)
The eBird Status and Trends range maps delineate the boundary of regions with non-zero relative abundance for a given species. We’ll start by aggregating the raster layers to a coarser resolution to speed up processing, then convert the boundaries of non-zero abundance regions to polygons. We’ll also convert the prediction areas to polygons so we can distinguish where the species is predicted to not occur from where it was not modeled.
# aggregate
<- aggregate(abd_season_proj, fact = 3)
abd_season_agg # raster to polygon, one season at a time
<- list()
range <- list()
pred_area for (s in names(abd_season_agg)) {
# range
<- rasterToPolygons(abd_season_agg[[s]],
range[[s]] fun = function(y) {y > 0},
digits = 6) %>%
st_as_sfc() %>%
# combine polygon pieces into a single multipolygon
st_set_precision(1e6) %>%
st_union() %>%
st_sf() %>%
# tag layers with season
mutate(season = s, layer = "range")
# prediction area
<- rasterToPolygons(abd_season_agg[[s]],
pred_area[[s]] fun = function(y) {!is.na(y)},
digits = 6) %>%
st_as_sfc() %>%
# combine polygon pieces into a single multipolygon
st_set_precision(1e6) %>%
st_union() %>%
st_sf() %>%
# tag layers with season
mutate(season = s, layer = "prediction_area")
}# combine the sf objects for all seasons
<- rbind(do.call(rbind, range), do.call(rbind, pred_area))
range row.names(range) <- NULL
print(range)
Converting from raster to polygons frequently yields tiny fragments of polygons and jagged polygon edges, which result in maps that aren’t very aesthetically pleasing. To address this, we’ll apply some algorithms from the smoothr
package to clean up the polygons.
# clean and smooth
<- (1.5 * prod(res(abd_season_agg)))
cell_area <- range %>%
range_smooth # drop fragment polygons smaller than 1.5 times the aggregated cell size
drop_crumbs(threshold = cell_area) %>%
# drop holes in polygons smaller than 1.5 times the aggregated cell size
fill_holes(threshold = cell_area) %>%
# smooth the polygon edges
smooth(method = "ksmooth", smoothness = 2)
# clip zeros to land border, range to buffered land to handle coastal species
<- split(range_smooth, range_smooth$layer)
range_split <- rbind(
range_smooth st_intersection(range_split$range, ne_land_buffer),
st_intersection(range_split$prediction_area, ne_land))
Finally, we can produce a seasonal range map in a similar fashion to the abundance map above.
# range map color palette
<- c(nonbreeding = "#1d6996",
range_palette prebreeding_migration = "#73af48",
breeding = "#cc503e",
postbreeding_migration = "#edad08",
migration = "#edad08",
year_round = "#6f4070")
# set up plot area
par(mar = c(0 , 0, 0, 0))
plot(ne_land, col = "#cfcfcf", border = NA,
xlim = c(ext@xmin, ext@xmax),
ylim = c(ext@ymin, ext@ymax))
# prediction region and explicit zeros
<- filter(range_smooth, layer == "prediction_area") %>%
annual_pred_area st_union()
plot(annual_pred_area, col = "#e6e6e6", border = NA, add = TRUE)
# lakes
plot(ne_lakes, col = "#ffffff", border = "#444444", lwd = 0.5, add = TRUE)
# land border
plot(ne_land, col = NA, border = "#444444", lwd = 0.5, add = TRUE)
# seasonal layer
for (s in intersect(season_order, unique(range_smooth$season))) {
# handle splitting of migration seasons into different colors
if (!split_migration && s %in% c("prebreeding_migration",
"postbreeding_migration")) {
<- "migration"
col_season else {
} <- s
col_season
}<- filter(range_smooth, season == s, layer == "range") %>%
rng_season st_geometry()
plot(rng_season, col = range_palette[col_season], border = NA, add = TRUE)
}# year round
if (show_yearround) {
# find common area between all seasons
<- filter(range_smooth, layer == "range")
range_combined <- range_combined[1, ]
range_yearround <- sf::st_geometry(range_combined)
range_combined for (i in 2:length(range_combined)) {
<- sf::st_intersection(range_yearround, range_combined[i])
range_yearround
}plot(st_geometry(range_yearround),
col = range_palette["year_round"], border = NA,
add = TRUE)
}# linework
plot(ne_rivers, col = "#ffffff", lwd = 0.75, add = TRUE)
plot(ne_state_lines, col = "#ffffff", lwd = 1.5, add = TRUE)
plot(ne_country_lines, col = "#ffffff", lwd = 2, add = TRUE)
# legend
<- rev(range_palette[legend_seasons])
rng_legend names(rng_legend) <- names(rng_legend) %>%
str_replace_all("_", " ") %>%
str_to_title()
legend("bottomleft", legend = names(rng_legend), fill = rng_legend)
title("Yellow-bellied Sapsucker Seasonal Range Map",
line = -1, cex.main = 1)