Exploring U.S. Census Data Across Seven MSAs with the IPUMS NHGIS API

import sys, subprocess, importlib.util, requests, io, zipfile, time
import pandas as pd
import numpy as np

for pkg in ["geopandas", "matplotlib", "folium", "mapclassify"]:
    if importlib.util.find_spec(pkg) is None:
        subprocess.check_call([sys.executable, "-m", "pip", "install", pkg])

import geopandas as gpd
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from matplotlib.ticker import FuncFormatter
import seaborn as sns
import warnings
warnings.filterwarnings("ignore")

API_KEY = "YOUR_IPUMS_API_KEY"
headers = {"Authorization": API_KEY}

url = "https://api.ipums.org/metadata/datasets"
params = {"collection": "nhgis", "version": 2}
response = requests.get(url, params=params, headers=headers)
response.raise_for_status()
data = response.json()

def get_msa_cbg(state_fips, county_fips, year=2020):
    url = f"https://www2.census.gov/geo/tiger/TIGER{year}/BG/tl_{year}_{state_fips}_bg.zip"
    gdf = gpd.read_file(url)
    return gdf[gdf["COUNTYFP"].isin(county_fips)].copy()

TCMA_county_fips = ["003", "019", "037", "053", "123", "139", "163"]
LA_county_fips = ["037", "059"]
Houston_county_fips = ["015", "039", "071", "157", "167", "201", "291", "339", "407", "473"]
Seattle_county_fips = ["033", "053", "061"]
Miami_county_fips = ["011", "086", "099"]
Atlanta_county_fips = ["013", "015", "035", "045", "057", "063", "067", "077", "085", "089", "097", "113", "117", "121", "135", "143", "149", "151", "159", "187", "199", "211", "217", "223", "227", "231", "247", "255", "297"]
Denver_county_fips = ["001", "005", "014", "019", "031", "035", "039", "047", "059", "093"]

TCMA_cbg = get_msa_cbg("27", TCMA_county_fips)
LA_cbg = get_msa_cbg("06", LA_county_fips)
Houston_cbg = get_msa_cbg("48", Houston_county_fips)
Seattle_cbg = get_msa_cbg("53", Seattle_county_fips)
Miami_cbg = get_msa_cbg("12", Miami_county_fips)
Atlanta_cbg = get_msa_cbg("13", Atlanta_county_fips)
Denver_cbg = get_msa_cbg("08", Denver_county_fips)

msa_dict = {
    "TCMA": TCMA_cbg,
    "Los Angeles": LA_cbg,
    "Houston": Houston_cbg,
    "Seattle": Seattle_cbg,
    "Miami": Miami_cbg,
    "Atlanta": Atlanta_cbg,
    "Denver": Denver_cbg
}

all_msa_cbg = gpd.GeoDataFrame(
    pd.concat(
        [gdf.to_crs(TCMA_cbg.crs).assign(MSA=name) for name, gdf in msa_dict.items()],
        ignore_index=True
    ),
    geometry="geometry",
    crs=TCMA_cbg.crs
)

all_msa_cbg.explore(
    column="MSA",
    categorical=True,
    tooltip=["MSA", "GEOID"],
    #tiles="CartoDB positron"
)

def summarize_msa_cbg(gdf, name="MSA", plot=True):
    print(f"{name}: {len(gdf):,} CBGs")
    print(f"Counties: {gdf['COUNTYFP'].nunique()}")
    if plot:
        ax = gdf.plot(figsize=(8, 8), color="lightblue", edgecolor="black", linewidth=0.2)
        ax.set_title(f"{name} Census Block Groups (n={len(gdf):,})")
        ax.axis("off")
    #return gdf[["STATEFP", "COUNTYFP", "GEOID"]]

summarize_msa_cbg(TCMA_cbg, "Twin Cities")
summarize_msa_cbg(LA_cbg, "Los Angeles")
summarize_msa_cbg(Houston_cbg, "Houston")
summarize_msa_cbg(Seattle_cbg, "Seattle")
summarize_msa_cbg(Miami_cbg, "Miami")
summarize_msa_cbg(Atlanta_cbg, "Atlanta")
summarize_msa_cbg(Denver_cbg, "Denver")

headers = {"Authorization": API_KEY}
dataset = "2017_2021_ACS5a" # 2021 ACS 5-year estimates

# Check table metadata
meta_url = f"https://api.ipums.org/metadata/datasets/{dataset}"
meta = requests.get(meta_url, headers=headers, params={"collection": "nhgis", "version": 2}).json()

tables = pd.DataFrame(meta["dataTables"])
tables[tables["description"].str.contains("Total Population|Median Household Income", case=False, na=False)]

extract = {
    "datasets": {
        dataset: {
            "dataTables": ["B01003", "B19013"], # "B01003" for total population, and "B19013" for median household income
            "geogLevels": ["blck_grp"] # Census block group (CBG) as the geographic level
        }
    },
    "geographicExtents": ["080", "120", "130", "270", "480", "530", "060"],  # the extent
    "dataFormat": "csv_no_header",
    "breakdownAndDataTypeLayout": "single_file",
    "description": "2021 ACS CBG population and household income"
}

url = "https://api.ipums.org/extracts/"
r = requests.post(url, headers={**headers, "Content-Type": "application/json"}, params={"collection": "nhgis", "version": 2}, json=extract)
result = r.json()
extract_number = result["number"]

print("Extract number:", extract_number)
print("Status:", result["status"])

status_url = f"https://api.ipums.org/extracts/{extract_number}"

while True:
    status = requests.get(status_url, headers=headers, params={"collection": "nhgis", "version": 2}).json()
    if status["status"] == "completed":
        break
    print(status["status"])
    time.sleep(5)

download_url = status["downloadLinks"]["tableData"]["url"]
zip_data = requests.get(download_url, headers=headers).content

with zipfile.ZipFile(io.BytesIO(zip_data)) as z:
    print(z.namelist())
    csv_file = [f for f in z.namelist() if f.endswith(".csv")][0]
    acs = pd.read_csv(z.open(csv_file))

acs.head(2)

msa_counties = {
    "TCMA": ("27", TCMA_county_fips),
    "LA": ("06", LA_county_fips),
    "Houston": ("48", Houston_county_fips),
    "Seattle": ("53", Seattle_county_fips),
    "Miami": ("12", Miami_county_fips),
    "Atlanta": ("13", Atlanta_county_fips),
    "Denver": ("08", Denver_county_fips)
}

msa_acs = {}

for name, (state, counties) in msa_counties.items():
    msa_acs[name] = acs[
        (acs["STATEA"].astype(str).str.zfill(2) == state) &
        (acs["COUNTYA"].astype(str).str.zfill(3).isin(counties))
    ].copy()

def link_acs_boundary(acs_df, boundary):
    acs_df = acs_df.copy()
    boundary = boundary.copy()
    acs_df["GEOID"] = acs_df["GEO_ID"].str[-12:] # match the GEOID and GEO_ID by their structures
    boundary["GEOID"] = boundary["GEOID"].astype(str)
    return boundary.merge(acs_df, on="GEOID", how="left")

msa_boundaries = {
    "TCMA": TCMA_cbg,
    "LA": LA_cbg,
    "Houston": Houston_cbg,
    "Seattle": Seattle_cbg,
    "Miami": Miami_cbg,
    "Atlanta": Atlanta_cbg,
    "Denver": Denver_cbg
}

msa_cbg_acs = {name: link_acs_boundary(msa_acs[name], boundary) for name, boundary in msa_boundaries.items()}

TCMA_cbg_acs = msa_cbg_acs["TCMA"]
LA_cbg_acs = msa_cbg_acs["LA"]
Houston_cbg_acs = msa_cbg_acs["Houston"]
Seattle_cbg_acs = msa_cbg_acs["Seattle"]
Miami_cbg_acs = msa_cbg_acs["Miami"]
Atlanta_cbg_acs = msa_cbg_acs["Atlanta"]
Denver_cbg_acs = msa_cbg_acs["Denver"]

for name, gdf in msa_cbg_acs.items():
    print(name, len(gdf), gdf["GEO_ID"].notna().sum())

def map_msa_socio(gdf, name, population_col="AON4E001", income_col="AOQIE001"):
    gdf = gdf.copy()
    gdf_proj = gdf.to_crs(gdf.estimate_utm_crs())
    gdf["pop_density"] = gdf[population_col] / (gdf_proj.geometry.area / 1e6) # calculate the population density

    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    gdf.plot(column="pop_density", cmap="YlOrRd", scheme="quantiles", k=5, legend=True, ax=axes[0], edgecolor="none")
    axes[0].set_title(f"{name} - Population Density")  # use pop density instead
    axes[0].axis("off")

    gdf.plot(column=income_col, cmap="YlGnBu", scheme="quantiles", k=5, legend=True, ax=axes[1], edgecolor="none")
    axes[1].set_title(f"{name} - Median Household Income")
    axes[1].axis("off")

    plt.tight_layout()
    plt.show()
    return gdf

for name, gdf in msa_cbg_acs.items():
    msa_cbg_acs[name] = map_msa_socio(gdf, name)

def violin_msa_socio(msa_cbg_acs, population_col="AON4E001", income_col="AOQIE001",
                     density_symlog=False, save_path=None):
    records = []
    for name, gdf in msa_cbg_acs.items():
        if gdf.empty:
            continue
        area = gdf.to_crs(gdf.estimate_utm_crs()).geometry.area / 1e6
        population = pd.to_numeric(gdf[population_col], errors="coerce")
        income = pd.to_numeric(gdf[income_col], errors="coerce")
        records.append(pd.DataFrame({
            "MSA": name,
            "pop_density": population.where(population >= 0) / area.where(area > 0),
            "income": income.where(income >= 0),
        }))

    df = pd.concat(records, ignore_index=True).replace([np.inf, -np.inf], np.nan)
    order = df.groupby("MSA")["pop_density"].median().sort_values(ascending=False).index

    with sns.axes_style("whitegrid"):
        fig, axes = plt.subplots(1, 2, figsize=(13, max(4, 0.45 * len(order) + 1.5)),
                                 sharey=True, layout="constrained")
        panels = [
            ("pop_density", "Population density", "People/km²", "#D96B6B"),
            ("income", "Median household income", "US dollars", "#55AFA5"),
        ]
        for ax, (column, title, xlabel, color) in zip(axes, panels):
            sns.violinplot(data=df, x=column, y="MSA", order=order, orient="h",
                           color=color, inner="quart", cut=0, linewidth=0.8, saturation=1, ax=ax)
            ax.set(title=title, xlabel=xlabel, ylabel="", xlim=(0, None))
            ax.set_axisbelow(True)
            ax.grid(axis="x", color="#E5E5E5", linewidth=0.7)
            ax.grid(axis="y", visible=False)
            ax.tick_params(axis="y", length=0)
            sns.despine(ax=ax, left=True)

        if density_symlog:
            axes[0].set_xscale("symlog", linthresh=100)
            axes[0].set_xlabel("People/km² (symlog; linear below 100)")

        axes[0].xaxis.set_major_formatter(FuncFormatter(lambda x, _: f"{x:,.0f}"))
        axes[1].xaxis.set_major_formatter(FuncFormatter(lambda x, _: f"${x:,.0f}"))

        if save_path:
            fig.savefig(save_path, dpi=600, bbox_inches="tight")
        plt.show()

    return df


socio_df = violin_msa_socio(msa_cbg_acs, density_symlog=True)

from matplotlib.patches import Patch

def map_joint_socio(msa_cbg_acs, population_col="AON4E001", income_col="AOQIE001",
                    within_msa=True, ncols=2, save_path=None):
    records, crs_by_msa = [], {}
    for name, gdf in msa_cbg_acs.items():
        if gdf.empty:
            continue
        crs_by_msa[name] = gdf.estimate_utm_crs()
        area = gdf.to_crs(crs_by_msa[name]).geometry.area / 1e6
        g = gdf.to_crs(4326)[["geometry"]].copy()
        population = pd.to_numeric(gdf[population_col], errors="coerce")
        income = pd.to_numeric(gdf[income_col], errors="coerce")
        g["MSA"] = name
        g["pop_density"] = population.where(population >= 0) / area.where(area > 0)
        g["income"] = income.where(income >= 0)
        records.append(g)

    maps = gpd.GeoDataFrame(pd.concat(records, ignore_index=True), crs=4326)
    cols = ["pop_density", "income"]
    maps[cols] = maps[cols].replace([np.inf, -np.inf], np.nan)
    valid = maps.dropna(subset=cols).copy()
    thresholds = valid.groupby("MSA")[cols].transform("median") if within_msa else valid[cols].median()
    high = valid[cols].ge(thresholds)

    labels = [
        "High density · Low income",
        "High density · High income",
        "Low density · Low income",
        "Low density · High income",
    ]
    colors = dict(zip(labels, ["#D96B6B", "#8064A2", "#D8D8D8", "#55AFA5"]))
    valid["profile"] = np.select(
        [high.pop_density & ~high.income, high.pop_density & high.income,
         ~high.pop_density & ~high.income],
        labels[:3], default=labels[3],
    )
    maps["profile"] = valid["profile"].reindex(maps.index).fillna("No data")
    colors["No data"] = "#F2F2F2"

    names = list(crs_by_msa)
    ncols = min(max(1, ncols), len(names))
    nrows = int(np.ceil(len(names) / ncols))
    fig, axes = plt.subplots(nrows, ncols, figsize=(6 * ncols, 5 * nrows), squeeze=False)

    for ax, name in zip(axes.flat, names):
        g = maps.loc[maps["MSA"].eq(name)].to_crs(crs_by_msa[name])
        g.plot(ax=ax, color=g["profile"].map(colors), edgecolor="white", linewidth=0.1)
        ax.set_title(name, fontsize=12, fontweight="bold")
        ax.set_axis_off()

    for ax in axes.flat[len(names):]:
        ax.set_axis_off()

    scope = "Within-MSA" if within_msa else "Pooled"
    fig.suptitle(f"Joint socioeconomic profiles — {scope.lower()} median thresholds", fontsize=14)
    fig.legend(handles=[Patch(facecolor=c, label=s) for s, c in colors.items()],
               loc="lower center", ncol=3, frameon=False)
    fig.tight_layout(rect=(0, 0.09, 1, 0.95))
    if save_path:
        fig.savefig(save_path, dpi=300, bbox_inches="tight")
    plt.show()

    return maps


profile_maps = map_joint_socio(msa_cbg_acs, within_msa=True)

def joint_socio_profiles(socio_df, within_msa=True, save_path=None):
    df = socio_df.copy()
    cols = ["pop_density", "income"]
    df[cols] = df[cols].apply(pd.to_numeric, errors="coerce")
    df[cols] = df[cols].replace([np.inf, -np.inf], np.nan)
    df = df.dropna(subset=["MSA"] + cols)
    df = df.loc[df[cols].ge(0).all(axis=1)].copy()

    thresholds = df.groupby("MSA")[cols].transform("median") if within_msa else df[cols].median()
    high = df[cols].ge(thresholds)

    labels = [
        "High density · Low income",
        "High density · High income",
        "Low density · Low income",
        "Low density · High income",
    ]
    df["profile"] = np.select(
        [high.pop_density & ~high.income, high.pop_density & high.income,
         ~high.pop_density & ~high.income],
        labels[:3], default=labels[3],
    )

    shares = pd.crosstab(df["MSA"], df["profile"], normalize="index")
    shares = shares.reindex(columns=labels, fill_value=0).mul(100)
    shares = shares.sort_values(labels[0], ascending=False)

    ax = shares.plot.barh(stacked=True, figsize=(11, max(4, 0.45 * len(shares))),
                          color=["#D96B6B", "#8064A2", "#D8D8D8", "#55AFA5"], width=0.75)
    for container in ax.containers:
        ax.bar_label(container, labels=[f"{v:.0f}%" if v >= 8 else "" for v in container.datavalues],
                     label_type="center", fontsize=9)

    ax.invert_yaxis()
    ax.set(xlim=(0, 100), xlabel="Share of valid CBGs (%)", ylabel="",
           title="Joint socioeconomic profiles")
    ax.legend(title=None, loc="upper center", bbox_to_anchor=(0.5, -0.16), ncol=2, frameon=False)
    ax.spines[["top", "right"]].set_visible(False)
    plt.tight_layout()
    if save_path:
        ax.figure.savefig(save_path, dpi=300, bbox_inches="tight")
    plt.show()

    return df, shares


profile_df, profile_shares = joint_socio_profiles(socio_df)