Philadelphia Shooting Victims Analysis

̥ # Introduction Gun violence is a significant issue in the city of Philadelphia. There were 771 shooting incidents in the city in 2024 per the Philadelphia Police Department. The goal of this project is to utilize the data provided by the Philadelphia Police Department via OpenDataPhilly on shooting incidents to build a picture of the fatal shooting incidents occurring in the city. Who is most likely to be the victim of a fatal incident (age, race, sex)? What conditions of an incident make it more likely to be fatal (wound type)? When and where are the most fatal incidents occurring (time of day/year, location)?

Understanding the trends of fatal gun violence incidents can assist in painting a clearer picture of the problem and point towards potential solutions. Understanding the time and location of the most dangerous shootings provides guidance for police on where to focus their efforts. Understanding who is most at risk of a fatal shooting guides social and health services towards populations that may benefit from interventions to assist shooting victims.

Data Overview

The data is collected by the Philadelphia Police Department and provided via OpenDataPhilly. Below, we will load the data and gather the available variables for answering our questions.

Looking at the available fields we have a rich potential set of data to answer our questions. The age, race, sex and latino fields can provide a profile of the most frequent shooting victims. Wound can inform which types of injuries are most often fatal. Location, time and date can be used to determine when and where fatal incidents occur most often. These fields will be compared to the fatal column in order to determine fatality incidence.

# Install necessary packages (run this only once)
# install.packages(c("tidyverse", "sf", "leaflet", "ggplot2", "viridis", "readr", "sp", "leaflet.extras", "readxl", "htmltools", "htmlwidgets", "shiny", "dplyr"))

#install.packages("RColorBrewer")

# Load the packages
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(sf)

## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE

library(leaflet)
library(ggplot2)
library(viridis)

## Loading required package: viridisLite

library(readr)
library(sp)
library(leaflet.extras)
library(htmltools)
library(htmlwidgets)
library(shiny)
library(dplyr)
library(RColorBrewer)
library(jsonlite)

## 
## Attaching package: 'jsonlite'
## 
## The following object is masked from 'package:shiny':
## 
##     validate
## 
## The following object is masked from 'package:purrr':
## 
##     flatten

library(tidyverse)
library(readxl)

shootings.raw = read_delim(
  "./philly_shootings.csv",
  delim=",",
  show_col_types = FALSE,
) %>% glimpse()

## Rows: 16,732
## Columns: 25
## $ the_geom             <chr> "0101000020E61000001049CE6BA7C652C0F0FF5490EA0344…
## $ the_geom_webmercator <chr> "0101000020110F0000790F72E295E45FC1A02823583D9452…
## $ objectid             <dbl> 634246, 634247, 634248, 634249, 634250, 634251, 6…
## $ year                 <dbl> 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2…
## $ dc_key               <dbl> 201502062409, 201503069342, 201512041905, 2015120…
## $ code                 <dbl> 3400, 400, 400, 400, 400, 400, 400, 400, 3400, 40…
## $ date_                <dttm> 2015-09-26, 2015-11-08, 2015-05-30, 2015-12-11, …
## $ time                 <time> 05:07:00, 12:50:00, 14:11:00, 18:00:00, 19:40:00…
## $ race                 <chr> "B", "A", "B", "B", "B", "B", "B", "B", "B", "B",…
## $ sex                  <chr> "M", "M", "M", "M", "F", "F", "F", "M", "F", "M",…
## $ age                  <dbl> 43, 36, 27, 23, 15, 21, 26, 33, 26, 43, 21, 35, 2…
## $ wound                <chr> "Chest", "Arm", "Stomach", "Leg", "Leg", "Leg", "…
## $ officer_involved     <chr> "N", "N", "N", "N", "N", "N", "N", "N", "N", "N",…
## $ offender_injured     <chr> "N", "N", "N", "N", "N", "N", "N", "N", "N", "N",…
## $ offender_deceased    <chr> "N", "N", "N", "N", "N", "N", "N", "N", "N", "N",…
## $ location             <chr> "700 BLOCK Adams Ave", "1700 BLOCK S 5th St", "69…
## $ latino               <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ point_x              <dbl> -75.10397, -75.15405, -75.23670, -75.22840, -75.1…
## $ point_y              <dbl> 40.03060, 39.92726, 39.91765, 39.91845, 40.04942,…
## $ dist                 <dbl> 2, 3, 12, 12, 14, 14, 14, 14, 14, 14, 14, 14, 14,…
## $ inside               <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ outside              <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ fatal                <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ lat                  <dbl> 40.03060, 39.92726, 39.91765, 39.91845, 40.04942,…
## $ lng                  <dbl> -75.10397, -75.15405, -75.23670, -75.22840, -75.1…

Exploratory Data Analysis

Handling NAs

First, we look at NAs to see where we may have missing data that needs to be cleaned up. Looking at the summary of NAs a group of 153 NAs shows up across multiple columns including fatal. We will investigate those to clean them up since they will prevent us from accurately summing the fatalities.

shootings.raw %>% summarise_all(~ sum(is.na(.)))

Based on this, it looks like certain columns, including fatal, are left blank in any incident that involves an officer. In these cases the only fatality data provided is offender deceased. However, this would leave out any other fatalities that occurred during the incident since not everyone injured in a shooting incident is necessarily an offender.

Critical data for our analysis is missing from all of these officer involved shootings - time, race, age, wound, latino and, most importantly, fatal. There are 40 fatal incidents with an officer involved for which there is incomplete data.

shootings.raw %>% filter(if_any(fatal, is.na))

shootings.raw %>% filter(officer_involved == 'Y')

shootings.officer = shootings.raw %>% filter(officer_involved == 'Y')

shootings.officer %>% filter(offender_deceased == 'Y')

We can make an adjusted table, adding a total_fatality column that includes offender_deceased. It will not be usable for the questions of the profile of the most frequent shooting victims and their wound types, or the time of day of incidents, since that data is missing, but it can be used for the geolocation and date questions.

This introduces some potential bias into our data because it is possible that we are excluding non-offender fatalities. It is unknown to us if those are excluded or just did not occur. However, it is more complete than not including these 40 fatalities at all.

This gives us a column, total_fatality, with no NAs that can be summed to a total of 3494 fatalities between 2015 and 2025.

shootings.adjusted = shootings.raw

shootings.adjusted$offender_fatal <- ifelse(shootings.raw$offender_deceased == "Y", 1, 0)

shootings.adjusted <- shootings.adjusted %>% 
                      mutate(fatal = replace_na(fatal, 0))

shootings.adjusted$total_fatality <- shootings.adjusted$offender_fatal + shootings.adjusted$fatal

shootings.adjusted

sum(shootings.adjusted$total_fatality)

## [1] 3496

Examining Data by Year

Looking at fatality totals by year we can see that partial data is included for 2025. This will obviously be skewed on an annual basis since 2025 is not complete yet. We can see here that fatal shooting incidents rose from 2015 to 2021, peaked in 2021, then fell from 2021-2024.

fatal.byyear = aggregate(shootings.adjusted['total_fatality'], by=shootings.adjusted['year'], sum)

fatal.15to24 <- fatal.byyear %>%
  filter(year != 2025)

ggplot(data = fatal.15to24, aes(x = year, y = total_fatality)) +
  ggtitle("Total Fatal Shooting Victims by Year 2015-2024") +
  geom_col() +
  scale_x_continuous(breaks = 2015:2024)

## Fatal Shooting Victim Profile

For this analysis we will be only looking at victims of criminal shootings, not officer involved shootings so we will only use the fatal column, not total_fatality. ṁ ### Age of Victims First let’s get the mean age of fatal shooting victim.

fatal.only <- shootings.adjusted %>%
  filter(fatal == 1)

fatal.only <- fatal.only[!is.na(fatal.only$age), ]

mean(fatal.only$age)

## [1] 30.72565

Now, we will look at the overall distribution of shooting victim age.

ggplot(data = shootings.adjusted, aes(x = age, y = fatal)) +
  ggtitle("Age Distribution of Fatal Shooting Victims") +
  geom_col()

## Warning: Removed 274 rows containing missing values or values outside the scale range
## (`geom_col()`).

## Counting shootings to which body parts were most fatal

fatal_df <- shootings.raw[shootings.raw$fatal == 1, ]

injury_counts <- table(fatal_df$wound)
injury_counts

## 
##       abdomen       Abdomen           Arm          back          Back 
##             1            86            13             2           116 
##      buttocks      Buttocks         chest         Chest          Foot 
##             1             7             5           463             2 
##         Groin          Hand          head          Head          HEAD 
##             3             2            15           755             3 
##           Hip           Leg         multi         MULTI      Multiple 
##             5            40             1             5          1474 
## Multiple/Head Multiple/HEAD MULTIPLE/HEAD          Neck      Shoulder 
##           314             1             1            53            22 
##       Stomach         TORSO       Unknown 
##            53             1            11

sorted_injury_counts <- sort(injury_counts,decreasing = TRUE)
print(sorted_injury_counts)

## 
##      Multiple          Head         Chest Multiple/Head          Back 
##          1474           755           463           314           116 
##       Abdomen          Neck       Stomach           Leg      Shoulder 
##            86            53            53            40            22 
##          head           Arm       Unknown      Buttocks         chest 
##            15            13            11             7             5 
##           Hip         MULTI         Groin          HEAD          back 
##             5             5             3             3             2 
##          Foot          Hand       abdomen      buttocks         multi 
##             2             2             1             1             1 
## Multiple/HEAD MULTIPLE/HEAD         TORSO 
##             1             1             1

Plotting it

barplot(sorted_injury_counts,las = 2,
        col = "steelblue",     # Bar color
        main = "Injury Types in Fatal Shootings",
        ylab = "Number of Fatal Cases",
        cex.names = 0.7)  

which age groups are mostly affected by fatal shootings

fatal_df <- shootings.raw[shootings.raw$fatal==1,]
head(fatal_df)

Sort by age

fatal_df$age_group <- cut(fatal_df$age,
                          breaks = seq(0, 100, by = 10),
                          right = FALSE,
                          labels = c("0-9", "10-19", "20-29", "30-39", 
                                     "40-49", "50-59", "60-69", "70-79", 
                                     "80-89", "90-99"))

age_counts <- table(fatal_df$age_group)
age_counts <- sort(age_counts, decreasing = TRUE)
print(age_counts)

## 
## 20-29 30-39 10-19 40-49 50-59 60-69   0-9 70-79 80-89 90-99 
##  1413   907   425   410   180    52    14    12     5     1

Plotting fatal shooting counts with respect to age

barplot(age_counts,
        main = "Fatal Shootings by Age Group",
        xlab = "Age Group",
        ylab = "Number of Fatal Shootings",
        col = "darkred",
        cex.names = 0.7)

Average Age over Time

First, we will calculate the mean age of fatal shooting victim by year.

age_avg <- fatal.only %>%
  group_by(year) %>%
  summarize(average_value = mean(age))

age_avg

First, we will calculate the mean age of fatal shooting victim by year.

age_avg <- fatal.only %>%
  group_by(year) %>%
  summarize(average_value = mean(age))

age_avg

Then, let’s visualize the average age by year for 2015 - 2014. We are excluding 2025 because it is a partial year of data.

ggplot(data = age_avg, aes(x = year, y = average_value)) +
  ggtitle("Average Age of Fatal Shooting Victims by Year 2015-2024") +
  geom_col() +
  scale_x_continuous(breaks = 2015:2024)

AK - Geospatial Analysis

Clean and prepare the data

# View structure
glimpse(shootings.adjusted)

## Rows: 16,732
## Columns: 27
## $ the_geom             <chr> "0101000020E61000001049CE6BA7C652C0F0FF5490EA0344…
## $ the_geom_webmercator <chr> "0101000020110F0000790F72E295E45FC1A02823583D9452…
## $ objectid             <dbl> 634246, 634247, 634248, 634249, 634250, 634251, 6…
## $ year                 <dbl> 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2…
## $ dc_key               <dbl> 201502062409, 201503069342, 201512041905, 2015120…
## $ code                 <dbl> 3400, 400, 400, 400, 400, 400, 400, 400, 3400, 40…
## $ date_                <dttm> 2015-09-26, 2015-11-08, 2015-05-30, 2015-12-11, …
## $ time                 <time> 05:07:00, 12:50:00, 14:11:00, 18:00:00, 19:40:00…
## $ race                 <chr> "B", "A", "B", "B", "B", "B", "B", "B", "B", "B",…
## $ sex                  <chr> "M", "M", "M", "M", "F", "F", "F", "M", "F", "M",…
## $ age                  <dbl> 43, 36, 27, 23, 15, 21, 26, 33, 26, 43, 21, 35, 2…
## $ wound                <chr> "Chest", "Arm", "Stomach", "Leg", "Leg", "Leg", "…
## $ officer_involved     <chr> "N", "N", "N", "N", "N", "N", "N", "N", "N", "N",…
## $ offender_injured     <chr> "N", "N", "N", "N", "N", "N", "N", "N", "N", "N",…
## $ offender_deceased    <chr> "N", "N", "N", "N", "N", "N", "N", "N", "N", "N",…
## $ location             <chr> "700 BLOCK Adams Ave", "1700 BLOCK S 5th St", "69…
## $ latino               <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ point_x              <dbl> -75.10397, -75.15405, -75.23670, -75.22840, -75.1…
## $ point_y              <dbl> 40.03060, 39.92726, 39.91765, 39.91845, 40.04942,…
## $ dist                 <dbl> 2, 3, 12, 12, 14, 14, 14, 14, 14, 14, 14, 14, 14,…
## $ inside               <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ outside              <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ fatal                <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ lat                  <dbl> 40.03060, 39.92726, 39.91765, 39.91845, 40.04942,…
## $ lng                  <dbl> -75.10397, -75.15405, -75.23670, -75.22840, -75.1…
## $ offender_fatal       <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ total_fatality       <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…

summary(shootings.adjusted)

##    the_geom         the_geom_webmercator    objectid           year     
##  Length:16732       Length:16732         Min.   :634246   Min.   :2015  
##  Class :character   Class :character     1st Qu.:638429   1st Qu.:2018  
##  Mode  :character   Mode  :character     Median :642612   Median :2020  
##                                          Mean   :642612   Mean   :2020  
##                                          3rd Qu.:646794   3rd Qu.:2022  
##                                          Max.   :650977   Max.   :2025  
##                                                                         
##      dc_key               code            date_                       
##  Min.   :1.502e+03   Min.   : 100.0   Min.   :2015-01-01 00:00:00.00  
##  1st Qu.:2.018e+11   1st Qu.: 300.0   1st Qu.:2018-04-02 00:00:00.00  
##  Median :2.020e+11   Median : 400.0   Median :2020-09-25 00:00:00.00  
##  Mean   :2.001e+11   Mean   : 428.5   Mean   :2020-05-14 23:53:22.39  
##  3rd Qu.:2.022e+11   3rd Qu.: 400.0   3rd Qu.:2022-06-27 06:00:00.00  
##  Max.   :2.025e+11   Max.   :3700.0   Max.   :2025-04-28 00:00:00.00  
##                      NA's   :167                                      
##      time              race               sex                 age        
##  Length:16732      Length:16732       Length:16732       Min.   :  0.00  
##  Class1:hms        Class :character   Class :character   1st Qu.: 21.00  
##  Class2:difftime   Mode  :character   Mode  :character   Median : 27.00  
##  Mode  :numeric                                          Mean   : 29.34  
##                                                          3rd Qu.: 35.00  
##                                                          Max.   :117.00  
##                                                          NA's   :274     
##     wound           officer_involved   offender_injured   offender_deceased 
##  Length:16732       Length:16732       Length:16732       Length:16732      
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##    location             latino         point_x          point_y     
##  Length:16732       Min.   :0.000   Min.   :-75.27   Min.   :39.88  
##  Class :character   1st Qu.:0.000   1st Qu.:-75.19   1st Qu.:39.97  
##  Mode  :character   Median :0.000   Median :-75.16   Median :39.99  
##                     Mean   :0.123   Mean   :-75.16   Mean   :39.99  
##                     3rd Qu.:0.000   3rd Qu.:-75.13   3rd Qu.:40.02  
##                     Max.   :1.000   Max.   :-74.96   Max.   :40.13  
##                     NA's   :153     NA's   :33       NA's   :33     
##       dist           inside           outside           fatal       
##  Min.   : 1.00   Min.   :0.00000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:15.00   1st Qu.:0.00000   1st Qu.:1.0000   1st Qu.:0.0000  
##  Median :22.00   Median :0.00000   Median :1.0000   Median :0.0000  
##  Mean   :20.91   Mean   :0.06333   Mean   :0.9367   Mean   :0.2066  
##  3rd Qu.:25.00   3rd Qu.:0.00000   3rd Qu.:1.0000   3rd Qu.:0.0000  
##  Max.   :77.00   Max.   :1.00000   Max.   :1.0000   Max.   :1.0000  
##  NA's   :2       NA's   :153       NA's   :153                      
##       lat             lng         offender_fatal     total_fatality  
##  Min.   :39.88   Min.   :-75.27   Min.   :0.000000   Min.   :0.0000  
##  1st Qu.:39.97   1st Qu.:-75.19   1st Qu.:0.000000   1st Qu.:0.0000  
##  Median :39.99   Median :-75.16   Median :0.000000   Median :0.0000  
##  Mean   :39.99   Mean   :-75.16   Mean   :0.002391   Mean   :0.2089  
##  3rd Qu.:40.02   3rd Qu.:-75.13   3rd Qu.:0.000000   3rd Qu.:0.0000  
##  Max.   :40.13   Max.   :-74.96   Max.   :1.000000   Max.   :1.0000  
##  NA's   :28      NA's   :28

shootings.spatial  <- shootings.adjusted %>%
  filter(!is.na(lat) & !is.na(lng)) %>%
  rename(latitude = lat, longitude = lng)

Convert to spatial data

shootings_sf <- st_as_sf(shootings.spatial , coords = c("longitude", "latitude"), crs = 4326)

Plotting Static Map

Visualizes Philadelphia shooting incidents using geospatial coordinates. The scatter plot reveals the geographic distribution, with red dots highlighting significant clustering of events in specific areas of the city.

ggplot(data = shootings_sf) +
  geom_sf(alpha = 0.3, color = "red") +
  theme_minimal() +
  labs(title = "Philadelphia Shooting Incidents",
       subtitle = "Based on geospatial coordinates",
       caption = "Source: philly_shootings.csv")

Interactive Map Using Leaflet

This interactive map displays shooting incidents in Philadelphia from 2015 to 2025. Each dot represents a single event and is color-coded by year. With all years selected, the visualization reveals the cumulative geographic distribution and extremely dense concentration of shootings in specific neighborhoods across the city. Futher incidents can be filtered based on year. Clicking a data point reveals a pop-up with specific details for that shooting, including its exact date, location, and fatality/injury information.

if (!"total_fatality" %in% names(shootings_sf) || !"year" %in% names(shootings_sf) || !"date_" %in% names(shootings_sf) || !"location" %in% names(shootings_sf) ) {
  stop("shootings_sf is missing one or more required columns: year, date_, location, total_fatality.")
}

if (!"unique_id" %in% names(shootings_sf)) {
  shootings_sf <- shootings_sf %>%
    dplyr::mutate(unique_id = paste0("marker-", dplyr::row_number(), "-", year))
}

# Create unique list of years
years <- if ("year" %in% names(shootings_sf)) sort(unique(shootings_sf$year)) else character(0)

# Create color palette
if (length(years) > 0) {
  if (length(years) <= 8) {
    year_palette_colors <- RColorBrewer::brewer.pal(max(3, length(years)), "Dark2")
    if (length(years) < 3) {
      year_palette_colors <- year_palette_colors[1:length(years)]
    }
  } else {
    year_palette_colors <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(min(length(years), 9), "Dark2"))(length(years))
  }
  year_colors <- leaflet::colorFactor(palette = year_palette_colors, domain = years)
} else {
  year_colors <- NULL
  warning("Year column not found or no unique years in shootings_sf. Map may not display year layers correctly.")
}

# Create base map centered on Philadelphia
center_lon <- -75.1652
center_lat <- 39.9526

m <- leaflet::leaflet(shootings_sf) %>%
  leaflet::setView(lng = center_lon, lat = center_lat, zoom = 11) %>%
  leaflet::addProviderTiles(leaflet::providers$CartoDB.Positron) %>%
  htmlwidgets::onRender("function(el, x) { el.classList.add('safegraph-map'); }")

# Add circle markers by year
if (length(years) > 0 && !is.null(year_colors)) {
  for (yr in years) {
    year_data <- shootings_sf %>% dplyr::filter(year == yr)
    if (nrow(year_data) > 0) {
      m <- m %>% leaflet::addCircleMarkers(
        data = year_data,
        radius = 5,
        color = ~year_colors(year),
        fillOpacity = 0.7,
        stroke = FALSE,
        group = as.character(yr),
        popup = ~paste(
          "<b>Year:</b>", htmltools::htmlEscape(year), "<br>",
          "<b>Date:</b>", htmltools::htmlEscape(date_), "<br>",
          "<b>Location:</b>", htmltools::htmlEscape(location), "<br>",
          "<b>Fatality:</b>", ifelse(total_fatality > 0, "Yes", "No"),
          "<b>Wound:</b>", htmltools::htmlEscape(wound)
        ),
        layerId = ~unique_id
      )
    }
  }
}

# Add layer control
if (length(years) > 0) {
  m <- m %>% leaflet::addLayersControl(
    overlayGroups = as.character(years),
    options = leaflet::layersControlOptions(collapsed = FALSE)
  )
}

## For button
year_groups_json_str <- jsonlite::toJSON(as.character(years)) 

select_toggle_js <- paste0(
  "function(el, x, data) {\n",
  "  setTimeout(function() {\n",
  "    var yearGroups = ", year_groups_json_str, ";\n",
  "    if (!yearGroups || yearGroups.length === 0) { console.warn('No year groups for toggle button.'); return; }\n",
  "    var layersControlContainer = el.querySelector('.leaflet-control-layers-overlays');\n",
  "    if (!layersControlContainer) { console.warn('Layers control container not found for toggle button.'); return; }\n",
  "\n",
  "    var toggleButton = L.DomUtil.create('button', 'leaflet-control-layers-selector safegraph-button');\n",
  "    toggleButton.style.marginBottom = '10px';\n",
  "\n",
  "    function getCurrentAllCheckedState() {\n",
  "      if (yearGroups.length === 0) return false;\n",
  "      for (var i = 0; i < yearGroups.length; i++) {\n",
  "        var groupStr = String(yearGroups[i]);\n",
  "        var checkbox = layersControlContainer.querySelector('input[type=\"checkbox\"][value=\"' + groupStr + '\"]');\n",
  "        if (!checkbox || !checkbox.checked) return false;\n",
  "      }\n",
  "      return true;\n",
  "    }\n",
  "\n",
  "    function updateButtonText() {\n",
  "      if (getCurrentAllCheckedState()) {\n",
  "        toggleButton.innerHTML = 'Deselect All Years';\n",
  "      } else {\n",
  "        toggleButton.innerHTML = 'Select All Years';\n",
  "      }\n",
  "    }\n",
  "\n",
  "    toggleButton.onclick = function() {\n",
  "      var currentlyAllChecked = getCurrentAllCheckedState();\n",
  "      var targetState = !currentlyAllChecked;\n",
  "      yearGroups.forEach(function(group) {\n",
  "        var groupStr = String(group);\n",
  "        var checkbox = layersControlContainer.querySelector('input[type=\"checkbox\"][value=\"' + groupStr + '\"]');\n",
  "        if (checkbox) {\n",
  "          if (checkbox.checked !== targetState) {\n",
  "            checkbox.click();\n",
  "          }\n",
  "        }\n",
  "      });\n",
  "      updateButtonText();\n",
  "    };\n",
  "\n",
  "    layersControlContainer.parentNode.insertBefore(toggleButton, layersControlContainer);\n",
  "    updateButtonText(); /* Set initial button text */\n",
  "  }, 250);\n",
  "}\n"
)

m <- m %>% htmlwidgets::onRender(select_toggle_js)

# Add legend
if (length(years) > 0 && !is.null(year_colors)) {
  m <- m %>% leaflet::addLegend(
    position = "bottomleft",
    pal = year_colors,
    values = shootings_sf$year,
    title = "Shooting Year",
    opacity = 1
  )
}

# CSS Styling using paste0
css_styling_js <- paste0(
  "function(el, x) {\n",
  "  var css = `\n", 
  "    .safegraph-map { border-radius: 8px; box-shadow: 0 2px 8px rgba(0,0,0,0.15); overflow: hidden; }\n",
  "    .safegraph-button { background-color: #f0f0f0; border: 1px solid #ccc; border-radius: 4px; padding: 4px 8px; font-size: 12px; cursor: pointer; transition: background-color 0.3s ease; font-family: 'Arial', sans-serif; color: #333; display: block; width: 100%; text-align: center; white-space: nowrap; box-sizing: border-box; }\n",
  "    .safegraph-button:hover { background-color: #e0e0e0; }\n",
  "    .leaflet-control-layers-list input[type='checkbox'] { margin-right: 4px; }\n",
  "    .leaflet-popup-content { font-size: 12px; line-height: 1.4; }\n",
  "    .leaflet-popup-content b { font-weight: 600; color: #2c3e50; }\n",
  "    .leaflet-control-layers { background-color: rgba(255,255,255,0.9); border-radius: 8px; box-shadow: 0 1px 4px rgba(0,0,0,0.1); padding: 6px 10px; font-size: 12px; color: #34495e; }\n",
  "    .leaflet-control-layers-title { font-weight: bold; margin-bottom: 6px; color: #2c3e50; }\n",
  "  `;\n", 
  "  var style = document.createElement('style');\n",
  "  style.type = 'text/css';\n",
  "  style.appendChild(document.createTextNode(css));\n",
  "  document.head.appendChild(style);\n",
  "  if (!el.classList.contains('safegraph-map')) { el.classList.add('safegraph-map'); }\n",
  "}\n"
)

m <- m %>% htmlwidgets::onRender(css_styling_js)

# Display map
m

Heatmap Visualization

Over the last decade, 16,704 total incidents have occurred across the city. To better visualize these widespread events, the data was organized into geographical clusters and displayed as a heatmaps for various types of incidents, revealing areas with the highest concentrations.

center_lon <- -75.1652
center_lat <- 39.9526

# Prepare base data for the map
shootings_sf_for_map <- shootings_sf %>%
  dplyr::mutate(
    incident_color = ifelse(total_fatality > 0, "red", "orange"),
    fatality_label = ifelse(total_fatality > 0, "Fatal", "Non-Fatal")
  )

# Create filtered datasets for new heatmap layers
officer_involved_fatal_sf <- shootings_sf_for_map %>%
  dplyr::filter(officer_involved == 'Y' & total_fatality > 0)

officer_involved_non_fatal_sf <- shootings_sf_for_map %>%
  dplyr::filter(officer_involved == 'Y' & total_fatality == 0)

offender_deceased_sf <- shootings_sf_for_map %>%
  dplyr::filter(offender_deceased == 'Y') 

all_officer_involved_sf <- shootings_sf_for_map %>%
  dplyr::filter(officer_involved == 'Y')
  
# Create the Leaflet map
map_enhanced_heatmaps <- leaflet::leaflet(data = shootings_sf_for_map) %>% # Base data
  leaflet::setView(lng = center_lon, lat = center_lat, zoom = 11) %>%
  leaflet::addProviderTiles(leaflet::providers$CartoDB.Positron, group = "Base Map") %>%
  
  # Clustered Markers (as before)
  leaflet::addAwesomeMarkers(
    lng = ~sf::st_coordinates(geometry)[,1], 
    lat = ~sf::st_coordinates(geometry)[,2],
    icon = ~leaflet::awesomeIcons(
      icon = 'info-circle', 
      library = 'glyphicon', 
      markerColor = incident_color, 
      iconColor = 'white'
    ),
    popup = ~paste(
      "<b>Year:</b>", htmltools::htmlEscape(year), "<br>",
      "<b>Date:</b>", htmltools::htmlEscape(date_), "<br>",
      "<b>Location:</b>", htmltools::htmlEscape(location), "<br>",
      "<b>Fatality:</b>", ifelse(total_fatality > 0, "Yes", "No"), "<br>",
      "<b>Wound:</b>", if ("wound" %in% names(.)) htmltools::htmlEscape(wound) else "N/A"
    ),
    clusterOptions = leaflet::markerClusterOptions(),
    group = "Shooting Incidents (Clusters)"
  ) %>%
  
  # Heatmap
  leaflet.extras::addHeatmap(
    lng = ~sf::st_coordinates(geometry)[,1], 
    lat = ~sf::st_coordinates(geometry)[,2],
    blur = 20, max = 0.05, radius = 15,
    group = "Heatmap (All Incidents)"
  ) %>%

  # Heatmap for All Fatal Incidents
  leaflet.extras::addHeatmap(
    data = dplyr::filter(shootings_sf_for_map, total_fatality > 0), 
    lng = ~sf::st_coordinates(geometry)[,1], 
    lat = ~sf::st_coordinates(geometry)[,2],
    blur = 20, max = 0.05, radius = 15,
    group = "Heatmap (All Fatal Incidents)" 
  ) %>%
  
  # Heatmap for Officer Involved & Fatal Incidents
  leaflet.extras::addHeatmap(
    data = officer_involved_fatal_sf,
    lng = ~sf::st_coordinates(geometry)[,1], 
    lat = ~sf::st_coordinates(geometry)[,2],
    blur = 20, max = 0.05, radius = 15,
    group = "Heatmap (Officer Involved & Fatal)"
  ) %>%
  
  # Heatmap for Officer Involved & Non-Fatal Incidents
  leaflet.extras::addHeatmap(
    data = officer_involved_non_fatal_sf,
    lng = ~sf::st_coordinates(geometry)[,1], 
    lat = ~sf::st_coordinates(geometry)[,2],
    blur = 20, max = 0.05, radius = 15, 
    group = "Heatmap (Officer Involved & Non-Fatal)"
  ) %>%
  
  # Heatmap for Offender Deceased Incidents
  leaflet.extras::addHeatmap(
    data = offender_deceased_sf,
    lng = ~sf::st_coordinates(geometry)[,1], 
    lat = ~sf::st_coordinates(geometry)[,2],
    blur = 20, max = 0.05, radius = 15,
    group = "Heatmap (Offender Deceased)"
  ) %>%

  # Heatmap for All Officer Involved Incidents
  leaflet.extras::addHeatmap(
    data = all_officer_involved_sf,
    lng = ~sf::st_coordinates(geometry)[,1], 
    lat = ~sf::st_coordinates(geometry)[,2],
    blur = 20, max = 0.05, radius = 15,
    group = "Heatmap (All Officer Involved)"
  ) %>%
  
  # Legend for Clustered Markers
  leaflet::addLegend(
    position = "bottomright", 
    colors = c("red", "orange"),
    labels = c("Fatal Incident", "Non-Fatal Incident"),
    title = "Incident Type (Markers)",
    opacity = 1,
    group = "Shooting Incidents (Clusters)" 
  ) %>%
  
  # Layers Control to switch between layers
  leaflet::addLayersControl(
    baseGroups = "Base Map",
    overlayGroups = c(
      "Shooting Incidents (Clusters)", 
      "Heatmap (All Incidents)", 
      "Heatmap (All Fatal Incidents)",
      "Heatmap (All Officer Involved)",
      "Heatmap (Officer Involved & Fatal)",
      "Heatmap (Officer Involved & Non-Fatal)", 
      "Heatmap (Offender Deceased)"
    ),
    options = leaflet::layersControlOptions(collapsed = FALSE) 
  ) %>%
  # Hide all heatmap layers by default so clusters are visible first
  leaflet::hideGroup(c(
      "Heatmap (All Incidents)", 
      "Heatmap (All Fatal Incidents)",
      "Heatmap (All Officer Involved)",
      "Heatmap (Officer Involved & Fatal)",
      "Heatmap (Officer Involved & Non-Fatal)",
      "Heatmap (Offender Deceased)"
  ))

# Display the map
map_enhanced_heatmaps