#!/usr/bin/env python3 """Reproduce the published analysis from the accompanying, checked source-value files. Usage: python reproduce-cold-storage-analysis.py [--data-dir PATH] Python 3.10+; standard library only. No network calls or external dependencies. This checks transcription-file consistency and calculations, not the accuracy of USDA's underlying estimates, individual facilities, or source pages after verification. """ from __future__ import annotations import argparse import csv import json import sys from decimal import Decimal, ROUND_HALF_UP from pathlib import Path from typing import Any STATE_BASE = "us-refrigerated-warehouse-capacity-by-state-2023-2025" HISTORY_FILE = "us-refrigerated-capacity-by-warehouse-type-2017-2025.csv" STOCK_FILE = "us-cold-storage-stocks-2026-07-31.csv" def read_csv(path: Path) -> list[dict[str, str]]: with path.open(encoding="utf-8-sig", newline="") as handle: reader = csv.DictReader(handle) rows = list(reader) if not rows or not reader.fieldnames: raise ValueError(f"Empty or headerless CSV: {path.name}") if any(None in row or any(v is None for v in row.values()) for row in rows): raise ValueError(f"Malformed CSV row: {path.name}") return rows def number(value: Any) -> Decimal | None: if value in (None, "", "D"): return None return Decimal(str(value)) def required(value: Any) -> Decimal: result = number(value) if result is None: raise ValueError(f"A disclosed numeric value was required, received {value!r}") return result def rounded(value: Decimal, places: int = 1) -> Decimal: return value.quantize(Decimal(1).scaleb(-places), rounding=ROUND_HALF_UP) def percent(numerator: Any, denominator: Any) -> Decimal | None: a, b = number(numerator), number(denominator) return a / b * 100 if a is not None and b is not None and b > 0 else None def expect(condition: bool, message: str) -> None: if not condition: raise ValueError(message) def check_value(actual: str, expected: Decimal | None, places: int, label: str) -> None: value = number(actual) if expected is None: expect(value is None, f"{label}: expected unavailable, found {actual}") else: expect(value == rounded(expected, places), f"{label}: {actual!r}, expected {rounded(expected, places)}") def output_number(value: Decimal | None, places: int = 1) -> float | None: return None if value is None else float(rounded(value, places)) def reproduce(data_dir: Path) -> dict[str, Any]: rows = read_csv(data_dir / f"{STATE_BASE}.csv") package = json.loads((data_dir / f"{STATE_BASE}.json").read_text(encoding="utf-8")) expect(len(rows) == 51, "Expected 50 states and one U.S. row") expect(len({r["state_abbr"] for r in rows}) == 51, "Duplicate state codes") by_state = {r["state_abbr"]: r for r in rows} national = by_state["US"] states = [r for r in rows if r["state_abbr"] != "US"] json_by_state = {r["state_abbr"]: r for r in package["rows"]} expect(set(json_by_state) == set(by_state), "CSV/JSON state sets differ") for row in rows: counterpart = json_by_state[row["state_abbr"]] expect(set(row) == set(counterpart), f"Columns differ: {row['state']}") for field, value in row.items(): other = counterpart[field] if value == "": expect(other is None, f"{row['state']}.{field}: empty CSV must map to JSON null") elif isinstance(other, (int, float)) and not isinstance(other, bool): expect(Decimal(value) == Decimal(str(other)), f"CSV/JSON numeric mismatch: {row['state']}.{field}") else: expect(value == other, f"CSV/JSON mismatch: {row['state']}.{field}") a = number(row["gross_refrigerated_2023_total"]) b = number(row["gross_refrigerated_2025_total"]) change = b - a if a is not None and b is not None else None check_value(row["change_2023_2025"], change, 0, row["state"] + " change") check_value(row["change_2023_2025_pct"], percent(change, a), 2, row["state"] + " change percent") check_value(row["freezer_share_of_gross_2025_pct"], percent(row["gross_freezer_2025_total"], b), 1, row["state"] + " freezer share") check_value(row["usable_share_of_gross_2025_pct"], percent(row["usable_refrigerated_2025_total"], b), 1, row["state"] + " usable share") if a is not None and b is not None: expect(row["identical_in_both_surveys"] == ("yes" if a == b else "no"), f"Identical flag mismatch: {row['state']}") for year in (2023, 2025): expect(int(row[f"warehouses_{year}_public"]) + int(row[f"warehouses_{year}_private"]) == int(row[f"warehouses_{year}_total"]), f"Location split mismatch: {row['state']} {year}") location_totals: dict[str, dict[str, int]] = {} for year in (2023, 2025): location_totals[str(year)] = {} for category in ("public", "private", "total"): key = f"warehouses_{year}_{category}" total = sum(int(row[key]) for row in states) expect(total == int(national[key]), f"National location total mismatch: {year} {category}") location_totals[str(year)][category] = total comparable = [r for r in states if (number(r["gross_refrigerated_2023_total"]) or 0) > 0 and (number(r["gross_refrigerated_2025_total"]) or 0) > 0] declines = [r["state"] for r in comparable if required(r["change_2023_2025"]) < 0] identical = [r["state"] for r in comparable if required(r["change_2023_2025"]) == 0] withheld = [r["state"] for r in states if r["gross_refrigerated_2025_total"] == "D"] zero = [r["state"] for r in states if r["gross_refrigerated_2025_total"] == "0"] expect((len(comparable), len(declines), len(identical), len(withheld), len(zero)) == (37, 11, 7, 8, 4), "Comparison-group counts do not match the page") component_states = [r for r in states if (number(r["gross_cooler_2025_total"]) or 0) > 0 and (number(r["gross_freezer_2025_total"]) or 0) > 0] freezer_ranking = sorted(component_states, key=lambda r: required(r["gross_freezer_2025_total"]) / required(r["gross_refrigerated_2025_total"]), reverse=True) expect(len(freezer_ranking) == 25 and freezer_ranking[3]["state_abbr"] == "NE", "Disclosed freezer-share ranking does not match the page") max_component_difference = max(abs(required(r["gross_cooler_2025_total"]) + required(r["gross_freezer_2025_total"]) - required(r["gross_refrigerated_2025_total"])) for r in component_states) expect(max_component_difference <= 1, "Unexpected state component discrepancy") national_differences = {str(year): int(required(national[f"gross_refrigerated_{year}_total"]) - required(national[f"gross_cooler_{year}_total"]) - required(national[f"gross_freezer_{year}_total"])) for year in (2023, 2025)} expect(national_differences == {"2023": 9, "2025": 11}, "National source discrepancies differ from the documented cells") ga_sc_change = required(by_state["GA"]["change_2023_2025"]) + required(by_state["SC"]["change_2023_2025"]) national_change = required(national["change_2023_2025"]) history = read_csv(data_dir / HISTORY_FILE) expect([int(r["survey_year"]) for r in history] == [2017, 2019, 2021, 2023, 2025], "Unexpected national series dates") baseline = history[0] for r in history: check_value(r["public_share_pct"], percent(r["gross_public"], r["gross_total"]), 1, r["survey_year"] + " public share") check_value(r["private_share_pct"], percent(r["gross_private_semiprivate"], r["gross_total"]), 1, r["survey_year"] + " private share") for field, column in [("index_2017_total", "gross_total"), ("index_2017_public", "gross_public"), ("index_2017_private", "gross_private_semiprivate")]: check_value(r[field], percent(r[column], baseline[column]), 1, r["survey_year"] + " " + field) stock_rows = read_csv(data_dir / STOCK_FILE) for r in stock_rows: check_value(r["public_share_pct"], percent(r["public_warehouse_stocks_1000_lb"], r["stocks_all_warehouses_1000_lb"]), 1, r["commodity"] + " public stock share") stock_by_name = {r["commodity"]: r for r in stock_rows} # Keep all commodity lines separate: child categories overlap parent totals. total_stock = next(r for r in stock_rows if r["stocks_all_warehouses_1000_lb"] == "8514843") freezer_stock = next(r for r in stock_rows if r["stocks_all_warehouses_1000_lb"] == "6602546") cooler_stock = next(r for r in stock_rows if r["stocks_all_warehouses_1000_lb"] == "1912297") expect(required(freezer_stock["stocks_all_warehouses_1000_lb"]) + required(cooler_stock["stocks_all_warehouses_1000_lb"]) == required(total_stock["stocks_all_warehouses_1000_lb"]), "Freezer/cooler stock total mismatch") last = history[-1] return { "status": "All consistency and arithmetic checks passed", "scope": "Reproduces the accompanying dated files; does not re-verify source pages or underlying survey estimates", "source_verification_date": "2026-09-12", "capacity_unit": "1,000 cubic feet", "state_count": len(states), "location_totals": location_totals, "comparable_positive_capacity_states": len(comparable), "declining_states": declines, "identical_positive_capacity_states": identical, "withheld_2025_gross_states": withheld, "published_zero_2025_gross_states": zero, "georgia_south_carolina_combined_change": int(ga_sc_change), "national_net_change": int(national_change), "georgia_south_carolina_share_of_net_change_pct": output_number(percent(ga_sc_change, national_change)), "national_capacity_change_2023_2025_pct": output_number(percent(national_change, national["gross_refrigerated_2023_total"])), "national_capacity_change_2017_2025_pct": output_number(percent(required(last["gross_total"])-required(baseline["gross_total"]), baseline["gross_total"])), "public_capacity_change_2017_2025_pct": output_number(percent(required(last["gross_public"])-required(baseline["gross_public"]), baseline["gross_public"])), "private_capacity_change_2017_2025_pct": output_number(percent(required(last["gross_private_semiprivate"])-required(baseline["gross_private_semiprivate"]), baseline["gross_private_semiprivate"])), "national_published_total_minus_component_sum": national_differences, "state_max_absolute_component_difference": int(max_component_difference), "disclosed_freezer_share_panel_size": len(freezer_ranking), "freezer_share_ranking": [{"rank": i, "state": r["state"], "freezer_share_pct": output_number(percent(r["gross_freezer_2025_total"], r["gross_refrigerated_2025_total"]))} for i, r in enumerate(freezer_ranking, 1)], "state_percentage_change_display": {r["state"]: output_number(percent(r["change_2023_2025"], r["gross_refrigerated_2023_total"])) for r in states}, "public_stock_share_july_2026_pct": output_number(percent(total_stock["public_warehouse_stocks_1000_lb"],total_stock["stocks_all_warehouses_1000_lb"])), "stock_units": "1,000 pounds", "stock_observation_date": "2026-07-31", "stock_release_date": "2026-08-24", "stock_revision_status": "Initial estimates in that release; not updated to subsequent revisions" } def main() -> int: parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) parser.add_argument("--data-dir", type=Path, default=Path(__file__).resolve().parent, help="Directory containing the four accompanying data files") args = parser.parse_args() try: print(json.dumps(reproduce(args.data_dir), ensure_ascii=False, indent=2)) return 0 except (OSError, ValueError, KeyError, TypeError, ArithmeticError) as error: print(f"Reproduction failed: {error}", file=sys.stderr) return 1 if __name__ == "__main__": raise SystemExit(main())