Shifra Dayak, shifrandayak@gmail.com, 202-215-7408
What the story’s about: This story is about trends in reported environmental spills from the 1960s until now in Connecticut. The story will examine how and why spills have changed over the years, and — zooming in on approximately the last 10 years — what places and people have been the most impacted by spills and why, as well as what implications those impacts and any follow-ups on the part of regulatory and response authorities have on environmental justice.
Why now: Environmental spills are fairly common and ongoing in Connecticut, with a recent oil spill closing down part of a major interstate highway, sewage spilling into waterways in the last year and more, so telling this story now can help contextualize those recent spills and inform audiences in Connecticut about how environmental hazards (or the lack thereof) in their communities fit into larger historical patterns/
Three sources:
Lori Saliby, the director of the state environment department’s emergency response and spill prevention division, for context on how spill reporting works and why we might be seeing some of the trends that appear in this dataset
State and local politicians in Connecticut that have been involved in environmental advocacy, some of whom are listed as part of the Center for Climate Integrity. In Hartford, which has seen the most spills per a chloropeth map at the end of this file, CCI lists at least six politicians
Residents in places that have experienced spills — what have their experiences been like? Do they feel that the cleanup and response process was adequate and took their community into account? Did the spill(s) impact their health? Ideally, at least a couple of these residents would be from places that have experienced high spill rates.
Potential impact: This story could help communities that have faced environmental hazards and injustice be brought to the forefront. If, through reporting, a lack of action on responding to these reported spills is exposed, the story could help remedy that and spark action from regulatory and response authorities.
Three things that surprised/interested me:
Reported spills peaked in 2001. In the analysis section of this markdown, there’s some exploration of what could have contributed to this, including tying in historical context about changing state hazardous waste regulations in late 2001 and 2002, as well as anxiety about anthrax during the 2001 attacks. I’d like to talk to some of the above sources to see if these theories are true, and if not, why else spills have declined steadily since then
While the top spill substances have relatively remained the same over the years, the share involving fuel oil has dropped significantly overtime. I’d like to know why that is
There’s a lot of variation between where the most spills have happened and where the highest rates of spills have been. Looking at just raw numbers shows patterns that line up with other environmental injustice maps, which point to central and southwestern Connecticut as facing lots of issues. But looking at per-1,000 rates shows us some random towns here and there that have seen disproportionate spills, which is surprising.
Context: The Connecticut Courant reported at the end of 2023 on sewage spills during the year and how they impacted waterways in the state. Most sewage-related spills are because of outdated sewer systems, according to the piece — an issue that the state’s environment department is working to remedy. That’s the most overarching spill-related reporting recently. Other stories have largely been about individual incidents and have mostly covered breaking news, rather than trends or implications of spills. But Connecticut Public reported last year on a new environmental justice screening tool for Connecticut that shows what communities and census tracts are most heavily impacted by environmental injustices. Notably, the map tool itself shows that the most hard-hit communities approximately line up with the parts of the state that have seen the most spills reported from 2010 until now.
Visuals:
Chart of total spills over time
Chart of spills broken down by type across the years
Maps of where the most spills have occurred in Connecticut (numbers-wise and per 1,000 residents)
Chart of top dischargers
Maps of biggest spills overlaid on income levels and non-white population percentage in different Census tracts
Behind-the-scenes video explaining my process for the code and visualizations
Other information to gather and questions left to ask:
What’s caused the decline in spill reports since 2001?
Are there any confounding factors that could have changed how spills have been reported throughout history? Most recently, the state’s energy and environment department reformed its spill regulations in March 2022, reducing the types of spills that need to be reported to the department. Did this and similar steps in the past impact reported spill numbers in any way?
Things to analyze and chart beyond this notebook
Are certain types of spills/substances more common in certain places?
What has the response record and rate for spills been like, and does it vary based on town or other factors?
How many spills have impacted waterways and what has the severity been like?
What are trends with corrective actions for spills? Outside of this dataset, have there been any legal responses or settlements related to any of these spills?
Mathematically, what’s the correlation between non-white population percentage and number of spills? What about the correlation between median household income and number of spills?
How does Connecticut compare to other states? Looking at recent years, is it doing better or worse than other parts of the nation in terms of spills?
Estimated delivery: Mid-July 2024
Here’s a step-by-step look at how I analyzed this data and gathered some major findings.
#Loading CSV and standardizing column header names
spills <- read_csv("spills.csv") |>
clean_names()## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
## dat <- vroom(...)
## problems(dat)## Rows: 198097 Columns: 30
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (26): Assigned to, Date Reported Time Reported, Release date and time, T...
## dbl (3): Year, Case No., Total Quantity Drums
## lgl (1): Responsible Party Accepts Responsibility (Y/N)
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Often, a small subset of companies is responsible for a lot of spills. I started by seeing if this was the case in Connecticut.
What entities have been at fault for spills between 1996 and 2022?
Before bringing this data into R, I ran it through OpenRefine to do some preliminary cleaning. The way that company/discharger names are stylized is very inconsistent, so I wanted a slightly more accurate count by combining some non-unique company names (e.g. EXXON MOBIL and exxon mobil). OpenRefine didn’t catch everything, so I wrote a couple lines of code to manually rename some items that seemed to be within the top 10 dischargers or so. I also standardized some discharger names that have changed over the years, such as companies that have now rebranded to become Eversource Energy. I also filtered out spills labeled as coming from certain dischargers or that didn’t have a discharger name attached to them — more of my reasoning is explained in comments below.
#Loading in and cleaning up
#Connecticut Light and Power does business as/is a subsidiary of Eversource, so I standardized all of those names
companies <- read_csv("companies_cleaned_spills.csv") |>
clean_names () |>
rename (responsible_party = responsibile_party_discharger) |>
mutate (responsible_party = case_when (
responsible_party == "cumberland farms" ~ "CUMBERLAND FARMS",
responsible_party == "Cumberland Farms" ~ "CUMBERLAND FARMS",
responsible_party == "Cumberland Farms Inc." ~ "CUMBERLAND FARMS",
responsible_party == "EVERSOURCE" ~ "EVERSOURCE ENERGY",
responsible_party == "Eversource" ~ "EVERSOURCE ENERGY",
responsible_party == "Eversource Energy" ~ "EVERSOURCE ENERGY",
responsible_party == "eversource" ~ "EVERSOURCE ENERGY",
responsible_party == "CONNECTICUT LIGHT & POWER" ~ "EVERSOURCE ENERGY",
responsible_party == "CT DOT" ~ "CT TRANSIT",
responsible_party == "ct dot" ~ "CT TRANSIT",
responsible_party == "DOT" ~ "CT TRANSIT",
responsible_party == "alliance energy" ~ "alliance global",
responsible_party == "Alliance Energy" ~ "alliance global",
responsible_party == "alliance" ~ "alliance global",
responsible_party == "ALLIANCE ENERGY" ~ "alliance global",
responsible_party == "CL AND P" ~ "EVERSOURCE ENERGY",
responsible_party == "cl and p" ~ "EVERSOURCE ENERGY",
responsible_party == "C. L. & P" ~ "EVERSOURCE ENERGY",
responsible_party == "CL&P" ~ "EVERSOURCE ENERGY",
responsible_party == "UNITED ILLUMINATING" ~ "United Illuminating",
responsible_party == "p&w" ~ "PRATT & WHITNEY",
responsible_party == "P&W" ~ "PRATT & WHITNEY",
responsible_party == "NORTHEAST UTILITIES" ~ "EVERSOURCE ENERGY",
responsible_party == "cwpm" ~ "CWPM",
responsible_party == "cwpm" ~ "CWPM",
responsible_party == "auto tlc" ~ "AUTOMATIC TLC",
responsible_party == "TLC" ~ "AUTOMATIC TLC",
TRUE ~ responsible_party
))## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
## dat <- vroom(...)
## problems(dat)## Rows: 198097 Columns: 30
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (26): Assigned to, Date Reported Time Reported, Release date and time, T...
## dbl (3): Year, Case No., Total Quantity Drums
## lgl (1): Responsible Party Accepts Responsibility (Y/N)
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#Grouping by dischargers
companies <- companies |>
group_by (responsible_party) |>
summarise (
number_of_spills = n()
) |>
arrange(desc(number_of_spills))
#Taking out spills that don't have a named discharger, since we're only looking at companies
companies_no_unknown <- companies |>
filter(!str_detect(responsible_party, regex('UNKNOWN|unk|N/A|same|as above|UNDETERMINED', ignore_case=TRUE)))
#Taking out spills attributed to motor vehicle accidents, since most don't involve a company
companies_no_unknown <- companies_no_unknown |>
filter(!str_detect(responsible_party, regex('M.V.A.|MVA|M.V.A|PRIVATE VEHICLE|MV ACCIDENT|PRIVATE AUTO|PRIVATE MV', ignore_case=TRUE)))
#Taking out spills attributed to "saa," since that broadly encompasses all state administrative agencies rather than one specific company
companies_no_unknown <- companies_no_unknown |>
filter(!str_detect(responsible_party, regex('saa', ignore_case=TRUE)))
#Exporting as a csv
companies_no_unknown <- companies_no_unknown |>
slice(1:10) |>
write_csv("top_company_dischargers.csv")Here’s a chart showing the top dischargers. Numbers might be slightly off, since we aggregated all subsidiaries and companies that used to go by different names but have since merged or rebranded. Taking that into account, it seems like Eversource and its subsidiaries are polluting/spilling far more than any other single entity in the state.
cat('<div class="flourish-embed flourish-chart" data-src="visualisation/18258842"><script src="https://public.flourish.studio/resources/embed.js"></script></div>')Eversource, Global Partners, United Illuminating, ExxonMobil and Automatic TLC are all energy or gas companies of some kind, making the top polluters overwhelmingly related to electricity, energy or oil and gas. Cumberland Farms is a gas station and convenience store, so depending on the nature of their spills, they could also be factored into this count. I’d say this is fairly unsurprising — a lot of spills tend to be oil and fuel-related.
I also noticed that lots of the spill incidents did not have a discharger name recorded. The code below shows us that more than 56% of total spills recorded between 1996 and 2022 aren’t attached to a specific responsible party. Talking to sources about the implications of this would be interesting — does this mean that responsible parties aren’t held accountable? Why are dischargers so often not recorded?
unknown_spills <- spills |>
filter(is.na(responsibile_party_discharger) | str_detect(responsibile_party_discharger, regex("UNKNOWN|UNK|unknown|Unknown|N/A|UNDETERMINED|unk.|UNKNOWN POLLUTER|UNK.|unknown at this time", ignore_case = TRUE)))
(nrow(unknown_spills) / nrow(spills)) *100## [1] 56.62024What parts of Connecticut have seen the most impacts from spills?
Let’s move on to grouping by place, to see if there are trends in where spills have been reported. To narrow down the scope, let’s look at roughly the last decade of recorded data (2010 and after).
#Filtering for 2010 and after
after_2010_spills <- spills |>
filter (year > 2009)
#Standardizing stylization of place names
after_2010_spills$town_of_release <- toupper(after_2010_spills$town_of_release)
#Cleaning up some messy/wrong names
after_2010_spills <- after_2010_spills |>
mutate(
town_of_release = case_when (
town_of_release == "GRONTON" ~ "GROTON",
town_of_release == "GROTOM" ~ "GROTON",
town_of_release == "MIFORD" ~ "MILFORD",
town_of_release == "MILLFORD" ~ "MILFORD",
town_of_release == "MONTOVILLE" ~ "MONTVILLE",
town_of_release == "OOLD SAYBROOK" ~ "OLD SAYBROOK",
town_of_release == "WATERTOWM" ~ "WATERTOWN",
town_of_release == "NORTHFORD" ~ "NORTH BRANFORD",
.default = town_of_release
)
)
#Grouping by place
spills_by_town <- after_2010_spills |>
group_by (town_of_release) |>
summarise (
reports_by_town = n()
) |>
arrange(desc(reports_by_town)) |>
slice(1:170)
#Exporting as CSV so we can map in Flourish
spills_by_town |>
write_csv("spills_by_town.csv")To clean up the data by town (in the chunk of code above), I cross-compared with name keys in a GeoJSON file of Connecticut municipalities, and took out any instances of towns that didn’t correspond to the shapefile data (or spills that had no town attached to them). I uploaded that GeoJSON data and the above data of spills by town to Flourish to create a custom map.
Here’s the resulting map of spills by town. More spills were reported from towns in central Connecticut and southwestern Connecticut.
cat('<div class="flourish-embed flourish-map" data-src="visualisation/18244615"><script src="https://public.flourish.studio/resources/embed.js"></script></div>')That’s just with raw numbers. To see what the impact on people is, we can further expand this map by normalizing the number of spills per 1,000 residents in each town. I used the same GeoJSON data, paired with data below.
#Reading in and joining a table that has populations for each town
town_populations <- read_csv("pop_towns2022.csv") |>
clean_names () |>
rename (town_of_release = town)## Rows: 169 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Town
## num (1): Population
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#Making joining easier by getting rid of the N/A row in the spills by town table and renaming
spills_by_town <- spills_by_town |>
arrange(town_of_release)
#Joining and calculating spills per 1000 people
spills_per_thousand <- spills_by_town |>
inner_join(town_populations, by=c("town_of_release")) |>
mutate (per_thousand_spill_rate = (reports_by_town/population)*1000) |>
arrange(desc(per_thousand_spill_rate))
#Exporting
spills_per_thousand |> write_csv("spills_per_thousand.csv")Here’s a map that shows those results charted.
cat ('<div class="flourish-embed flourish-map" data-src="visualisation/18244185"><script src="https://public.flourish.studio/resources/embed.js"></script></div>')The patterns look different — the southwestern and central towns that saw the largest numbers of spills didn’t necessarily see the highest spill rates. The highest rates were spread out across the state, with Franklin, Stafford and Bloomfield at the top of the list. This would be another interesting thing to ask sources about. Why are there so many spills in these places, relative to the number of people who live there? Some higher rates make sense — Groton, for instance, is home to a large U.S. Navy submarine contracting company. But others, like Franklin (a small rural town) might point to less obvious answers.
I was curious about tying this into environmental justice, so I decided to contrast impacts on places with the makeup of those places.
What’s the demographic makeup of parts of the state that have experienced spills?
For this, I’m going to look at how non-white percentage and median household income might correlate with the amount of spills by town. I’m using the Census API to load in data on these variables, starting with non-white percentage (which I had to calculate manually, since the Census Bureau does not provide an aggregate measure of non-white population).
census_vars <- load_variables(2021, "acs5", cache=TRUE)
#Loading in data for white population
ct_demographics <- get_acs(geography = "tract",
variables = c(whitepop ="B02001_002"),
state="CT",
year= 2021)## Getting data from the 2017-2021 5-year ACS#Loading in data for total population of each census tract
ct_population <- get_acs(geography = "tract",
variables = c(totalpop = "B01003_001"),
state = "CT",
year = 2021)## Getting data from the 2017-2021 5-year ACS#Joining the two tables, cleaning up column names, and removing any columns we don't need
ct_demographics <- ct_demographics |>
inner_join(ct_population, by=c("GEOID")) |>
clean_names() |>
select(-moe_x, -name_y, -moe_y, -variable_x, -variable_y) |>
rename (census_tract = name_x,
white_pop = estimate_x,
total_pop = estimate_y
)
#Calculating the total non-white population percentage
ct_demographics <- ct_demographics |>
mutate (
white_pct = (white_pop/total_pop)*100
) |>
mutate(
non_white_pct = (100-white_pct)
)
#Exporting
ct_demographics |>
write_csv("ct_demographics.csv")I want to overlay those non-white percentages for each tract with some of the largest spills that have happened in Connecticut, to see if there are connections. After looking at the data arranged by spill, I decided the best way to define “top” was anything recorded as being 10,000 gallons or more. Since the demographic data is from the most recent five-year Census, I’m going to filter the spill data by that time period as well, so the timeframes match up. Then, I’m going to geocode the addresses of the top spills from that period so we can see them on a map. Since the addresses are entered in inconsistent formats, I’m just going to look up each address and manually enter the latitude and longitude values into the dataset.
#Selecting the years we want to look at
top_2017_2021_spills <- spills |>
filter(year %in% c(2017, 2018, 2019, 2020, 2021)) |>
arrange(total_quantity_gallons)
#Standardizing numbers in the quantity column
top_2017_2021_spills <- top_2017_2021_spills |>
mutate(total_quantity_gallons = gsub("[+\\-~]", "", total_quantity_gallons)) |>
mutate(total_quantity_gallons = as.numeric(total_quantity_gallons))## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `total_quantity_gallons = as.numeric(total_quantity_gallons)`.
## Caused by warning:
## ! NAs introduced by coercion#Figuring out the top spills
top_2017_2021_spills <- top_2017_2021_spills |>
filter(total_quantity_gallons > 9000) |>
arrange(desc(total_quantity_gallons))
#Taking out one inconsistently-reported spill
top_2017_2021_spills <- top_2017_2021_spills |>
filter(!str_detect(release_substance, "CHILLER WATER 1500 GALLONS OF NUTRAL PH WATER"))
#Adding lat and long values
top_2017_2021_spills <- top_2017_2021_spills |>
mutate (latitude = case_when (
location_of_reported_release == "1082 north benson road, McInness Road runs parralel to stream" ~ "41.16168860651233",
location_of_reported_release == "263 farmington ave" ~ "41.732834",
location_of_reported_release == "200 bloomfield ave" ~ "41.79582769811641",
location_of_reported_release == "11 GLASTONBURY AVE" ~ "41.665117241187346",
location_of_reported_release == "Dayville 56 Alexander Parkway" ~ "41.872053930599705",
location_of_reported_release == "55 SHELLAND ST" ~ "41.22265023721389",
location_of_reported_release == "239 West Service Rd." ~ "41.79738307648949",
location_of_reported_release == "6900 MAIN ST" ~ "41.25207965073338",
location_of_reported_release == "6 HILLCREST DR" ~ "41.559327103632825",
location_of_reported_release == "400 main street" ~ "41.74789958168573",
location_of_reported_release == "1254 Laurel Ridge Drive" ~ "41.59475217004565",
location_of_reported_release == "206 garfield street" ~ "41.693631187132574",
location_of_reported_release == "6 howard ave" ~ "41.16300933956072",
location_of_reported_release == "2 Forbes Street" ~ "41.77641331928542",
location_of_reported_release == "149 west ville ave" ~ "41.3967166877614",
location_of_reported_release == "94 Baldwin Ave." ~ "41.535345140683674",
location_of_reported_release == "Rt. 189 and Rt. 44 - University of Hartford - Reich Family Pavilion" ~ "41.80195983049168"
)) |>
mutate(longitude = case_when (
location_of_reported_release == "1082 north benson road, McInness Road runs parralel to stream" ~ "-73.25377220245423",
location_of_reported_release == "263 farmington ave" ~ "-72.791318",
location_of_reported_release == "200 bloomfield ave" ~ "-72.71465723160469",
location_of_reported_release == "11 GLASTONBURY AVE" ~ "-72.63791540397118",
location_of_reported_release == "Dayville 56 Alexander Parkway" ~ "-71.89576719880588",
location_of_reported_release == "55 SHELLAND ST" ~ "-73.098632516012",
location_of_reported_release == "239 West Service Rd." ~ "72.65604183843769",
location_of_reported_release == "6900 MAIN ST" ~ "-73.09598977724993",
location_of_reported_release == "6 HILLCREST DR" ~ "72.86271589083927",
location_of_reported_release == "400 main street" ~ "-72.63383413280191",
location_of_reported_release == "1254 Laurel Ridge Drive" ~ "-72.48578458186292",
location_of_reported_release == "206 garfield street" ~ "-72.732589905816",
location_of_reported_release == "6 howard ave" ~ "-73.20855214061626",
location_of_reported_release == "2 Forbes Street" ~ "-72.60675027697646",
location_of_reported_release == "149 west ville ave" ~ "-73.47991147515009",
location_of_reported_release == "94 Baldwin Ave." ~ "-73.03074004630756",
location_of_reported_release == "Rt. 189 and Rt. 44 - University of Hartford - Reich Family Pavilion" ~ "-72.71653747512859"
))
#Exporting
top_2017_2021_spills |>
write_csv("top_spills.csv")Here’s the resulting map. We can see that most of the largest spills seem to be clustered in and around areas that also have higher non-white population percentages. Fewer areas with smaller non-white populations contain those same clusters of big spills.
cat('<div class="flourish-embed flourish-map" data-src="visualisation/18287818"><script src="https://public.flourish.studio/resources/embed.js"></script></div>')Note that this may not be 100% accurate — the way these spills are reported is messy. As shown in the code above, I removed one spill because the amount noted in the release_substance column was very different than the amount entered in the quantity column. It’s possible there were other badly-entered spills, but I’m creating this map under the assumption that most reported values are accurate and do genuinely represent the biggest spills during the five-year Census period.
What about when we compare those same major spills to socioeconomic data? Here’s the same process repeated, but with median household income this time instead of non-white population.
#Loading in data and selecting relevant columns
ct_income <- get_acs(geography = "tract",
variables = c(median_inc ="B19013_001"),
state="CT",
year= 2021) |>
clean_names() |>
rename (median_income = estimate) |>
select(-moe, -variable)## Getting data from the 2017-2021 5-year ACS#Exporting
ct_income |>
write_csv("ct_median_income.csv")After using the same process of creating a custom map in Flourish with Census tract boundaries, here’s the result.
cat('<div class="flourish-embed flourish-map" data-src="visualisation/18288615"><script src="https://public.flourish.studio/resources/embed.js"></script></div>')The correlation doesn’t appear as strong here, because the major spills are spread across tracts with varying household incomes. But there are definitely still some notable trends. Namely, the richest part of the state (the southwest corner) doesn’t house as many spills as areas like Hartford, which have median household incomes on the lower side.
Why can’t we calculate actual correlation to confirm the relationship between race and spills, or income and spills? That’s because the Census Bureau doesn’t provide data on those variables at a municipality level. The spill data from Connecticut’s environment department, on the other hand, is all recorded on a municipality level, so there’s no way to compare the two.
We could mutate the dataset to assign each of the spills to its respective county and then compare that to county-level Census data. But Connecticut only has eight counties, and my evaluation of Census tracts showed that there’s significant variation of non-white population and median household income within these counties. So aggregating data by county would have provided a look at the data that’s far too broad to see any nuances.
With more time, I could clean the formatting of the addresses and use a geocoding API to assign each spill to its respective Census tract, and then calculate correlation with tract-level data.
To end, I wanted to get an overview of spills generally.
What have spill numbers looked like between 1996 and 2022?
#Taking out 2024 because not enough data has been recorded yet, and then grouping by year to calculate how many spills reported per year
spills_by_year <- spills |>
filter (year != "2024") |>
group_by(year) |>
summarise (
spill_number = n()
)#Exporting the new grouped table
spills_by_year |> write_csv("spills_by_year.csv")Charting this shows us that spills peaked in 2001 and have been on a relatively steady decline since then, despite some minor variation from year to year.
cat('<div class="flourish-embed flourish-chart" data-src="visualisation/18240536"><script src="https://public.flourish.studio/resources/embed.js"></script></div>')I was also curious what the biggest types of spills recorded in this dataset have been, and how they may have changed over the years.
What are the trends in top substances spilled over the years?
spills_by_type <- spills |>
group_by(release_substance) |>
summarise (
substance_type = n()
) |>
arrange(desc(substance_type))Grouping by release substance (above) confirms that those same three substances have made up the most spills across all years in this dataset. Now, let’s break this down further and group by substance type and year, to see how some of the top substances spilled may vary across this dataset. One possible hole: the way that substances are reported in this dataset is very non-standardized, so these numbers could be undercounting the actual number of spills in each substance category. For instance, most fuel oil spills are written as “#2 FUEL OIL,” but some may be written as “#2 fuel oil” and thus not counted. Cleaning the data and then redoing this count will make it more accurate, but this is just a rough estimate.
#Fuel oil
fuel_oil <- spills |>
filter(release_substance == "#2 FUEL OIL") |>
group_by (year) |>
summarise (
fuel_oil = n()
) |>
arrange(year)
#Antifreeze
antifreeze <- spills |>
filter(release_substance == "ANTIFREEZE") |>
group_by (year) |>
summarise (
antifreeze = n()
) |>
arrange(year)
#Transformer oil
transformer_oil <- spills |>
filter(release_substance == "TRANSFORMER OIL") |>
group_by (year) |>
summarise (
transformer_oil = n()
) |>
arrange(year)
#Gasoline
gasoline <- spills |>
filter(release_substance == "GASOLINE") |>
group_by (year) |>
summarise (
gasoline = n()
) |>
arrange(year)
#Diesel fuel
diesel_fuel <- spills |>
filter(release_substance == "DIESEL FUEL") |>
group_by (year) |>
summarise (
diesel_fuel = n()
) |>
arrange(year)
#Joining all of these
substance_dfs <- list(antifreeze, diesel_fuel, fuel_oil, gasoline, transformer_oil)
substances_across_years <-reduce(substance_dfs, left_join, by = "year")
#Calculating the numbers as a percentage of all spills
substances_pct_total <- substances_across_years |>
inner_join(spills_by_year, by=c("year")) |>
rename (total_spills = spill_number) |>
mutate (antifreeze = (antifreeze/total_spills)*100,
diesel_fuel = (diesel_fuel/total_spills)*100,
fuel_oil = (fuel_oil/total_spills)*100,
gasoline = (gasoline/total_spills)*100,
transformer_oil = (transformer_oil/total_spills)*100
)
#Renaming column names for more clarity
substances_pct_total <- substances_pct_total |>
rename (
antifreeze_pct = antifreeze,
diesel_fuel_pct = diesel_fuel,
fuel_oil_pct = fuel_oil,
gasoline_pct = gasoline,
transformer_oil_pct = transformer_oil
)
#Exporting
substances_pct_total |>
write_csv("substances_pct_total.csv")When we chart by substance by year, we can see that patterns have remained relatively the same — except for fuel oil, which no longer makes up as large a share of total reported spills as it used to. There were some years where antifreeze and transformer oil also peaked, making up the biggest share of spills. Possible reasons for these variations could be another interesting thing to ask sources about.
cat('<div class="flourish-embed flourish-chart" data-src="visualisation/18242246"><script src="https://public.flourish.studio/resources/embed.js"></script></div>')Also, if we look at 2001, fuel oil related spills are at their second-highest. But the other top substances don’t peak — so these top substances aren’t the only thing causing the massive peak of reports in 2001, evidently. So what is causing the peak?
What could have contributed to the peak of reported spills in 2001?
One big factor that could have caused the peak in that year is the 2001 anthrax scare.
As compared to the year immediately following (2002), a bigger chunk of all spills reported in 2001 seem to be related to anthrax or anthrax scares, with keywords like “powder” and “white substance.” See below for the code I used to figure that out:
#Filtering for year, grouping by type of substance, filtering for relevant release substance, then calculating the percentage of those that make up the total reports in 2001
spill_types_2001 <- spills |>
filter (year == "2001") |>
group_by(release_substance) |>
summarise (
reports_by_substance = n()
) |>
arrange(desc(reports_by_substance))
spill_types_2001 |>
mutate(total_reports = sum(reports_by_substance)) |>
filter(str_detect(release_substance, "WHITE | SUSPICIOUS | ANTHRAX")) |>
mutate(total_anthrax_related_pct = (sum(reports_by_substance))/total_reports*100)## # A tibble: 111 × 4
## release_substance reports_by_substance total_reports total_anthrax_relate…¹
## <chr> <int> <int> <dbl>
## 1 WHITE POWDER 305 10421 5.58
## 2 UNKNOWN WHITE POWD… 117 10421 5.58
## 3 WHITE POWDER SUBST… 13 10421 5.58
## 4 WHITE POWDERY SUBS… 9 10421 5.58
## 5 SUSPECTED ANTHRAX 7 10421 5.58
## 6 UNKNOWN WHITE SUBS… 5 10421 5.58
## 7 WHITE POWDER IN EN… 4 10421 5.58
## 8 WHITE SUBSTANCE 4 10421 5.58
## 9 ? OF ANTHRAX 3 10421 5.58
## 10 ENVELOPE WITH WHIT… 3 10421 5.58
## # ℹ 101 more rows
## # ℹ abbreviated name: ¹total_anthrax_related_pct#Doing the same thing for 2002, to see if anthrax was less of a concern in that year
spill_types_2002 <- spills |>
filter (year == "2002") |>
group_by(release_substance) |>
summarise (
reports_by_substance = n()
) |>
arrange(desc(reports_by_substance))
spill_types_2002 |>
mutate(total_reports = sum(reports_by_substance)) |>
filter(str_detect(release_substance, "WHITE | SUSPICIOUS | ANTHRAX")) |>
mutate(total_anthrax_related_pct = (sum(reports_by_substance))/total_reports*100)## # A tibble: 27 × 4
## release_substance reports_by_substance total_reports total_anthrax_relate…¹
## <chr> <int> <int> <dbl>
## 1 WHITE POWDER 53 8812 1.28
## 2 UNKNOWN WHITE POWD… 30 8812 1.28
## 3 WHITE POWDERY SUBS… 3 8812 1.28
## 4 UNKNOWN WHITE POWD… 2 8812 1.28
## 5 WHITE POWDER IN AN… 2 8812 1.28
## 6 WHITE POWDER SUBST… 2 8812 1.28
## 7 GASOLINE & WHITE S… 1 8812 1.28
## 8 SUBSTANCE TURNED O… 1 8812 1.28
## 9 UNKNOWN GREYISH WH… 1 8812 1.28
## 10 UNKNOWN SUSPICIOUS… 1 8812 1.28
## # ℹ 17 more rows
## # ℹ abbreviated name: ¹total_anthrax_related_pctPer the code results, anthrax-related reports in 2001 made up more than 5 percent of total reports that year, whereas anthrax-related reports in 2002 made up just over 1 percent of total reports that year, showing a significant difference and lending support to the hypothesis that anthrax reports drove at least part of the massive peak in 2001 reports.