Spatial Data & Mapping¶
PyPUMS can attach Census cartographic boundary shapefiles to any Census
query, returning a GeoDataFrame that is ready for mapping, spatial joins, and
geospatial analysis. Shapefile downloads are handled by
pygris, with automatic local caching so
files are only downloaded once.
Requires the spatial extras
Geometry features depend on geopandas and pygris (plus their dependencies). Install them with:
The geometry=True flag¶
The fastest way to get spatial data is to pass geometry=True to any of the
main data retrieval functions. PyPUMS will automatically download the
corresponding cartographic boundary shapefile (via pygris), merge it with the
tabular data on the GEOID column, and return a GeoDataFrame. Downloaded
shapefiles are cached locally so subsequent calls are fast.
Tip
Notice the geometry column — that is what makes this a
GeoDataFrame instead of a plain DataFrame. Each row carries its
polygon boundary, ready for mapping or spatial analysis.
GEOID geometry … variable value 0 09 MULTIPOLYGON (((-72.22593 41.29384, -72.22523 … … P1_001N 3605944 1 10 MULTIPOLYGON (((-75.56752 39.5102, -75.56477 3… … P1_001N 989948 2 11 POLYGON ((-77.11976 38.93434, -77.11253 38.940… … P1_001N 689545 3 12 MULTIPOLYGON (((-80.17628 25.52505, -80.17395 … … P1_001N 21538187 4 13 MULTIPOLYGON (((-81.27939 31.30792, -81.27716 … … P1_001N 10711908
[5 rows x 5 columns]
Interactive preview — 2020 Census population by state
Every geometry=True call returns a GeoDataFrame you can map directly.
Here is the decennial population data from the tab above rendered as a
choropleth:
Coordinate reference system¶
All geometry returned by PyPUMS is in NAD83 (EPSG:4269), which is the
native CRS of the Census Bureau’s cartographic boundary files. You can verify this on any
GeoDataFrame:
Reprojecting¶
You will often need to reproject to a different CRS depending on your use case.
Most web mapping libraries (Leaflet, Mapbox, Folium) expect WGS 84:
For area or distance calculations, use a projected CRS appropriate to your region. For example, EPSG:2229 (NAD83 / California zone 5, in feet):
Tip
The EPSG code you pick should match the part of the country you are studying. epsg.io is a handy search tool for finding the right projection.
Shapefile resolution¶
Cartographic boundary files come in three resolutions. The default
is 500k, which strikes a good balance between detail and file size.
| Resolution | Description | Best for |
|---|---|---|
500k |
1:500,000 (default) | Most maps and analysis |
5m |
1:5,000,000 | National-level overview maps |
20m |
1:20,000,000 | Thumbnail or small inset maps |
You control the resolution through the attach_geometry() function (see below)
or when calling it indirectly via geometry=True. The geometry=True flag
always uses 500k.
Supported geographies¶
The following geography levels have matching cartographic boundary shapefiles:
| Geography | Requires state? |
Notes |
|---|---|---|
state |
No | All 50 states + DC + territories |
county |
No | All US counties |
tract |
Yes | Census tracts (sub-county) |
block group |
Yes | Block groups (sub-tract) |
place |
Yes | Census-designated and incorporated places |
congressional district |
No | US House districts |
zcta |
No | ZIP Code Tabulation Areas |
puma |
Yes | Public Use Microdata Areas |
cbsa |
No | Core-Based Statistical Areas (metros) |
csa |
No | Combined Statistical Areas |
Sub-state geographies require the state parameter
For tract, block group, place, and puma, the Census Bureau
publishes shapefiles per state. PyPUMS will raise a ValueError if you
omit the state parameter for these geography levels.
attach_geometry() standalone function¶
If you already have a DataFrame with a GEOID column (for example, from a
cached query or an external source), you can attach geometry after the fact
with attach_geometry():
GEOID NAME ... estimate moe
0 17031010100 Census Tract 101; Cook County; Illinois … 69460 21834 1 17031010201 Census Tract 102.01; Cook County; Illinois … 49639 24247 2 17031010202 Census Tract 102.02; Cook County; Illinois … 55119 15618 3 17031010300 Census Tract 103; Cook County; Illinois … 65871 14559 4 17031010400 Census Tract 104; Cook County; Illinois … 49017 8306
[5 rows x 5 columns]
[5 rows x 6 columns]
Before and after
Notice how the first print() shows a plain DataFrame with five
columns. After attach_geometry(), the result gains a geometry column
and becomes a GeoDataFrame — ready for mapping or spatial analysis.
Signature:
attach_geometry(
df, # DataFrame with a GEOID column
geography, # geography level name
state=None, # state FIPS or abbreviation
year=2023, # data year
resolution="500k", # "500k", "5m", or "20m"
cache=True, # cache shapefiles locally (via pygris)
) -> GeoDataFrame
Plotting with Altair¶
Altair works naturally with GeoDataFrame objects
returned by PyPUMS. Use mark_geoshape() to create choropleth maps:
import altair as alt
chart = alt.Chart(gdf).mark_geoshape().encode(
color=alt.Color("estimate:Q", title="Median Household Income"),
tooltip=["NAME:N", "estimate:Q"],
).project("albersUsa").properties(
width=600, height=400,
title="Median Household Income by State",
)
chart
Interactive preview — basic choropleth
For more control, use a color scale and configure the legend:
chart = alt.Chart(gdf).mark_geoshape(stroke="white", strokeWidth=0.3).encode(
color=alt.Color(
"estimate:Q",
scale=alt.Scale(scheme="yellowgreenblue"),
legend=alt.Legend(title="Median Income ($)"),
),
tooltip=["NAME:N", alt.Tooltip("estimate:Q", format="$,")],
).project("albersUsa").properties(width=600, height=400)
chart
Interactive preview — styled choropleth with custom color scale
Spatial joins¶
Combine Census data with non-Census spatial data using geopandas.sjoin().
This is useful when you have points (e.g., business locations, crime incidents)
and want to know which Census geography they fall into.
import geopandas as gpd
import pypums
# Census tracts with population data.
tracts = pypums.get_acs(
geography="tract",
variables="B01001_001",
state="TX",
county="201",
year=2023,
geometry=True,
)
# Your own point data (e.g., store locations).
stores = gpd.read_file("stores.geojson")
# Make sure both layers share the same CRS.
stores = stores.to_crs(tracts.crs)
# Spatial join: assign each store to its tract.
joined = gpd.sjoin(stores, tracts, how="left", predicate="within")
print(joined.head())
store_name geometry index_right GEOID NAME variable estimate moe
0 Downtown Grocery POINT (-95.362 29.764) 42 48201312100 Census Tract 3121, Harris County, Texas B01001_001 5241 498
1 Heights Market POINT (-95.392 29.793) 18 48201311300 Census Tract 3113, Harris County, Texas B01001_001 3812 341
2 Midtown Express POINT (-95.383 29.738) 56 48201313200 Census Tract 3132, Harris County, Texas B01001_001 4105 512
3 Galleria Deli POINT (-95.461 29.739) 73 48201410300 Census Tract 4103, Harris County, Texas B01001_001 6723 620
4 Memorial Pantry POINT (-95.503 29.779) 81 48201430100 Census Tract 4301, Harris County, Texas B01001_001 3290 402
Tip
After the join, each store row now carries the tract’s GEOID, NAME,
and population estimate. You can use this to answer questions like
“which stores are in the most populated tracts?” or aggregate store
counts per tract.
Note
sjoin() requires both GeoDataFrames to be in the same CRS. Always
call .to_crs() on one of them before joining if they differ.
Dot density mapping¶
as_dot_density() converts polygon data into random points, where each
dot represents a fixed number of people (or any other count). This is a
popular technique for visualizing racial or ethnic composition at fine
spatial scales.
GEOID ... B03002_012M
0 06001400100 … 107 1 06001400200 … 95 2 06001400300 … 305 3 06001400400 … 280 4 06001400500 … 236
[5 rows x 11 columns]
geometry value
0 POINT (-122.22019 37.86309) White 1 POINT (-122.24366 37.86566) White 2 POINT (-122.21322 37.858) White 3 POINT (-122.22063 37.86959) White 4 POINT (-122.24579 37.84953) White 5 POINT (-122.25481 37.84437) White 6 POINT (-122.25605 37.83734) White 7 POINT (-122.24739 37.84439) White 8 POINT (-122.24799 37.84545) White 9 POINT (-122.25714 37.84095) White
# Extract point coordinates for Altair.
dots = dots.copy()
dots["lon"] = dots.geometry.x
dots["lat"] = dots.geometry.y
# Plot.
import altair as alt
chart = alt.Chart(dots).mark_circle(size=1, opacity=0.6).encode(
longitude="lon:Q",
latitude="lat:Q",
color=alt.Color(
"value:N",
scale=alt.Scale(
domain=["White", "Black", "Asian", "Hispanic"],
range=["#1f77b4", "#2ca02c", "#d62728", "#ff7f0e"],
),
),
tooltip=["value:N"],
).project("mercator").properties(width=600, height=600)
chart
Interactive preview — dot density map of Alameda County, CA (sample data)
Each dot represents 500 people. The spatial distribution of dots reveals neighborhood-level patterns of racial and ethnic composition.
Signature:
as_dot_density(
gdf, # GeoDataFrame with polygon geometries
values, # dict of {column_name: label}
dots_per_value=100, # data units per dot
seed=None, # random seed for reproducibility
) -> GeoDataFrame # point GeoDataFrame with "geometry" and "value" columns
Population-weighted interpolation¶
When you need to transfer data from one set of boundaries to another (for
example, from 2010 tracts to 2020 tracts, or from tracts to school districts),
use interpolate_pw(). It distributes values proportionally based on
population weights.
from pypums.spatial import interpolate_pw
# Source: income data on 2010 tracts.
from_gdf = ... # GeoDataFrame with "median_income" and "POP" columns
# Target: 2020 tract boundaries.
to_gdf = ... # GeoDataFrame with 2020 tract polygons
result = interpolate_pw(
from_gdf,
to_gdf,
value_col="median_income",
weight_col="POP",
extensive=True,
)
print(result.head())
GEOID NAME median_income geometry
0 17031010100 Census Tract 101, Cook County, Illinois 58421.3 POLYGON ((-87.63 41.90, -87.63 41.90, -87.62...
1 17031010201 Census Tract 102.01, Cook County, Illinois 73812.7 POLYGON ((-87.63 41.91, -87.63 41.91, -87.62...
2 17031010202 Census Tract 102.02, Cook County, Illinois 45129.5 POLYGON ((-87.64 41.91, -87.64 41.91, -87.63...
3 17031010300 Census Tract 103, Cook County, Illinois 36291.0 POLYGON ((-87.64 41.92, -87.64 41.92, -87.63...
4 17031010400 Census Tract 104, Cook County, Illinois 30187.4 POLYGON ((-87.65 41.92, -87.65 41.91, -87.64...
Info
The result is a GeoDataFrame with the target boundaries and the interpolated values. Population-weighted interpolation ensures that areas with more people contribute more to the transferred value than sparsely populated areas.
Signature:
interpolate_pw(
from_gdf, # source GeoDataFrame
to_gdf, # target GeoDataFrame
value_col, # column name in from_gdf to interpolate
weight_col="POP", # column name in from_gdf with population weights
extensive=True, # True for counts (sum), False for rates (weighted avg)
) -> GeoDataFrame
| Parameter | extensive=True |
extensive=False |
|---|---|---|
| Use for | Counts, totals | Rates, medians, averages |
| Method | Proportional distribution | Population-weighted average |
| Example | Total population | Median household income |
Memory considerations¶
Spatial queries download shapefiles and store full polygon geometries in memory. A few tips for working with large datasets:
Keep memory usage in check
- Use a coarser resolution (
5mor20m) if you only need a national overview and do not need precise boundaries. - Filter by state whenever possible. Tract-level shapefiles for a single state are much smaller than a nationwide download.
-
Drop the geometry column when you no longer need it:
-
Save to file and work from disk instead of re-downloading:
-
Shapefile caching is automatic. PyPUMS uses pygris with caching enabled, so shapefiles are downloaded once and reused from a local cache directory (
~/.cache/pygris/on Linux,~/Library/Caches/pygris/on macOS).
Troubleshooting¶
ImportError: geopandas is required for spatial operations
: Install the spatial extra: uv add "pypums[spatial]". This pulls in
geopandas, pygris, shapely, and pyproj.
ValueError: geography='tract' requires a state parameter
: Sub-state geographies (tract, block group, place, puma) need a
state argument. Pass a FIPS code, abbreviation, or full name.
Geometry column is all None
: The Census Bureau may not have shapefiles for the geography level and
year you requested. Try a different year or a broader geography (e.g.,
county instead of block group).
CRS mismatch when combining DataFrames : All PyPUMS geometry is returned in EPSG:4269 (NAD83). If you are combining with data in a different CRS, reproject first:
```python
gdf = gdf.to_crs(epsg=4326) # convert to WGS84
print(gdf.crs)
```
```
EPSG:4326
```
MemoryError with tract-level national data
: Nationwide tract shapefiles are very large. Filter by state, use a
coarser resolution (20m), or work state-by-state and concatenate
after dropping unneeded columns.
Dot-density map is slow or crashes
: as_dot_density() generates one point per dots_per_value people. For
large areas, increase dots_per_value (e.g., dots_per_value=100
instead of dots_per_value=1) or subset to a single county before
generating dots.
See Also¶
- ACS Data — Using
get_acs(geometry=True)to retrieve spatial ACS data - Geography & FIPS — Supported geographies and FIPS code structure