exposure_to_lead_in_schools_nys.qmd

---
title: "NYS School Water Quality: Exposure to Lead"
author: "Hasim Engin"
date: "February 2, 2024"
execute:
  echo: true
  warning: false
  eval: false
bibliography: lead-references.bib
---

```{r, context="setup", include=FALSE}


rm(list=ls())

list.of.packages <- c("tidyverse","scales", "sf", "tmap","leaflet","raster","shiny","tidycensus","DT","dataset","shinyWidgets")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
#import libraries
library(tidyverse)
library(scales)
library (sf)
library(tmap)
library(leaflet)
library(tools)
library(raster)
library(shiny)
library(tidycensus)
library(DT)
library(datasets)
library(shinyWidgets)
```

```{r, context="data", echo=FALSE}

# this changes prepairing all the datasets for following sections
# quarto require all inputs in shiny application to be in the same chink and 
# the chunk should be context="data"
# this section will not be visible to user
# However each part of the process will be viisble to user in their respective section


##setwd("D:/Grid3/school_lead_data")
#school_lead_data<-"Lead_Testing_in_School_Drinking_Water_Sampling_and_Results_Compliance_Year_2016.csv"
get_data<-"https://raw.githubusercontent.com/ciesin-geospatial/TOPSTSCHOOL-module-1-water/main/Lead_Testing_in_School_Drinking_Water_Sampling_and_Results_Compliance_Year_2016.csv"
school_lead_df<-read_csv(url(get_data))

# extract xy coordinates
school_lead_df<-school_lead_df %>% mutate(location_temp=str_extract(Location,"(?<=\\().*(?=\\))")) %>% 
  separate(location_temp,c("lat","long"),sep=",")


# remove space from fields names
names(school_lead_df) <- names(school_lead_df) %>% stringr::str_replace_all("\\s","_")


# extract year fro date field. we are using Date_Results_Updated for the date
school_lead_df<-school_lead_df %>%  mutate(year=format(as.Date(Date_Results_Updated, format="%d/%m/%Y"),"%Y"))

school_lead_df<-school_lead_df %>%  mutate(lead_summary_by_school=case_when(
                                          Number_of_Outlets_above_15_ppb ==0~ "has lead < 15ppb",
                                          Number_of_Outlets_above_15_ppb >0~ "has lead >15ppb",
                                          TRUE ~ "no data"))

school_locations<-school_lead_df %>% filter(!is.na(lat)) %>% 
  # convert the dataset to spatial dataset by using WGS84 projection
  st_as_sf(coords = c("long", "lat"),crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") 

get_census_data<- function(state_abr){
  census_api_key("94dee3107dc86d32c8039896b244eeb4aad72eb7")
  county_census<-get_decennial(
  geography = "county", 
  variables = c("total_pop"="P1_001N","total_white"="P1_003N", "total_black"="P1_004N",
                "total_asian"="P1_006N","total_hispanic"="P2_002N"),
  year = 2020,
  state=state_abr,
  output = "wide") %>% dplyr::select(GEOID,total_pop,total_white,total_black, total_asian,total_hispanic)
# get school age (5-17) population estimates from ACS 
county_acs<- get_acs(
  geography = "county", 
  variables = c("white_5_9"="B01001A_004","white_10_14"="B01001A_005","white_15_17"="B01001A_006",
                "black_5_9"="B01001B_004","black_10_14"="B01001B_005","black_15_17"="B01001B_006",
                "asian_5_9"="B01001D_004","asian_10_14"="B01001D_005","asian_15_17"="B01001D_006",
                "hispanic_5_9"="B01001I_004","hispanic_10_14"="B01001I_005","hispanic_15_17"="B01001I_006"),
  year = 2020,
  state=state_abr,
  output = "wide",
  geometry = TRUE) %>% mutate(white_5_17=white_5_9E+white_10_14E+white_15_17E,
             black_5_17=black_5_9E+black_10_14E+black_15_17E,
             asian_5_17=asian_5_9E+asian_10_14E+asian_15_17E,
             hispanic_5_17=hispanic_5_9E+hispanic_10_14E+hispanic_15_17E)

county_population<-county_acs %>% left_join(county_census, by="GEOID")%>%  
  separate(NAME, c("county","state"),",")
return (county_population)}


### ===== selection by dropdowns ===========###

state_list<-datasets::state.abb
list_of_counties<-school_lead_df %>% pull(County) %>% unique()
list_of_counties<-append(list_of_counties,"All Counties")
field_list<-c("total_white","total_black","total_asian","total_hispanic",
              "white_5_17","black_5_17","hispanic_5_17","asian_5_17")
classification_methods<-c("sd", "equal", "pretty", "quantile", "jenks")
lead_field_list<-c("school_count","outlets_under_15_ppb","outlets_above_15_ppb")
lead_data_all_fields<-names(school_lead_df ) %>% unique()


select_state<-column(3, selectInput("state", label = "Select a state",
              choices = state_list, selected = "NY"))
select_county<-column(3, selectInput("county", label = "Select a county",
              choices = list_of_counties, selected = "Bronx"))
select_fields1<-column(3, selectInput("field1", label = "Select a field for map (left)",
              choices = field_list, selected = "total_white"))
select_fields2<-column(3, selectInput("field2", label = "Select a field for maps (right)",
              choices = field_list, selected = "white_5_17"))
select_classification_method<-column(3, selectInput("classification", label = "Classification method",
              choices = classification_methods, selected = "jenks"))
select_fields3<-column(3, selectInput("field3", label = "Select a field",
              choices = lead_field_list, selected = "school_count"))
select_fields4<-  column(4, selectInput("fields", label = "Select fields",
              choices = lead_data_all_fields, selected = "County",multiple = TRUE))

theme<-theme(axis.text=element_text(size=12, face="bold"),
             axis.title=element_text(size=13, face="bold"),
              plot.title=element_text(size=13, face="bold"),
              panel.background = element_rect(fill = "lightblue",colour = "grey")
             )

```

# Lead Testing in NYS School Drinking Water

## Overview

This use case explores the risks of exposure to lead via drinking water. Lead contamination is a serious issue that poses severe health risks and requires remedial action. In this lesson, we will analyze data on lead levels in [NY State schools](https://health.data.ny.gov/Health/Lead-Testing-in-School-Drinking-Water-Sampling-and/rkyy-fsv9/data) collected from 2016, 2017, 2018, and 2019 and compare it with population characteristics at the county level in New York State (NYS) to understand its impact. No “safe” levels of lead have been established, but we will discuss what level of lead can be detected. The lesson examines the sources of lead exposure and its adverse effects. The review also discusses the importance of data transparency, public participation, visualizing contamination as maps/graphs, and estimating population risks to address water quality issues.

## Learning Objectives

-   Learn about lead as a drinking water contaminant and its health impacts.
-   Visualize lead contamination data geographically and statistically.
-   Discuss the importance of safe drinking water in the context of lead contamination.
-   Understand and apply open data principles in the context of lead contamination data.
-   Explore different data classification methods (equal, quantile, natural breaks, standard deviation).

## Introduction

Access to clean and safe drinking water is significant to ensure public health [@NIEHS]. Drinking water contaminants may have both short-term and long-term negative health impacts. One such contaminant that can have detrimental effects is lead, which is particularly harmful to a child's development [@levallois2018]. Children's central nervous systems and cognitive function have been linked to harm from lead exposure, even at low levels [@lanphear2005]. Therefore, it is essential to address lead-contaminated water in schools and homes.

Lead exposure has no safe threshold; therefore, safety depends on ensuring that lead levels in water are below the legal thresholds set by the World Health Organization (WHO) and the US Environmental Protection Agency (EPA) [@EPA2024]. According to guidelines issued by the WHO, lead concentrations in drinking water should not exceed 10 ppb (Parts Per Billion) [@WHO2022], and 15 ppb is the action level set by the EPA [@EPA2024].

It is crucial to highlight that the degree of lead able to be detected depends upon many variables [@schock1990]. The main source of lead in drinking water is the corrosion of lead-containing plumbing materials [@EPA2024]. Older plumbing systems with lead pipes and solder may leak lead into drinking water, particularly in regions with acidic water [@EPA2024]. Furthermore, runoff from contaminated soil and industrial discharges can also introduce lead into water systems. Developing effective preventive and remedial strategies requires understanding the factors that contribute to lead contamination.

## LEAD IN NEW YORK SCHOOLS

New York State (NYS) Lead Testing in School Drinking Water dataset shows the school drinking water lead sampling and results information reported by each NYS public school and Boards of Cooperative Educational Services (BOCES) [(NYS Department of Health)](https://health.data.ny.gov/Health/Lead-Testing-in-School-Drinking-Water-Sampling-and/rkyy-fsv9/data). More information on the NYS dataset sampling is available [here](https://www.health.ny.gov/environmental/water/drinking/lead/lead_testing_of_school_drinking_water.htm).

Analysis of the dataset reveals that as of 2022, 1,864 schools had lead outlets testing higher than 15 ppb [@NYS2016]. While 527 schools finished their remediation, 1,851 schools reported taking remedial action. There are now 12 schools with outlets exceeding 15 ppb in operation, indicating possible continuous exposure. However, there are gaps in following up and documenting the corrective measures. More transparency is necessary for schools with high exposure to lead to address the hazards of lead contamination and implement improved repeated testing protocols. The New York State Department of Health provides guidelines, rules, and resources on lead testing and remediation in schools [@DOH2023]. However, it seems that there is currently a lack of financial and technical support for schools to handle lead hazards.

::: {.callout-tip style="color: #7d2748;"}
### Lead in the News

Thousands of people were exposed to dangerously high lead levels in their drinking water when the [Flint water crisis](https://www.nrdc.org/stories/flint-water-crisis-everything-you-need-know) broke out in 2014. A study by Virginia Tech researchers, through their resident-organized sampling to testing data of 252 homes, revealed that lead levels in the city had increased [@NRDC]. Over 17% of samples tested higher than the federal "action level" of 15 ppb, which calls for the need for corrective action. More than 40% had lead readings higher than 5 ppb, which the researchers deemed indicative of a "very serious" issue.

Even years after the crisis began, elevated lead levels remained in Flint's schools. An article by [The New York Times](https://www.nytimes.com/2019/11/06/us/politics/flint-michigan-schools.html) discusses how, in 2019, drinking water samples from 30 Flint school buildings still exhibit excessive lead levels. The elevated levels demonstrate remaining problems with a prolonged impact on children's health and development. Schools have an obligation to supply their pupils with clean drinking water. The Flint water crisis brought to light the long-term consequences of prolonged exposure to lead, especially for vulnerable groups such as children.

Lead exposure is cumulative, as noted in the article by the [Tampa Bay Times](https://projects.tampabay.com/projects/2018/investigations/school-lead/hillsborough-disclosure/). The duration of lead exposure in Flint has lasting impacts on the public, not only affecting physical health but also leading to psychological consequences of communities not being able to trust their drinking water [@brooks2021]. Rather than waiting for concerns to arise, schools can detect contamination issues early and take corrective action by implementing a lead testing program. Better learning outcomes are made possible by shielding children from lead exposure. The Flint water crisis made clear how crucial it is to conduct proactive lead testing, monitor the situation, and take prompt, corrective action. It also made clear how important it is to be transparent, involve the community, and address barriers to clean drinking water access.

In 2021, the Biden-Harris administration announced an ambitious [Lead Pipe and Paint Action Plan](https://www.whitehouse.gov/briefing-room/statements-releases/2021/12/16/fact-sheet-the-biden-harris-lead-pipe-and-paint-action-plan/). This comprehensive \$15 billion effort intends to promptly replace all lead service lines and pipes that are contaminating drinking water systems across the country. The plan has a provision of providing a lead remediation grant of \$9 billion to disadvantaged communities through the Water Infrastructure Improvements for the Nation Act (WIIN) program, including for schools and childcare centers at EPA.
:::

### Read data

To work with the NYS data, first we will read the NYS school lead testing results from 2016 to 2019. The dataset is hosted on a GitHub repository and we will read the dataset by using the dataset URL.

```{r, eval=FALSE}
# Dataset url on GitHub repository.
data_url<-"https://raw.githubusercontent.com/renastechschool/Python_tutorials/main/Lead_Testing_in_School_Drinking_Water_Sampling_and_Results_Compliance_Year_2016_formated.csv?token=GHSAT0AAAAAACNH7S3BJGTQXNGH4UPQCJI6ZNQB3VA"
# Read dataset. Input data wrapped by url mothod. This allows to read data from a url.
school_lead_df<-read_csv(url(data_url))

```

### Preparing the lead dataset for the analysis

All datasets require some pre-cleaning and formatting. In the section below, we will format field names. R does not like field names with spaces, so we need to convert space to an underscore "\_". Also, we need to extract the year from the date field for the next step of our work.

```{r, eval=FALSE}
# There are empty spaces in the field names which R does not like.
# Replace empty space with "_".
names(school_lead_df) <- names(school_lead_df) %>% stringr::str_replace_all("\\s","_")


# Extract the year from the date field. We are using Date_Results_Updated for the date.
school_lead_df<-school_lead_df %>%  mutate(year=format(as.Date(Date_Results_Updated, format="%d/%m/%Y"),"%Y"))

# Data reports lead level by outlets if a outlet lead level is above or under 15ppb.
# The line below categorized the schools if they have any outlet above or under 15ppb.
school_lead_df<-school_lead_df %>%  mutate(lead_summary_by_school=case_when(
                                          Number_of_Outlets_above_15_ppb ==0~ "has lead < 15ppb",
                                          Number_of_Outlets_above_15_ppb >0~ "has lead >15ppb",
                                          TRUE ~ "no data"))

```

### Get familiar with the dataset

Getting familiar with the dataset is the first step of an analysis. To understand the attributes, we will query data by geographic region (county) and different attributes (fields).

The code below creates a [Shiny app](https://shiny.rstudio.com/) that allows users to select a county and specific fields from a data frame (school_lead_df), and then it displays the corresponding data table based on the user's selections.

::: column-margin
::: {.callout-tip style="color: #5a7a2b;"}
## Coding Review

A [Shiny app](https://shiny.rstudio.com/) is a web application framework for R programming language that allows you to create interactive web applications directly from R code. It's part of the RStudio ecosystem and is widely used for creating interactive data visualizations, dashboards, and web-based tools without needing to know HTML, CSS, or JavaScript.
:::
:::

```{r, context="render"}


#| panel: fill
fluidPage(
      fluidRow(style = "padding-bottom: 30px;background-color:#f1f2f3;",
        select_county,select_fields4),
      fluidRow(column(12, DT::dataTableOutput("table"))))
```

```{r, context="server"}


#| context: server
output$table<- DT::renderDataTable({
if (input$county=="All Counties")
    {school_lead_df %>% dplyr::select(input$fields) }

  else(school_lead_df %>% filter(County==input$county)%>% dplyr::select(input$fields))
})
```

<br>

### Converting from tabular data to geospatial data

A dataset needs to have a geometry attribute to plot the data on a map or to conduct different spatial analyses. The NYS dataset has xy coordinates of schools. The xy coordinates will allow us to convert the tabular dataset to a spatial dataset.

We first need to address that the xy coordinates are not properly formatted. The coordinates are currently stored with school addresses, for example: 31-02 67 AVENUE Queens, NY 11364(40.74779141700003, -73.74551716499997). We need to extract the coordinates (40.74779141700003, -73.74551716499997) and store the value on each side of the comma as a separate field. The first number refers to the y coordinate (latitude), and the second number refers to the x coordinate (longitude).

While converting the data, we also need to know the projection of xy coordinates. XY coordinates can be in different projection systems. Projection information is typically stored in the metadata of a dataset. However, in the NYS dataset, there is not any metadata attached to the dataset.

The most commonly used geographic coordinate system is the [WORLD GEODETIC SYSTEM 1984 (WGS 84)](https://earth-info.nga.mil/index.php?dir=wgs84&action=wgs84). We will use the WGS84 projection to convert the NYS dataset to spatial data.

<br>

```{r ,context: "data"}

# In the datasets xy coordinates of schools are combined with school addresses 
# and recorded under the location field. 
# For example: 231-02 67 AVENUEQueens, NY 11364(40.74779141700003, -73.74551716499997)
# The xy coordinates, (40.74779141700003, -73.74551716499997), need to be extracted, 
# sparated by a ",", and saved under different fields names.
# First, we will extract the xy coordinates and save under "lat","long" fields.
school_lead_df<-school_lead_df %>% mutate(location_temp=str_extract(Location,"(?<=\\().*(?=\\))")) %>% 
  separate(location_temp,c("lat","long"),sep=",")

# Next, we will create geometry attributes by using lat long fields and the geographic coordinate system.
# This step allows us to map the school locations and do geospatial anaylysis.
school_locations<-school_lead_df %>%  filter(!is.na(lat)) %>% 
  # Convert the dataset to spatial dataset by using WGS84 projection.
  st_as_sf(coords = c("long", "lat"),crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") 

```

### Mapping the dataset

We will now map the school locations in New York State from the dataset.

```{r ,context: "render"}

#|panel: fill
fluidPage(
      fluidRow(style = "padding-bottom: 30px;background-color:#f1f2f3;",
               select_county),
      fluidRow(column(12,leafletOutput("map"))))

```

<br>

```{r ,context: "server"}
#| context: server

school_loc<-reactive({school_locations %>% filter(County==input$county)})

output$map<-renderLeaflet({

  
#  if (input$county=="All Counties"){school_loc <-school_locations
#    tmap<-tm_basemap(leaflet::providers$Esri.WorldImagery)+
#   tm_shape(school_loc, name="school locations") +
#   tm_dots(id="",col="County",palette="magma",
#         popup.vars=c("School name: "="School" ),
#         legend.show = FALSE)
# tmap_leaflet(tmap)
#  }
#   
#   else{school_loc <-school_locations %>% filter(County==input$county)
  
  tmap<-tm_basemap(leaflet::providers$Esri.WorldImagery)+
  tm_shape(school_loc(), name="school locations") +
  tm_dots(id="",col="County",palette="magma",
        popup.vars=c("School name: "="School" ),
        legend.show = FALSE)
tmap_leaflet(tmap)


})

```

<br>

### Obtain population data from the US Census Bureau

We will next pull population data from the US Census Bureau by using the Census Bureau's Application Programming Interface (API).

```{r , eval=FALSE }
#| context: "render"
#|panel: fill

fluidPage(
      fluidRow(style = "padding-bottom: 20px;background-color:#f1f2f3;",
               select_state,select_fields1, select_fields2, select_classification_method),

      fluidRow(column(6, plotOutput("plot1")),
                column(6, plotOutput("plot2"))),
      fluidRow(column(6, leafletOutput("map1")),
                column(6,leafletOutput("map2"))))

```

```{r , eval=FALSE}
#| context: server

# get total population from census 2020
county_census<-reactive({
  req(input$state)
  get_census_data(input$state) })


output$plot1<-renderPlot({

  county_census()%>%st_drop_geometry() %>%
    dplyr::select(contains("total")) %>% gather(key,value) %>%
    group_by(key) %>% summarise(Count=sum(value)) %>%
    filter(key!="total_pop") %>%
    mutate(percent = prop.table(Count),  prc=paste0( "%",round(percent*100, 1)," \n(", comma(Count/1000),"K)")) %>%
     ggplot(aes(x=key,y=Count,fill=key))+
    geom_bar(stat="identity")+
    geom_text(aes(label = prc),size=5, hjust=0.7,
    position=position_stack(0.9))+guides(fill=FALSE)+
     labs(x="",y=" ", fill="", title = "Total Population by Race")+
    coord_flip()+theme
})

output$plot2<-renderPlot({

    county_census() %>%st_drop_geometry() %>%
    dplyr::select("white_5_17","black_5_17","hispanic_5_17","asian_5_17") %>%
    gather(key,value) %>%
    group_by(key) %>% summarise(Count=sum(value)) %>%
    mutate(percent = prop.table(Count),  prc=paste0( "%",round(percent*100, 1)," \n(", comma(Count/1000),"K)")) %>%
     ggplot(aes(x=key,y=Count,fill=key))+
    geom_bar(stat="identity")+
    geom_text(aes(label = prc),size=5, hjust=0.7,
    position=position_stack(0.9))+guides(fill=FALSE)+
     labs(x="",y=" ", fill="", title = "Total Populaiton age 5-17 by Race")+
    coord_flip()+theme
})


output$map1<-renderLeaflet({

tmap<-tm_basemap(leaflet::providers$Esri.WorldImagery)+
  tmap::tm_shape(county_census(), name="NYS counties") +
  tm_polygons(col=input$field1,style=input$classification, palette="RdYlGn")
tmap_leaflet(tmap)

})

output$map2<-renderLeaflet({

tmap<-tm_basemap(leaflet::providers$Esri.WorldImagery)+
  tmap::tm_shape(county_census(), name="NYS counties") +
  tm_polygons(col=input$field2,style=input$classification, n=6,palette="RdYlGn")
tmap_leaflet(tmap)

})



```

<br>

### Explore your data spatially

```{r }
#| contex: render
#| panel: fill
fluidPage(
      fluidRow(style = "padding-bottom: 20px;background-color:#f1f2f3;",
              select_state,select_county),

      fluidRow(column(5,"", plotOutput("plot")),
                column(7,"", leafletOutput("map"))))
```

```{r }
#| context: server

counties_bry <- reactive({req(input$state)
                              get_census_data(input$state) })

school_loc <- reactive({school_locations %>%filter(County==input$county)})


output$plot<-renderPlot({
  if (input$county=="All Counties"){school_lead_df_<-school_lead_df}
  else{school_lead_df_<-school_lead_df %>% filter(County==input$county)}
  school_lead_df_ %>%
  group_by(lead_summary_by_school) %>% summarize(Count=n()) %>%
    mutate(percent = prop.table(Count),  prc=paste0( "   %",round(percent*100, 1)," (", comma(Count),")")) %>%
  ggplot(aes(x=lead_summary_by_school,
                      y=Count,fill=lead_summary_by_school))+
    geom_bar(stat="identity")+
  geom_text(aes(label = prc),size=5,
    position=position_stack(0.9))+guides(fill=FALSE)+
    scale_fill_manual(values=c("green3","red","grey"))+
     labs(x="",y="Count of School", fill="", title = "Count of schools based on lead status",
         caption="has lead < 15ppb: None of the outlets have lead over 15ppb\nhas lead >15ppb: At least a outlet has lead over 15ppb")+theme
})


output$map<-renderLeaflet({
tmap<-tm_basemap(leaflet::providers$Esri.WorldImagery)+
  tmap::tm_shape(counties_bry() , name="NYS counties") +tm_polygons("black",alpha=0,popup.vars=c("County Name :  "="county"))+
  tm_shape(counties_bry() ,name="NYS counties ") +  tm_borders("black", lwd = 2) +
  tm_shape(school_loc(), name="School locations") +
tm_dots(id="something",col="lead_summary_by_school",palette=c("has lead < 15ppb"="green","has lead < 15ppb"="red"),
        popup.vars=c("School: "="School" ),
        legend.show = FALSE)
tmap_leaflet(tmap)

})


```

<br>

------------------------------------------------------------------------

### Aggregate school locations in county boundary and compare with county population

Lastly, we will aggregate points into county boundaries.

```{r }
#| contex: render
#| panel: fill
 fluidPage(

      fluidRow(style = "padding-bottom: 20px;background-color:#f1f2f3;",

               select_state,select_fields1,select_fields3, select_classification_method),

      fluidRow(column(6,"", leafletOutput("map1")),

               column(6,"", leafletOutput("map2"))))
```

```{r}
#| context: server

# get total population from census 2020

county_bry<-reactive({
  req(input$state)

get_census_data(input$state)})

output$map1<-renderLeaflet({

 county_bry<- st_transform(county_bry(), crs(school_locations))

  school_locations_join<-st_join(school_locations,county_bry,join = st_intersects) %>%

  group_by(GEOID,county) %>% summarise(school_count=n(),

                               outlets_under_15_ppb=sum(Number_of_Outlets_under_15_ppb, na.rm =TRUE),

                               outlets_above_15_ppb=sum(Number_of_Outlets_above_15_ppb, na.rm=TRUE)) %>%

  dplyr::select(GEOID,county,school_count,outlets_under_15_ppb,outlets_above_15_ppb) %>% st_drop_geometry()

  county_with_school<-county_bry%>% left_join(school_locations_join, by="GEOID")

tmap<-tm_basemap(leaflet::providers$Esri.WorldImagery)+

  tmap::tm_shape(county_with_school, name="county boundaries") +

  tm_polygons(col=input$field1,style=input$classification, n=6,palette="YlOrRd")

tmap_leaflet(tmap)

})

output$map2<-renderLeaflet({

 county_bry<- st_transform(county_bry(), crs(school_locations))

  school_locations_join<-st_join(school_locations,county_bry,join = st_intersects) %>%

  group_by(GEOID,county) %>% summarise(school_count=n(),

                               outlets_under_15_ppb=sum(Number_of_Outlets_under_15_ppb, na.rm =TRUE),

                               outlets_above_15_ppb=sum(Number_of_Outlets_above_15_ppb, na.rm=TRUE)) %>%

  dplyr::select(GEOID,county,school_count,outlets_under_15_ppb,outlets_above_15_ppb) %>% st_drop_geometry()

  county_with_school<-county_bry%>% left_join(school_locations_join, by="GEOID")

tmap<-tm_basemap(leaflet::providers$Esri.WorldImagery)+

  tmap::tm_shape(county_with_school, name="county boundaries") +

  tm_polygons(col=input$field3,style=input$classification, n=6,palette="YlOrRd")

tmap_leaflet(tmap)

})


```

### ENSURING SAFE DRINKING WATER

Access to safe drinking water is a fundamental human right and a pillar of public health [@li2023]. However, the cases discussed underscore the critical need to address the risk of lead contamination. From New York to Flint, Michigan, these examples serve as reminders of the need for action, supported by transparent and accessible data.

A balance is needed between short and long-term action - promptly addressing risks in the short term by closing down contaminated water sources, improving monitoring, and transparency in reporting water quality issues [@EPA2024]. In the long term, managing the costs of major infrastructure changes required to prevent such crises and exposures, such as replacing lead service lines and school plumbing [@EPA2024]. Community collaboration, precise surveillance databases, and flexible adaptation strategies are critical in ensuring access to safe drinking water.

Maintaining transparency and exchanging information with communities regarding the outcomes of lead testing and mitigation strategies is essential; empowering the public with knowledge and access to data can catalyze action and foster leadership within communities to address these challenges. Educating individuals about water quality issues, both globally and locally, is paramount.

## In this lesson, you learned...

Congratulations! Now you should be able to:

-   Read into R to analyze a dataset.
-   Convert survey dataset to spatial data using xy coordinates.
-   Plot locations on a map.
-   Obtain population data from US Census Bureau by using census API.
-   Aggregate a dataset to a boundary.
-   Create your own map.

## Explore the data

<iframe src="https://topstschool.shinyapps.io/nyc-lead" width="100%" height="800" title="iframe test" data-external="1">

</iframe>