Below is a choropleth map showing the number of aggravated assaults in the Nashville-Metropolitan area in summer 2026, based on data from the Metro Nashville Police Department Incidents API. It initially included each area’s ZCTA and the number of assaults in that area. The map now includes the population and assault rate per 1,000 residents within each area.
I used Microsoft Copilot to generate some suggestions on improving the initial script for the map. Some of the suggestions included classification breaks instead of a continuous color scale, restricting the map so that it only included ZIP codes in Davidson County, and mapping the assault rates of each area instead of assault counts. While mapping the assault rates was an excellent suggestion on paper, some of the low-population areas with higher assault counts looked very skewed on the map’s legend, particularly ZCTA 37213, with a population of 40 and 9 reported aggravated assaults (that map showed that area as bright red whereas other areas were mostly white). Instead of having the map represent the assault rate for each area, I had Copilot help me write R code to include the population and assault rate while having the map show the number of assaults.
These statistics can give more insight and context for each area’s assault count. For instance, ZCTA 37210’s assault rate is rather high considering its assault count relative to its population. The previous map would have shown just the assault count, and one would have to research the population of that area to better understand why there are so many assaults in that area. By including both the population and assault rate, we can have a much better explanation as to why they have high assault counts.
Below is the code used to make the map. Copilot first told me to add the tidycensus package, which essentially acts as a bridge between the crime data and Census demographic data. Next, it added three additional steps after Step 13. In those steps, the code downloaded ZCTA population data, joined that data to the ZCTA layer, and then calculated the assault rates. Finally, it updated the popup content to include the population data and assault rates.
############################################################
# 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)
library(tidycensus)
# 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)
############################################################
# 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: 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 13: Join Incident Counts to ZCTAs
#
# Attach the incident totals to the corresponding
# ZCTA polygons.
############################################################
Nashville_ZCTAs <- Nashville_ZCTAs |>
left_join(
ZCTACounts,
by = "ZCTA5CE20"
)
############################################################
# Step 13A: Download ZCTA Population Data
#
# Census table B01003 contains total population.
############################################################
ZCTA_Population <- get_acs(
geography = "zcta",
variables = "B01003_001",
year = 2020
) |>
select(
ZCTA5CE20 = GEOID,
Population = estimate
)
############################################################
# Step 13B: Join Population Data
############################################################
Nashville_ZCTAs <- Nashville_ZCTAs |>
left_join(
ZCTA_Population,
by = "ZCTA5CE20"
)
############################################################
# Step 13C: Calculate Assault Rates
#
# Rate = assaults per 1,000 residents
############################################################
Nashville_ZCTAs <- Nashville_ZCTAs |>
mutate(
AssaultRate = (Incidents / Population) * 1000
)
############################################################
# Step 14: 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$Incidents
)
############################################################
# Step 15: Create Popup Content
############################################################
ZCTAPopup <- ~paste0(
"<strong>ZCTA:</strong> ", ZCTA5CE20,
"<br><strong>Population:</strong> ",
format(Population, big.mark = ","),
"<br><strong>Aggravated Assaults:</strong> ",
Incidents,
"<br><strong>Assault Rate:</strong> ",
round(AssaultRate, 1),
" per 1,000 residents"
)
############################################################
# Step 16: Create the Choropleth Map
#
# Polygon Color = Number of Incidents
#
# This map highlights which ZIP code areas have
# the greatest concentrations of aggravated assaults.
############################################################
ChoroplethMap <- leaflet(Nashville_ZCTAs) |>
addProviderTiles("Esri.WorldStreetMap") |>
addPolygons(
fillColor = ~pal(Incidents),
fillOpacity = 0.7,
color = "black",
weight = 1,
popup = ZCTAPopup
) |>
addLegend(
position = "bottomright",
pal = pal,
values = ~Incidents,
title = "Aggravated Assaults"
)
ChoroplethMap