grosbeak.mapping

0.1 Setup

Code

library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

Code

library(ebirdst)

Please cite the eBird Status and Trends data using: 
  Fink, D., T. Auer, A. Johnston, M. Strimas-Mackey, S. Ligocki, O. Robinson, 
  W. Hochachka, L. Jaromczyk, C. Crowley, K. Dunham, A. Stillman, I. Davies, 
  A. Rodewald, V. Ruiz-Gutierrez, C. Wood. 2023.
  eBird Status and Trends, Data Version: 2022; Released: 2023. Cornell Lab of
  Ornithology, Ithaca, New York. https://doi.org/10.2173/ebirdst.2022

This version of the package provides access to the 2022 version of the eBird
Status and Trends Data Products. Access to the 2022 data will be provided 
until November 2024 when it will be replaced by the 2023 data. At that 
point, you will be required to update this R package and transition to using 
the new data.

Code

library(fields)

Warning: package 'fields' was built under R version 4.3.3

Loading required package: spam

Warning: package 'spam' was built under R version 4.3.3

Spam version 2.11-0 (2024-10-03) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction 
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

Attaching package: 'spam'

The following objects are masked from 'package:base':

    backsolve, forwardsolve

Loading required package: viridisLite

Try help(fields) to get started.

Code

library(ggplot2)
library(lubridate)

Attaching package: 'lubridate'

The following objects are masked from 'package:base':

    date, intersect, setdiff, union

Code

library(rnaturalearth)
library(sf)

Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

Code

library(terra)

terra 1.7.55

Attaching package: 'terra'

The following object is masked from 'package:fields':

    describe

Code

library(tidyr)

Attaching package: 'tidyr'

The following object is masked from 'package:terra':

    extract

Code

library(tidyterra)

Warning: package 'tidyterra' was built under R version 4.3.3

Attaching package: 'tidyterra'

The following object is masked from 'package:stats':

    filter

Code

library(RColorBrewer)
#extract <- terra::extract
#devtools::install_github("ropensci/rnaturalearthhires")

0.2 Download data

Code

#set_ebirdst_access_key("d3q8inorcsog", overwrite = TRUE)
#download the Marsh Wren data
#path <- ebirdst_download_status("Rose-breasted Grosbeak", download_ranges = T)
#path <- ebirdst_download_status("Black-headed Grosbeak", download_ranges = T)
#above only needs to be run once
#if the data package has already been downloaded just do this:
#path <- get_species_path("Marsh Wren")

#read in sample sheet with sample sex info
samps<-read.csv("~/Desktop/grosbeak.data.csv")
samps<-samps[samps$passed.genomic.filtering == "TRUE",] #retain only samples that passed filtering
table(samps$state)

    Colorado       Kansas South Dakota 
           8            9          121 

Code

table(samps$site)

 0  1  2  3  4  5  6  7  8  9 10 11 12 
 8 19  3 16 10  4  8  9 10  9 10 23  9 

#pull in seasonal abundance data for RBGB and check it out

Code

#print the dates for the breeding season for our two focal species
ebirdst_runs[ebirdst_runs$scientific_name == "Pheucticus melanocephalus" | ebirdst_runs$scientific_name == "Pheucticus ludovicianus",]

# A tibble: 2 × 28
  species_code scientific_name          common_name is_resident breeding_quality
  <chr>        <chr>                    <chr>       <lgl>       <chr>           
1 bkhgro       Pheucticus melanocephal… Black-head… FALSE       3               
2 robgro       Pheucticus ludovicianus  Rose-breas… FALSE       3               
# ℹ 23 more variables: breeding_start <date>, breeding_end <date>,
#   nonbreeding_quality <chr>, nonbreeding_start <date>,
#   nonbreeding_end <date>, postbreeding_migration_quality <chr>,
#   postbreeding_migration_start <date>, postbreeding_migration_end <date>,
#   prebreeding_migration_quality <chr>, prebreeding_migration_start <date>,
#   prebreeding_migration_end <date>, resident_quality <chr>,
#   resident_start <date>, resident_end <date>, has_trends <lgl>, …

Code

# load seasonal mean relative abundance at low res
abd_seasonal <- load_raster("Rose-breasted Grosbeak", 
                            product = "abundance", 
                            period = "seasonal",
                            metric = "mean",
                            resolution = "3km")
#extract just the breeding season relative abundance
abd_breeding <- abd_seasonal[["breeding"]]
plot(abd_breeding, axes = FALSE)

Code

# boundaries of states in the united states
region_boundary1 <- ne_states(iso_a2 = "US", returnclass = "sf") %>%
  filter(iso_a2 == "US", !postal %in% c("AK", "HI")) %>%
  transmute(state = iso_3166_2)

region_boundary <- ne_countries(continent = "North America")
#plot(region_boundary)

region_boundary_proj <- st_transform(region_boundary, st_crs(abd_breeding))
# crop and mask to boundary
abd_breeding_mask <- crop(abd_breeding, region_boundary_proj) |> mask(region_boundary_proj)
# map the cropped data
plot(abd_breeding_mask, axes = FALSE)

Code

# find the centroid of the region
region_centroid <- region_boundary |> 
  st_geometry() |> 
  st_transform(crs = 4326) |> 
  st_centroid() |> 
  st_coordinates() |> 
  round(1)

# define projection centered on Kansas
crs_laea <- paste0("+proj=laea +lat_0=", "38.5", " +lon_0=", "-98.4")

# transform to the custom projection using nearest neighbor resampling
abd_breeding_laea <- project(abd_breeding_mask, crs_laea) |> trim()

|---------|---------|---------|---------|
=========================================
                                          

Code

# map the cropped and projected data
plot(abd_breeding_laea, axes = FALSE, breakby = "cases")

#repeat for BHGB and check it out

Code

# load seasonal mean relative abundance at low res
BHabd_seasonal <- load_raster("Black-headed Grosbeak", 
                            product = "abundance", 
                            period = "seasonal",
                            metric = "mean",
                            resolution = "3km")
#extract just the breeding season relative abundance
BHabd_breeding <- BHabd_seasonal[["breeding"]]
plot(BHabd_breeding, axes = FALSE)

Code

# boundaries of states in the united states
region_boundary1 <- ne_states(iso_a2 = "US", returnclass = "sf") %>%
  filter(iso_a2 == "US", !postal %in% c("AK", "HI")) %>%
  transmute(state = iso_3166_2)

region_boundary <- ne_countries(continent = "North America")
#plot(region_boundary)

BHregion_boundary_proj <- st_transform(region_boundary, st_crs(BHabd_breeding))
# crop and mask to boundary
BHabd_breeding_mask <- crop(BHabd_breeding, BHregion_boundary_proj) |> mask(BHregion_boundary_proj)
# map the cropped data
plot(BHabd_breeding_mask, axes = FALSE)

Code

# find the centroid of the region
region_centroid <- region_boundary |> 
  st_geometry() |> 
  st_transform(crs = 4326) |> 
  st_centroid() |> 
  st_coordinates() |> 
  round(1)

# define projection centered on Kansas
crs_laea <- paste0("+proj=laea +lat_0=", "38.5", " +lon_0=", "-98.4")

# transform to the custom projection using nearest neighbor resampling
BHabd_breeding_laea <- project(BHabd_breeding_mask, crs_laea) |> trim()

|---------|---------|---------|---------|
=========================================
                                          

Code

# map the cropped and projected data
plot(BHabd_breeding_laea, axes = FALSE, breakby = "cases")

1 plot them together

Code

# transform to the custom projection using nearest neighbor resampling
BHabd_breeding_laea <- project(BHabd_breeding_mask, crs_laea) |> trim()

|---------|---------|---------|---------|
=========================================
                                          

Code

#remove 0 cells from BH raster which would overwrite the presence data from the rose-breasted
y <- classify(BHabd_breeding_laea, cbind(-Inf, 0, NA))

#remove 0 cells from RB raster 
abd_breeding_laea <- classify(abd_breeding_laea, cbind(-Inf, 0, NA))

# define the map plotting extent with the region boundary polygon
region_boundary_laea <- region_boundary |> 
  st_geometry() |> 
  st_transform(crs_laea)
# 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)
# 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)

#color-code separately
plot(region_boundary_laea, border = "black")
#plot(countries, col = "grey90", border = "grey90", add = TRUE)
# add relative abundance
plot(abd_breeding_laea,
     breaks = breaks, col = c("#e6e6e6",brewer.pal(9, 'Reds'),"#450009"), 
     maxcell = ncell(abd_breeding_laea),
     legend = FALSE, add = TRUE)
plot(y, breaks = breaks, col = c("#e6e6e6",brewer.pal(9, 'Oranges'),"#682003"), maxcell = ncell(abd_breeding_laea),legend = FALSE, add=TRUE)

Code

#color schemes shown below
# blues: "#e6e6e6" "#F7FBFF" "#DEEBF7" "#C6DBEF" "#9ECAE1" "#6BAED6" "#4292C6" "#2171B5" "#08519C" "#08306B" "#02214f"
# reds: "#e6e6e6" "#FFF5F0" "#FEE0D2" "#FCBBA1" "#FC9272" "#FB6A4A" "#EF3B2C" "#CB181D" "#A50F15" "#67000D" "#450009"
#oranges: "#FFF5EB" "#FEE6CE" "#FDD0A2" "#FDAE6B" "#FD8D3C" "#F16913" "#D94801" "#A63603" "#7F2704" "#682003"

#project the sampling points
prove<- vect(samps,geom=c("decimallongitude", "decimallatitude"),crs="+proj=longlat")
proven <- project(prove, crs_laea)

###save plot
## 1. Open a pdf file
#pdf("~/Desktop/grosbeak.breeding.distribution.map.sampling.test.pdf", width=15, height=20) 
## 2. Create a plot
#color-code separately
plot(region_boundary_laea, border = "white")
plot(countries, col = "grey90", border = NA, add = TRUE)
# add relative abundance
plot(abd_breeding_laea, breaks = breaks, col = c("#e6e6e6",brewer.pal(9, 'Reds'),"#450009"),
     maxcell = ncell(abd_breeding_laea), legend = FALSE, add = TRUE)
plot(y, breaks = breaks, col = c("#e6e6e6",brewer.pal(9, 'Oranges'),"#682003"),
     maxcell = ncell(abd_breeding_laea),legend = FALSE, add=TRUE)
plot(states, col = NA, border = "grey50", add = TRUE) #add state outline
plot(countries[[4]], border = "grey50", add=TRUE) #add canada outline
plot(countries[[28]], border = "grey50", add=TRUE) #add mexico outline
#add sampling points
plot(proven, add=TRUE, col="black", bg="white", pch=21, cex=1)

Code

## 3. Close the pdf file
#dev.off()

#test with yellow instead of orange
###save plot
## 1. Open a pdf file
#pdf("~/Desktop/grosbeak.breeding.distribution.map.sampling.test2.pdf", width=15, height=20) 
## 2. Create a plot
#color-code separately
plot(region_boundary_laea, border = "white")
plot(countries, col = "grey90", border = NA, add = TRUE)
# add relative abundance
plot(abd_breeding_laea, breaks = breaks, col = c("#e6e6e6",brewer.pal(9, 'Reds'),"#450009"),
     maxcell = ncell(abd_breeding_laea), legend = FALSE, add = TRUE)
plot(y, breaks = breaks,
     col = c("#e6e6e6","#fffdd0","#fffbae","#fff86b","#fff64e",
             "#fff319","#f5e900","#e5da00","#d5cb00","#b4ab00","#878000"),
     maxcell = ncell(abd_breeding_laea),legend = FALSE, add=TRUE)
plot(states, col = NA, border = "grey50", add = TRUE) #add state outline
plot(countries[[4]], border = "grey50", add=TRUE) #add canada outline
plot(countries[[28]], border = "grey50", add=TRUE) #add mexico outline
#add sampling points
plot(proven, add=TRUE, col="black", bg="white", pch=21, cex=1)

Code

## 3. Close the pdf file
#dev.off()

1.1 repeat for only South Dakota

Code

# set boundary
region_boundary <- ne_states(iso_a2 = "US") |> 
  filter(name == "South Dakota")

#transform
region_boundary_proj <- st_transform(region_boundary, st_crs(abd_breeding))
# crop and mask to boundary
abd_breeding_mask <- crop(abd_breeding, region_boundary_proj) |> mask(region_boundary_proj)
# map the cropped data
plot(abd_breeding_mask, axes = FALSE)

Code

# 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) |> trim()

# map the cropped and projected data for Rose-breasted Grosbeak
plot(abd_breeding_laea, axes = FALSE, breakby = "cases")

Code

#repeat for black-headed grosbeak
BHregion_boundary_proj <- st_transform(region_boundary, st_crs(BHabd_breeding))
# crop and mask to boundary
BHabd_breeding_mask <- crop(BHabd_breeding, BHregion_boundary_proj) |> mask(BHregion_boundary_proj)
# map the cropped data
plot(BHabd_breeding_mask, axes = FALSE)

Code

# transform to the custom projection using nearest neighbor resampling
BHabd_breeding_laea <- project(BHabd_breeding_mask, crs_laea) |> trim()

# map the cropped and projected data
plot(BHabd_breeding_laea, axes = FALSE, breakby = "cases")

Code

### plot them together

#remove 0 cells from BH raster which would overwrite the presence data from the rose-breasted
y <- classify(BHabd_breeding_laea, cbind(-Inf, 0, NA))

#remove 0 cells from RB raster 
abd_breeding_laea <- classify(abd_breeding_laea, cbind(-Inf, 0, NA))

# define the map plotting extent with the region boundary polygon
region_boundary_laea <- region_boundary |> 
  st_geometry() |> 
  st_transform(crs_laea)
# 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)
# 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)

#color-code separately
plot(region_boundary_laea, border = "black", lwd=5)
#plot(countries, col = "grey90", border = "grey90", add = TRUE)
# add relative abundance
plot(abd_breeding_laea,
     breaks = breaks, col = c("#e6e6e6",brewer.pal(9, 'Reds'),"#450009"), 
     maxcell = ncell(abd_breeding_laea),
     legend = FALSE,add=TRUE)
plot(y, breaks = breaks, col = c("#e6e6e6","#fffdd0","#fffbae","#fff86b","#fff64e",
             "#fff319","#f5e900","#e5da00","#d5cb00","#b4ab00","#878000"),
     maxcell = ncell(abd_breeding_laea),legend = FALSE, add=TRUE)
plot(region_boundary_laea, border = "black", lwd=5,add=TRUE)

Code

#reproject the sampling points and add them
prove<- vect(samps,geom=c("decimallongitude", "decimallatitude"),crs="+proj=longlat")
proven <- project(prove, crs_laea)

## 1. Open a pdf file
#pdf("~/Desktop/grosbeak.rad/SD.grosbeak.breeding.distribution.map.sampling.pdf", width=6, height=4) 
## 2. Create a plot
plot(region_boundary_laea, border = "black", lwd=5)
plot(abd_breeding_laea,
     breaks = breaks, col = c("#e6e6e6",brewer.pal(9, 'Reds'),"#450009"), 
     maxcell = ncell(abd_breeding_laea),
     legend = FALSE,add=TRUE)
plot(y, breaks = breaks, col = c("#e6e6e6","#fffdd0","#fffbae","#fff86b","#fff64e",
             "#fff319","#f5e900","#e5da00","#d5cb00","#b4ab00","#878000"),
     maxcell = ncell(abd_breeding_laea),legend = FALSE, add=TRUE)
plot(region_boundary_laea, border = "black", lwd=5,add=TRUE)
#add sampling points
plot(proven, add=TRUE, col="black", bg="white", pch=21, cex=1)

Code

## 3. close plot
#dev.off()

## 1. Open a pdf file
#pdf("~/Desktop/SD.grosbeak.map.sampling.pdf", width=6, height=4) 
## 2. Create a plot
plot(region_boundary_laea, border = "black", lwd=5)
#plot(abd_breeding_laea,breaks = breaks, col = c("#e6e6e6",brewer.pal(9, 'Reds'),"#450009"), maxcell = ncell(abd_breeding_laea),legend = FALSE,add=TRUE)
#plot(y, breaks = breaks, col = c("#e6e6e6",brewer.pal(9, 'Blues'),"#02214f"), maxcell = ncell(abd_breeding_laea),legend = FALSE, add=TRUE)
#add sampling points
plot(proven, add=TRUE, col="black", bg="white", pch=21, cex=1)

Code

## 3. close plot
#dev.off()

1.2 repeat in reverse order, to assess whether one raster is blocking the other

Code

### plot them together

#remove 0 cells from BH raster which would overwrite the presence data from the rose-breasted
y <- classify(BHabd_breeding_laea, cbind(-Inf, 0, NA))

#remove 0 cells from RB raster 
abd_breeding_laea <- classify(abd_breeding_laea, cbind(-Inf, 0, NA))

# define the map plotting extent with the region boundary polygon
region_boundary_laea <- region_boundary |> 
  st_geometry() |> 
  st_transform(crs_laea)
# 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)
# 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)

#color-code separately
plot(region_boundary_laea, border = "black", lwd=5)
#plot(countries, col = "grey90", border = "grey90", add = TRUE)
# add relative abundance
plot(y, breaks = breaks, col = c("#e6e6e6","#fffdd0","#fffbae","#fff86b","#fff64e",
             "#fff319","#f5e900","#e5da00","#d5cb00","#b4ab00","#878000"),
     maxcell = ncell(abd_breeding_laea),legend = FALSE, add=TRUE)
plot(abd_breeding_laea,
     breaks = breaks, col = c("#e6e6e6",brewer.pal(9, 'Reds'),"#450009"), 
     maxcell = ncell(abd_breeding_laea),
     legend = FALSE,add=TRUE)
plot(region_boundary_laea, border = "black", lwd=5,add=TRUE)

Code

#reproject the sampling points and add them
prove<- vect(samps,geom=c("decimallongitude", "decimallatitude"),crs="+proj=longlat")
proven <- project(prove, crs_laea)

## 1. Open a pdf file
#pdf("~/Desktop/grosbeak.rad/SD.grosbeak.breeding.distribution.map.sampling.pdf", width=6, height=4) 
## 2. Create a plot
plot(region_boundary_laea, border = "black", lwd=5)
plot(y, breaks = breaks, col = c("#e6e6e6","#fffdd0","#fffbae","#fff86b","#fff64e",
             "#fff319","#f5e900","#e5da00","#d5cb00","#b4ab00","#878000"),
     maxcell = ncell(abd_breeding_laea),legend = FALSE, add=TRUE)
plot(abd_breeding_laea,
     breaks = breaks, col = c("#e6e6e6",brewer.pal(9, 'Reds'),"#450009"), 
     maxcell = ncell(abd_breeding_laea),
     legend = FALSE,add=TRUE)
plot(region_boundary_laea, border = "black", lwd=5,add=TRUE)
#add sampling points
plot(proven, add=TRUE, col="black", bg="white", pch=21, cex=1)

Code

## 3. close plot
#dev.off()

This second approach actually made the hybrid zone appear even further west, which is not accurate according to our phenotypic and genomic data. So in the paper, we will use the first approach, which is more accurate according to our transect.