The original choropleth map showed aggravated assault incidents in Nashville’s ZIP Code Tabulation Areas (ZCTAs). To make comparisons more meaningful, I changed the map to display assault rates per 1,000 residents instead of raw incident counts. This shift helps viewers identify areas with higher assault rates while considering population differences, emphasizing the map’s focus on analytical insights.
To create this modification, I used the tidycensus
package’s get_acs() function to retrieve population
estimates from the U.S. Census Bureau’s American Community Survey (ACS).
I requested variable B01003_001, representing the total
population by ZCTA. The population data were then joined to the crime
counts using a common ZCTA identifier, enabling the calculation of
assault rates:
(Incidents / Population) * 1000
The choropleth map was subsequently adjusted so that the colors of the polygons, the map legend, and the pop-up information were based on the assault rate rather than the total number of incidents. This change helps viewers interpret the data more accurately by highlighting Nashville-area ZCTAs with relatively high aggravated assault rates instead of merely identifying areas with the largest number of reported assaults.
During verification, I reviewed population and assault rates for each ZCTA and identified ZCTA 37213 as an outlier. Although its population estimate of 40 residents was confirmed, its assault rate was abnormally high, distorting the color scale. Removing ZCTA 37213 improved visual differentiation among the remaining ZCTAs and more clearly revealed geographic patterns in assault rates.
Artificial intelligence (AI) was used for brainstorming and coding assistance throughout this project. AI suggested several potential map enhancements, including displaying assault rates per capita, adding additional information to map popups, and creating a toggle-able layer for individual crime incidents. The most valuable suggestion was incorporating ACS population data to calculate aggravated assault rates per 1,000 residents, which significantly improved the map’s analytical value and enabled fairer comparisons among Nashville-area ZCTAs. One suggestion I chose not to implement was adding a second layer to display individual crime incidents. While this modification could have been helpful, I felt that the rate-based choropleth provided greater analytical value while remaining suitable for an introductory mapping assignment.
############################################################
# Lesson 2: Mapping Crime by Area with a Choropleth Map
#
# This script creates a choropleth map showing how
# aggravated assaults are distributed across ZIP Code
# Tabulation Areas (ZCTAs).
#
# Instead of mapping individual incidents, this lesson:
#
# 1. Assigns each incident to a ZCTA.
# 2. Counts incidents within each ZCTA.
# 3. Colors ZCTAs based on incident totals.
#
# Darker shades represent more incidents.
############################################################
############################################################
# Step 1: Install Required Packages
#
# Check whether the required packages are installed.
# Install any that are missing.
############################################################
required_packages <- c(
"jsonlite",
"tidyverse",
"sf",
"leaflet",
"tigris",
"tidycensus"
)
installed_packages <- rownames(installed.packages())
for (pkg in required_packages) {
if (!pkg %in% installed_packages) {
install.packages(pkg)
}
}
############################################################
# Step 2: Load Required Packages
#
# jsonlite = download data from the API
# tidyverse = clean and summarize data
# sf = perform spatial analysis
# leaflet = create interactive maps
# tigris = download Census boundaries
############################################################
library(jsonlite)
library(tidyverse)
library(sf)
library(leaflet)
# Store Census downloads locally so they do not need
# to be downloaded every time the script is run.
options(tigris_use_cache = TRUE)
library(tigris)
library(tidycensus)
############################################################
# Step 3: Download Aggravated Assault Data
#
# Retrieve Summer 2026 aggravated assault incidents
# (NIBRS code 13A) from Nashville's crime API.
############################################################
base_url <- paste0(
"https://services2.arcgis.com/HdTo6HJqh92wn4D8/",
"arcgis/rest/services/",
"Metro_Nashville_Police_Department_Incidents_view/",
"FeatureServer/0/query"
)
query <- paste(
"Incident_Occurred >= DATE '2026-06-01'",
"AND Incident_Occurred < DATE '2026-09-01'",
"AND Offense_NIBRS = '13A'"
)
crime_url <- paste0(
base_url,
"?where=",
URLencode(query, reserved = TRUE),
"&outFields=*",
"&f=json"
)
CrimeData <- fromJSON(
crime_url
)$features$attributes
############################################################
# Step 4: Prepare the Data
#
# Convert dates into a readable format and remove
# records that are missing coordinates.
############################################################
CrimeData <- CrimeData |>
mutate(
Incident_Occurred = as.POSIXct(
Incident_Occurred / 1000,
origin = "1970-01-01",
tz = "America/Chicago"
)
) |>
filter(
!is.na(Longitude),
!is.na(Latitude)
)
cat(
"Number of aggravated assault incidents:",
nrow(CrimeData),
"\n"
)
############################################################
# Step 5: Convert Incidents into a Spatial Layer
#
# Unlike Lesson 1, this lesson requires spatial
# analysis. The sf package converts latitude and
# longitude values into geographic points that can
# participate in GIS operations.
############################################################
CrimeData_sf <- st_as_sf(
CrimeData,
coords = c("Longitude", "Latitude"),
crs = 4326,
remove = FALSE
)
############################################################
# Step 6: Download ZCTA Boundaries
#
# Download ZIP Code Tabulation Area (ZCTA)
# boundaries from the U.S. Census Bureau.
############################################################
TN_ZCTAs <- zctas(
cb = TRUE,
year = 2020
)
############################################################
# Step 7: Keep Tennessee ZCTAs
#
# Tennessee ZIP codes generally begin with "37".
# Keep only ZCTAs that are likely to be located
# within Tennessee.
############################################################
TN_ZCTAs <- TN_ZCTAs |>
filter(
startsWith(ZCTA5CE20, "37")
)
############################################################
# Step 8: Match Coordinate Systems
#
# Spatial layers must use the same coordinate
# reference system (CRS) before they can be
# combined through a spatial join.
############################################################
TN_ZCTAs <- st_transform(
TN_ZCTAs,
st_crs(CrimeData_sf)
)
############################################################
# Step 9: Assign Each Incident to a ZCTA
#
# This spatial join determines which ZCTA contains
# each crime incident.
#
# Conceptually:
#
# Point
# ↓
# Falls inside
# ↓
# Polygon
############################################################
CrimeWithZCTA <- st_join(
CrimeData_sf,
TN_ZCTAs |>
select(ZCTA5CE20)
)
############################################################
# Step 10: Verify the Spatial Join
#
# Confirm that incidents were successfully assigned
# to ZCTAs.
############################################################
CrimeWithZCTA |>
st_drop_geometry() |>
count(is.na(ZCTA5CE20))
############################################################
# Step 11: Count Incidents by ZCTA
#
# Aggregate individual incidents into area-level
# totals.
#
# This transforms:
#
# One row = One incident
#
# into:
#
# One row = One ZCTA
############################################################
ZCTACounts <- CrimeWithZCTA |>
st_drop_geometry() |>
count(
ZCTA5CE20,
name = "Incidents"
) |>
arrange(desc(Incidents))
ZCTACounts
############################################################
# Step 12: Download ZCTA Population Data
#
# Retrieve total population estimates for each ZCTA
# from the American Community Survey (ACS).
############################################################
ZCTAPopulation <- get_acs(
geography = "zcta",
variables = "B01003_001",
year = 2020
) |>
select(
ZCTA5CE20 = GEOID,
Population = estimate
)
head(ZCTAPopulation)
############################################################
# Step 13: Calculate Assault Rates
#
# Join population data to incident counts and
# calculate aggravated assault rates per
# 1,000 residents.
############################################################
ZCTACounts <- ZCTACounts |>
left_join(
ZCTAPopulation,
by = "ZCTA5CE20"
) |>
mutate(
AssaultRate = (Incidents / Population) * 1000
)
head(ZCTACounts)
############################################################
# Step 14: Keep Only ZCTAs Found in the Data
#
# Retain only the ZCTAs containing one or more
# aggravated assault incidents.
############################################################
Nashville_ZCTAs <- TN_ZCTAs |>
filter(
!is.na(ZCTA5CE20),
ZCTA5CE20 %in% unique(ZCTACounts$ZCTA5CE20)
)
############################################################
# Step 15: Join Incident Counts to ZCTAs
#
# Attach the incident totals to the corresponding
# ZCTA polygons.
############################################################
Nashville_ZCTAs <- Nashville_ZCTAs |>
left_join(
ZCTACounts,
by = "ZCTA5CE20"
)
# Deleting the weird ZIP code
Nashville_ZCTAs <- Nashville_ZCTAs %>%
filter(NAME20 != 37213)
############################################################
# Step 16: Create a Choropleth Color Palette
#
# Choropleth maps use color to represent numerical
# values.
#
# Lighter colors = fewer incidents
# Darker colors = more incidents
############################################################
pal <- colorNumeric(
palette = "Reds",
domain = Nashville_ZCTAs$AssaultRate
)
############################################################
# Step 17: Create Popup Content
#
# Display information about each ZCTA when users
# click a polygon.
############################################################
ZCTAPopup <- ~paste0(
"<strong>ZCTA:</strong> ",
ZCTA5CE20,
"<br><strong>Aggravated Assaults:</strong> ",
Incidents,
"<br><strong>Population:</strong> ",
format(Population, big.mark = ","),
"<br><strong>Assault Rate:</strong> ",
round(AssaultRate, 2),
" per 1,000 residents"
)
############################################################
# Step 18: Create the Choropleth Map
#
# Polygon Color = Number of Incidents
#
# This map highlights which ZIP code areas have
# the greatest concentrations of aggravated assaults.
############################################################
# Nashville_ZCTAs <- Nashville_ZCTAs %>%
# filter(NAME20 != 37213)
ChoroplethMap <- leaflet(Nashville_ZCTAs) |>
addProviderTiles("Esri.WorldStreetMap") |>
addPolygons(
fillColor = ~pal(AssaultRate),
fillOpacity = 0.7,
color = "black",
weight = 1,
popup = ZCTAPopup
) |>
addLegend(
position = "bottomright",
pal = pal,
values = ~AssaultRate,
title = "Assaults per 1,000 Residents"
)
ChoroplethMap