diff --git a/changelog.d/us-ct-planning-region-cbsa.fixed.md b/changelog.d/us-ct-planning-region-cbsa.fixed.md new file mode 100644 index 0000000000..36005b116d --- /dev/null +++ b/changelog.d/us-ct-planning-region-cbsa.fixed.md @@ -0,0 +1 @@ +Connecticut blocks in the US block ladder now carry their CBSA. The OMB 2023 delineations list Connecticut by its nine 2022 planning regions, while 2020 blocks carry the eight old counties, so every CT block had `cbsa_code` 0. The builder now maps each CT block to its planning region through its town, using the 2020 BAF county-subdivision layer and Census's CT county-to-county-subdivision crosswalk. The join refuses a delineated county-equivalent that no populated block reaches, and the builder's load-back check refuses a ladder that leaves a CBSA-delineated state entirely without a CBSA. diff --git a/packages/microcosm-build/src/microcosm/build/us_runtime/block_ladder_sources.py b/packages/microcosm-build/src/microcosm/build/us_runtime/block_ladder_sources.py index 000813d96b..8b4c21baf6 100644 --- a/packages/microcosm-build/src/microcosm/build/us_runtime/block_ladder_sources.py +++ b/packages/microcosm-build/src/microcosm/build/us_runtime/block_ladder_sources.py @@ -22,16 +22,32 @@ state against the summary-level 040 state row). - OMB CBSA delineations (``list1_2023.xlsx``): county → CBSA rows keyed by the FIPS state + county code columns. +- 2020 BAF county-subdivision layer (``_MCD``): pipe-delimited + ``BLOCKID|COUNTYFP|COUSUBFP``; 2020 blocks nest in county subdivisions. +- Census Connecticut crosswalk (``ct_cou_to_cousub_crosswalk.txt``): + pipe-delimited, a quoted multi-line header, one row per county subdivision + giving its 2020 county and its 2022 planning region (the county-equivalent + that replaced it), then a glossary after a blank line. + +Connecticut is the one state whose county-equivalents differ between the +two vintages joined here. Its 2020 blocks carry the eight old counties +(09001-09015), while the OMB 2023 delineations list its nine 2022 planning +regions (09110-09190). A CT block reaches its CBSA through its town, because +planning regions are unions of whole towns. """ from __future__ import annotations +import csv import json -from collections.abc import Iterable, Mapping -from typing import Any +from collections.abc import Iterable, Iterator, Mapping +from typing import TYPE_CHECKING, Any import numpy as np +if TYPE_CHECKING: + from microcosm.build.us_runtime.geography_ladder import UsBlockLadder + #: The 50 states plus DC: (state FIPS, USPS abbreviation, P.L. 94-171 #: directory name). Puerto Rico and the island territories are outside the #: US artifact's state spine. @@ -102,6 +118,30 @@ #: Microcosm at-large convention ``state_fips * 100 + 00``. _CD_DELEGATE = "98" +#: Connecticut's state FIPS: the state whose OMB 2023 county-equivalents +#: (2022 planning regions) are not its 2020 block counties. +CT_STATE_FIPS = "09" +#: County-subdivision code for areas with no town ("County subdivisions not +#: defined": the coastal water of Fairfield, Middlesex, New Haven and New +#: London counties). Fairfield's is split between two planning regions, so the +#: code does not identify one; no populated 2020 block carries it. +_COUSUB_NOT_DEFINED = "00000" +_CT_CROSSWALK_COLUMNS = ( + "STATEFP", + "OLD_COUNTYFP", + "OLD_COUNTY_NAMELSAD", + "NEW_COUNTYFP", + "NEW_COUNTY_NAMELSAD", + "COUSUBFP", + "OLD_COUSUB_GEOID", + "NEW_COUSUB_GEOID", + "COUSUB_NAMELSAD", + "COUSUBNS", + "COUSUB_LSAD", + "COUSUB_FUNCSTAT", + "COUSUB_CLASSFP", +) + def parse_national_cd_bef(lines: Iterable[str]) -> dict[int, int]: """Parse the national CD BEF into block geoid → CD geoid (SSDD). @@ -224,6 +264,59 @@ def parse_baf_place_file(lines: Iterable[str], *, label: str) -> dict[int, int]: return result +def parse_baf_county_subdivision_file( + lines: Iterable[str], *, label: str +) -> dict[int, str]: + """Parse a BAF MCD file into block geoid → 10-digit county-subdivision geoid. + + The geoid is state + 2020 county + county-subdivision FIPS, the key the + Census Connecticut crosswalk uses for a town's 2020 identity. A row whose + ``COUNTYFP`` disagrees with the block's own county prefix is refused. + """ + + iterator = iter(lines) + header = _required_header(iterator, source=label) + if [part.strip().upper() for part in header.split("|")] != [ + "BLOCKID", + "COUNTYFP", + "COUSUBFP", + ]: + raise ValueError( + f"{label} header must be 'BLOCKID|COUNTYFP|COUSUBFP', got {header!r}." + ) + result: dict[int, str] = {} + for line_number, line in enumerate(iterator, start=2): + stripped = line.rstrip("\n") + if not stripped.strip(): + continue + parts = stripped.split("|") + if len(parts) != 3: + raise ValueError( + f"{label} line {line_number} must have three fields, got {stripped!r}." + ) + block_raw, county_raw, cousub_raw = (part.strip() for part in parts) + block = _block_geoid(block_raw, source=f"{label} line {line_number}") + if not (county_raw.isdigit() and len(county_raw) == 3): + raise ValueError( + f"{label} line {line_number} has COUNTYFP {county_raw!r}; " + "expected 3 digits." + ) + if not (cousub_raw.isdigit() and len(cousub_raw) == 5): + raise ValueError( + f"{label} line {line_number} has COUSUBFP {cousub_raw!r}; " + "expected 5 digits." + ) + if block_raw[2:5] != county_raw: + raise ValueError( + f"{label} line {line_number} puts block {block_raw} in county " + f"{county_raw}, but the block geoid's county is {block_raw[2:5]}." + ) + if block in result: + raise ValueError(f"{label} assigns block {block_raw} more than once.") + result[block] = block_raw[:5] + cousub_raw + return result + + def parse_pl_geo_blocks(lines: Iterable[str], *, state_fips: str) -> dict[int, int]: """Parse a P.L. 94-171 geo header into populated block geoid → POP100. @@ -339,6 +432,181 @@ def parse_cbsa_delineations( return result +def parse_ct_planning_region_crosswalk(lines: Iterable[str]) -> dict[str, str]: + """Parse the Census CT crosswalk into town → 2022 planning region. + + Keys are 2020 county-subdivision geoids (``09`` + old county + town, + the value :func:`parse_baf_county_subdivision_file` yields); values are + the 5-digit planning-region county-equivalent FIPS the OMB 2023 + delineations use. The table ends at the first blank row (a glossary + follows). "County subdivisions not defined" rows are skipped: Fairfield + County's is listed under two planning regions, so the code names no + single region. A populated block there would fail the block join instead. + """ + + reader = csv.reader(_csv_lines(lines), delimiter="|") + header: list[str] | None = None + for row in reader: + if any(cell.strip() for cell in row): + # "STATEFP\n(INCITS38)" → "STATEFP": the parenthesized standard + # is a second header line inside the quoted cell. + header = [_crosswalk_column_name(cell) for cell in row] + break + if header is None: + raise ValueError("CT planning-region crosswalk is empty.") + if tuple(header) != _CT_CROSSWALK_COLUMNS: + raise ValueError( + "CT planning-region crosswalk header must be " + f"{list(_CT_CROSSWALK_COLUMNS)}, got {header}." + ) + column = {name: index for index, name in enumerate(header)} + result: dict[str, str] = {} + for row in reader: + if not any(cell.strip() for cell in row): + break + where = f"CT planning-region crosswalk line {reader.line_num}" + if len(row) != len(header): + raise ValueError( + f"{where} must have {len(header)} fields, got {len(row)}: {row!r}." + ) + cells = {name: row[index].strip() for name, index in column.items()} + state = cells["STATEFP"] + old_county = cells["OLD_COUNTYFP"] + new_county = cells["NEW_COUNTYFP"] + cousub = cells["COUSUBFP"] + if state != CT_STATE_FIPS: + raise ValueError( + f"{where} has STATEFP {state!r}; expected {CT_STATE_FIPS!r}." + ) + for name, value, width in ( + ("OLD_COUNTYFP", old_county, 3), + ("NEW_COUNTYFP", new_county, 3), + ("COUSUBFP", cousub, 5), + ): + if not (value.isdigit() and len(value) == width): + raise ValueError( + f"{where} has {name} {value!r}; expected {width} digits." + ) + old_geoid = state + old_county + cousub + new_geoid = state + new_county + cousub + if (cells["OLD_COUSUB_GEOID"], cells["NEW_COUSUB_GEOID"]) != ( + old_geoid, + new_geoid, + ): + raise ValueError( + f"{where} geoids {cells['OLD_COUSUB_GEOID']!r}/" + f"{cells['NEW_COUSUB_GEOID']!r} disagree with its code columns " + f"({old_geoid}/{new_geoid})." + ) + if cousub == _COUSUB_NOT_DEFINED: + continue + region = state + new_county + existing = result.get(old_geoid) + if existing is not None and existing != region: + raise ValueError( + f"CT planning-region crosswalk puts town {old_geoid} in both " + f"{existing} and {region}." + ) + result[old_geoid] = region + if not result: + raise ValueError("CT planning-region crosswalk contained no town rows.") + return result + + +def ct_planning_region_by_block( + blocks: Iterable[int], + *, + cousub_by_block: Mapping[int, str], + planning_region_by_cousub: Mapping[str, str], +) -> dict[int, str]: + """Map each Connecticut block to its 2022 planning region, via its town. + + ``blocks`` are the populated blocks to map; blocks outside Connecticut + are ignored. 2020 blocks nest in towns and planning regions are unions of + whole towns, so the mapping is exact. A CT block without a town, or in a + town the crosswalk does not list, is a source defect and raises. + """ + + result: dict[int, str] = {} + missing_town: list[int] = [] + unmapped_town: set[str] = set() + ct_state = int(CT_STATE_FIPS) + for block in blocks: + if block // 10**13 != ct_state: + continue + cousub = cousub_by_block.get(block) + if cousub is None: + missing_town.append(block) + continue + region = planning_region_by_cousub.get(cousub) + if region is None: + unmapped_town.add(cousub) + continue + result[block] = region + if missing_town: + examples = [f"{block:015d}" for block in sorted(missing_town)[:5]] + raise ValueError( + f"{len(missing_town)} Connecticut block(s) have no county " + f"subdivision in the BAF MCD layer; examples: {examples}." + ) + if unmapped_town: + raise ValueError( + "Connecticut county subdivision(s) absent from the planning-region " + f"crosswalk: {sorted(unmapped_town)[:10]}." + ) + return result + + +def cbsa_delineated_states( + cbsa_by_county: Mapping[str, int], states: Iterable[str] +) -> list[str]: + """Return the ``states`` (2-digit FIPS) with any CBSA-delineated territory.""" + + delineated = {county[:2] for county in cbsa_by_county} + return sorted(state for state in set(states) if state in delineated) + + +def us_block_ladder_cbsa_coverage_failures(ladder: UsBlockLadder) -> list[str]: + """Return one failure per CBSA-delineated state whose blocks all lack a CBSA. + + Every state with OMB-delineated territory has populated blocks inside a + CBSA, so a state whose every block carries ``cbsa_code == 0`` is a join + defect (Connecticut's planning regions against 2020 counties was one). + The states come from the artifact's ``cbsa_delineated_states`` record. + Artifacts built before that record fall back to every state they + contain: the OMB 2023 delineations cover territory in all 50 states and + DC. For a ladder this module's join just built the check cannot fail, + since the join refuses an unreached delineated county-equivalent. Its + job is to catch a ladder built elsewhere or before that refusal existed. + """ + + block_state = ladder.block_geoid // 10**13 + present = {f"{state:02d}" for state in np.unique(block_state).tolist()} + recorded = ladder.metadata.get("cbsa_delineated_states") + if recorded is None: + delineated = sorted(present) + elif isinstance(recorded, list) and all( + isinstance(state, str) and state.isdigit() and len(state) == 2 + for state in recorded + ): + delineated = sorted(set(recorded) & present) + else: + return [ + "metadata cbsa_delineated_states must be a list of 2-digit state " + f"FIPS strings, got {recorded!r}" + ] + failures: list[str] = [] + for state in delineated: + in_state = block_state == int(state) + if not (ladder.cbsa_code[in_state] > 0).any(): + failures.append( + f"state {state}: all {int(in_state.sum()):,} populated blocks " + "have cbsa_code 0, although the CBSA delineation covers " + "territory there" + ) + return failures + + def assemble_us_block_ladder( *, block_population: Mapping[int, int], @@ -348,6 +616,7 @@ def assemble_us_block_ladder( place_by_block: Mapping[int, int], cbsa_by_county: Mapping[str, int], metadata: Mapping[str, Any], + cbsa_county_by_block: Mapping[int, str] | None = None, ) -> dict[str, np.ndarray]: """Join the parsed sources into the ladder artifact's NPZ payload. @@ -355,6 +624,14 @@ def assemble_us_block_ladder( block the CD BEF does not cover is a source defect, not a skippable row. SLD and place maps may legitimately not cover a block (states without a layer); absent entries mean unassigned. + + A block's CBSA is looked up by its county-equivalent in the delineation + vintage: its 2020 county, unless ``cbsa_county_by_block`` names another + (Connecticut's planning regions, from :func:`ct_planning_region_by_block`). + A state with any such entry must have one for every populated block. + Every delineated county-equivalent in a built state must be reached by at + least one populated block; one that is not means the delineation and the + blocks disagree on county-equivalents, and the join refuses. """ blocks = np.asarray(sorted(block_population), dtype=np.int64) @@ -378,10 +655,30 @@ def assemble_us_block_ladder( place = np.asarray( [place_by_block.get(block, 0) for block in blocks.tolist()], dtype=np.int32 ) - cbsa = np.asarray( - [cbsa_by_county.get(f"{block:015d}"[:5], 0) for block in blocks.tolist()], - dtype=np.int32, + counties, county_index = np.unique( + _cbsa_county_equivalents(blocks, cbsa_county_by_block or {}), + return_inverse=True, + ) + county_codes = [f"{county:05d}" for county in counties.tolist()] + reached = set(county_codes) + built_states = {county[:2] for county in reached} + unreached = sorted( + county + for county in cbsa_by_county + if county[:2] in built_states and county not in reached ) + if unreached: + raise ValueError( + f"{len(unreached)} CBSA-delineated county-equivalent(s) in the " + "built states match no populated block's county-equivalent; " + f"examples: {unreached[:10]}. The delineation and the blocks use " + "different county-equivalents (Connecticut's 2022 planning " + "regions against 2020 counties is the known case: pass " + "cbsa_county_by_block from ct_planning_region_by_block)." + ) + cbsa = np.asarray( + [cbsa_by_county.get(county, 0) for county in county_codes], dtype=np.int32 + )[county_index] return { "block_geoid": blocks, "population": population, @@ -394,6 +691,63 @@ def assemble_us_block_ladder( } +def _cbsa_county_equivalents( + blocks: np.ndarray, cbsa_county_by_block: Mapping[int, str] +) -> np.ndarray: + """Each block's delineation county-equivalent as a 5-digit FIPS int. + + ``blocks`` is sorted. Entries for blocks outside it (unpopulated blocks) + are validated and otherwise ignored. + """ + + counties = blocks // 10**10 + if not cbsa_county_by_block: + return counties + override_block = np.fromiter( + cbsa_county_by_block, dtype=np.int64, count=len(cbsa_county_by_block) + ) + override_county = np.empty(len(override_block), dtype=np.int64) + for index, (block, county) in enumerate(cbsa_county_by_block.items()): + if not (isinstance(county, str) and county.isdigit() and len(county) == 5): + raise ValueError( + f"cbsa_county_by_block maps block {block:015d} to {county!r}; " + "expected a 5-digit county-equivalent FIPS." + ) + if int(county[:2]) != block // 10**13: + raise ValueError( + f"cbsa_county_by_block maps block {block:015d} to " + f"county-equivalent {county} in another state." + ) + override_county[index] = int(county) + position = np.searchsorted(blocks, override_block) + populated = position < len(blocks) + populated[populated] = blocks[position[populated]] == override_block[populated] + counties[position[populated]] = override_county[populated] + remapped = np.zeros(len(blocks), dtype=bool) + remapped[position[populated]] = True + missing = ~remapped & np.isin(blocks // 10**13, np.unique(override_block // 10**13)) + if missing.any(): + examples = [f"{block:015d}" for block in blocks[missing][:5].tolist()] + raise ValueError( + f"{int(missing.sum())} populated block(s) lack a county-equivalent " + "in cbsa_county_by_block although other blocks in their state have " + f"one; examples: {examples}." + ) + return counties + + +def _csv_lines(lines: Iterable[str]) -> Iterator[str]: + """Lines with a leading byte-order mark removed; it would hide the quoting.""" + for index, line in enumerate(lines): + yield line.removeprefix("\ufeff") if index == 0 else line + + +def _crosswalk_column_name(cell: str) -> str: + """Header cell → column name, dropping a parenthesized standard suffix.""" + tokens = cell.split("(", 1)[0].split() + return tokens[0] if tokens else "" + + def _five_digit_code(value: str) -> int | None: """Normalize a spreadsheet cell to a 5-digit code; None if it is not one.""" if value.endswith(".0"): @@ -418,11 +772,17 @@ def _block_geoid(value: str, *, source: str) -> int: __all__ = [ + "CT_STATE_FIPS", "US_STATES", "assemble_us_block_ladder", + "cbsa_delineated_states", + "ct_planning_region_by_block", + "parse_baf_county_subdivision_file", "parse_baf_district_file", "parse_baf_place_file", "parse_cbsa_delineations", + "parse_ct_planning_region_crosswalk", "parse_national_cd_bef", "parse_pl_geo_blocks", + "us_block_ladder_cbsa_coverage_failures", ] diff --git a/packages/microcosm-build/tests/engine_free/us/test_us_block_ladder_sources.py b/packages/microcosm-build/tests/engine_free/us/test_us_block_ladder_sources.py index ea638558b7..7e1c251da2 100644 --- a/packages/microcosm-build/tests/engine_free/us/test_us_block_ladder_sources.py +++ b/packages/microcosm-build/tests/engine_free/us/test_us_block_ladder_sources.py @@ -8,11 +8,16 @@ from microcosm.build.us_runtime.block_ladder_sources import ( US_STATES, assemble_us_block_ladder, + cbsa_delineated_states, + ct_planning_region_by_block, + parse_baf_county_subdivision_file, parse_baf_district_file, parse_baf_place_file, parse_cbsa_delineations, + parse_ct_planning_region_crosswalk, parse_national_cd_bef, parse_pl_geo_blocks, + us_block_ladder_cbsa_coverage_failures, ) @@ -279,3 +284,579 @@ def test_assemble_refuses_populated_block_without_a_district() -> None: cbsa_by_county={}, metadata={}, ) + + +# -------------------------------------------------------------------------- +# Connecticut: 2022 planning regions against 2020 block counties +# -------------------------------------------------------------------------- + +# Rows copied from Census's ct_cou_to_cousub_crosswalk.txt (BOM already +# stripped, as the builder reads it): a quoted two-line header, town rows, +# "County subdivisions not defined" rows listed under several planning +# regions, and a glossary after blank lines. +_CT_CROSSWALK = ( + '"STATEFP\n(INCITS38)"|"OLD_COUNTYFP\n(INCITS31)"|OLD_COUNTY_NAMELSAD|' + '"NEW_COUNTYFP\n(INCITS31)"|NEW_COUNTY_NAMELSAD|COUSUBFP|OLD_COUSUB_GEOID|' + 'NEW_COUSUB_GEOID|COUSUB_NAMELSAD|"COUSUBNS\n(INCITS446)"|COUSUB_LSAD|' + "COUSUB_FUNCSTAT|COUSUB_CLASSFP\n" + "09|001|Fairfield County|120|Greater Bridgeport Planning Region|08070|" + "0900108070|0912008070|Bridgeport town|00213396|43|C|T5\n" + "09|001|Fairfield County|120|Greater Bridgeport Planning Region|00000|" + "0900100000|0912000000|County subdivisions not defined|00000000|00|F|Z9\n" + "09|001|Fairfield County|140|Naugatuck Valley Planning Region|68170|" + "0900168170|0914068170|Shelton town|00213504|43|C|T5\n" + "09|001|Fairfield County|190|Western Connecticut Planning Region|00000|" + "0900100000|0919000000|County subdivisions not defined|00000000|00|F|Z9\n" + "09|001|Fairfield County|190|Western Connecticut Planning Region|73070|" + "0900173070|0919073070|Stamford town|00213511|43|C|T5\n" + "09|009|New Haven County|170|South Central Connecticut Planning Region|52070|" + "0900952070|0917052070|New Haven town|00213471|43|C|T5\n" + "\n" + "\n" + "GLOSSARY\n" + "STATEFP = State FIPS Code / ANSI INCITS 38 Code\n" +) + +# Real 2020 blocks (BAF MCD layer): Bridgeport, Shelton, Stamford (all in +# 2020 Fairfield County 09001) and New Haven (09009). +_BRIDGEPORT_BLOCK = 90010701001000 +_SHELTON_BLOCK = 90011101001000 +_STAMFORD_BLOCK = 90010201011000 +_NEW_HAVEN_BLOCK = 90091401011000 +_CT_MCD = [ + "BLOCKID|COUNTYFP|COUSUBFP\n", + "090010701001000|001|08070\n", + "090011101001000|001|68170\n", + "090010201011000|001|73070\n", + "090091401011000|009|52070\n", +] +# OMB 2023 delineation rows for these planning regions. +_CT_CBSA_BY_COUNTY = { + "09120": 14860, # Greater Bridgeport → Bridgeport-Stamford-Danbury + "09140": 47930, # Naugatuck Valley → Waterbury-Shelton + "09170": 35300, # South Central → New Haven + "09190": 14860, # Western → Bridgeport-Stamford-Danbury +} +_DE_BLOCK = 100010401001000 + + +def _ladder_metadata(**extra: object) -> dict: + return { + "schema_version": 1, + "kind": "us_block_ladder", + "block_vintage": "2020_tabulation_blocks", + "sampling_basis": "population", + "layers": { + layer: {"vintage": "test", "source": "test source"} + for layer in ("congressional_district", "sldu", "sldl", "place", "cbsa") + }, + **extra, + } + + +def _assemble( + block_population: dict[int, int], + *, + cbsa_by_county: dict[str, int], + cbsa_county_by_block: dict[int, str] | None = None, + metadata: dict | None = None, +) -> dict[str, np.ndarray]: + return assemble_us_block_ladder( + block_population=block_population, + cd_by_block={block: (block // 10**13) * 100 + 1 for block in block_population}, + sldu_by_block={}, + sldl_by_block={}, + place_by_block={}, + cbsa_by_county=cbsa_by_county, + metadata=metadata if metadata is not None else _ladder_metadata(), + cbsa_county_by_block=cbsa_county_by_block, + ) + + +def _ct_planning_regions() -> dict[int, str]: + return ct_planning_region_by_block( + [_BRIDGEPORT_BLOCK, _SHELTON_BLOCK, _STAMFORD_BLOCK, _NEW_HAVEN_BLOCK], + cousub_by_block=parse_baf_county_subdivision_file( + _CT_MCD, label="BlockAssign_ST09_CT_MCD.txt" + ), + planning_region_by_cousub=parse_ct_planning_region_crosswalk( + _CT_CROSSWALK.splitlines(keepends=True) + ), + ) + + +def test_parse_ct_planning_region_crosswalk_reads_the_census_layout() -> None: + result = parse_ct_planning_region_crosswalk(_CT_CROSSWALK.splitlines(keepends=True)) + + # Towns only: the not-defined rows name no single planning region. + assert result == { + "0900108070": "09120", + "0900168170": "09140", + "0900173070": "09190", + "0900952070": "09170", + } + # Line endings, or a BOM the caller did not strip, do not change the parse. + assert parse_ct_planning_region_crosswalk(_CT_CROSSWALK.splitlines()) == result + assert ( + parse_ct_planning_region_crosswalk( + ("\ufeff" + _CT_CROSSWALK).splitlines(keepends=True) + ) + == result + ) + + +@pytest.mark.parametrize( + ("old", "new", "match"), + [ + ("|08070|0900108070|", "|08070|0900108071|", "disagree with its code"), + ("09|009|New Haven", "10|009|New Haven", "STATEFP '10'"), + ("|52070|0900952070|", "|5207|0900952070|", "COUSUBFP '5207'"), + ("Shelton town|00213504|43|C|T5", "Shelton town|00213504|43|C", "fields"), + ], +) +def test_parse_ct_planning_region_crosswalk_refuses_malformed_rows( + old: str, new: str, match: str +) -> None: + assert _CT_CROSSWALK.count(old) == 1 + with pytest.raises(ValueError, match=match): + parse_ct_planning_region_crosswalk( + _CT_CROSSWALK.replace(old, new).splitlines(keepends=True) + ) + + +def test_parse_ct_planning_region_crosswalk_refuses_a_town_in_two_regions() -> None: + conflicting = _CT_CROSSWALK.replace( + "\n\n\nGLOSSARY", + "\n09|001|Fairfield County|190|Western Connecticut Planning Region|08070|" + "0900108070|0919008070|Bridgeport town|00213396|43|C|T5\n\n\nGLOSSARY", + ) + with pytest.raises(ValueError, match="both 09120 and 09190"): + parse_ct_planning_region_crosswalk(conflicting.splitlines(keepends=True)) + + +def test_parse_ct_planning_region_crosswalk_refuses_a_table_without_towns() -> None: + header, rest = _CT_CROSSWALK.split("COUSUB_CLASSFP\n", 1) + undefined_only = "".join( + line + "\n" for line in rest.split("\n") if "|00000|" in line + ) + with pytest.raises(ValueError, match="no town rows"): + parse_ct_planning_region_crosswalk( + (header + "COUSUB_CLASSFP\n" + undefined_only + "\nGLOSSARY\n").splitlines( + keepends=True + ) + ) + + +def test_parse_ct_planning_region_crosswalk_refuses_a_changed_header() -> None: + with pytest.raises(ValueError, match="header must be"): + parse_ct_planning_region_crosswalk( + _CT_CROSSWALK.replace("NEW_COUNTY_NAMELSAD", "NEW_NAME").splitlines( + keepends=True + ) + ) + + +def test_parse_baf_county_subdivision_file_keys_towns_by_2020_county() -> None: + result = parse_baf_county_subdivision_file(_CT_MCD, label="MCD") + + assert result == { + _BRIDGEPORT_BLOCK: "0900108070", + _SHELTON_BLOCK: "0900168170", + _STAMFORD_BLOCK: "0900173070", + _NEW_HAVEN_BLOCK: "0900952070", + } + + +@pytest.mark.parametrize( + ("row", "match"), + [ + ("090010701001000|009|08070\n", "block geoid's county is 001"), + ("090010701001000|001|8070\n", "COUSUBFP '8070'"), + ("090010701001000|001\n", "three fields"), + ], +) +def test_parse_baf_county_subdivision_file_refuses_bad_rows( + row: str, match: str +) -> None: + with pytest.raises(ValueError, match=match): + parse_baf_county_subdivision_file([_CT_MCD[0], row], label="MCD") + + +def test_parse_baf_county_subdivision_file_refuses_repeats_and_other_headers() -> None: + with pytest.raises(ValueError, match="more than once"): + parse_baf_county_subdivision_file([*_CT_MCD, _CT_MCD[1]], label="MCD") + with pytest.raises(ValueError, match="header must be"): + parse_baf_county_subdivision_file( + ["BLOCKID|COUSUBFP\n", "090010701001000|08070\n"], label="MCD" + ) + + +def test_ct_planning_region_by_block_follows_each_town() -> None: + # Fairfield County (09001) splits across three planning regions, so the + # 2020 county alone cannot name a block's county-equivalent. + assert _ct_planning_regions() == { + _BRIDGEPORT_BLOCK: "09120", + _SHELTON_BLOCK: "09140", + _STAMFORD_BLOCK: "09190", + _NEW_HAVEN_BLOCK: "09170", + } + + +def test_ct_planning_region_by_block_refuses_blocks_it_cannot_place() -> None: + crosswalk = parse_ct_planning_region_crosswalk( + _CT_CROSSWALK.splitlines(keepends=True) + ) + with pytest.raises(ValueError, match="no county subdivision"): + ct_planning_region_by_block( + [_BRIDGEPORT_BLOCK, _DE_BLOCK], + cousub_by_block={}, + planning_region_by_cousub=crosswalk, + ) + # A populated block in an undefined (water) area names no planning region. + with pytest.raises(ValueError, match="absent from the planning-region"): + ct_planning_region_by_block( + [_BRIDGEPORT_BLOCK], + cousub_by_block={_BRIDGEPORT_BLOCK: "0900100000"}, + planning_region_by_cousub=crosswalk, + ) + # Non-CT blocks are not its business. + assert ( + ct_planning_region_by_block( + [_DE_BLOCK], cousub_by_block={}, planning_region_by_cousub=crosswalk + ) + == {} + ) + + +def test_assemble_looks_up_connecticut_cbsa_by_planning_region(tmp_path) -> None: + from microcosm.build.us_runtime import load_us_block_ladder + + cbsa_by_county = {**_CT_CBSA_BY_COUNTY, "10001": 20100} + blocks = { + _BRIDGEPORT_BLOCK: 40, + _SHELTON_BLOCK: 30, + _STAMFORD_BLOCK: 20, + _NEW_HAVEN_BLOCK: 10, + _DE_BLOCK: 5, + } + metadata = _ladder_metadata( + cbsa_delineated_states=cbsa_delineated_states(cbsa_by_county, ["09", "10"]) + ) + payload = _assemble( + blocks, + cbsa_by_county=cbsa_by_county, + cbsa_county_by_block=_ct_planning_regions(), + metadata=metadata, + ) + path = tmp_path / "ladder.npz" + np.savez_compressed(path, **payload) + ladder = load_us_block_ladder(path) + + cbsa = dict( + zip(ladder.block_geoid.tolist(), ladder.cbsa_code.tolist(), strict=True) + ) + assert cbsa == { + _BRIDGEPORT_BLOCK: 14860, + _SHELTON_BLOCK: 47930, + _STAMFORD_BLOCK: 14860, + _NEW_HAVEN_BLOCK: 35300, + _DE_BLOCK: 20100, + } + assert us_block_ladder_cbsa_coverage_failures(ladder) == [] + + +def test_assemble_refuses_ct_blocks_keyed_by_2020_county() -> None: + # The pre-fix join: CT blocks looked up by 09001/09009 reach none of the + # planning-region rows, which used to leave every CT block at CBSA 0. + with pytest.raises(ValueError, match=r"4 CBSA-delineated county-equivalent"): + _assemble( + {_BRIDGEPORT_BLOCK: 40, _NEW_HAVEN_BLOCK: 10}, + cbsa_by_county=_CT_CBSA_BY_COUNTY, + ) + + +def test_assemble_refuses_a_state_only_partly_remapped() -> None: + with pytest.raises(ValueError, match="lack a county-equivalent"): + _assemble( + {_BRIDGEPORT_BLOCK: 40, _NEW_HAVEN_BLOCK: 10}, + cbsa_by_county={"09120": 14860}, + cbsa_county_by_block={_BRIDGEPORT_BLOCK: "09120"}, + ) + + +@pytest.mark.parametrize( + ("county", "match"), + [("10001", "in another state"), ("0912", "5-digit"), (9120, "5-digit")], +) +def test_assemble_refuses_malformed_county_equivalents( + county: object, match: str +) -> None: + with pytest.raises(ValueError, match=match): + _assemble( + {_BRIDGEPORT_BLOCK: 40}, + cbsa_by_county={}, + cbsa_county_by_block={_BRIDGEPORT_BLOCK: county}, + ) + + +def test_assemble_ignores_remap_entries_for_unpopulated_blocks() -> None: + # The MCD layer covers every block, populated or not; entries for blocks + # outside the ladder (before, between and after its rows) change nothing. + regions = { + **_ct_planning_regions(), + 90010101011000: "09190", + 90010900001000: "09120", + 99999999999999: "09170", + } + payload = _assemble( + {_BRIDGEPORT_BLOCK: 40, _NEW_HAVEN_BLOCK: 10}, + cbsa_by_county={"09120": 14860, "09170": 35300}, + cbsa_county_by_block=regions, + ) + + assert payload["block_geoid"].tolist() == [_BRIDGEPORT_BLOCK, _NEW_HAVEN_BLOCK] + assert payload["cbsa_code"].tolist() == [14860, 35300] + + +def test_assemble_applies_the_whole_state_rule_to_unpopulated_entries() -> None: + # A remap that names only unpopulated CT blocks still claims the state. + with pytest.raises(ValueError, match="2 populated block"): + _assemble( + {_BRIDGEPORT_BLOCK: 40, _NEW_HAVEN_BLOCK: 10}, + cbsa_by_county={}, + cbsa_county_by_block={90010101011000: "09190"}, + ) + + +def test_assemble_ignores_delineation_rows_outside_the_built_states() -> None: + # A smoke build of Delaware alone must not trip over Connecticut rows. + payload = _assemble({_DE_BLOCK: 5}, cbsa_by_county=_CT_CBSA_BY_COUNTY) + + assert payload["cbsa_code"].tolist() == [0] + + +def _write_loaded_ladder(tmp_path, payload: dict[str, np.ndarray]): + from microcosm.build.us_runtime import load_us_block_ladder + + path = tmp_path / "ladder.npz" + np.savez_compressed(path, **payload) + return load_us_block_ladder(path) + + +def _legacy_ct_payload(**metadata_extra: object) -> dict[str, np.ndarray]: + """A ladder shaped like the pre-fix national build: CT all at CBSA 0.""" + payload = _assemble( + {_BRIDGEPORT_BLOCK: 40, _NEW_HAVEN_BLOCK: 10, _DE_BLOCK: 5}, + cbsa_by_county={"10001": 20100}, + metadata=_ladder_metadata(**metadata_extra), + ) + assert payload["cbsa_code"].tolist() == [0, 0, 20100] + return payload + + +def test_cbsa_coverage_check_flags_the_legacy_connecticut_ladder(tmp_path) -> None: + # No cbsa_delineated_states record (artifacts built before it): every + # state present is taken as delineated, so all-zero CT fails. + ladder = _write_loaded_ladder(tmp_path, _legacy_ct_payload()) + + failures = us_block_ladder_cbsa_coverage_failures(ladder) + + assert len(failures) == 1 + assert failures[0].startswith("state 09: all 2 populated blocks") + + +def test_cbsa_coverage_check_follows_the_recorded_states(tmp_path) -> None: + recorded = _write_loaded_ladder( + tmp_path, _legacy_ct_payload(cbsa_delineated_states=["09", "10"]) + ) + assert [ + failure[:8] for failure in us_block_ladder_cbsa_coverage_failures(recorded) + ] == ["state 09"] + # A state the delineation does not cover may legitimately be all zero. + uncovered = _write_loaded_ladder( + tmp_path, _legacy_ct_payload(cbsa_delineated_states=["10"]) + ) + assert us_block_ladder_cbsa_coverage_failures(uncovered) == [] + malformed = _write_loaded_ladder( + tmp_path, _legacy_ct_payload(cbsa_delineated_states="09") + ) + assert "must be a list" in us_block_ladder_cbsa_coverage_failures(malformed)[0] + + +def test_cbsa_coverage_check_covers_every_recorded_state(tmp_path) -> None: + # CT is fine and Delaware, second in order, is all zero. + payload = _assemble( + {_BRIDGEPORT_BLOCK: 40, _DE_BLOCK: 5}, + cbsa_by_county={"09120": 14860}, + cbsa_county_by_block={_BRIDGEPORT_BLOCK: "09120"}, + metadata=_ladder_metadata(cbsa_delineated_states=["09", "10", "36"]), + ) + failures = us_block_ladder_cbsa_coverage_failures( + _write_loaded_ladder(tmp_path, payload) + ) + + # New York is recorded but has no blocks in this ladder: not a failure. + assert [failure[:8] for failure in failures] == ["state 10"] + + +@pytest.mark.parametrize("recorded", [["9"], [9], ["090"], [None]]) +def test_cbsa_coverage_check_refuses_malformed_state_entries( + tmp_path, recorded: list +) -> None: + ladder = _write_loaded_ladder( + tmp_path, _legacy_ct_payload(cbsa_delineated_states=recorded) + ) + + assert "must be a list" in us_block_ladder_cbsa_coverage_failures(ladder)[0] + + +def test_cbsa_delineated_states_keeps_built_states_with_delineated_territory() -> None: + assert cbsa_delineated_states( + {"09120": 14860, "10001": 20100, "72127": 41980}, ["09", "10", "11"] + ) == ["09", "10"] + + +# -------------------------------------------------------------------------- +# Invariants of the CBSA join, for all inputs +# -------------------------------------------------------------------------- + +_STATES = ("09", "10", "36") + + +def _block_strategy(st): + return st.builds( + lambda state, county, tract_block: int( + f"{state}{county:03d}{tract_block:010d}" + ), + st.sampled_from(_STATES), + st.integers(1, 4), + st.integers(1, 10**10 - 1), + ) + + +def test_cbsa_join_invariants_hold_for_generated_ladders(tmp_path_factory) -> None: + """For any blocks, delineation and CT remap: + + 1. the join refuses exactly when a delineated county-equivalent in a + built state is reached by no populated block; + 2. otherwise each block's cbsa_code is the delineation's code for its + county-equivalent (0 when undelineated); + 3. every other array is independent of the CBSA inputs; + 4. a ladder the join accepts passes the loader-level coverage check. + """ + pytest.importorskip("hypothesis") + from hypothesis import given, settings + from hypothesis import strategies as st + + from microcosm.build.us_runtime import load_us_block_ladder + + @settings(max_examples=200, deadline=None) + @given( + blocks=st.dictionaries( + _block_strategy(st), st.integers(1, 1000), min_size=1, max_size=12 + ), + regions=st.lists(st.sampled_from(["09110", "09120", "09190"]), max_size=12), + delineation=st.dictionaries( + st.sampled_from( + [ + "09001", + "09002", + "09110", + "09120", + "09190", + "10001", + "10003", + "36001", + "36004", + ] + ), + st.integers(10000, 99999), + max_size=6, + ), + remap_ct=st.booleans(), + ) + def check(blocks, regions, delineation, remap_ct) -> None: + ordered = sorted(blocks) + ct_blocks = [block for block in ordered if f"{block:015d}"[:2] == "09"] + remap = ( + { + block: regions[index % len(regions)] + for index, block in enumerate(ct_blocks) + } + if remap_ct and regions + else {} + ) + key = {block: remap.get(block, f"{block:015d}"[:5]) for block in ordered} + built = {f"{block:015d}"[:2] for block in ordered} + unreached = { + county + for county in delineation + if county[:2] in built and county not in set(key.values()) + } + if unreached: + with pytest.raises(ValueError, match="CBSA-delineated county-equivalent"): + _assemble( + blocks, cbsa_by_county=delineation, cbsa_county_by_block=remap + ) + return + metadata = _ladder_metadata( + cbsa_delineated_states=cbsa_delineated_states(delineation, built) + ) + payload = _assemble( + blocks, + cbsa_by_county=delineation, + cbsa_county_by_block=remap, + metadata=metadata, + ) + assert payload["cbsa_code"].tolist() == [ + delineation.get(key[block], 0) for block in ordered + ] + baseline = _assemble(blocks, cbsa_by_county={}, metadata=metadata) + for name, values in payload.items(): + if name != "cbsa_code": + np.testing.assert_array_equal(values, baseline[name]) + path = tmp_path_factory.mktemp("ladder") / "ladder.npz" + np.savez_compressed(path, **payload) + assert us_block_ladder_cbsa_coverage_failures(load_us_block_ladder(path)) == [] + + check() + + +def test_ct_crosswalk_round_trips_generated_towns() -> None: + """Any town table written in the Census layout parses back to itself.""" + pytest.importorskip("hypothesis") + from hypothesis import given, settings + from hypothesis import strategies as st + + header = _CT_CROSSWALK.split("COUSUB_CLASSFP\n", 1)[0] + "COUSUB_CLASSFP\n" + + @settings(max_examples=200, deadline=None) + @given( + towns=st.dictionaries( + st.tuples(st.integers(1, 15), st.integers(1, 99999)), + st.integers(110, 190), + min_size=1, + max_size=20, + ), + undefined=st.lists( + st.tuples(st.integers(1, 15), st.integers(110, 190)), max_size=4 + ), + ) + def check(towns, undefined) -> None: + rows = [ + f"09|{old:03d}|X County|{new:03d}|Y Planning Region|{cousub:05d}|" + f"09{old:03d}{cousub:05d}|09{new:03d}{cousub:05d}|Z town|0|43|A|T1\n" + for (old, cousub), new in towns.items() + ] + [ + f"09|{old:03d}|X County|{new:03d}|Y Planning Region|00000|" + f"09{old:03d}00000|09{new:03d}00000|County subdivisions not defined|" + "00000000|00|F|Z9\n" + for old, new in undefined + ] + text = header + "".join(rows) + "\n\nGLOSSARY\n" + assert parse_ct_planning_region_crosswalk(text.splitlines(keepends=True)) == { + f"09{old:03d}{cousub:05d}": f"09{new:03d}" + for (old, cousub), new in towns.items() + } + + check() diff --git a/packages/microcosm-build/tests/engine_free/us/test_us_block_ladder_tool.py b/packages/microcosm-build/tests/engine_free/us/test_us_block_ladder_tool.py new file mode 100644 index 0000000000..d52e517331 --- /dev/null +++ b/packages/microcosm-build/tests/engine_free/us/test_us_block_ladder_tool.py @@ -0,0 +1,265 @@ +"""End-to-end tests for tools/build_us_block_ladder_artifact.py on a local cache. + +Every source the builder downloads is written into a temporary cache under the +file name the builder caches it as, so ``main()`` runs offline; any network +request fails the test. +""" + +from __future__ import annotations + +import importlib.util +import json +import urllib.request +import zipfile +from pathlib import Path + +import numpy as np +import pytest + +from microcosm.build.us_runtime import load_us_block_ladder +from microcosm.build.us_runtime.block_ladder_sources import ( + us_block_ladder_cbsa_coverage_failures, +) +from test_support.paths import paths_for + +_TEST_PATHS = paths_for("microcosm-build") + +# Real 2020 blocks (BAF MCD layer): Bridgeport, Shelton, Stamford (all in 2020 +# Fairfield County 09001), New Haven (09009), and one Delaware block. +_BLOCKS = { + "090010701001000": ("08070", "04", 40), + "090011101001000": ("68170", "03", 30), + "090010201011000": ("73070", "04", 20), + "090091401011000": ("52070", "03", 10), +} +_DE_BLOCK = "100010401001000" + +# Rows copied from Census's ct_cou_to_cousub_crosswalk.txt, BOM included. +_CT_CROSSWALK = ( + "" + '"STATEFP\n(INCITS38)"|"OLD_COUNTYFP\n(INCITS31)"|OLD_COUNTY_NAMELSAD|' + '"NEW_COUNTYFP\n(INCITS31)"|NEW_COUNTY_NAMELSAD|COUSUBFP|OLD_COUSUB_GEOID|' + 'NEW_COUSUB_GEOID|COUSUB_NAMELSAD|"COUSUBNS\n(INCITS446)"|COUSUB_LSAD|' + "COUSUB_FUNCSTAT|COUSUB_CLASSFP\n" + "09|001|Fairfield County|120|Greater Bridgeport Planning Region|08070|" + "0900108070|0912008070|Bridgeport town|00213396|43|C|T5\n" + "09|001|Fairfield County|120|Greater Bridgeport Planning Region|00000|" + "0900100000|0912000000|County subdivisions not defined|00000000|00|F|Z9\n" + "09|001|Fairfield County|140|Naugatuck Valley Planning Region|68170|" + "0900168170|0914068170|Shelton town|00213504|43|C|T5\n" + "09|001|Fairfield County|190|Western Connecticut Planning Region|73070|" + "0900173070|0919073070|Stamford town|00213511|43|C|T5\n" + "09|009|New Haven County|170|South Central Connecticut Planning Region|52070|" + "0900952070|0917052070|New Haven town|00213471|43|C|T5\n" + "\n" + "\n" + "GLOSSARY\n" + "STATEFP = State FIPS Code / ANSI INCITS 38 Code\n" +) +# (CBSA Code, CBSA Title, FIPS State Code, FIPS County Code), as in +# list1_2023.xlsx. Puerto Rico is outside the state spine and never built. +_DELINEATION_ROWS = ( + ("14860", "Bridgeport-Stamford-Danbury, CT", "09", "120"), + ("47930", "Waterbury-Shelton, CT", "09", "140"), + ("35300", "New Haven, CT", "09", "170"), + ("14860", "Bridgeport-Stamford-Danbury, CT", "09", "190"), + ("20100", "Dover, DE", "10", "001"), + ("41980", "San Juan-Bayamón-Caguas, PR", "72", "127"), +) + + +def _load_tool_module(): + path = _TEST_PATHS.repository / "tools" / "build_us_block_ladder_artifact.py" + spec = importlib.util.spec_from_file_location( + "build_us_block_ladder_artifact", path + ) + module = importlib.util.module_from_spec(spec) + assert spec.loader is not None + spec.loader.exec_module(module) + return module + + +def _pl_line(summary_level: str, geocode: str, population: int) -> str: + fields = [""] * 97 + fields[0], fields[2], fields[9], fields[90] = ( + "PLST", + summary_level, + geocode, + str(population), + ) + return "|".join(fields) + "\n" + + +def _write_zip(path: Path, members: dict[str, str]) -> None: + with zipfile.ZipFile(path, "w") as archive: + for name, text in members.items(): + archive.writestr(name, text) + + +def _write_cache(cache: Path, *, with_ct_mcd: bool = True) -> None: + import openpyxl + + cache.mkdir() + _write_zip( + cache / "cd119.zip", + { + "NationalCD119.txt": "GEOID,CDFP\n" + + "".join(f"{block},{cd}\n" for block, (_, cd, _) in _BLOCKS.items()) + + f"{_DE_BLOCK},00\n" + }, + ) + workbook = openpyxl.Workbook() + sheet = workbook.active + sheet.append(("List 1. CORE BASED STATISTICAL AREAS (CBSAs)",)) + sheet.append(("CBSA Code", "CBSA Title", "FIPS State Code", "FIPS County Code")) + for row in _DELINEATION_ROWS: + sheet.append(row) + workbook.save(cache / "list1_2023.xlsx") + ct_population = sum(population for _, _, population in _BLOCKS.values()) + _write_zip( + cache / "ct2020.pl.zip", + { + "ctgeo2020.pl": _pl_line("040", "09", ct_population) + + "".join( + _pl_line("750", block, population) + for block, (_, _, population) in _BLOCKS.items() + ) + }, + ) + _write_zip( + cache / "de2020.pl.zip", + {"degeo2020.pl": _pl_line("040", "10", 5) + _pl_line("750", _DE_BLOCK, 5)}, + ) + ct_baf = { + "BlockAssign_ST09_CT_SLDU.txt": "BLOCKID|DISTRICT\n" + + "".join(f"{block}|001\n" for block in _BLOCKS), + "BlockAssign_ST09_CT_SLDL.txt": "BLOCKID|DISTRICT\n" + + "".join(f"{block}|120\n" for block in _BLOCKS), + "BlockAssign_ST09_CT_INCPLACE_CDP.txt": "BLOCKID|PLACEFP\n" + + "".join(f"{block}|\n" for block in _BLOCKS), + } + if with_ct_mcd: + ct_baf["BlockAssign_ST09_CT_MCD.txt"] = "BLOCKID|COUNTYFP|COUSUBFP\n" + "".join( + f"{block}|{block[2:5]}|{cousub}\n" + for block, (cousub, _, _) in _BLOCKS.items() + ) + _write_zip(cache / "BlockAssign_ST09_CT.zip", ct_baf) + _write_zip( + cache / "BlockAssign_ST10_DE.zip", + {"BlockAssign_ST10_DE_SLDU.txt": f"BLOCKID|DISTRICT\n{_DE_BLOCK}|001\n"}, + ) + (cache / "ct_cou_to_cousub_crosswalk.txt").write_text( + _CT_CROSSWALK, encoding="utf-8" + ) + + +@pytest.fixture(autouse=True) +def _offline(monkeypatch): + def refuse(*args, **kwargs): + raise AssertionError(f"unexpected download: {args[0]!r}") + + monkeypatch.setattr(urllib.request, "urlopen", refuse) + + +def test_builder_gives_connecticut_blocks_their_planning_region_cbsa( + tmp_path, capsys +) -> None: + tool = _load_tool_module() + cache = tmp_path / "cache" + _write_cache(cache) + out = tmp_path / "ladder.npz" + + tool.main(["--out", str(out), "--cache-dir", str(cache), "--states", "09,10"]) + + ladder = load_us_block_ladder(out) + cbsa = dict( + zip( + (f"{block:015d}" for block in ladder.block_geoid.tolist()), + ladder.cbsa_code.tolist(), + strict=True, + ) + ) + assert cbsa == { + "090010201011000": 14860, # Stamford → Western → Bridgeport-Stamford + "090010701001000": 14860, # Bridgeport → Greater Bridgeport + "090011101001000": 47930, # Shelton → Naugatuck Valley → Waterbury + "090091401011000": 35300, # New Haven → South Central → New Haven + _DE_BLOCK: 20100, + } + assert us_block_ladder_cbsa_coverage_failures(ladder) == [] + metadata = ladder.metadata + assert metadata["cbsa_delineated_states"] == ["09", "10"] + assert metadata["layers"]["cbsa"]["vintage"] == "omb_2023_delineations" + assert ( + metadata["layers"]["cbsa"]["ct_planning_region_crosswalk_url"] + == tool.CT_PLANNING_REGION_CROSSWALK_URL + ) + crosswalk_source = metadata["source_files"]["ct_planning_region_crosswalk"] + assert crosswalk_source["path"] == str(cache / "ct_cou_to_cousub_crosswalk.txt") + assert len(crosswalk_source["sha256"]) == 64 + summary = json.loads(out.with_suffix(".summary.json").read_text()) + assert summary["cbsa_population_share_by_state"] == {"09": 1.0, "10": 1.0} + assert summary["cbsa_delineated_states"] == ["09", "10"] + capsys.readouterr() + + +def test_builder_without_connecticut_never_reads_the_crosswalk( + tmp_path, capsys +) -> None: + tool = _load_tool_module() + cache = tmp_path / "cache" + _write_cache(cache) + # Absent from the cache: reading it would attempt a (refused) download. + (cache / "ct_cou_to_cousub_crosswalk.txt").unlink() + out = tmp_path / "ladder.npz" + + tool.main(["--out", str(out), "--cache-dir", str(cache), "--states", "10"]) + + ladder = load_us_block_ladder(out) + assert ladder.cbsa_code.tolist() == [20100] + assert ladder.metadata["cbsa_delineated_states"] == ["10"] + assert "ct_planning_region_crosswalk" not in ladder.metadata["source_files"] + capsys.readouterr() + + +def test_builder_refuses_connecticut_without_its_county_subdivision_layer( + tmp_path, +) -> None: + tool = _load_tool_module() + cache = tmp_path / "cache" + _write_cache(cache, with_ct_mcd=False) + + with pytest.raises(SystemExit, match="BlockAssign_ST09_CT_MCD.txt"): + tool.main( + [ + "--out", + str(tmp_path / "ladder.npz"), + "--cache-dir", + str(cache), + "--states", + "09,10", + ] + ) + assert not (tmp_path / "ladder.npz").exists() + + +def test_builder_output_matches_the_pre_fix_join_outside_connecticut( + tmp_path, capsys +) -> None: + """Differential: outside CT the join equals the 2020-county lookup.""" + tool = _load_tool_module() + cache = tmp_path / "cache" + _write_cache(cache) + out = tmp_path / "ladder.npz" + tool.main(["--out", str(out), "--cache-dir", str(cache), "--states", "09,10"]) + capsys.readouterr() + + ladder = load_us_block_ladder(out) + by_county = { + f"{state}{county}": int(code) for code, _, state, county in _DELINEATION_ROWS + } + outside = ladder.block_geoid // 10**13 != 9 + expected = np.asarray( + [by_county.get(f"{block:015d}"[:5], 0) for block in ladder.block_geoid[outside]] + ) + np.testing.assert_array_equal(ladder.cbsa_code[outside], expected) diff --git a/tools/build_us_block_ladder_artifact.py b/tools/build_us_block_ladder_artifact.py index 1627059a46..6fc08c2f5d 100644 --- a/tools/build_us_block_ladder_artifact.py +++ b/tools/build_us_block_ladder_artifact.py @@ -2,15 +2,22 @@ Downloads (with a local cache) the 2020 P.L. 94-171 geographic headers (block populations), the 119th Congressional District block equivalency -file, the 2020 Block Assignment Files (SLDU/SLDL/INCPLACE_CDP), and the OMB -CBSA delineation workbook; joins them at 2020-tabulation-block grain; and -writes one national NPZ artifact whose embedded metadata records a vintage -and source per derived layer (``vintage_policy: error`` — the loader refuses -an artifact missing any of them). No per-area files, the standing rule. +file, the 2020 Block Assignment Files (SLDU/SLDL/INCPLACE_CDP, plus the MCD +layer for Connecticut), the OMB CBSA delineation workbook, and Census's +Connecticut county-to-county-subdivision crosswalk; joins them at +2020-tabulation-block grain; and writes one national NPZ artifact whose +embedded metadata records a vintage and source per derived layer +(``vintage_policy: error`` — the loader refuses an artifact missing any of +them). No per-area files, the standing rule. + +Connecticut blocks reach their CBSA through their town's 2022 planning +region, the county-equivalent the OMB 2023 delineations list the state by +(see ``block_ladder_sources``). The artifact is self-checked by loading it back through ``microcosm.build.us_runtime.load_us_block_ladder`` before the summary is -written, so a published ladder is by construction a loadable ladder. +written, so a published ladder is by construction a loadable ladder; the +loaded ladder must also leave no CBSA-delineated state without a CBSA. Example: uv run python tools/build_us_block_ladder_artifact.py \ @@ -37,13 +44,19 @@ from microcosm.build.us_runtime import load_us_block_ladder from microcosm.build.us_runtime.block_ladder_sources import ( + CT_STATE_FIPS, US_STATES, assemble_us_block_ladder, + cbsa_delineated_states, + ct_planning_region_by_block, + parse_baf_county_subdivision_file, parse_baf_district_file, parse_baf_place_file, parse_cbsa_delineations, + parse_ct_planning_region_crosswalk, parse_national_cd_bef, parse_pl_geo_blocks, + us_block_ladder_cbsa_coverage_failures, ) CD119_BEF_URL = ( @@ -63,6 +76,10 @@ "https://www2.census.gov/programs-surveys/metro-micro/geographies/" "reference-files/2023/delineation-files/list1_2023.xlsx" ) +CT_PLANNING_REGION_CROSSWALK_URL = ( + "https://www2.census.gov/geo/docs/reference/ct_change/" + "ct_cou_to_cousub_crosswalk.txt" +) LAYER_VINTAGES = { "congressional_district": { @@ -95,9 +112,12 @@ "vintage": "omb_2023_delineations", "source": ( "OMB Bulletin No. 23-01 CBSA delineation file (list1_2023.xlsx), " - "county to CBSA" + "county-equivalent to CBSA; Connecticut blocks by 2022 planning " + "region via the 2020 BAF MCD layer and the Census CT " + "county-to-county-subdivision crosswalk" ), "url": CBSA_DELINEATIONS_URL, + "ct_planning_region_crosswalk_url": CT_PLANNING_REGION_CROSSWALK_URL, }, } @@ -220,6 +240,7 @@ def main(argv: list[str] | None = None) -> None: sldu_by_block: dict[int, str] = {} sldl_by_block: dict[int, str] = {} place_by_block: dict[int, int] = {} + ct_cousub_by_block: dict[int, str] = {} missing_layers: dict[str, list[str]] = {"SLDU": [], "SLDL": [], "INCPLACE_CDP": []} for fips, usps, dirname in states: _log(f" state {fips} {usps}") @@ -261,6 +282,36 @@ def main(argv: list[str] | None = None) -> None: target.update(parse_baf_place_file(lines, label=member)) else: target.update(parse_baf_district_file(lines, label=member)) + if fips == CT_STATE_FIPS: + # Connecticut's CBSAs are delineated by planning region, which a + # block reaches only through its town (see block_ladder_sources). + member = f"BlockAssign_ST{fips}_{usps}_MCD.txt" + lines = _optional_baf_layer(baf_zip, member) + if lines is None: + raise SystemExit( + f"{baf_zip} has no {member}; Connecticut blocks need their " + "county subdivision to reach a CBSA." + ) + ct_cousub_by_block = parse_baf_county_subdivision_file(lines, label=member) + + cbsa_county_by_block: dict[int, str] = {} + if any(fips == CT_STATE_FIPS for fips, _, _ in states): + crosswalk = _download(CT_PLANNING_REGION_CROSSWALK_URL, args.cache_dir) + source_files["ct_planning_region_crosswalk"] = { + "path": str(crosswalk), + "sha256": _sha256(crosswalk), + } + cbsa_county_by_block = ct_planning_region_by_block( + block_population, + cousub_by_block=ct_cousub_by_block, + planning_region_by_cousub=parse_ct_planning_region_crosswalk( + crosswalk.read_text(encoding="utf-8-sig").splitlines(keepends=True) + ), + ) + _log( + f" Connecticut: {len(cbsa_county_by_block):,} blocks in " + f"{len(set(cbsa_county_by_block.values()))} planning regions" + ) metadata = { "schema_version": 1, @@ -272,6 +323,11 @@ def main(argv: list[str] | None = None) -> None: "missing_baf_layers": { layer: sorted(values) for layer, values in missing_layers.items() }, + # States the delineation covers territory in: each must come out with + # a CBSA somewhere (us_block_ladder_cbsa_coverage_failures). + "cbsa_delineated_states": cbsa_delineated_states( + cbsa_by_county, (fips for fips, _, _ in states) + ), "source_files": source_files, } _log(f" assembling {len(block_population):,} populated blocks") @@ -283,11 +339,20 @@ def main(argv: list[str] | None = None) -> None: place_by_block=place_by_block, cbsa_by_county=cbsa_by_county, metadata=metadata, + cbsa_county_by_block=cbsa_county_by_block, ) args.out.parent.mkdir(parents=True, exist_ok=True) np.savez_compressed(args.out, **payload) ladder = load_us_block_ladder(args.out) # self-check: must load cleanly + coverage_failures = us_block_ladder_cbsa_coverage_failures(ladder) + if coverage_failures: + raise SystemExit( + f"{args.out} failed the CBSA coverage check:\n " + + "\n ".join(coverage_failures) + ) + block_state = ladder.block_geoid // 10**13 + in_cbsa = ladder.cbsa_code > 0 summary = { "output": str(args.out.resolve()), "output_sha256": _sha256(args.out), @@ -300,6 +365,14 @@ def main(argv: list[str] | None = None) -> None: "layer_vintages": ladder.layer_vintages, "states": [fips for fips, _, _ in states], "missing_baf_layers": metadata["missing_baf_layers"], + "cbsa_delineated_states": metadata["cbsa_delineated_states"], + "cbsa_population_share_by_state": { + f"{state:02d}": float( + ladder.population[(block_state == state) & in_cbsa].sum() + / ladder.population[block_state == state].sum() + ) + for state in np.unique(block_state).tolist() + }, "source_files": source_files, } summary_path = (