In TrashBirdEcology/bbsRDM: Calculate regime detection metrics using bird communities

Setup and Overview

This document briefly reviews the functionality of bbsRDM. Although this package can be used to calculate and visualize BBS data using time series, the example at hand presents an application to large spatial transects across the continental United States.

knitr::opts_chunk$set(
  collapse = TRUE,
  # comment = "#>", 
  eval=TRUE, 
  echo=FALSE, 
  # cache.path = "cache/",
  cache=TRUE,
  warning=FALSE, 
  message=FALSE,
  out.width="65%", 
  fig.show=""
)

options(knitr.graphics.auto_pdf = TRUE)

# load("runthrough.RData")

First, install the bbsRDM and regimeDetectionMeasures packages from GitHub.

devtools::install_github("trashbirdecology/bbsRDM",
                         force = FALSE, dependencies = FALSE)
devtools::install_github("trashbirdecology/regimedetectionmeasures",
                         force = FALSE, dependencies = FALSE)

library(bbsRDM)
library(regimeDetectionMeasures)
library(sp)
library(raster)
library(feather)
library(here)
library(ggplot2)

ggplot2::theme_set(theme_bw())

Next, create diretories for storing the BBS and results locally (if directories exist, this BBS data will not be downloaded (unless you specify 'downloadBBSData=TRUE'), however, the results will be overwritten!).

# a. Create a directory to store and/or load the BBS data as feathers
bbsDir <- here::here("bbs_raw_data")
# If bbsDir is NOT empty then we will NOT download any data.
if (length(list.files(bbsDir, pattern = "*.feather")) > 0) {
downloadBBSData = FALSE
} else
(
{
dir.create(bbsDir)
downloadBBSData = TRUE
}
) # If this returns a warning, proceed with caution as directory already exists, 
     # results WILL be written over
# b. Create a directory to store and/or load the BBS data as feathers
resultsDir <- here::here("myResults")
dir.create(resultsDir)
# c. Create directory for storing early warning signal results
dir.create(here::here("myResults/ews"))
# d. Create directory for storing distance travelled results
dir.create( here::here("myResults/distances"))

Retrieve the BBS Data and Save Locally

If necessary, download all the state data. Downloading the entire BBS data takes a while, so only run if you have not recently downloaded the BBS data. If you only need a subset of the data, please see the documentation for bbsAssistant::get_regions and bbsRDM::getDataBBS.

# a. Load the regional .txt file from Patuxent's FTP server 
# (you must be connected to the internet to perform this step)
regions <- bbsAssistant::get_regions()

# b. Create a series or one filenames for states, regions
regionFileName <- regions$zipFileName %>% na.omit()

# c.  Download and unzip the BBS data.
if(downloadBBSData==TRUE){
for(i in 1:length(regionFileName)){
        bbsData <-  importDataBBS(
            # arguments for getDataBBS()
            file = regionFileName[i],
            dir =  "ftp://ftpext.usgs.gov/pub/er/md/laurel/BBS/DataFiles/States/",
            year = NULL,
            aou = NULL,
            countrynum = NULL,
            states = NULL,
            #  arguments for getRouteInfo():
            routesFile = "routes.zip",
            routesDir =  "ftp://ftpext.usgs.gov/pub/er/md/laurel/BBS/DataFiles/",
            RouteTypeID = 1,
            # one or more of c(1,2,3)
            Stratum = NULL,
            BCR = NULL
        )


# d. Save the unzipped files to disk.
birdsToFeathers(dataIn  = bbsData,
                newDir  = bbsDir,
                filename = regionFileName[i])
# e. Clear large object from memory
rm(bbsData)}}

Create a Spatial Sampling Grid Across North America

Next we build a spatial sampling grid for aligning BBS data to regularly spaced cells. This is important for spatial interpretation of the regime detection metric results.

# Define the grid's cell size (lat, long; unit:degrees)
        ## 1 deg latitude ~= 69 miles
        ## 1 deg longitude ~= 55 miles
cs <-
    c(0.5, 0.5)  # default is cell size 0.5 deg lat x 0.5 deg long

# Create the grid
routes_gridList <- createSamplingGrid(cs = cs)

Now we load in the BBS data from the feathers we created and align with the sampling grid. This requires a bit of memory, proceed with caution.

feathers <- NULL
featherNames <- list.files(bbsDir, pattern = ".feather")
featherNames <- str_c("/", featherNames) #add separator
for (i in 1:length(featherNames)) {
  feather <- NULL
  feather <- loadBirdFeathers(newDir  = bbsDir,
                              filename = featherNames[i]) 

  feather <- feather %>%
    dplyr::rename(lat = latitude,
                  long = longitude) %>%
    left_join(routes_gridList$routes_grid, 
              by = c("countrynum", "statenum", "route", "lat", "long"))

  feathers <- rbind(feathers, feather)
  rm(feather)
 }

Subset the BBS data (optional)

By species (using AOU numbers)

Although subsetting the species is optional, this package contains features for subsetting by AOU code, functional groups, or by spatial location (e.g. remove all Montana observaitons).

Subset species according to AOU species codes (i.e. by family, genera, etc..)
For this example we will remove shorebirds, wading birds, and waterfowl (i.e., AOU species' codes 0000:2880). *See R/subsetByAOU.R source code or documentation for options (see: subset.by)

# Subset the species
feathers <- subsetByAOU(myData = feathers, subset.by= 'remove.shoreWaderFowl')

Subset species according to functional traits (or body mass).

Note: eval = FALSE...change to true if you wish to evaluate in rmd knitting.

# Create a single list of mass and functional traits. 
funMass <-
    funcMass(dataWD = here::here("data/"),
             fxn = TRUE, # get functional trait data?
             mass = TRUE) # get body mass data?
# Combine the functional trqits and/or body mass
bbsData <-
    mergeFunMassBBS(bbsData = feathers[1:1e3,], funMass = funMass)

Calculate regime detection metrics across space or time

First, define the parameters required to calcualte the metrics.

# Which metrics do you want to calculate?
metrics.to.calc <- c("distances", "ews")

# If calculating "EWSs, you can calculate select metrics. 
## Default = all early-warning signals, FI, and VI
to.calc = c("EWS", "FI", "VI")

# Choose spatial or temporal analysis
direction <-
    "South-North" # choose one of : 'South-North', 'East-West', or 'temporal'

# Choose the fill value as the unobserved value. 
fill = 0

# Minimum number of sites within each transect
min.samp.sites = 8

# Minimum number of BBS routes within each moving window 
min.window.dat = 3

# Which Equation of Fisher Information to use (default = 7.12)
fi.equation = "7.12"

# By what % of the entire data should the window move? 
winMove = 0.25

# Define some filtering and labeling parameters based on 
# direction of spatial analysis  
if (direction == "South-North") {
    dir.use =  unique(feathers$colID) %>% na.omit(colID) %>% sort()}
if (direction == "East-West") {
    dir.use = unique(feathers$rowID) %>% na.omit(rowID) %>% sort()}

Define the years we want to analyze

For this (spatial) example, we will analyze only every tenth year

# Get all possible years
years.use = unique(feathers$year)
# Keep only the years which are divisible by T
T = 10
years.use  <- years.use[which(years.use %% T == 0 & years.use > 1975)] %>% sort()

Calculate Regime Detection Measures Across Spatial Transects

This section will loop through years.use and dir.use, running each BBS route (temporal analysis) or spatial transect by year (spatial analysis) at a time. Results are saved in directories created in section below. Please note: depending on the # of years and # spatial transects, this could take a while.

######################################
### NOTE: change eval=TRUE if the ####
### results are not yet stored in file.
######################################
for (j in 1:(length(dir.use))) {
    # For east-west analysis
        if (direction == "East-West"){
            birdsData <- feathers %>%
                filter(rowID == dir.use[j]) %>%
                mutate(direction = direction,
                       dirID = dir.use[j])
    }
    # For south-north analysis
        if (direction == "South-North"){
            birdsData <- feathers %>%
            filter(colID == dir.use[j]) %>%
            mutate(direction = direction,
                   dirID = dir.use[j])
    }

      if (nrow(birdsData) < min.samp.sites) {
            next(print(paste0("Not enough data to analyze. 
                              Skipping j-loop ", dir.use[j])))
    }


    # Analyze the subset of data
    for (i in 1:length(years.use)){
        # a. Subset the data 
        birdData <- birdsData %>%
            filter(year == years.use[i]) %>%
            dplyr::rename(variable = aou,
                          value = stoptotal)

        if (nrow(birdData) == 0){
            next
        }

        # b. Munge the data further
        birdData <- mungeSubsetData(birdData)

        # c. Calculate metrics
        calculateMetrics(dataIn = birdData, metrics.to.calc, 
                         direction = direction,  yearInd = years.use[i])
        message(paste0("End i-loop (years) ", i, " of ",  
                       length(years.use)))
    }  
    message(paste0("End j-loop (transects) ", j, " of ",  
                   length(dir.use)))
}

Import and munge the results to prepare for visualization

First, use the function 'importResults' to import and combine the results as created in the previous code chunk. The following chunk chunk will import the EWS results and the distance results separately, combining each into their own data frames.

#  Import EWS results
results_ews <-
    importResults(resultsDir = resultsDir, myPattern = 'ews',
                  subset.by = direction) %>%
# assign the end of the window as the cellID
    mutate(cellID = cellID_max)     
# glimpse(results_ews)

#  Import distance results
results_dist <-
    importResults(resultsDir = resultsDir, myPattern = 'distances', subset.by = direction)
# glimpse(results_dist)

Next, get the results to align with our sampling grid for visualizing results across space.

#  Get the spatial sampling grid coordinates
coords_grd <-
    cbind(routes_gridList$sp_grd@data,
          coordinates(routes_gridList$sp_grd)) %>%
    rename(lat = s2,
           long = s1,
           cellID  = id)

# Join coords_grd with results
# note: a full join will likely produce many cells with NO results data..
# but NO lat or long should == NA!
results_dist <-
    full_join(coords_grd,
              results_dist) %>%
    na.omit(metricType)

results_ews <-
    full_join(coords_grd,
              results_ews) %>%
    na.omit(metricType) %>%
    dplyr::select(-cellID_min,-cellID_max, -winStart  , -winStop)

# Set coordinate system and projection of both results
coordinates(results_dist) <-
    coordinates(results_ews) <- c("long", "lat")
sp::proj4string(results_dist) <-
    sp::proj4string(results_ews) <-
    sp::CRS("+proj=longlat +datum=WGS84")

Visualize Results: Temporally

First, specify plotting parameters (below). We can visualize either the distance results (results_dist) or the early-warning signal results (results_ews). For this example we will visualize the distance results.

# Define the results we want to visualize
plotResults <- results_dist
# Sort the years
year.ind <- unique(plotResults@data$year) %>% sort()
# Create a label for plotting, depending on direction
sortVar.lab <-
    ifelse(unique(plotResults@data$direction)=="South-North",
           "latitude",
           "longitude")

# Specify the transect # we want to see
dirID.ind <- 13
metric.ind <- "s"
# build plot
(p.s <- sort.year.line(plotResults, metric.ind, year.ind, dirID.ind, dirInd,  scale = T, center = T))

Plot the distance travelled, s for transect # 13 (Figure \@ref(fig:plots)). This measure indicates how quickly the multivariable system is moving through phase space. Notice that the absolute value of s increases over time. This is unsurprising given that the number of species in BBS censuses has increased over time. The importance of this phenomenon should be evaluated at the local scale.

# Specify the transect # we want to see
dirID.ind <- 13
metric.ind <- "dsdt"
(p.dsdt <- sort.year.line(plotResults, metric.ind, year.ind, dirID.ind, dirInd,  scale = T, center = T))

Plot the velocity of the distance travelled, dsdt for transect # 13 (Figure \@ref(fig:plotdsdt)). Rapid changes in the velocity of the distance travelled indicates rapid species turnover, or large changes in the abundances of a few species. The ecological significance of rapid change in community abundances, or among the state variables used to calculate distance travelled and velocity, should be evaluated at local and sub-regional scales.

Visualize Results: Spatially

library(viridis)
temp <- as.data.frame(plotResults) %>%
    filter(
        year %in% c(1980, 1990, 2000, 2010),
           metricType == "s") %>%
    group_by(dirID, year, analysis) %>%
    mutate(scaledMetricValue.bycol = base::scale(metricValue, center =T, scale = T)) %>%
    ungroup() %>%
    na.omit(metricValue)

ggplot(temp) + geom_raster(aes(x = long, y = lat, fill = scaledMetricValue.bycol)) +
    coord_fixed(ratio = 1) +
    scale_fill_viridis(direction = -1) +
    facet_wrap(~year, strip.position = "top", nrow=2) +
    theme(legend.position = "none")

Plot the distance travelled across the entire range (Figure \@ref(fig:plotsus)).

library(viridis)
temp <- as_tibble(plotResults) %>%
    filter(year %in% c(1980, 1990, 2000, 2010),
           metricType == "dsdt") %>%
    group_by(dirID, year) %>%
    mutate(logmetric = log(metricValue)) %>% 
    ungroup()


ggplot(temp) + geom_raster(aes(x = long, y = lat, fill = logmetric)) +
    coord_fixed(ratio = 1) +
    scale_fill_viridis(direction = -1) +
    facet_wrap(~year, strip.position = "top", nrow=2) +
    theme(legend.position = "none")

Plot the velocity of the distance travelled across the entire range (Figure \@ref(fig:plotdsdtus)).

t <- as_tibble(plotResults@data) %>%
    filter(year %in% c(1980, 1990, 2000, 2010),
           metricType == "s") %>%
    group_by(dirID, year, analysis) %>%
    mutate(scaledMetricValue.bycol = base::scale(metricValue, center =T, scale = T)) %>%
    ungroup() %>%
    group_by(metricType, cellID) %>%
    arrange(year) %>%
    mutate(dScaledMetricValue.bycol = scaledMetricValue.bycol - lag(scaledMetricValue.bycol))

ggplot() + geom_raster(data = t, aes(x = long, y = lat, fill = dScaledMetricValue.bycol)) +
    coord_fixed(ratio = 1) +
    scale_fill_viridis(direction = -1) +
    facet_wrap(~year,  strip.position = "top") +
    theme(legend.position = "none")+
    geom_point(data = milBases)

# save.image("runthrough.RData")

Military Bases

Use the function getMilBases() to retrieve MIRTA military base point shapefiles:

milBases <- getMilBases() # default arguments will download the MIRTA FY 2018 points

sp::proj4string(milBases) <- 
    sp::CRS("+proj=longlat +datum=WGS84")# Set projection to WGS84 lat long friendly

Next, we need to assign the row and column ID numbers from our sampling grid to each base location, such that we can visualize the results of the regime detection measures within our spatial sampling grid.

# Extract the CELL ID in which each military base falls.
milBases.df <- milBases %>% as.data.frame() %>% 
    rename(long = coords.x1, lat = coords.x2) 
milBases.df$cellID <-  milBases %over%  routes_gridList$sp_grd 

# remove bases outside our sampling grid
milBases.df <- milBases.df %>% filter(!is.na(cellID))

ggplot(milBases.df, aes(long, lat, color = COMPONENT))+ geom_point()+
    ggtitle("Military Bases in the Continental U.S.")

Plot the location of military bases (Figure \@ref(fig:plotmilbases)).

TrashBirdEcology/bbsRDM documentation built on July 21, 2019, 2:18 a.m.