How to plot ZIP codes on a map in R with zipcodeR

Mapping ZIP code data in R is a two-layer problem: you need the data (your values, keyed by ZIP) and the geometry (polygons describing where each ZIP sits on the map). zipcodeR is great for the data layer — it ships with demographic, geographic, and centroid attributes for every US ZIP — but it is not a geometry source. For the polygons you want tigris, which pulls ZIP Code Tabulation Areas (ZCTAs) directly from the U.S. Census Bureau.

This post walks through both a static choropleth (ggplot2) and an interactive map (leaflet).

Setup

install.packages(c("zipcodeR", "tigris", "sf", "ggplot2", "dplyr"))

library(zipcodeR)
library(tigris)
library(sf)
library(ggplot2)
library(dplyr)

options(tigris_use_cache = TRUE)  # cache the Census downloads locally

For the example I will use median household income from the zipcodeR built-in database, mapped over New Jersey ZCTAs. Substitute your own dataset for the value column.

1. Pull the polygon geometry

tigris::zctas() downloads the ZCTA polygons. The year argument controls which Census vintage you get; the state argument is the simplest way to filter:

nj_zctas <- zctas(year = 2020, state = "NJ", cb = TRUE)

cb = TRUE requests the cartographic-boundary version, which is generalized and renders much faster than the full TIGER/Line geometries. For analytical work where polygon area matters (e.g., density calculations) use cb = FALSE for the full-precision boundaries.

The result is an sf data frame with one row per ZCTA. The 5-digit ZIP code is in the ZCTA5CE20 column for the 2020 vintage (older vintages use ZCTA5CE10).

2. Get your data layer

For the example, pull median household income from zipcodeR’s zip_code_db:

nj_income <- zip_code_db %>%
  filter(state == "NJ") %>%
  select(zipcode, median_household_income, population)

For your own data this is wherever your value column lives — an analytics extract, a CSV, a database query. The only requirement is that one of the columns is the 5-digit ZIP, character type with leading zeros preserved.

3. Join the data onto the polygons

nj_map_data <- nj_zctas %>%
  left_join(nj_income, by = c("ZCTA5CE20" = "zipcode"))

A few rows may end up with NA values in your data column if your dataset doesn’t cover every ZCTA. That’s fine — ggplot will render those polygons in the missing-value color (gray by default).

4. Render the static choropleth

ggplot(nj_map_data) +
  geom_sf(aes(fill = median_household_income), color = "white", size = 0.05) +
  scale_fill_viridis_c(
    option = "magma",
    trans = "log10",
    labels = scales::dollar_format(),
    name = "Median household\nincome"
  ) +
  labs(
    title = "Median household income by ZIP code",
    subtitle = "New Jersey, 2020",
    caption = "Source: zipcodeR (US Census ACS). Geometry: tigris (US Census ZCTAs)."
  ) +
  theme_void() +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    legend.position = "right"
  )

A few render decisions worth flagging:

• Color scale. viridis is the right default — it is perceptually uniform and color-blind friendly. The magma option works well for income/density data because the dark end reads as “low” intuitively.
• Log transform. Income, population, and density data span orders of magnitude; a log transform makes the choropleth legible. Otherwise a few extreme values dominate the color scale.
• theme_void(). Strips the axis grid for a clean cartographic look. For thematic maps you almost never want the grid.

5. Interactive map with leaflet

For a web-deliverable interactive map, swap ggplot2 for leaflet:

library(leaflet)

# Project ZCTAs to WGS84 for leaflet
nj_map_data_wgs <- st_transform(nj_map_data, 4326)

pal <- colorNumeric(
  palette = "magma",
  domain = nj_map_data_wgs$median_household_income,
  na.color = "#cccccc"
)

leaflet(nj_map_data_wgs) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    fillColor = ~pal(median_household_income),
    fillOpacity = 0.75,
    color = "white",
    weight = 0.5,
    label = ~paste0(ZCTA5CE20, ": $", scales::comma(median_household_income))
  ) %>%
  addLegend(
    pal = pal,
    values = ~median_household_income,
    title = "Median income"
  )

The result is a pan-and-zoom web map with hover labels. Drop it into an R Markdown / Quarto document or a Shiny app and it just works.

ZCTA vs USPS ZIP — important caveat

ZCTAs (Census ZIP Code Tabulation Areas) are areal polygons. USPS ZIPs are operational route designations. They are usually equivalent, but not always:

• USPS adds and changes ZIPs continuously; ZCTAs are updated only with each decennial Census.
• USPS PO Box-only ZIPs and some unique ZIPs (large corporate campuses, military addresses) have no corresponding ZCTA polygon.
• In a few cases, multiple USPS ZIPs are aggregated into one ZCTA for cartographic purposes.

For statistical work — almost all academic, public health, and policy mapping — ZCTAs are correct and you should use them. For operational work (mail delivery, logistics) USPS ZIPs are correct and ZCTAs are an approximation.

Related work

• zipcodeR project page — full feature list, FAQ, and CRAN/GitHub links.
• How to calculate distance between ZIP codes in R — proximity calculations and radius searches.
• How to assign ZIP codes to geographic regions in R — county/state/MSA rollups and custom region crosswalks.
• New Jersey Population Density Map (3D) — tract-level rendering of NJ population density using rayshader, built with the same toolchain.

zipcodeR is free and open source under the GPL-3 license. tigris is published by Kyle Walker (Texas Christian University) under the MIT license. Both are CRAN-stable and part of the standard R geospatial toolkit.