Geography

School district, county, and college-location crosswalks used by every later step

html`<style>
  .cell-output-display svg text { font-size: 15px !important; }
  .cell-output-display table { font-size: 16px; }
</style>`

Overview

This script builds the geographic crosswalk that every later step in this project joins against.

• Objective: tie every Census block in the Wasatch Front region (Box Elder, Davis, Weber, Salt Lake, and Utah counties) to the school district and county it falls in, and to the K-12 school(s) and college(s) located in it — so that PopulationSim’s synthetic population (3_synthetic.qmd) and real enrollment data (2_enrollment.qmd) can both be rolled up to district/county level as exact joins, with no area-weighting or estimation involved.
• What it pulls: county, school district, Travel Analysis Zone (TAZ), and Census block boundaries; K-12 school and college locations from both UGRC/SGID (the state-preferred source) and NCES (kept alongside, not as a replacement, so schools and colleges can also be joined back to federal CCD/IPEDS data by ID); and the WFRC/MAG travel demand model boundary, used later to handle Box Elder County’s partial coverage in that model. TAZ boundaries come from WFRC’s WasatchFrontTAZ feature service rather than the TDM-GEO-AnalysisZones repo directly, so the pull stays current without tracking that repo’s releases — TAZID is what is kept as a reference layer. The synthetic join uses MAZ, not TAZ. MAZ boundaries come from a draft microzone layer supplied for Utah (_data/raw/geometry/microzones_draft_*.zip), not yet a published WFRC layer. Block boundaries are still cached and checked, but no crosswalk uses them.
• Output: two crosswalks, one per use. _data/intermediate/crosswalk_school.csv maps each K-12 school (UGRC, NCES public, NCES private) and each MAZ to its district. _data/intermediate/crosswalk_college.csv maps each college (UGRC, NCES) and each MAZ to its county. MAZ rows appear in both files, so PopulationSim output rolls up through MAZ rather than through blocks. Section 10 checks coverage and maps the study area (Figure 1).

Setup

import sqlalchemy
import pandas as pd
import geopandas as gpd
import requests
from pygris import blocks, counties
from pygris.helpers import validate_state, validate_county
from pathlib import Path

paths = {
    "raw": Path("_data/raw/"),
    "intermediate": Path("_data/intermediate/"),
    "output": Path("_outputs/"),
}

constants = {
    "state": "UT",
    "counties": {"Box Elder", "Davis", "Weber", "Salt Lake", "Utah"},
    "crs": "EPSG:26912",
    "data_year": 2023,  # NCES/EDGE school year
}

# Public, read-only credential documented by UGRC for the Open SGID: https://gis.utah.gov/sgid/open-sgid/
OPENSGID_URI = "postgresql+psycopg://agrc:agrc@opensgid.ugrc.utah.gov:5432/opensgid"
opensgid_engine = sqlalchemy.create_engine(OPENSGID_URI)

# Derive FIPS codes via pygris instead of hardcoding them
state_fips = validate_state(constants["state"])
county_fips = [
    validate_county(state_fips, county)
    for county in constants["counties"]
]
fips_codes = [f"{state_fips}{county}" for county in county_fips]
fips_codes

Using FIPS code '49' for input 'UT'
Using FIPS code '057' for input 'Weber'
Using FIPS code '035' for input 'Salt Lake'
Using FIPS code '011' for input 'Davis'
Using FIPS code '003' for input 'Box Elder'
Using FIPS code '049' for input 'Utah'

['49057', '49035', '49011', '49003', '49049']

# Derive the NCES EDGE/CCD year-code (2023 -> "2324") so a future vintage swap below is a single
# constants change, not a hunt through every URL in this file.
year_code = f"{str(constants['data_year'])[-2:]}{str(constants['data_year'] + 1)[-2:]}"
year_code

'2324'

1. County boundaries

county_boundary_path = paths["raw"] / "geometry" / "county_boundary.gpkg"

# 1. Download only if not already cached locally
if not county_boundary_path.exists():
    all_counties = counties(state=state_fips, year=2023, cache=True).to_crs(constants["crs"])
    all_counties.to_file(county_boundary_path)

# 2. Always read from local storage -- full state, not just our 5 counties. Every later section
#    tags its schools/colleges/districts against all_counties so those statewide .gpkg caches carry
#    a correct county_fips everywhere, not just inside our study area; bound_counties (below) is the
#    5-county-filtered version used to restrict the crosswalk itself.
all_counties = gpd.read_file(county_boundary_path)

county_list_sql = ", ".join(f"'{fips}'" for fips in county_fips)

bound_counties = gpd.read_file(
    county_boundary_path, sql=f"SELECT * FROM county_boundary WHERE COUNTYFP IN ({county_list_sql})"
)

bound_counties.explore()

2. School district boundaries

district_boundary_path = paths["raw"] / "geometry" / "school_district_boundary.gpkg"

# 1. Download only if not already cached locally
if not district_boundary_path.exists():
    all_districts = gpd.read_postgis(
        "select * from boundaries.school_district_boundaries", opensgid_engine, geom_col="shape"
    ).to_crs(constants["crs"])

    all_districts["county_fips"] = gpd.sjoin(
        all_districts.set_geometry(all_districts.centroid),
        all_counties[["COUNTYFP", "geometry"]],
        predicate="within",
        how="left",
    )["COUNTYFP"]

    all_districts.to_file(district_boundary_path)

# 2. Read + filter in one step, pushed down to SQLite
county_list_sql = ", ".join(f"'{fips}'" for fips in county_fips)

bound_districts = gpd.read_file(
    district_boundary_path,
    sql=f"SELECT * FROM school_district_boundary WHERE county_fips IN ({county_list_sql})",
)

bound_districts.explore()

3. Travel Analysis Zone (TAZ) boundaries

taz_boundary_path = paths["raw"] / "geometry" / "taz_boundary.gpkg"

# 1. Download only if not already cached locally
if not taz_boundary_path.exists():
    taz_url = (
        "https://services1.arcgis.com/taguadKoI1XFwivx/ArcGIS/rest/services/"
        "WasatchFrontTAZ/FeatureServer/0/query"
    )
    taz_pages = []
    offset = 0
    while True:
        page = gpd.read_file(
            f"{taz_url}?where=1=1&outFields=*&f=geojson"
            f"&orderByFields=FID&resultOffset={offset}&resultRecordCount=2000"
        )
        if page.empty:
            break
        taz_pages.append(page)
        offset += len(page)
        if len(page) < 2000:
            break
    all_taz = pd.concat(taz_pages, ignore_index=True).to_crs(constants["crs"])

    count_resp = requests.get(f"{taz_url}?where=1=1&returnCountOnly=true&f=json").json()
    assert len(all_taz) == count_resp["count"], (len(all_taz), count_resp["count"])

    all_taz["county_fips"] = all_taz["CO_FIPS"].astype(int).astype(str).str.zfill(3)

    all_taz.to_file(taz_boundary_path)

# 2. Read + filter in one step, pushed down to SQLite
county_list_sql = ", ".join(f"'{fips}'" for fips in county_fips)

bound_taz = gpd.read_file(
    taz_boundary_path, sql=f"SELECT * FROM taz_boundary WHERE county_fips IN ({county_list_sql})"
)

bound_taz.explore()

4. Micro Analysis Zone (MAZ) boundaries

maz_boundary_path = paths["raw"] / "geometry" / "maz_boundary.gpkg"

# 1. Build the cache once from the local zip (no download -- it's a draft Andy supplied)
if not maz_boundary_path.exists():
    all_maz = gpd.read_file(
        "zip://_data/raw/geometry/microzones_draft_20260610AndyaddedmazidforUtah.zip"
    ).to_crs(constants["crs"])
    all_maz = all_maz[all_maz.geometry.notna()].copy()
    all_maz = all_maz.set_geometry(all_maz.geometry.buffer(0))

    all_maz["county_fips"] = gpd.sjoin(
        all_maz.set_geometry(all_maz.centroid),
        all_counties[["COUNTYFP", "geometry"]],
        predicate="within",
        how="left",
    )["COUNTYFP"]

    all_maz.to_file(maz_boundary_path)

# 2. Read + filter in one step, pushed down to SQLite
bound_maz = gpd.read_file(
    maz_boundary_path, sql=f"SELECT * FROM maz_boundary WHERE county_fips IN ({county_list_sql})"
)

bound_maz.explore()

5. Census blocks

block_boundary_path = paths["raw"] / "geometry" / "block_boundary.gpkg"
block_counties_path = paths["raw"] / "geometry" / "block_boundary.counties.txt"

# Auto-invalidate the cache if constants["counties"] has changed since it was built -- pygris'
# blocks() call below is scoped to the CURRENT county list at fetch time (unlike every other
# cached layer in this file, which fetches statewide and filters at read-time), so a stale cache
# here would silently miss any newly-added county's blocks instead of erroring.
current_counties = sorted(county_fips)
cached_counties = (
    block_counties_path.read_text().splitlines() if block_counties_path.exists() else []
)
if block_boundary_path.exists() and current_counties != cached_counties:
    block_boundary_path.unlink()

# 1. Download only if not already cached locally
if not block_boundary_path.exists():
    all_blocks = blocks(state=state_fips, county=county_fips, year=2023, cache=True).to_crs(
        constants["crs"]
    )

    all_blocks["district_name"] = gpd.sjoin(
        all_blocks.set_geometry(all_blocks.centroid),
        bound_districts[["name", "geometry"]],
        predicate="within",
        how="left",
    )["name"]

    all_blocks.to_file(block_boundary_path)
    block_counties_path.write_text("\n".join(current_counties))

# 2. Read + filter in one step, pushed down to SQLite
county_list_sql = ", ".join(f"'{fips}'" for fips in county_fips)

bound_blocks = gpd.read_file(
    block_boundary_path, sql=f"SELECT * FROM block_boundary  WHERE COUNTYFP20 IN({county_list_sql})"
)

bound_blocks.explore()

6. Pre-K to 12 school locations

UGRC Pre-K to 12

ugrc_school_path = paths["raw"] / "geometry" / "ugrc_school.gpkg"

# Current Open SGID layer. Kept for reference only -- the analysis uses the BY2023 snapshot below.
# 1. Download only if not already cached locally
if not ugrc_school_path.exists():
    all_ugrc_schools = gpd.read_postgis(
        "select * from society.schools_prek_to_12", opensgid_engine, geom_col="shape"
    ).to_crs(constants["crs"])

    all_ugrc_schools["county_fips"] = gpd.sjoin(
        all_ugrc_schools, all_counties[["COUNTYFP", "geometry"]], predicate="within", how="left"
    )["COUNTYFP"]

    all_ugrc_schools.to_file(ugrc_school_path)

# 2. Read + filter in one step, pushed down to SQLite
county_list_sql = ", ".join(f"'{fips}'" for fips in county_fips)

ugrc_schools_current = gpd.read_file(
    ugrc_school_path,
    sql=f"SELECT * FROM ugrc_school WHERE county_fips IN ({county_list_sql})",
)

UGRC Pre-K to 12, BY2023 snapshot (used)

# Copy of UGRC Data/BY2023 from https://github.com/WFRCAnalytics/TDM-INP-K-12-Enrollment, the
# snapshot behind the model's K-12 enrollment input. Its SchoolYear field is 2025 (pulled March 2025).
ugrc_school_by2023_shp = paths["raw"] / "ugrc_by2023" / "Schools_PreKto12.shp"
ugrc_school_by2023_path = paths["intermediate"] / "ugrc_school_by2023.gpkg"

# Shapefile names are cut to 10 characters, so restore the Open SGID names that downstream code uses
truncated_names = {
    "SchoolNumb": "schoolnumber",
    "CharterSch": "charterschool",
    "PrivateSch": "privateschool",
    "OnlineScho": "onlineschool",
    "SchoolLeve": "schoollevel",
    "FederalLoc": "federallocalecode",
    "PhoneNumbe": "phonenumber",
    "AmericanIn": "americanindian",
    "MultipleRa": "multiplerace",
    "PacificIsl": "pacificislander",
    "Economical": "economicallydisadvantaged",
    "EnglishLea": "englishlearner",
    "StudentDis": "studentdisability",
}
all_ugrc_schools_by2023 = (
    gpd.read_file(ugrc_school_by2023_shp)
    .rename(columns=truncated_names)
    .rename(columns=str.lower)
    .to_crs(constants["crs"])
)
assert set(all_ugrc_schools_by2023.columns) <= set(ugrc_schools_current.columns), (
    "BY2023 columns do not match the Open SGID schema"
)
all_ugrc_schools_by2023["county_fips"] = gpd.sjoin(
    all_ugrc_schools_by2023, all_counties[["COUNTYFP", "geometry"]], predicate="within", how="left"
)["COUNTYFP"]
all_ugrc_schools_by2023.to_file(ugrc_school_by2023_path)

ugrc_schools = all_ugrc_schools_by2023[all_ugrc_schools_by2023["county_fips"].isin(county_fips)]

ugrc_schools.explore()

NCES public schools

nces_public_school_path = paths["raw"] / "geometry" / "nces_public_school.gpkg"

if not nces_public_school_path.exists():
    all_nces_public_schools = gpd.read_file(
        "https://nces.ed.gov/opengis/rest/services/K12_School_Locations/"
        f"EDGE_ADMINDATA_PUBLICSCH_{year_code}/MapServer/1/query?where=LSTATE='UT'&outFields=*&f=geojson"
    ).to_crs(constants["crs"])

    all_nces_public_schools["county_fips"] = gpd.sjoin(
        all_nces_public_schools,
        all_counties[["COUNTYFP", "geometry"]],
        predicate="within",
        how="left",
    )["COUNTYFP"]

    all_nces_public_schools.to_file(nces_public_school_path)

nces_public_schools = gpd.read_file(
    nces_public_school_path,
    sql=f"SELECT * FROM nces_public_school WHERE county_fips IN ({county_list_sql})",
)

nces_public_schools.explore()

# TOTAL is documented by NCES as "Total Students (Includes AE)" -- verified 2026-10-01 this has no
# effect on current Utah data (named adult-ed CCD entries all report 0). Re-check if this ever
# trips, since it would mean enrollment figures downstream include adult-ed headcounts alongside K-12.
adult_ed_total = nces_public_schools.query(
    "SCH_NAME.str.contains('Adult', case=False, na=False)"
)["TOTAL"].sum()
assert adult_ed_total == 0, "Adult-ed schools now report nonzero TOTAL -- revisit NCES enrollment sourcing"

NCES private school

nces_private_school_path = paths["raw"] / "geometry" / "nces_private_school.gpkg"

if not nces_private_school_path.exists():
    # EDGE_GEOCODE_PRIVATESCH is only published every other school year (confirmed 2026-10-01:
    # years 15_16, 1718, 1920, 2122, 2324 exist; 1819/2021/2223 don't) -- step back one cycle if
    # this year's layer isn't published yet.
    private_year_code = year_code
    private_url = (
        "https://nces.ed.gov/opengis/rest/services/K12_School_Locations/"
        f"EDGE_GEOCODE_PRIVATESCH_{private_year_code}/MapServer/0/query?where=STATE='UT'&outFields=*&f=geojson"
    )

    if not requests.get(private_url).json().get("features"):
        fallback_year = constants["data_year"] - 1
        private_year_code = f"{str(fallback_year)[-2:]}{str(fallback_year + 1)[-2:]}"
        print(f"NOTE: {year_code} private-school layer not published, using {private_year_code}")
        private_url = (
            "https://nces.ed.gov/opengis/rest/services/K12_School_Locations/"
            f"EDGE_GEOCODE_PRIVATESCH_{private_year_code}/MapServer/0/query?where=STATE='UT'&outFields=*&f=geojson"
        )

    all_nces_private_schools = gpd.read_file(private_url).to_crs(constants["crs"])

    all_nces_private_schools["county_fips"] = gpd.sjoin(
        all_nces_private_schools,
        all_counties[["COUNTYFP", "geometry"]],
        predicate="within",
        how="left",
    )["COUNTYFP"]

    all_nces_private_schools.to_file(nces_private_school_path)

nces_private_schools = gpd.read_file(
    nces_private_school_path,
    sql=f"SELECT * FROM nces_private_school WHERE county_fips IN ({county_list_sql})",
)

nces_private_schools.explore()

7. College campus locations

UGRC higher education

ugrc_college_path = paths["raw"] / "geometry" / "ugrc_college.gpkg"

# 1. Download only if not already cached locally
if not ugrc_college_path.exists():
    all_ugrc_colleges = gpd.read_postgis(
        "select * from society.schools_higher_education", opensgid_engine, geom_col="shape"
    ).to_crs(constants["crs"])

    all_ugrc_colleges["county_fips"] = gpd.sjoin(
        all_ugrc_colleges, all_counties[["COUNTYFP", "geometry"]], predicate="within", how="left"
    )["COUNTYFP"]

    all_ugrc_colleges.to_file(ugrc_college_path)

# 2. Read + filter in one step, pushed down to SQLite
county_list_sql = ", ".join(f"'{fips}'" for fips in county_fips)

ugrc_colleges = gpd.read_file(
    ugrc_college_path,
    sql=f"SELECT * FROM ugrc_college WHERE county_fips IN ({county_list_sql})",
)

ugrc_colleges.explore()

NCES postsecondary (colleges)

nces_college_path = paths["raw"] / "geometry" / "nces_college.gpkg"

if not nces_college_path.exists():
    all_nces_colleges = gpd.read_file(
        "https://nces.ed.gov/opengis/rest/services/Postsecondary_School_Locations/"
        f"EDGE_GEOCODE_POSTSECONDARYSCH_{year_code}/MapServer/0/query?where=STATE='UT'&outFields=*&f=geojson"
    ).to_crs(constants["crs"])

    all_nces_colleges["county_fips"] = gpd.sjoin(
        all_nces_colleges, all_counties[["COUNTYFP", "geometry"]], predicate="within", how="left"
    )["COUNTYFP"]

    all_nces_colleges.to_file(nces_college_path)

nces_colleges = gpd.read_file(
    nces_college_path,
    sql=f"SELECT * FROM nces_college WHERE county_fips IN ({county_list_sql})",
)

nces_colleges.explore()

8. Model / MPO region boundary

model_boundary_path = paths["raw"] / "geometry" / "model_boundary.gpkg"

# 1. Download only if not already cached locally
if not model_boundary_path.exists():
    model_boundary = gpd.read_file(
        "https://services1.arcgis.com/taguadKoI1XFwivx/ArcGIS/rest/services/"
        "WFRC_Administrative_and_Planning_Area_Boundaries/FeatureServer/4/query"
        "?where=1=1&outFields=*&f=geojson"
    ).to_crs(constants["crs"])
    model_boundary.to_file(model_boundary_path)

# 2. Always read from local storage
bound_model = gpd.read_file(model_boundary_path)

bound_model.explore()

9. Generate crosswalks

district_polys = bound_districts[["name", "geometry"]]
county_polys = bound_counties[["COUNTYFP", "geometry"]]

# Point sources join directly to polygons, so no blocks are involved
point_sources = [
    ("school", "UGRC", ugrc_schools, "schoolid"),
    ("school", "NCES public", nces_public_schools, "NCESSCH"),
    ("school", "NCES private", nces_private_schools, "PPIN"),
    ("college", "UGRC", ugrc_colleges, "xid"),
    ("college", "NCES", nces_colleges, "UNITID"),
]

rows = []
for unit, source, pts, id_col in point_sources:
    point_district = gpd.sjoin(pts[[id_col, "geometry"]], district_polys, predicate="within", how="left")
    point_county = gpd.sjoin(pts[[id_col, "geometry"]], county_polys, predicate="within", how="left")
    assert len(point_district) == len(pts) == len(point_county), f"{source} {unit}: duplicate matches"
    rows.append(pd.DataFrame({
        "Unit": unit,
        "Source": source,
        "ID": pts[id_col].astype(str).values,
        "District Name": point_district["name"].values,
        "COUNTYFP": point_county["COUNTYFP"].values,
    }))

# MAZ rows: district and county from the MAZ centroid, so MAZs with no blocks are covered
maz_points = bound_maz.set_geometry(bound_maz.centroid)[["maz_id", "geometry"]]
maz_district = gpd.sjoin(maz_points, district_polys, predicate="within", how="left")
maz_county = gpd.sjoin(maz_points, county_polys, predicate="within", how="left")
assert len(maz_district) == len(maz_points) == len(maz_county), "MAZ centroid matches more than one polygon"
rows.append(pd.DataFrame({
    "Unit": "maz",
    "Source": "PopulationSim",
    "ID": maz_points["maz_id"].astype(str).values,
    "District Name": maz_district["name"].values,
    "COUNTYFP": maz_county["COUNTYFP"].values,
}))

long = pd.concat(rows, ignore_index=True)
school_crosswalk = long[long["Unit"].isin(["school", "maz"])].copy()
college_crosswalk = long[long["Unit"].isin(["college", "maz"])].copy()

for name, df in [("school", school_crosswalk), ("college", college_crosswalk)]:
    print(name, "rows:", len(df), "duplicate (Source, ID):", df.duplicated(["Source", "ID"]).sum())

school_crosswalk.to_csv(paths["intermediate"] / "crosswalk_school.csv", index=False)
college_crosswalk.to_csv(paths["intermediate"] / "crosswalk_college.csv", index=False)

school rows: 26331 duplicate (Source, ID): 11
college rows: 24620 duplicate (Source, ID): 0

# Not strictly needed for a one-shot Quarto render (the process exit reclaims the connection
# either way), but matters in an interactive session where re-running Setup without restarting
# the kernel would otherwise orphan the previous engine's pooled connection to Open SGID.
opensgid_engine.dispose()

10. Check the crosswalk

Every later join depends on these match rates, so coverage is checked here before anything downstream uses the crosswalk.

The map is drawn in OJS, and the coverage summary is printed by Python. Python first writes the map layers as GeoJSON (WGS84, which is what browsers expect) so OJS can read them.

# OJS reads GeoJSON, not GeoDataFrames, so export each map layer in WGS84 for the {ojs} map below
map_layers = {
    "counties": bound_counties,
    "districts": bound_districts,
    "schools": ugrc_schools,
    "colleges": ugrc_colleges,
}
for layer_name, layer in map_layers.items():
    layer.to_crs("EPSG:4326").to_file(
        paths["intermediate"] / f"map_{layer_name}.geojson", driver="GeoJSON"
    )

study_counties = FileAttachment("_data/intermediate/map_counties.geojson").json()
study_districts = FileAttachment("_data/intermediate/map_districts.geojson").json()
study_schools = FileAttachment("_data/intermediate/map_schools.geojson").json()
study_colleges = FileAttachment("_data/intermediate/map_colleges.geojson").json()

Plot.plot({
  width: 640,
  height: 720,
  projection: "mercator",
  style: {background: "transparent", color: "inherit"},
  marks: [
    Plot.geo(study_districts, {fill: "#eeeeee", stroke: "#bbbbbb", strokeWidth: 0.5}),
    Plot.geo(study_counties, {stroke: "black", strokeWidth: 1.2}),
    Plot.geo(study_schools, {r: 1.8, fill: "#1f78b4", stroke: "none"}),
    Plot.geo(study_colleges, {r: 4, fill: "#e31a1c", stroke: "none"})
  ]
})

(a)

(b)

(c)

(d)

(e)

Figure 1: Study area. Gray fill is school district boundaries, black lines are county boundaries, blue dots are K-12 schools (UGRC BY2023 snapshot, SchoolYear 2025) and red dots are colleges (UGRC Open SGID).

for name in ["school", "college"]:
    df = pd.read_csv(paths["intermediate"] / f"crosswalk_{name}.csv", dtype=str)
    cov = df.groupby(["Unit", "Source"]).agg(
        rows=("ID", "size"),
        with_district=("District Name", "count"),
        with_county=("COUNTYFP", "count"),
    )
    print(f"--- crosswalk_{name} ---")
    print(cov.to_string())

--- crosswalk_school ---
                       rows  with_district  with_county
Unit   Source                                          
maz    PopulationSim  24531          24531        24531
school NCES private     105            105          105
       NCES public      742            742          742
       UGRC             953            953          953
--- crosswalk_college ---
                        rows  with_district  with_county
Unit    Source                                          
college NCES              47             47           47
        UGRC              42             42           42
maz     PopulationSim  24531          24531        24531

Crosswalk coverage by unit and source. Rows are school and college points plus MAZs; MAZ rows appear in both files.