NYC Food Access: SNAP Retailer Analysis

Analyze SNAP retailer density by income level across NYC ZIP codes.

Introduction

30 Second Summary

Not every neighborhood has equal access to grocery stores. In cities like New York, the number of stores that accept food assistance benefits can vary dramatically from one ZIP code to the next.

In this project, you will use real government data to answer a specific question: do lower-income NYC ZIP codes have fewer SNAP-authorized food retailers per resident than higher-income ones? You will load data from two sources, clean it, join it on a shared geographic key, compute a normalized metric, and produce a visual comparison with a one-sentence finding.

What You'll Build

You can open your finished notebook, run every cell top to bottom, and walk someone through a data-driven finding about food access inequality in New York City. Every number traces back to a government source, and every decision you made (why per-capita, why you dropped certain rows, why you split at the median) has a reason you can explain in an interview.

By the end of this project, you'll have:

  • A Jupyter notebook that pulls live data from the USDA and the Census Bureau, cleans it, and merges it on a shared geographic key so you can trace every number back to its origin.
  • A per-capita metric that normalizes retailer counts by population, letting you make fair comparisons across ZIP codes of different sizes instead of falling into the "more people equals more stores" trap.
  • A scatter plot and bar chart built with matplotlib that make your finding visible at a glance, ready to drop into a presentation or portfolio.
  • Secret Mission: Break down your analysis by store type to discover whether the overall pattern hides a more interesting story about supermarkets vs. convenience stores.

Are there any prerequisites?

You need basic Python syntax (variables, loops, functions). No prior experience with pandas, numpy, or matplotlib is required. This project teaches you those libraries as you go.

Before We Start

Every data analysis project starts with a question. Without a clear question, you end up exploring aimlessly and producing charts that look impressive but say nothing specific.

In this step, you will set up everything you need to start analyzing real government data. You will install the Python libraries, get access to the Census Bureau API, and download the SNAP retailer dataset so that Step 2 can begin with loading and exploring data immediately.

In this step, get ready to:
  • Install pandas, numpy, matplotlib, and requests.
  • Get a free Census API key for pulling income and population data.
  • Download the SNAP Retailer Locator CSV from the USDA.
Create your notebook and frame the research question

Good analysis starts with a specific, falsifiable question. Writing it down first keeps you honest. Every decision you make later (which columns to keep, how to normalize, where to split) flows from this question.

  • Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
  • Type VS Code and press Enter to open it.
  • Open your project folder by clicking File in the top menu bar, then Open Folder.
  • Select your project folder (the one with the data/ subfolder already inside it) and click Open.
  • Create a new Jupyter notebook by pressing Cmd+Shift+P (macOS) or Ctrl+Shift+P (Windows) to open the Command Palette.
  • Type Jupyter: Create New Blank Notebook and select it.
  • Save the notebook immediately as analysis.ipynb in your project folder root (not inside data/).

Your first cell should be a Markdown cell that states the question driving this entire project.

  • Click on the first cell in the notebook.
  • Change the cell type from Code to Markdown using the dropdown at the bottom-right of the cell (or press M while the cell is selected but not in edit mode).
  • Type your research question. A strong version: "Do NYC ZIP codes with below-median income have fewer SNAP retailers per 1,000 residents than ZIP codes with above-median income?"
  • Run the cell by pressing Shift+Enter. You should see your question rendered as formatted text.

Why write the question first?

The research question is your compass. When you later face decisions (should I keep convenience stores in the count? inner join or left join?), you can ask: which choice answers my question more honestly? Without this, data analysis drifts into aimless exploration.

Install the Python libraries

You need four libraries for this project. pandas handles tabular data. numpy provides numerical operations. matplotlib creates charts. requests makes HTTP calls to the Census API.

  • Open a terminal inside VS Code by pressing Ctrl+` (backtick key, below Escape).
  • Install all four libraries by running this command:
pip install pandas numpy matplotlib requests

What does this command do?

This tells pip (Python's package manager) to download and install four packages from the Python Package Index. You only need to do this once per environment.

Now add your imports cell to the notebook so the libraries are available in every future cell.

  • Click below your Markdown cell to create a new code cell in analysis.ipynb.
  • Type the following imports into the cell:
# Cell 1: Imports
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import requests

What does this code do?

  • import pandas as pd loads the pandas library and gives it the shorthand pd so you can write pd.read_csv() instead of pandas.read_csv().
  • import numpy as np loads numpy for numerical operations like computing medians.
  • import matplotlib.pyplot as plt loads the plotting interface you will use to create scatter plots and bar charts.
  • import requests loads the HTTP library you will use to call the Census API.
  • Run the cell by pressing Shift+Enter. If it completes with no error, the imports succeeded.

Before you move on, do you think the following verification command will succeed or produce an error?

  • Confirm all four libraries are importable by running this in your VS Code terminal:
python -c "import pandas; import numpy; import matplotlib; import requests; print('All good')"

You should see All good printed with no errors.

Seeing ModuleNotFoundError?

  • Make sure your pip install command completed without errors. Scroll up in the terminal to check for red text.
  • If your system has both Python 2 and Python 3, try pip3 install pandas numpy matplotlib requests and then use python3 in the verification command.
  • help me fix a ModuleNotFoundError when importing pandas or numpy
Get your Census API key and download the SNAP data

You need two external data sources for this project. The Census Bureau provides income and population numbers by ZIP code. The USDA provides the list of every store in the country authorized to accept SNAP benefits.

  • Open your browser and go to the Census API key signup page:
https://api.census.gov/data/key_signup.html
  • Check the I agree to the terms of service checkbox.
  • Enter your name and email address, then submit the form.
  • Check your email inbox. You will receive an email with your API key code in the message. Save this key somewhere safe (a notes app or text file on your computer). You will paste it into your notebook in a later step.

What is this API key for?

The Census Bureau requires an API key for all queries to its data API. The key is free and does not expire. It identifies your requests so the Census Bureau can enforce rate limits (500 queries per day without a key, unlimited with one).

Use a .com, .net, .org, .gov, or .edu email

The Census Bureau requires that the email you use to request the key has one of these domain extensions. A personal Gmail or Outlook address works fine.

  • Now navigate to the SNAP Retailer Locator data page:
https://www.fns.usda.gov/snap/retailer-locator
  • Scroll down to the SNAP Retailer Locator Data section.
  • Download the current SNAP retailer CSV file (the one with latitude and longitude coordinates, NOT the historical data zip file).
  • Save the downloaded CSV file into your project's data/ folder.
  • Rename the file to snap_retailers.csv so it matches the code you will write later.

What is in this dataset?

The SNAP Retailer Locator CSV is a national dataset of every store currently authorized to accept SNAP benefits. It includes store name, store type, address, ZIP code, and latitude/longitude coordinates. You will filter it down to NYC in the next step.

  • Create an empty output/ folder in your project root for saving figures later. Run this in your terminal:
mkdir output
  • Confirm your project folder has the right structure by running:
ls

You should see analysis.ipynb, data/, and output/ listed.

Don't see the CSV in your data folder?

  • Check your browser's Downloads folder. The file may still be there instead of in data/. Move it manually.
  • Make sure you downloaded the current data CSV, not the historical zip file. The current file is a plain .csv file.
  • help me find and move the SNAP retailer CSV into my project folder

✔️ Awesome, I've got everything!

Great. Double check that you have saved analysis.ipynb with the imports cell, and that snap_retailers.csv is in your data/ folder.

ⓧ I'd like to double check the full code

Your analysis.ipynb should have two cells so far:

  • Cell 1 (Markdown): Your research question about SNAP retailer density and income in NYC.
  • Cell 2 (Code): The imports block shown below.
# Cell 1: Imports
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import requests

Your environment is ready and your data is in place. Next up, you will load the SNAP retailer CSV into a DataFrame, explore what is inside, and count how many retailers exist per NYC ZIP code.

Load and Explore SNAP Retailer Data

Your environment is set up with all four libraries installed and your SNAP retailer CSV saved locally. Now it is time to actually look at the data.

Before you can compute anything meaningful, you need to understand what columns exist in the dataset, filter it down to NYC, and count how many retailers appear per ZIP code. This gives you the numerator for the per-capita calculation you will build in a later step.

In this step, get ready to:
  • Load the SNAP retailer CSV and identify the columns you need.
  • Filter the dataset to NYC ZIP codes only.
  • Count retailers per ZIP and inspect the store type distribution.
Load the CSV and inspect it

The SNAP Retailer Locator CSV is a national dataset containing every store authorized to accept SNAP benefits. It includes store name, type, address, state, ZIP code, and coordinates.

The exact column names depend on the version of the file you downloaded. Part of this step is you figuring out what those names are by running .info() on the loaded DataFrame.

  • In your notebook, create a new cell below your imports (Cell 2) and load the CSV by running this code:
# Cell 2: Load SNAP retailer data
snap = pd.read_csv('data/snap_retailers.csv')
print(snap.info())
print(snap.head())

What does this code do?

  • pd.read_csv() reads your CSV file into a pandas DataFrame called snap.
  • .info() prints a summary showing every column name, its data type, and how many non-null values it has.
  • .head() shows the first five rows so you can see actual data values.
  • Run the cell.

You should see a table of column names with their dtypes (like object, int64, float64) and a preview of the first few rows. Look for the column that contains the ZIP code, the column that contains the state, and the column that contains the store type.

Seeing a FileNotFoundError?

Make sure your CSV is saved at data/snap_retailers.csv relative to where your notebook lives. If your file has a different name, update the path in pd.read_csv() to match.

help me fix a FileNotFoundError when loading my SNAP CSV in pandas

Now that you can see the column names, assign them to variables so the rest of your code stays readable even if column names change.

  • Create a new cell (Cell 3) and assign your column name variables by running this code:
# Cell 3: Identify your columns
ZIP_COL = 'Zip5'        # ADJUST THIS based on your .info() output
STATE_COL = 'State'     # ADJUST THIS
TYPE_COL = 'Store_Type' # ADJUST THIS

What does this code do?

You are creating three variables that store the exact column names from YOUR data. The values shown here (Zip5, State, Store_Type) are based on typical USDA formatting.

Check your .info() output from the previous cell. If your column names differ, change the strings in quotes to match exactly what you see.

  • Compare the variable values to the column names from your .info() output and adjust if needed.
  • Run the cell.

Not sure which column is the ZIP?

Look at your .head() output. The ZIP column will contain five-digit numbers (like 10001 or 11201). The state column will contain two-letter abbreviations (like NY). The store type column will contain descriptive text (like names of store categories).

help me identify which columns in my SNAP data correspond to ZIP code, state, and store type

Filter to NYC and count by ZIP

The CSV contains retailers from across the entire United States. You only care about NYC, so you need to filter by both state and ZIP code prefix.

  • Create a new cell (Cell 4) and filter the data to NYC by running this code:
# Cell 4: Filter to NYC
nyc_zips_mask = snap[ZIP_COL].astype(str).str[:3].isin(
    ['100', '101', '102', '103', '104', '110', '111', '112', '113', '114', '116']
)
state_mask = snap[STATE_COL].str.upper().isin(['NY', 'NEW YORK'])
snap_nyc = snap[nyc_zips_mask & state_mask].copy()
print(f"NYC SNAP retailers: {len(snap_nyc)}")

What does this code do?

NYC ZIP codes start with specific three-digit prefixes: 100-104 for Manhattan and the Bronx, 110-114 and 116 for Brooklyn, Queens, and Staten Island.

  • nyc_zips_mask converts the ZIP column to strings, takes the first three characters, and checks if they match any NYC prefix.
  • state_mask ensures the row is actually in New York (not a ZIP collision with another state).
  • The & combines both conditions so only rows matching BOTH pass through.
  • .copy() creates an independent DataFrame so edits to snap_nyc do not accidentally modify the original.
  • Run the cell.

You should see a count printed showing how many SNAP retailers fall within NYC ZIP codes. Expect a number in the thousands.

Seeing 0 NYC SNAP retailers?

This usually means your column names do not match. Go back to Cell 3 and verify that ZIP_COL and STATE_COL exactly match the column headers from your .info() output (including capitalisation).

help me debug why my NYC filter returns zero rows

Now you know how many total retailers exist in NYC. Next, group them by ZIP to see how retailers are distributed.

  • Create a new cell (Cell 5) and count retailers per ZIP by running this code:
# Cell 5: Count retailers per ZIP
zip_counts = snap_nyc.groupby(ZIP_COL).size().reset_index(name='retailer_count')
print(f"Number of unique NYC ZIPs with SNAP retailers: {len(zip_counts)}")
print(zip_counts.sort_values('retailer_count', ascending=False).head(10))

What does this code do?

  • .groupby(ZIP_COL).size() groups all rows by ZIP code and counts how many retailers fall in each one.
  • .reset_index(name='retailer_count') converts the result from a Series back into a clean DataFrame with a named column.
  • The sorted .head(10) shows you which ZIP codes have the most SNAP retailers.

Before you run it, think about what you expect. Do you think all NYC ZIP codes will have roughly the same number of retailers, or will some have dramatically more?

  • Run the cell.

You should see a list of ZIP codes (five-digit strings starting with 100, 110, 112, 103, 104, or 116) alongside their retailer counts. Most NYC ZIPs will have between 10 and 200+ retailers.

Seeing unexpected ZIP codes or very few results?

If you see ZIPs that do not start with NYC prefixes, your filter in Cell 4 may not have applied correctly. Re-run Cell 4 first, then re-run Cell 5. Also confirm you ran the cells in order (Jupyter cells must execute sequentially to share variables).

help me fix my groupby results showing wrong ZIP codes

Inspect your counts and make a decision about store types

Some ZIP codes probably have surprisingly high retailer counts. Before moving on, check what kinds of stores are driving those numbers.

  • Create a new cell (Cell 6) and inspect the store type distribution by running this code:
# Cell 6: Inspect store types
print(snap_nyc[TYPE_COL].value_counts())

What does this code do?

value_counts() tallies how many times each unique store type appears in your NYC-filtered data. This reveals whether most SNAP retailers are supermarkets, convenience stores, gas stations, or something else entirely.

  • Run the cell.

You will likely notice that convenience stores and small groceries vastly outnumber supermarkets. For now, keep all store types included in your analysis.

Why keep all store types for now?

Lumping all types together gives you the broadest view of SNAP access. But this decision has consequences. A ZIP code with 150 convenience stores and 0 supermarkets tells a very different story than one with 30 supermarkets.

Note this decision. You will revisit it in the Secret Mission, where disaggregating by store type may reveal a completely different pattern.

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook with all six cells.

ⓧ I'd like to double check the full code

Your notebook should contain these cells in order. Confirm each one is present and matches what you have.

Cell 1: Imports (from Step 1).

# Cell 1: Imports
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import requests

Cell 2: Load SNAP retailer data.

# Cell 2: Load SNAP retailer data
snap = pd.read_csv('data/snap_retailers.csv')
print(snap.info())
print(snap.head())

Cell 3: Identify your columns.

# Cell 3: Identify your columns
ZIP_COL = 'Zip5'        # ADJUST THIS based on your .info() output
STATE_COL = 'State'     # ADJUST THIS
TYPE_COL = 'Store_Type' # ADJUST THIS

Cell 4: Filter to NYC.

# Cell 4: Filter to NYC
nyc_zips_mask = snap[ZIP_COL].astype(str).str[:3].isin(
    ['100', '101', '102', '103', '104', '110', '111', '112', '113', '114', '116']
)
state_mask = snap[STATE_COL].str.upper().isin(['NY', 'NEW YORK'])
snap_nyc = snap[nyc_zips_mask & state_mask].copy()
print(f"NYC SNAP retailers: {len(snap_nyc)}")

Cell 5: Count retailers per ZIP.

# Cell 5: Count retailers per ZIP
zip_counts = snap_nyc.groupby(ZIP_COL).size().reset_index(name='retailer_count')
print(f"Number of unique NYC ZIPs with SNAP retailers: {len(zip_counts)}")
print(zip_counts.sort_values('retailer_count', ascending=False).head(10))

Cell 6: Inspect store types.

# Cell 6: Inspect store types
print(snap_nyc[TYPE_COL].value_counts())

You now have the numerator for your analysis: retailer counts per ZIP code. Next up, you will pull population and income data from the Census API so you can normalize these counts and split by income level.

Pull Census Data and Merge

You now have a count of SNAP retailers per ZIP code for NYC. But raw counts on their own don't answer your research question.

A ZIP code with 200 retailers and 500,000 residents has very different food access than one with 200 retailers and 20,000 residents. You need population and income data to normalize your counts and split by income level. The Census Bureau provides both through a free API, keyed to geographic units called ZCTAs (ZIP Code Tabulation Areas) that align closely with the ZIP codes in your SNAP data.

In this step, get ready to:
  • Query the Census ACS 5-Year API for median income and population by ZCTA.
  • Clean and filter the Census data to NYC.
  • Merge your retailer counts with Census demographics on a shared geographic key.
Query the Census API for income and population

The Census ACS 5-Year API returns data for every ZCTA in the country. You need two variables: B19013_001E (median household income) and B01003_001E (total population).

  • Create a new cell in your notebook (Cell 7) and paste the following code:
# Cell 7: Query Census API for income and population by ZCTA
CENSUS_API_KEY = 'your-api-key-here'

census_url = (
    'https://api.census.gov/data/2024/acs/acs5'
    '?get=NAME,B19013_001E,B01003_001E'
    '&for=zip%20code%20tabulation%20area:*'
    f'&key={CENSUS_API_KEY}'
)
response = requests.get(census_url)
census_raw = response.json()

census_df = pd.DataFrame(census_raw[1:], columns=census_raw[0])
print(census_df.head())

What does this code do?

This cell builds a URL that asks the Census Bureau for income and population data for every ZCTA in the US. The response comes back as a JSON list of lists, where the first list is column headers and everything after is data.

You convert that raw JSON into a DataFrame so you can manipulate it with pandas. Replace 'your-api-key-here' with the actual API key you received from the Census Bureau in Step 1.

  • Replace 'your-api-key-here' with the Census API key you saved in Step 1.
  • Run the cell.

You should see a table with columns like NAME, B19013_001E, B01003_001E, and zip code tabulation area. The values are all strings at this point.

Getting a 204 or error response?

  • Double-check that you replaced the placeholder with your actual API key (no extra spaces or quotes around it).
  • Make sure the key is inside the single quotes: CENSUS_API_KEY = 'abc123...'.
  • If you get a KeyError on response.json(), print response.status_code to see what the Census server returned.

Help me debug my Census API call.

Filter Census data to NYC ZCTAs and clean

The API returns every ZCTA in the country. You need to rename the columns to something readable, convert strings to numbers, and filter down to NYC.

  • Create a new cell (Cell 8) and paste this code:
# Cell 8: Clean Census data
census_df = census_df.rename(columns={
    'B19013_001E': 'median_income',
    'B01003_001E': 'population',
    'zip code tabulation area': 'zcta'
})

census_df['median_income'] = pd.to_numeric(census_df['median_income'], errors='coerce')
census_df['population'] = pd.to_numeric(census_df['population'], errors='coerce')

census_df = census_df[census_df['median_income'] > 0]
census_df = census_df[census_df['population'] > 0]

nyc_zcta_mask = census_df['zcta'].str[:3].isin(
    ['100', '101', '102', '103', '104', '110', '111', '112', '113', '114', '116']
)
census_nyc = census_df[nyc_zcta_mask].copy()
print(f"NYC ZCTAs with valid income data: {len(census_nyc)}")

What does this code do?

First, you rename the cryptic Census variable codes to human-readable names. Then you convert the income and population columns from strings to numbers (the Census API returns everything as text).

The filter census_df['median_income'] > 0 drops ZCTAs where the Census suppressed the estimate (it uses special negative sentinel values like -666666666 for areas with too few households). Finally, you apply the same NYC ZIP prefix filter from Step 2.

  • Run the cell.

You should see a count of NYC ZCTAs with valid income data. This number tells you how many geographic units you have to work with for your analysis.

Seeing 0 NYC ZCTAs?

  • Check that Cell 7 ran successfully first. If census_df is empty, the rename and filter will produce nothing.
  • Verify the ZCTA column name matches what the Census returned. Print census_df.columns to see the actual names before renaming.

Help me figure out why my Census data filtering returned zero rows.

Merge SNAP counts with Census data

You now have two DataFrames: zip_counts (retailer counts by ZIP from your SNAP data) and census_nyc (income and population by ZCTA from the Census). The merge connects them on the shared geographic key so each row has retailer count, population, AND income.

What is an inner join?

An inner join keeps only rows that have a match in BOTH DataFrames. If a ZIP code appears in your SNAP data but not in Census data (or vice versa), that row gets dropped. This means you lose some ZIPs, but every remaining row has complete data for analysis.

The alternative (a left join) would keep all ZIPs from one side and fill the other columns with NaN. For this analysis, complete rows are more useful than maximum coverage.

  • Create a new cell (Cell 9) and paste this code:
# Cell 9: Merge SNAP counts with Census data
zip_counts[ZIP_COL] = zip_counts[ZIP_COL].astype(str).str.zfill(5)
census_nyc['zcta'] = census_nyc['zcta'].astype(str).str.zfill(5)

merged = census_nyc.merge(
    zip_counts,
    left_on='zcta',
    right_on=ZIP_COL,
    how='inner'
)
print(f"Merged rows: {len(merged)}")
print(f"Missing values:\n{merged.isna().sum()}")

What does this code do?

Before merging, you pad both ZIP columns to exactly 5 digits with str.zfill(5). This prevents mismatches where one side stores "10001" and the other stores "10001" with a leading space or different type.

The how='inner' parameter means you only keep ZCTAs that exist in both datasets. You will lose some rows here. That is a deliberate tradeoff: every row in merged has a retailer count, a population number, and an income estimate.

Before you run this cell, think: how many rows do you expect? You had about 170-200 unique NYC ZIPs with SNAP retailers, and a similar number of Census ZCTAs. Will the inner join keep all of them, or will some get dropped?

  • Run the cell.

You should see the number of merged rows (likely fewer than both input DataFrames) and a check of missing values across all columns. If a column shows zero NaN values, your inner join worked cleanly.

Seeing 0 merged rows?

  • The most common cause is mismatched types. Print zip_counts[ZIP_COL].dtype and census_nyc['zcta'].dtype to confirm both are strings.
  • Print a few values from each column side by side. Look for leading/trailing spaces or different digit counts.
  • Make sure you ran Cells 7 and 8 before Cell 9. The census_nyc variable must already exist.

Help me debug why my merge returned zero rows.

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook with all nine cells.

ⓧ I'd like to double check the full code

Your notebook should now have nine cells. Scroll through and confirm each is present in this order:

  • Cell 1: Imports (pandas, numpy, matplotlib, requests).
  • Cell 2: Load SNAP retailer data from CSV.
  • Cell 3: Assign column name variables (ZIP_COL, STATE_COL, TYPE_COL).
  • Cell 4: Filter to NYC using ZIP prefix and state masks.
  • Cell 5: Group by ZIP and count retailers.
  • Cell 6: Inspect store type distribution.
  • Cell 7: Query Census API for income and population by ZCTA.
  • Cell 8: Rename columns, convert to numeric, drop invalid rows, filter to NYC ZCTAs.
  • Cell 9: Pad ZIP codes to 5 digits and merge with an inner join.

You now have a single DataFrame where every row is a NYC ZCTA with a retailer count, a population number, and a median income estimate. Next, you will compute a per-capita metric and split by income level to answer your research question.

Compute Per-Capita Metric and Compare Income Groups

Your merged DataFrame now connects retailer counts to population and income data for each NYC ZIP code. But you cannot just compare raw retailer counts between ZIP codes.

A ZIP code with 80,000 residents will naturally have more stores than one with 15,000, regardless of income level. To answer your research question fairly, you need to normalize by population, then split by income and compare.

In this step, get ready to:
  • Compute a per-capita metric that accounts for population differences.
  • Split NYC ZCTAs into below-median and above-median income groups.
  • Produce a specific numeric finding that answers your research question.
Compute retailers per 1,000 residents

The formula is straightforward: divide the retailer count by the population, then multiply by 1,000. This gives you a rate that is comparable across ZIP codes of different sizes.

  • Create a new cell (Cell 10) in your notebook and add this code to compute the per-capita metric:
# Cell 10: Compute retailers per 1,000 residents
merged['retailers_per_1000'] = (merged['retailer_count'] / merged['population']) * 1000
print(merged[['zcta', 'retailer_count', 'population', 'median_income', 'retailers_per_1000']].head(10))

What does this code do?

The formula (retailer_count / population) * 1000 converts raw counts into a rate per 1,000 people. A ZIP with 50 retailers and 25,000 residents gets a score of 2.0. A ZIP with 50 retailers and 50,000 residents gets 1.0.

This makes the two ZIPs comparable, even though they have the same raw count. Without this normalization, you would incorrectly conclude they have equal access.

  • Run Cell 10. You should see a table with five columns: zcta, retailer_count, population, median_income, and retailers_per_1000. The new column should show decimal values, typically between 0.5 and 10.0.

Seeing infinity or NaN values?

If any rows have a population of zero, dividing by zero produces infinity. Check with merged[merged['population'] == 0]. If you find any, go back to Cell 8 and confirm your filter census_df[census_df['population'] > 0] is in place.

Help me debug NaN or infinity values in my retailers_per_1000 column

Classify ZCTAs by income

You now have a comparable rate for every ZIP code. The next question is: do lower-income areas have a different rate than higher-income areas? To answer this, you need a dividing line between "low" and "high" income.

You will use the median of all ZCTA median incomes as your split point. Every ZCTA below that value goes into one group. Every ZCTA at or above it goes into the other.

  • Create a new cell (Cell 11) and add this code to classify and compare:
# Cell 11: Classify by income and compare
citywide_median = merged['median_income'].median()
print(f"Citywide median of ZCTA median incomes: ${citywide_median:,.0f}")

merged['income_group'] = np.where(
    merged['median_income'] < citywide_median,
    'Below Median',
    'Above Median'
)

comparison = merged.groupby('income_group')['retailers_per_1000'].agg(['mean', 'median', 'count'])
print(comparison)

What does this code do?

  • citywide_median computes the median of all ZCTA median household incomes. This is your split point between "low" and "high" income areas.
  • np.where() creates a new column that labels each ZCTA as either 'Below Median' or 'Above Median' based on whether its income falls below or at/above the citywide median.
  • .agg(['mean', 'median', 'count']) computes three summary statistics for each income group: the average retailer density, the median density, and how many ZCTAs are in each group.
  • Run Cell 11. You should see the citywide median income printed, followed by a comparison table with two rows ('Above Median' and 'Below Median') and three columns ('mean', 'median', 'count').
  • Check the 'count' column. Each group should have at least 5 ZCTAs. If one group has fewer than 5, something went wrong in your filtering or merge.

Why use the median as the split point?

The median splits your ZCTAs into two roughly equal-sized groups. This avoids the problem of picking an arbitrary income threshold (like $50,000) that might put most ZCTAs on one side. With a median split, you always get a balanced comparison.

One group has fewer than 5 ZCTAs?

This usually means your merge in Cell 9 lost too many rows. Run merged.shape to check your total row count. If you have fewer than 15 rows total, revisit your ZIP filtering logic in Cell 4 and Cell 8 to confirm the ZIP prefix lists match.

Help me debug why my income groups have too few ZCTAs

Assess your finding

The comparison table gives you a quick picture, but you need crisp, quotable numbers. This cell isolates the two group means and computes the percentage difference between them.

  • Create a new cell (Cell 12) and add this code to compute your final numbers:
# Cell 12: State your finding
below = merged[merged['income_group'] == 'Below Median']['retailers_per_1000'].mean()
above = merged[merged['income_group'] == 'Above Median']['retailers_per_1000'].mean()
pct_diff = ((below - above) / above) * 100

print(f"\nBelow-median income ZCTAs: {below:.2f} SNAP retailers per 1,000 residents")
print(f"Above-median income ZCTAs: {above:.2f} SNAP retailers per 1,000 residents")
print(f"Difference: {pct_diff:+.1f}%")

What does this code do?

  • below and above extract the mean retailers-per-1000 for each income group separately.
  • pct_diff computes the percentage difference relative to the above-median group. A positive value means below-median ZCTAs have MORE retailers per capita. A negative value means they have FEWER.

Before you run this cell, think about what you expect. Do you think lower-income ZIP codes will have more or fewer SNAP retailers per 1,000 residents than higher-income ones?

  • Run Cell 12. You should see two specific numbers and a percentage difference printed.

Look at your output. You now have a concrete, numeric answer to your research question. Was the direction what you expected?

  • Create a new Markdown cell below Cell 12. Write a one-sentence finding using your actual numbers. For example: "NYC ZIP codes with below-median household income have X.XX SNAP retailers per 1,000 residents compared to Y.YY in above-median ZIP codes, a Z% difference."

This is your resume bullet point

A specific, quantified finding is far more powerful than "analyzed SNAP data." In an interview, you can say: "I found that low-income NYC ZIP codes have X% more/fewer SNAP retailers per capita than high-income ones." That specificity shows you can extract insights, not just run code.

Both groups show identical values?

If below and above are identical or pct_diff is exactly 0.0%, check that your income_group column actually has both values. Run merged['income_group'].value_counts() to confirm both groups are populated.

Help me debug why my income group comparison shows no difference

✔️ Awesome, I've got everything!

Great. Make sure you have saved your notebook with all cells through Cell 12 plus your Markdown finding.

ⓧ I'd like to double check the full code

Scroll through your notebook and confirm you have these three new cells in order, each producing output when run:

  • Cell 10 computes retailers_per_1000 and prints a 10-row preview with five columns.
  • Cell 11 computes citywide_median, creates the income_group column, and prints a comparison table with mean, median, and count.
  • Cell 12 extracts below, above, and pct_diff, then prints the two group means and percentage difference.
  • A Markdown cell below Cell 12 with your one-sentence finding using real numbers from your output.

You now have a quantified answer to your research question. Numbers in a print statement are convincing to you, but not to anyone else. Next up, you will turn your finding into clear visualizations that make the pattern obvious at a glance.

Visualize Your Finding

You have a specific numeric finding about SNAP retailer density across income groups. But numbers in a print() statement are not convincing to anyone who glances at your work.

A well-labeled chart lets someone see the pattern instantly. In this step, you will learn matplotlib basics by turning your numeric finding into two publication-ready visualizations.

In this step, get ready to:
  • Create a scatter plot showing income vs. retailer density for every ZCTA.
  • Build a bar chart comparing the two income group means.
  • Save your visualizations and write a short interpretation of what you see.
Create a scatter plot of income vs. retailer density

A scatter plot shows the relationship between two continuous variables. Each dot represents one ZCTA, letting you see whether the pattern is clean or noisy.

  • Create a new code cell in your notebook (Cell 13) and add the following scatter plot code:
# Cell 13: Scatter plot
fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(merged['median_income'], merged['retailers_per_1000'], alpha=0.6, edgecolors='black', linewidths=0.5)
ax.axvline(x=citywide_median, color='red', linestyle='--', label=f'Median Income (${citywide_median:,.0f})')
ax.set_xlabel('Median Household Income ($)')
ax.set_ylabel('SNAP Retailers per 1,000 Residents')
ax.set_title('SNAP Retailer Density vs. Income in NYC ZIP Codes')
ax.legend()
plt.tight_layout()
plt.savefig('output/scatter_income_vs_density.png', dpi=150)
plt.show()

What does this code do?

  • plt.subplots(figsize=(10, 6)) creates a figure that is 10 inches wide and 6 inches tall.
  • ax.scatter() plots each ZCTA as a semi-transparent dot with a thin black edge so overlapping points are visible.
  • ax.axvline() draws a red dashed vertical line at the citywide median income. This is your visual split point between the two income groups.
  • plt.savefig() saves the figure to your output/ folder before plt.show() displays it inline.
  • Run the cell.

You should see a scatter plot with dots spread across the chart. A red dashed vertical line marks the median income split point. Dots to the left of the line are your below-median ZCTAs. Dots to the right are above-median.

Seeing a blank plot or no output?

  • Make sure you ran all previous cells (1-12) first. The variables merged and citywide_median must already exist in memory.
  • If you see a FileNotFoundError on the savefig line, confirm your output/ folder exists in the same directory as your notebook.

Help me debug my matplotlib scatter plot

Create a bar chart comparing group means

The scatter shows every data point, which is great for seeing the spread. A bar chart distills the comparison into two clean numbers that anyone can read at a glance.

  • Create a new code cell (Cell 14) and add this bar chart code:
# Cell 14: Bar chart
fig, ax = plt.subplots(figsize=(8, 5))
groups = ['Below Median', 'Above Median']
means = [below, above]
bars = ax.bar(groups, means, color=['#d62728', '#2ca02c'], edgecolor='black')

for bar, val in zip(bars, means):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.05,
            f'{val:.2f}', ha='center', va='bottom', fontweight='bold')

ax.set_ylabel('SNAP Retailers per 1,000 Residents')
ax.set_title('SNAP Retailer Density by Income Group (NYC)')
ax.set_ylim(0, max(means) * 1.3)
plt.tight_layout()
plt.savefig('output/bar_income_comparison.png', dpi=150)
plt.show()

What does this code do?

  • ax.bar() draws two bars. Red for below-median, green for above-median income ZCTAs.
  • The for loop places the exact numeric value above each bar so the reader never has to guess heights from the y-axis.
  • ax.set_ylim(0, max(means) * 1.3) adds 30% headroom above the tallest bar so the text labels are not cut off.
  • Run the cell.

You should see two bars side by side with bold numeric labels above each one. The difference between the two values is your finding, visible at a glance.

Numeric labels overlapping or cut off?

  • If labels are clipped at the top, increase the multiplier in set_ylim from 1.3 to 1.5 and re-run.
  • If you see a NameError for below or above, make sure Cell 12 ran successfully first.

Help me fix my bar chart labels

Verify your saved files and write an interpretation

Both cells saved PNG files to your output/ folder automatically. Time to confirm they landed and reflect on what the visualizations tell you.

  • Before you check, predict: how many files do you expect to see in output/?
  • Confirm your saved images by running this in a new cell:
import os
print(os.listdir('output/'))

You should see scatter_income_vs_density.png and bar_income_comparison.png listed in the output.

  • Now add a new Markdown cell below your last code cell. To do this, click the + Markdown button in the notebook toolbar (or click the cell type dropdown and select Markdown).
  • Write 2-3 sentences interpreting what you see in your two charts. Consider these questions:
  • Does the scatter plot show a clean linear relationship, or is it noisy with lots of spread?
  • Are there outlier ZCTAs that break the pattern? (Dots far from the cluster.)
  • What might explain those outliers? (A tiny ZCTA with one store, or a commercial district with few residents?)

Why write an interpretation?

Charts without interpretation are decoration. An interviewer will ask you what the chart means, not just how you made it. Practicing that narration now means you can explain your work later without scrambling.

  • Run the Markdown cell (Shift+Enter) to render it.

Output folder is empty?

  • If os.listdir('output/') returns an empty list, scroll back up to Cells 13 and 14 and re-run them. The plt.savefig() call only executes when the cell runs.
  • Make sure your output/ folder sits in the same directory as your notebook file, not inside data/.

Help me figure out why my plots are not saving

✔️ Awesome, I've got everything!

Great. Make sure you have saved your notebook with all cells intact.

ⓧ I'd like to double check the full code

Your notebook now has 14 code cells plus one Markdown interpretation cell. Scroll through and confirm these two new cells are present in order:

  • Cell 13 starts with # Cell 13: Scatter plot and ends with plt.show().
  • Cell 14 starts with # Cell 14: Bar chart and ends with plt.show().
  • A Markdown cell with your 2-3 sentence interpretation follows Cell 14.

You now have two presentation-ready charts and a written interpretation of your finding. Your analysis is complete. Next up, you will push further by breaking down the data by store type in the Secret Mission.

Secret mission

Disaggregate by Store Type

Your overall finding treats all SNAP retailers the same. A Whole Foods and a gas station both count as one retailer. What if you break the analysis apart by store type? Low-income areas might have MORE total SNAP retailers but FEWER supermarkets. Find out whether the headline number masks a more interesting story underneath.

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely locally, so there are no ongoing costs. Your Census API key is free and does not expire. Decide whether to keep your files for future reference or delete them entirely.

Resources you used:

  • analysis.ipynb notebook with all analysis cells.
  • data/snap_retailers.csv downloaded from USDA FNS.
  • output/ folder containing generated PNG charts.
  • Python libraries: pandas, numpy, matplotlib, requests installed via pip.

Keep everything running

No action needed. Choose this if you plan to revisit your analysis, extend it with new data, or use it as a reference for future projects.

  • Your notebook and data files stay where they are.
  • Your Census API key remains active and free.
  • The installed Python libraries remain available for other projects.

Pause - I'll come back to this later

Close VS Code to free up memory. Your files and libraries stay on disk so you can pick up where you left off.

  • Close the VS Code window with your project.
  • Your notebook, data, and output files remain in your project folder.
  • When you return, open the project folder in VS Code and your notebook will be ready.

Delete - I don't want to use this again

Remove all project files and optionally uninstall the Python libraries.

  • Close VS Code.
  • Delete your project folder (the one containing analysis.ipynb) by running the command below for your operating system.

On macOS/Linux, open Terminal and run:

cd ~/Desktop
rm -rf snap-analysis

On Windows, open PowerShell and run:

cd ~\Desktop
Remove-Item -Recurse -Force snap-analysis

Replace snap-analysis with the actual name of your project folder if you named it differently.

  • Optionally, uninstall the Python libraries by running:
pip uninstall pandas numpy matplotlib requests -y

Should I uninstall the libraries?

Only uninstall if you are sure you will not use pandas, numpy, matplotlib, or requests in other projects. These are common data science libraries you will likely need again.

Nice Work!

Nice Work!

You just built a complete data analysis from scratch using real government data. You can now open your notebook, walk someone through every cell, and state a specific numeric finding about food access in NYC backed by your own work.

You've learned how to:

  • Load real government data from two different sources (USDA FNS and the Census Bureau) and join them on a shared geographic key using pandas merge operations.
  • Compute a per-capita normalized metric that avoids the naive "more people = more stores" fallacy, and split your data by income level to produce a specific numeric finding.
  • Create publication-ready matplotlib visualizations including scatter plots and bar charts with proper labels, titles, and annotations.
  • Secret Mission: You disaggregated your analysis by store type and discovered whether the overall pattern holds or breaks apart by category.

Ready to quiz yourself?