example.qmd

---
title: "Example application: establishing a baseline"
format: html
editor: source

execute:
  echo: true
---

The purpose of mapping underserved communities is to monitor and report on TBEP's contribution to increasing equitable and fair access to environmental benefits in Tampa Bay. In line with the goals of the Justice40 Initiative, TBEP's Equity Strategy has a goal to ensure that at least 40% of the benefits from BIL-funded activities are distributed to underserved communities.

Recognizably, there are a myriad of benefits that flow from the diverse conservation, restoration, research, and outreach activities that TBEP produces, spanning the many provisioning, supporting, regulating, and cultural services provided by nature. While we support ongoing efforts to develop additional methodologies to quantify these unique services provided by different activities, TBEP considers the following factors when defining and measuring the flow of benefits from our BIL-funded activities for each fiscal year:

-   The number of BIL-funded projects that benefit underserved communities
-   The amount of BIL funding/investment dollars in projects that benefit underserved communities

There are two restoration themes that BIL funds will contribute to:

1.  Building Neighborhood Resilience through Tidal Tributary Restoration
2.  Nitrogen Management through Resilient Wastewater Infrastructure

Projects within these two themes provide both local and downstream benefits to communities reliant on the ecosystem services provided by Tampa Bay's extensive network of tributaries, canals, and other waterways. Therefore, for a BIL-funded project to benefit an underserved community (and thus contribute toward our equity target), it must meet the following criteria:

-   The project is located within one mile of an underserved community; *OR*
-   The project is located in a drainage basin that contains at least one underserved community, *AND*
-   The underserved community is located *downstream* of the project.

Consideration of these metrics relative to all BIL-funded projects will inform TBEP's progress toward the 40% target. Notably, these metrics are based on projects with physical locations and represent indicators of distributive equity only. Additional metrics will be developed in the future that can apply to other types of projects (including research, outreach, and education) and measure progress toward procedural and recognitional equity.

Below, we have provided an example of how we can utilize the maps of underserved communities to estimate progress toward achieving the 40% distributive equity target by generating a baseline of the equitable nature of TBEP habitat projects over the previous 5 years.

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

```{r}
#| eval: true
library(sf)
library(leaflet)
library(dplyr)
library(RColorBrewer)
library(ggplot2)
library(gridExtra)
library(reshape2)
library(leafsync)
```

Load the map of underserved communities and drainage basins that will guide BIL-funded projects.

```{r}
load(file = 'data/underserved_tract.RData')
load(file = 'data/underserved_dbasins.RData')

# project to EPSG 6443 (Florida West) and keep only those meeting underserved definition
underserved_tract <- underserved_tract %>%
  filter(underserved == "Yes") %>%
  st_transform(crs = 6443)
underserved_dbasins <- underserved_dbasins %>%
  filter(pct_under > 0) %>%
  st_transform(crs = 6443)
```

## Baseline equity stats for 2017-2021

The EPA's National Estuary Program Online Reporting Tool (NEPORT) is a database that NEP staff use for reporting habitat projects occurring within each NEP on an annual basis. TBEP's NEPORT database contains point locations of all habitat restoration projects in the Tampa Bay watershed since 2006. It also includes a number of project characteristics, such as the name and description of the project, the type of habitat and activities associated with it, lead and other partners, primary funding source, and total project costs. We will use this data to establish a baseline of benefits to underserved communities from past TBEP programs.

Load the NEPORT data.

```{r}
load(file = 'data/tbepNEPORT.RData')

tbepNEPORT <- st_transform(tbepNEPORT, crs = 6443)
```

Filter the dataset. For this baseline, we only want habitat projects from 2017-2021 (five years before BIL), and only those projects in which TBEP is the lead partner, primary funder, and/or grant administrator.

```{r}
baseline5yr <- tbepNEPORT %>%
  filter(Year > 2016 & Year < 2022,
         grepl('TBEP|Tampa Bay Estuary Program', Lead_Imple) |
         grepl('TBEP|Tampa Bay Estuary Program|TBERF|Tampa Bay Environmental Restoration Fund|Bay Mini-Grant|RESTORE', Main_Fundi)) %>%
  st_zm()
baseline5yr_wgs84 <- st_transform(baseline5yr, crs = 4326)

leaflet(baseline5yr_wgs84) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addCircleMarkers(
    fillColor = ~"blue",
    weight = 0.5,
    opacity = 1,
    radius = 4,
    color = "black",
    fillOpacity = 0.7,
    popup = ~paste("ID:", ID),
    group = "TBEP Habitat Projects (2017-2021)"
  ) %>%
  addLayersControl(
    overlayGroups = c("TBEP Habitat Projects (2017-2021)"),
    options = layersControlOptions(collapsed = FALSE)
  )
```

Create a 1 mile buffer around the habitat projects.

```{r}
baseline5yr_buff <- st_buffer(baseline5yr, dist = 5280)
```

Identify the projects that are located within 1 mile of underserved communities.

```{r}
# projects within 1 mile of underserved communities
baseline5yr_buff_tractint <- st_intersection(baseline5yr_buff, underserved_tract)

summaries <- as.data.frame(baseline5yr_buff_tractint) %>%
  group_by(ID) %>%
  summarise(communities_1mile = n())

# merge with NEPORT dataset
baseline5yr <- left_join(baseline5yr, summaries, by = 'ID') %>%
  mutate(communities_1mile = coalesce(communities_1mile, 0))
```

Of the projects over 1 mile from underserved communities, identify which ones are in a drainage basin containing underserved communities.

```{r}
# projects within 1 mile of underserved communities
baseline5yr_dbint <- st_intersection(baseline5yr, underserved_dbasins) %>%
  filter(communities_1mile == 0)

summariesdb <- as.data.frame(baseline5yr_dbint)
```

To identify which of these projects have downstream benefits to underserved communities, we need to visualize the various tributaries, canals, and other waterbodies flowing throughout the watershed. The Florida National Hydrography Dataset's Flowlines (24k) layer is available from the [Florida Department of Environmental Protection](https://geodata.dep.state.fl.us/datasets/FDEP::florida-national-hydrography-dataset-nhd-flowlines-24k/about){target="_blank"}. Because this is a large shapefile, we have already saved a clipped version of it. However, you can download and read the layer on your own using
<https://ca.dep.state.fl.us/arcgis/rest/services/OpenData/NHD/MapServer/4/query?outFields=*&where=1%3D1&f=geojson>.

```{r}
#| cache: true
load(file = 'data/tb_flowlines.RData')
tb_flowlines <- st_zm(tb_flowlines)
tb_flowlines_wgs84 <- st_transform(tb_flowlines, crs = 4326)

leaflet(tb_flowlines_wgs84) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolylines(
    color = ~"blue",
    weight = 1,
    opacity = 0.7
  )
```

View the projects in relation to the flowlines. Adding in the underserved communities and drainage basin boundaries will help you determine if there are downstream benefits to these communities for each project.

```{r}
#| cache: true

# Transform all data to WGS84 for leaflet
baseline5yr_dbint_wgs84 <- st_transform(baseline5yr_dbint, crs = 4326)
underserved_dbasins_wgs84 <- st_transform(underserved_dbasins, crs = 4326)
underserved_tract_wgs84 <- st_transform(underserved_tract, crs = 4326)
tb_flowlines_wgs84 <- st_transform(tb_flowlines, crs = 4326)

# Create combined leaflet map with multiple layers
leaflet() %>%
  addProviderTiles("CartoDB.Positron") %>%
  # Add baseline5yr_dbint layer
  addCircleMarkers(
    data = baseline5yr_dbint_wgs84,
    fillColor = "yellow",
    weight = 0.5,
    opacity = 1,
    radius = 4,
    color = "black",
    fillOpacity = 0.7,
    popup = ~paste("ID:", ID),
    group = "Projects in DBs with underserved communities"
  ) %>%
  # Add underserved_dbasins layer
  addPolygons(
    data = underserved_dbasins_wgs84,
    fillColor = "black",
    weight = 0.5,
    opacity = 1,
    color = "black",
    fillOpacity = 0,
    group = "DBs with underserved communities"
  ) %>%
  # Add underserved_tract layer
  addPolygons(
    data = underserved_tract_wgs84,
    fillColor = "red",
    weight = 0.5,
    opacity = 1,
    color = "darkred",
    fillOpacity = 0.7,
    popup = ~paste("ID:", ID),
    group = "Underserved communities"
  ) %>%
  # Add tb_flowlines layer
  addPolylines(
    data = tb_flowlines_wgs84,
    color = "blue",
    weight = 1,
    opacity = 0.7,
    group = "Flowlines"
  ) %>%
  # Add layer control
  addLayersControl(
    overlayGroups = c("Projects in DBs with underserved communities",
                      "DBs with underserved communities",
                      "Underserved communities",
                      "Flowlines"),
    options = layersControlOptions(collapsed = FALSE)
  )
```

There are 6 projects occurring within a drainage basin containing underserved communities that are more than 1 mile away from the nearest community. We can investigate each one relative to the flowlines to determine if there is an underserved community downstream that could benefit from the project.

First, we can see a point in the midde of Hillsborough Bay. This point is actually 3 projects (IDs 226, 264 and 265) on Fantasy Island. Since this is offshore on an uninhabited island, we will not consider these projects as benefiting downstream communities.

Second, we can see two adjacent projects in eastern Clearwater (IDs 229 and 273). We can see that the projects are adjacent to the main tributary of the drainage basin, but the underserved community is located upstream of the projects. This projects will therefore not contribute direct benefits to the nearby underserved community.

Finally, there is a project in South St. Petersburg (ID 270) adjacent to Lake Maggiore. In this case, we can see that the benefits provided to the lake can be expected to flow downstream to a couple underserved communities north of the lake.

Keep only the project with expected benefits to downstream underserved communities and include the number of communities benefiting (2).

```{r}
summariesdb <- summariesdb %>%
  filter(ID == 270) %>%
  mutate(communities_downstream = 2)

# merge with NEPORT dataset
baseline5yr <- left_join(baseline5yr, summariesdb, by = 'ID') %>%
  mutate(communities_downstream = coalesce(communities_downstream, 0),
         communities_1mile = coalesce(communities_1mile.x, 0),
         underserved_communities = ifelse(communities_1mile > 0 | communities_downstream > 0, "Yes", "No"),
         Year = Year.x,
         Total_Project_Costs = Total_Proj.x,
         EPA320_Costs = EPA_Sectio.x) %>%
  rowwise() %>%
  select(ID, Year, Total_Project_Costs, EPA320_Costs, communities_1mile, communities_downstream, underserved_communities)
```

Get annual summaries for our measures of distributive equity.

```{r}
equitystats_annual <- as.data.frame(baseline5yr) %>%
  group_by(Year, underserved_communities) %>%
  summarise(Projects_N = n(),
            Project_Costs = sum(Total_Project_Costs),
            Section320_Funds = sum(EPA320_Costs))

totalstats_annual <- equitystats_annual %>%
  group_by(Year) %>%
  summarize(All_Projects = sum(Projects_N),
            All_Project_Costs = sum(Project_Costs),
            All_Section320_Funds = sum(Section320_Funds))
# by year
baselinestats_annual <- left_join(equitystats_annual, totalstats_annual, by = 'Year') %>%
  mutate(Pct_Projects = Projects_N/All_Projects * 100,
         Pct_Project_Costs = Project_Costs/All_Project_Costs * 100,
         Pct_Section320_Funds = coalesce(Section320_Funds/All_Section320_Funds * 100, 0),) %>%
  filter(underserved_communities == "Yes")

# over entire period
baselinestats_5yr <- baselinestats_annual %>%
  group_by(underserved_communities) %>%
  summarize(Projects_N = sum(Projects_N),
            Project_Costs = sum(Project_Costs),
            Section320_Funds = sum(Section320_Funds),
            All_Projects = sum(All_Projects),
            All_Project_Costs = sum(All_Project_Costs),
            All_Section320_Funds = sum(All_Section320_Funds)) %>%
  mutate(Pct_Projects = Projects_N/All_Projects * 100,
         Pct_Project_Costs = Project_Costs/All_Project_Costs * 100,
         Pct_Section320_Funds = Section320_Funds/All_Section320_Funds * 100)
```

If we look at the table for the entire 5 year period, we can see that 30% of habitat projects provided benefits to underserved communities. Additionally, 34% of all investment dollars and 3% of Section 320 funds went to projects benefiting underserved commmunities.

We can visualize the annual stats with ggplot, relative to our 40% target.

```{r}
plot1 <- ggplot(baselinestats_annual, aes(x = as.factor(Year), y = Pct_Projects)) +
  geom_bar(stat="identity") +
  xlab("Year") +
  ylab("Habitat Projects (%)") + 
  geom_hline(yintercept = 40, linetype = "dashed", color = "red") +
  ylim(0,100)
  
plot2 <- ggplot(baselinestats_annual, aes(x = as.factor(Year), y = Pct_Section320_Funds)) +
  geom_bar(stat="identity") +
  xlab("Year") +
  ylab("Section 320 Funds (%)") +
  geom_hline(yintercept = 40, linetype = "dashed", color = "red") +
  ylim(0,100)

plot3 <- ggplot(baselinestats_annual, aes(x = as.factor(Year), y = Pct_Project_Costs)) +
  geom_bar(stat="identity") +
  xlab("Year") +
  ylab("Total Project Costs (%)") +
  geom_hline(yintercept = 40, linetype = "dashed", color = "red") +
  ylim(0,100)

grid.arrange(plot1, plot2, plot3, ncol = 3, nrow = 1)

```

### Mapping benefits to communities

We can modify the steps above to identify the communities that are benefiting from TBEP activities. In this section, we will identify these communities and characterize the benefits flowing to them.

Load the underserved and overburdened later and intersect with the buffered project locations.

```{r}
load(file = 'data/under_over.RData')

under_over <- under_over %>%
  st_transform(crs = 6443)

# communities within 1 mile of habitat projects
under_over_1mile <- st_intersection(under_over, baseline5yr_buff)
```

We also need to manually add in the project that was over 1 mile but still had downstream benefits to two underserved communities, which we identified previously. You can see the IDs for that project (270) and the two communities (12103020600, 12103020500) in the map created earlier.

```{r}
project <- as.data.frame(baseline5yr_buff) %>%
  filter(ID == 270) %>%
  rename("ID.1" = "ID")

community1 <- as.data.frame(under_over) %>%
  filter(ID == 12103020600)

community2 <- as.data.frame(under_over) %>%
  filter(ID == 12103020500)

community1.1 <- cbind(community1, project)
community2.1 <- cbind(community2, project)

# communities downstream of habitat projects
under_over_downstream <- rbind(community1.1, community2.1) %>%
  mutate(Shape.x.x = NULL) %>%
  rename("Shape.x.x" = "geometry")

# merge datasets
under_over_benefits <- rbind(as.data.frame(under_over_1mile), under_over_downstream)

```

Create a table characterizing the underserved tracts that are benefiting from habitat projects for each year and merge with the polygon layer.

```{r}
summaries2 <- under_over_benefits %>%
  group_by(ID, Year, Restoratio) %>%
  summarise(Projects_N = n(),
            Project_Costs = sum(Total_Proj))

under_over_projects <- full_join(under_over, summaries2, by = 'ID') %>%
  rename("Restoration_Technique" = "Restoratio") %>%
  mutate(Restoration_Technique = coalesce(Restoration_Technique, "None"),
         Projects_N = coalesce(Projects_N, 0),
         Project_Costs = coalesce(Project_Costs, 0))
```

#### Distribution of benefits to underserved communities

Map the number of projects and the investment dollars spent on projects benefiting each underserved community during 2017-2021. You can see that past projects have provided more benefits to communities of south St. Petersburg and southeastern Tampa, with those in southeast and northern Tampa receiving greater investment. Such maps can help us identify where we can direct future efforts to more equitably distribute the benefits of TBEP activities across the watershed.

```{r}
under_over_summary1 <- under_over_projects %>%
  group_by(ID) %>%
  summarise(Projects_N = sum(Projects_N),
            Project_Costs = sum(Project_Costs))

# Transform to WGS84 for leaflet
under_over_summary1_wgs84 <- st_transform(under_over_summary1, crs = 4326)

# Create color palettes for the two maps
pal_projects <- colorNumeric(palette = "viridis", domain = under_over_summary1_wgs84$Projects_N)
pal_costs <- colorNumeric(palette = "plasma", domain = under_over_summary1_wgs84$Project_Costs)

# Create the projects map
projectmap <- leaflet(under_over_summary1_wgs84) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(
    fillColor = ~pal_projects(Projects_N),
    weight = 0.5,
    opacity = 1,
    color = "black",
    fillOpacity = 0.7,
    popup = ~paste("ID:", ID, "<br>Projects:", Projects_N)
  ) %>%
  addLegend("bottomright", 
            pal = pal_projects, 
            values = ~Projects_N,
            title = "No. of Projects",
            opacity = 1)

# Create the investment map
investmentmap <- leaflet(under_over_summary1_wgs84) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(
    fillColor = ~pal_costs(Project_Costs),
    weight = 0.5,
    opacity = 1,
    color = "black",
    fillOpacity = 0.7,
    popup = ~paste("ID:", ID, "<br>Investment:", scales::dollar(Project_Costs))
  ) %>%
  addLegend("bottomright", 
            pal = pal_costs, 
            values = ~Project_Costs,
            title = "Investment Dollars",
            opacity = 1)

# Display both maps side by side
sync(projectmap, investmentmap)
```

#### Environmental burdens facing these underserved communities

We can also look at the burdens facing the underserved communities that have benefited from these past projects. For example, the graph generated below shows how frequently TBEP habitat projects have benefited underserved communities facing pollution and other environmental justice issues. We can see that the majority of underserved communities benefiting from past TBEP projects have a high density of leaking underground storage tanks, are in close proximity to brownfields, and have higher risks of flooding under climate change.

```{r}
under_over_burdens <- under_over_benefits %>%
  select(ID, thresholdEJ_agloss, thresholdEJ_floodr, thresholdEJ_greens, thresholdEJ_pm2.5, thresholdEJ_trafic, thresholdEJ_wastew, 
         thresholdEJ_hwaste, thresholdEJ_ugtank, thresholdEJ_sfunds, thresholdEJ_redlin, thresholdEJ_phmine, thresholdEJ_bfield) %>%
  melt(id.vars = "ID", variable.name = "Burden", value.name = "Threshold_Met") %>%
  mutate(Threshold_Met = coalesce(Threshold_Met,0)) %>%
  distinct() %>%
  group_by(Burden) %>%
  summarise(Communities_N = sum(Threshold_Met)) %>%
  mutate(Burden = c("Climate Change: Ag Loss Risk",
         "Climate Change: Flood Risk",
         "Nature Deprivation",
         "Air Pollution: PM2.5",
         "Air Pollution: Traffic",
         "Water Pollution: Wasetwater",
         "Pollution: Hazardous Waste",
         "Pollution: Underground Tanks",
         "Pollution: Superfund Sites",
         "Legacy Effects: Redlining",
         "Pollution: Phosphate Mines",
         "Pollution: Brownfield Sites"))
  
ggplot(under_over_burdens, aes(x = reorder(Burden, -Communities_N), y = Communities_N)) +
  geom_bar(aes(fill = Burden), color = "black", stat = "identity") +
  xlab("\nEnvironmental Burden") +
  ylab("Number of Underserved Communities Benefiting from Projects\n") +
  theme(legend.position = "none") +
  coord_flip() + 
  scale_x_discrete(limits = rev) +
  scale_y_continuous(position = "right")
```

#### Types of activities undertaken in these underserved communities

You can make additional maps to characterize communities based on the restoration activities. For example, we can map the restoration techniques of projects that benefited underserved communities. The maps generated below show the underserved communities benefiting from planting (green), oyster gardening (black), shoreline enhancement/stabilization (brown), and invasive species control/removal (blue) activities.

```{r}
under_over_techniques <- full_join(under_over, summaries2, by = 'ID') %>%
  rename("Restoration_Technique" = "Restoratio") %>%
  mutate(Restoration_Technique = coalesce(Restoration_Technique, "None"),
         Projects_N = coalesce(Projects_N, 0),
         Project_Costs = coalesce(Project_Costs, 0)) %>%
  group_by(ID) %>%
  summarise(ShorelineEnhance = sum(Restoration_Technique == "Shoreline Stabilization/Enhancement"),
            VegetationBuffer = sum(Restoration_Technique == "Vegetation Buffer"),
            Planting = sum(Restoration_Technique == "Planting"),
            ReefConstruction = sum(Restoration_Technique == "Reef Construction - Natural Materials"),
            OysterGardening = sum(Restoration_Technique == "Oyster Gardening"),
            DebrisRemoval = sum(Restoration_Technique == "Debris Removal"),
            InvasivesControl = sum(Restoration_Technique == "Invasives Control/Removal - Vegetation"),
            Other = sum(Restoration_Technique == "Other"))

# Transform to WGS84
under_over_techniques_wgs84 <- st_transform(under_over_techniques, crs = 4326)

# Create individual leaflet maps for each technique
map_planting <- under_over_techniques_wgs84 %>%
  filter(Planting > 0) %>%
  leaflet() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(
    fillColor = "green",
    weight = 0.5,
    opacity = 1,
    color = "black",
    fillOpacity = 0.7,
    popup = ~paste("ID:", ID, "<br>Planting Projects:", Planting),
    group = "Planting"
  ) %>%
  addLayersControl(
    overlayGroups = c("Planting"),
    options = layersControlOptions(collapsed = FALSE)
  )

map_oyster <- under_over_techniques_wgs84 %>%
  filter(OysterGardening > 0) %>%
  leaflet() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(
    fillColor = "black",
    weight = 0.5,
    opacity = 1,
    color = "black",
    fillOpacity = 0.7,
    popup = ~paste("ID:", ID, "<br>Oyster Gardening Projects:", OysterGardening),
    group = "Oyster Gardening"
  ) %>%
  addLayersControl(
    overlayGroups = c("Oyster Gardening"),
    options = layersControlOptions(collapsed = FALSE)
  )

map_shore <- under_over_techniques_wgs84 %>%
  filter(ShorelineEnhance > 0) %>%
  leaflet() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(
    fillColor = "brown",
    weight = 0.5,
    opacity = 1,
    color = "black",
    fillOpacity = 0.7,
    popup = ~paste("ID:", ID, "<br>Shoreline Enhancement Projects:", ShorelineEnhance),
    group = "Shoreline Enhancement"
  ) %>%
  addLayersControl(
    overlayGroups = c("Shoreline Enhancement"),
    options = layersControlOptions(collapsed = FALSE)
  )

map_invasives <- under_over_techniques_wgs84 %>%
  filter(InvasivesControl > 0) %>%
  leaflet() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(
    fillColor = "blue",
    weight = 0.5,
    opacity = 1,
    color = "black",
    fillOpacity = 0.7,
    popup = ~paste("ID:", ID, "<br>Invasives Control Projects:", InvasivesControl),
    group = "Invasives Control"
  ) %>%
  addLayersControl(
    overlayGroups = c("Invasives Control"),
    options = layersControlOptions(collapsed = FALSE)
  )

sync(map_planting, map_oyster, map_shore, map_invasives)
```