In jacob-ogre/ch.first.pass: Analysis of critical habitat changes 2001-2011

library(ch.first.pass)
library(dplyr)
library(plotly)
library(stringr)
library(viridis)
data("all_sp")
all_sp <- dplyr::select(all_sp, -from_to, -Count)
# vignette: >
#   %\VignetteIndexEntry{Vignette Title}
#   %\VignetteEngine{knitr::rmarkdown}
#   %\VignetteEncoding{UTF-8}

Introduction

Critical habitat occupies a contentious position in endangered species policy (see, e.g., James and Ward 2016). Although much has been written about the intent, scope, and effectiveness of critical habitat designations (Gibbs and Currie 2012; Mullen, Peterson and Todd 2013; Nelson et al. 2015; Taylor, Suckling and Rachlinski 2005), very little is known about its current ecological condition. This knowledge gap makes it difficult, if not impossible, to understand whether the U.S. Fish and Wildlife Service (FWS) is adequately protecting critical habitat or the extent to which critical habitat is furthering recovery. In this working document, we test the hypothesis that remotely-sensed data--particularly the National Land Cover Database (NLCD)--can help close this pervasive knowledge gap and provide a glimpse into the condition of critical habitat. Using a set of 42 ESA-listed species with critical habitat designated before 2000, we calculate the change in acreage of potentially suitable land cover types during the period 2001-2011. We expect that high rates of land cover change and declining species status signal inadequate protection of critical habitat, and may reveal instances of "destruction or adverse modification." Conversely, we expect that low rates of land cover conversion and high species recovery indicate appropriate levels of protection. Our preliminary results indicate that our 42 listed species have not witnessed significant losses of critical habitat during the study period.

Methods

To monitor land cover change during the period 2001-2011, we eliminated all listed species with critical habitat designated after 2000. We also eliminated fishes and species where habitat disturbance was not the primary threat to their recovery. These included wolves--threatened by hunting--and condors--threatened by lead--among several others. Ultimately, we were left with 42 species with critical habitat designated before 2000. We used only publically available data: critical habitat (provided by the Fish and Wildlife Service through its Environmental Conservation Online System [ECOS]) and the National Land Cover Database (provided by the Multi-Resolution Land Characteristics Consortium [MRLC]). Within the NLCD, we used several products, including: 2001 Land Cover (2011 edition), 2001 Percent Developed Imperviousness (2011 editon), 2011 Land Cover, 2011 Percent Developed Imperviousness, and 2001-2011 Land Cover From To Change Index. We made a geodatabase of critical habitat with a seperate feature class for each of our 42 case study species. Using ArcGIS Model Builder, we clipped each of the NLCD products to each of the feature classes (i.e., 5 rasters for each of our case study species). For the land cover products, we used arcpy scripts to calculate the percent cover (2001 and 2011), the total acreage change, and the percent acreage change of each land cover type for each feature class. For the percent developed imperviousness, we used arcpy scripts to to calculate the total number of acres with > 50% developed imperviousness for 2001 and 2011 for each of our 42 species, using those numbers to then calculate total and percent acreage change in imperviousness.

Mount Graham red squirrel habitats in 2001 (top) and in 2011 (bottom) after the large wildfires of 2004.

Green is evergreen forest; tan and whitish are shrub/scrub and grassland.

---

Least Bell's Vireo habitat in CH in 2001 (top), 2011 (middle), and the

change in habitats from 2001-2011.

Reds are classified as development of varying intensity; shades of brown and gray are 'natural' habitat.

Species selection

TBD

NLCD

TBD

Results

TBD

Species

no_chg <- filter(all_sp, from == to)
change <- filter(all_sp, from != to)

chg_by_sp <- tapply(change$Acres,
                    INDEX = change$Species,
                    FUN = sum, na.rm = TRUE)
sme_by_sp <- tapply(no_chg$Acres,
                    INDEX = no_chg$Species,
                    FUN = sum, na.rm = TRUE)

# And set up the data frame...
by_spp <- data.frame(Species = names(chg_by_sp),
                     chg_by_sp = as.vector(chg_by_sp))
tmp_sp <- data.frame(Species = names(sme_by_sp),
                     sme_by_sp = as.vector(sme_by_sp))
by_spp <- full_join(by_spp, tmp_sp, by = "Species")
by_spp$sp_total <- by_spp$chg_by_sp + by_spp$sme_by_sp
by_spp$chg_pct <- (by_spp$chg_by_sp / by_spp$sp_total) * 100
by_spp <- arrange(by_spp, desc(chg_pct))

The percentage of critical habitat that changed from one NLCD class to another ranged from r round(min(by_spp$chg_pct, na.rm = TRUE), digits = 1) to r round(max(by_spp$chg_pct, na.rm = TRUE), digits = 1) per species. The average amount of change was r round(mean(by_spp$chg_pct, na.rm = TRUE), digits = 1)% (median = r round(median(by_spp$chg_pct, na.rm = TRUE), digits = 1)%).

plot_a <- plot_ly(data = by_spp,
                  x = by_spp$Species,
                  y = by_spp$chg_pct, 
                  type = "bar",
                  marker = list(color = substr(viridis(1), 0, 7))) %>%
          layout(xaxis = list(title = "", tickangle = 60), 
                 yaxis = list(title = "% CH changed"), 
                 margin = list(b = 160))

plot_b <- plot_ly(by_spp,
                  x = by_spp$Species,
                  y = by_spp$chg_by_sp, 
                  type = "bar",
                  marker = list(color = substr(viridis(1), 0, 7))) %>%
          layout(xaxis = list(title = "", tickangle = 60), 
                 yaxis = list(title = "CH changed (acres)"), 
                 margin = list(b = 160))

# subplot(plot_a, plot_b, margin = 0.05, nrows = 2) %>% 
#     layout(title = "CH changes by percent (top) and area (bottom)", 
#            showlegend = FALSE)
htmltools::tagList(list(plot_a, plot_b))

The proportions of habitat transitions varies by species. For example, the Mt. Graham red squirrel had the largest proportion of its habitat change, primarily because of a [[year]] fire

make_heatmap(all_sp, "MountGrahamRedSquirrel", "MountGrahamRedSquirrel")

Whooping crane, which saw 3%, or >11,000 acres, of habitat change within critical habitat:

make_heatmap(all_sp, "WhoopingCrane", "WhoopingCrane")

Habitats

We next asked if there were systematic habitat transitions across species. MORE

all_hab <- tapply(all_sp$Acres,
                  INDEX = all_sp$from,
                  FUN = sum, na.rm = TRUE)
chg_fr_hb <- tapply(change$Acres,
                    INDEX = change$from,
                    FUN = function(x) -sum(x, na.rm = TRUE))
chg_to_hb <- tapply(change$Acres,
                    INDEX = change$to,
                    FUN = sum, na.rm = TRUE)

# And set up the data frame...
tot_hab <- data.frame(habitat = names(all_hab),
                      tot_hab = as.vector(all_hab))
by_hab <- data.frame(habitat = names(chg_to_hb),
                     chg_to_hb = as.vector(chg_to_hb))
tmp_hb <- data.frame(habitat = names(chg_fr_hb),
                     chg_fr_hb = as.vector(chg_fr_hb))
by_hab <- full_join(by_hab, tmp_hb, by = "habitat")
by_hab <- full_join(by_hab, tot_hab, by = "habitat")
by_hab$net_chg <- by_hab$chg_fr_hb + by_hab$chg_to_hb
by_hab$chg_pct <- (by_hab$net_chg / by_hab$tot_hab) * 100
by_hab <- arrange(by_hab, desc(chg_pct))

plot_a <- plot_ly(by_hab,
                  x = by_hab$habitat,
                  y = by_hab$chg_pct, 
                  type = "bar",
                  marker = list(color = substr(viridis(1), 0, 7))) %>%
          layout(xaxis = list(title = "", tickangle = 60), 
                 yaxis = list(title = "% CH changed"), 
                 margin = list(b = 160))

plot_b <- plot_ly(by_hab,
                  x = by_hab$habitat,
                  y = by_hab$net_chg, 
                  type = "bar",
                  marker = list(color = substr(viridis(1), 0, 7))) %>%
          layout(xaxis = list(title = "", tickangle = 60), 
                 yaxis = list(title = "CH changed (acres)"), 
                 margin = list(b = 160))

# subplot(plot_a, plot_b, margin = 0.05, nrows = 2) %>% 
#     layout(title = "CH changes by percent (top) and area (bottom)", 
#            showlegend = FALSE)
htmltools::tagList(list(plot_a, plot_b))

To see which habitat was converted to which habitat, we can create a heatmap. Note that the acres are log10-transformed because the area of shrub/scrub is so large relative to other areas:

make_heatmap(all_sp, "", "All species")

Discussion

TBD

Appendix

To facilitate species-by-species evaluation of habitat changes, we provide 'from:to' heatmaps for each of the 42 species in the dataset:

plots <- list()
un_sp <- unique(all_sp$Species)
for(i in 1:length(un_sp)) {
  plots[[i]] <- make_heatmap(all_sp,
                             species = un_sp[i],
                             title = un_sp[i],
                             height = 600)
}
htmltools::tagList(plots)

jacob-ogre/ch.first.pass documentation built on May 18, 2019, 8 a.m.