Mapping underserved communities

Below, we have provided instructions for replicating our methodology for mapping underserved communities in Tampa Bay. To view instructions for downloading the necessary source data from EJScreen, see Getting demographic data.

Load the required R packages (install first as needed).

library(sf)
library(mapview)
library(dplyr)
library(RColorBrewer)

Load the census tract data and map only the spatial data.

load(file = 'data/tb_tract.RData')

tb_tract %>%
  select(-everything()) %>% 
  mapview(layer.name = "Census tracts")

You can see that some tracts are only representative of the bay. We can clean up this data by retaining only the tracts in which the total population recorded in the latest American Community Survey (“ACSTOTPOP” attribute field) was above zero. This will remove tracts in which no people reside (e.g., large waterbodies, parks, or other natural areas).

tb_tract <- tb_tract %>%
  filter(ACSTOTPOP > 0)

Projects funded by BIL

The 2021 Bipartisan Infrastructure Law (BIL), implemented through the EPA and NEPs, provides approximately $909,800 in annual funding to each NEP over the 2022-2026 period. In this section, we outline our approach for ensuring that projects using BIL funding will sustain and increase investments in underserved communities, and the benefits that flow to them.

The EPA uses the EJScreen Supplemental Demographic Index to identify disadvantaged communities at the block group level. The Index is based on an average of the following five demographic variables:

  • Percent classified as low income
  • Percent unemployed
  • Percent linguistically isolated (i.e., limited English speaking)
  • Percent with less than a high school education
  • Percent with low life expectancy

Based on the average percent across these variables, the EPA recommends flagging communities that fall within the 80th percentile (or higher) nationally as potentially disadvantaged. However, this index represents a relatively exclusive indicator of disadvantage. Some communities may face significant challenges but only meet the national threshold in two out of these five variables, ultimately falling below the EPA’s Supplemental Demographic Index threshold.

TBEP considers a more inclusive definition of underserved communities, in which a census tract must meet or exceed the 80th percentile in at least two of the five demographic screening variables. This approach increases the number of communities that could be considered underserved while reducing the number of errors in identification from meeting just one threshold (e.g., a community of predominantly wealthy retirees/unemployed residents). Aggregation to the census tract level also permits linking this layer with other socioeconomic and environmental data provided by the federal government at this same community delineation.

Run the code below to (1) count the number of demographic thresholds met in each tract, and (2) identify which tracts will be classified as “underserved” communities more broadly.

underserved_tract <- tb_tract %>%
  mutate(threshold_income = ifelse(P_LWINCPCT >= 80, 1, 0),
         threshold_unempl = ifelse(P_UNEMPPCT >= 80, 1, 0),
         threshold_lingui = ifelse(P_LNGISPCT >= 80, 1, 0),
         threshold_educat = ifelse(P_LESHSPCT >= 80, 1, 0),
         threshold_lifexp = ifelse(P_LIFEEXPCT >= 80, 1, 0)) %>%
  rowwise() %>%
  select(matches('^threshold|^ID')) %>% 
  mutate(threshold_N = sum(threshold_income,threshold_unempl,threshold_lingui,threshold_educat,threshold_lifexp, na.rm = TRUE)) %>%
  mutate(underserved = ifelse(threshold_N > 1, "Yes", "No"))

View by Census Tract

View a map showing the number of thresholds met per census tract (you may adapt to a different color scale of your choice).

underserved_tract %>% 
  select(threshold_N) %>% 
  mapview(zcol = "threshold_N", col.regions = brewer.pal(6, "Reds"), layer.name = "No. of Thresholds Met")

View the tracts that meet our more inclusive definition of underserved communities. The areas in red are those that rank in the 80th percentile (or greater) nationally in 2 or more of the demographic screening variables. They will serve as priority areas for increasing the equitable distribution of benefits from TBEP’s non-BIL funded environmental programs.

underserved_tract %>% 
  select(underserved) %>% 
  mapview(zcol = "underserved", col.regions = list("gray","red"), layer.name = "Underserved Communities")

You can save this final data as an RData object for future use.

# save the layer as an RData object
save(underserved_tract, file = 'data/underserved_tract.RData')

View by Drainage Basin

While the census tract delineations provide a practical map for identifying target neighborhoods and stakeholders that may be working for or in these underserved communities, it is not well-aligned with geographic delineations that are most informative for planning conservation and and restoration projects across the watershed.

For environmental planning purposes, TBEP will thus characterize unique water body assessment units according to the proportion of each unit containing underserved census tracts. This approach serves as a bridge between social and ecological units relevant to different stakeholders across Tampa Bay.

The Florida Department of Environmental Protection provides the Waterbody IDs (WBIDs) dataset, available here, which includes polygons delineating the drainage basins surrounding water body assessment units.

Load the drainage basin data.

dbasins <- st_read('https://ca.dep.state.fl.us/arcgis/rest/services/OpenData/WBIDS/MapServer/0/query?outFields=*&where=1%3D1&f=geojson')
Reading layer `OGRGeoJSON' from data source 
  `https://ca.dep.state.fl.us/arcgis/rest/services/OpenData/WBIDS/MapServer/0/query?outFields=*&where=1%3D1&f=geojson' 
  using driver `GeoJSON'
Simple feature collection with 6788 features and 17 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -87.6349 ymin: 24.39631 xmax: -79.97431 ymax: 31.00065
Geodetic CRS:  WGS 84

Intersect the layer with the Tampa Bay watershed.

load(file = 'data/tbshed.RData')

sf_use_s2(FALSE)

tb_dbasins <- dbasins %>% 
  st_intersection(tbshed)

Since we will be doing area calculations, we need to set a projected coordinate system for the drainage basin and census tract layers. The most appropriate CRS for Tampa Bay is EPSG 6443.

tb_dbasins <- st_transform(tb_dbasins, crs = 6443)
underserved_tract <- st_transform(underserved_tract, crs = 6443)

Keep only the land area within each drainage basin and census tract. We can exclude water areas by intersecting the layers with Florida’s shoreline (available here from the Florida Fish and Wildlife Conservation Commission).

flshore <- st_read('https://atoll.floridamarine.org/arcgis/rest/services/FWC_GIS/OpenData_Shoreline/MapServer/0/query?outFields=*&where=1%3D1&f=geojson')
Reading layer `OGRGeoJSON' from data source 
  `https://atoll.floridamarine.org/arcgis/rest/services/FWC_GIS/OpenData_Shoreline/MapServer/0/query?outFields=*&where=1%3D1&f=geojson' 
  using driver `GeoJSON'
Simple feature collection with 15470 features and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -87.63477 ymin: 24.52114 xmax: -80.03118 ymax: 31.00091
Geodetic CRS:  WGS 84
flshore <- st_transform(flshore, crs = 6443)

tb_dbasins <- tb_dbasins %>% 
  st_intersection(flshore)

underserved_tract <- underserved_tract %>% 
  st_intersection(flshore) %>%
  filter(underserved == "Yes")

Calculate (1) the total land area of each drainage basin in the watershed and (2) the total area of underserved tracts within each drainage basin. Units are in feet for this projected coordinate system.

# total basin area
basinareas <- tb_dbasins %>%
  mutate(area_ft = st_area(tb_dbasins)) %>%
  group_by(WBID) %>%
  summarise(db_area_ft = sum(area_ft))

# underserved area within basins
basin_int <- st_intersection(basinareas, underserved_tract)

basin_underareas <- basin_int %>%
  mutate(area_int_ft = st_area(basin_int)) %>%
  group_by(WBID) %>%
  summarise(under_area_ft = sum(area_int_ft)) %>%
  as.data.frame()

# join area estimates
areasjoined <- left_join(basinareas, basin_underareas, by = 'WBID') %>%
  rowwise() %>%
  select(matches('^WBID|^db|^under')) %>%
  as.data.frame()

tb_dbasins_under <- left_join(tb_dbasins, areasjoined, by = 'WBID')

Calculate the proportion of each drainage basin containing underserved communities.

underserved_dbasins <- tb_dbasins_under %>%
  mutate(under_area_ft = ifelse(is.na(under_area_ft), 0, under_area_ft)) %>%
  mutate(pct_under = under_area_ft/db_area_ft * 100) %>%
  mutate(pct_under = as.numeric(pct_under)) %>% 
  select(HUC, WBID, WATERBODY_NAME, pct_under)

Create a map showing priority drainage basins based on the presence of underserved communities. Hover over the drainage basins to view the name of the water body assessment unit.

mapviewOptions("basemaps.color.shuffle" = FALSE)

underserved_dbasins %>% 
  select(WATERBODY_NAME, pct_under) %>% 
  mapview(zcol = "pct_under", col.regions = brewer.pal(6, "Reds"), layer.name = "Underserved Tracts (% of DB)", label = "WATERBODY_NAME")

Exclude drainage basins with no underserved communities present.

mapviewOptions("basemaps.color.shuffle" = FALSE)

underserved_dbasins %>% 
  select(WATERBODY_NAME, pct_under) %>% 
  filter(pct_under > 0) %>%
  mapview(zcol = "pct_under", col.regions = brewer.pal(6, "Reds"), layer.name = "Underserved Tracts (% of DB)", label = "WATERBODY_NAME")

Save this final data as an RData object for future use.

# save the layer as an RData object
save(underserved_dbasins, file = 'data/underserved_dbasins.RData')