diff --git a/AGENTS.md b/AGENTS.md index 0fa4c1f..e2c3bef 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -4,66 +4,90 @@ This project downloads population and geography data from the Census Bureau API ## Pipeline overview -The full pipeline is driven by `scripts/reproduce.sh`, which sources `scripts/pipeline_config.sh` for all configuration. Steps run in order: - -1. **`scripts/setup.sh`** — scaffolds the directory tree (`census_raw/`, `census_geographies/`, `study_area_sources/`, `study_areas/`, `outputs/`, etc.) -2. **`pipeline/download_population_tables.py`** — downloads decennial census race/ethnicity counts (TOTPOP, WHITE, BLACK, POC, etc.) via Census API; uses IPUMS/NHGIS extracts for 1980 and 1990 -3. **`pipeline/download_geographies.py`** — downloads TIGER/Line shapefiles (2000–2020 via Census API; 1980/1990 via IPUMS NHGIS) -4. **`pipeline/build_census_geographies.py`** — joins population tables to shapefiles, producing one attributed shapefile per year/level in `census_geographies/` -5. **`scripts/build_study_areas.sh`** → **`pipeline/build_study_areas.py`** — builds study area boundary polygons (e.g. CBSA outlines from county-component `.xls` files) into `study_areas/definitions/` -6. **`scripts/overlaps.sh`** → **`pipeline/overlaps.py`** — clips census geography shapefiles to each study area boundary (parallelized over years); outputs clipped shapefiles to `study_areas//` and coverage stats to `outputs//coverage_stats.csv` -7. **`pipeline/gen_duals.py`** — builds the dual adjacency graph from each clipped shapefile; contracts zero-population nodes and ensures full connectivity; outputs `*_orig.json` and `*_connected.json` alongside each shapefile -8. **`pipeline/calculate_metrics.py`** — computes ~80 segregation metrics per study area / year from each connected graph JSON; outputs one CSV row per area; errors logged to `outputs//metric_failures.csv` -9. **`pipeline/generate_figures.py`** — reads aggregated metric CSVs and produces publication figures +The full pipeline is driven by `scripts/reproduce.sh`. Configuration lives in `pipeline/config.yaml` and is loaded by `pipeline/config.py`. Steps run in order: + +1. **`scripts/setup.sh`** — scaffolds the directory tree +2. **`pipeline/download/download_population_tables.py`** — downloads decennial census race/ethnicity counts (TOTPOP, WHITE, BLACK, POC, etc.) via Census API; uses IPUMS/NHGIS extracts for 1980 and 1990 +3. **`pipeline/download/download_geographies.py`** — downloads TIGER/Line shapefiles (2000–2020 via Census API; 1980/1990 via IPUMS NHGIS) +4. **`pipeline/preprocessing/census_geographies.py`** — joins population tables to shapefiles, producing one attributed shapefile per state/year/level in `data/processed/census_geographies/` +5. **`pipeline/preprocessing/study_areas.py`** — builds study area boundary polygons (e.g. CBSA outlines from county-component `.xls` files) into `data/processed/study_area_definitions/` +6. **`pipeline/preprocessing/overlaps.py`** — clips census geography shapefiles to each study area boundary; outputs clipped shapefiles to `data/processed/clipped_geographies/` +7. **`pipeline/graphs.py`** — builds the dual adjacency graph from each clipped shapefile; drops zero-population nodes and ensures full connectivity; outputs `*_connected.json` files to `data/processed/dual_graphs/` +8. **`pipeline/metrics.py`** — computes ~80 segregation metrics per study area / year from each connected graph JSON; outputs one CSV row per area; errors logged to `outputs//metric_failures.csv` +9. **`pipeline/visualization/generate_figures.py`** — reads the metrics CSV and produces publication figures ## Configuration -All pipeline behavior is controlled by environment variables (with defaults in `scripts/pipeline_config.sh`): +All pipeline behavior is controlled by `pipeline/config.yaml`: -| Variable | Default | Options | +| Key | Example | Options | |---|---|---| -| `STUDY_AREA_TYPE` | `cbsa` | `cbsa`, `county` | -| `CENSUS_GEOGRAPHY_TYPE` | `tracts` | `tracts`, `block_groups`, `blocks`, `counties` | -| `CENSUS_GEOGRAPHY_YEARS` | `2020 2010 2000 1990 1980` | space-separated year list | -| `STUDY_AREA_VINTAGE` | `2020` | year | -| `RUN_NAME` | `_in_` | string | -| `RUN_OUTPUT_DIR` | `outputs/` | path | +| `study_area_type` | `cbsa` | `cbsa`, `max_city`, `max_county`, `county` | +| `census_geography_type` | `tracts` | `tracts`, `block_groups`, `blocks`, `counties` | +| `census_geography_years` | `[2020, 2010, 2000]` | list of years | +| `study_area_vintage` | `2020` | year | -For `STUDY_AREA_TYPE=cbsa`, a delineation file matching `list1_*_.xls` must exist in `study_area_sources/`. A Census API key and IPUMS API key are required for downloads. +For `study_area_type: cbsa`, a delineation file matching `list1_*.xls` must exist in `data/raw/study_area_sources/`. A Census API key and IPUMS API key are required for downloads. ## Important directories -| Directory | Contents | Notes | -|---|---|---| -| `study_area_sources/` | CBSA delineation `.xls` files | Required input; not generated | -| `census_raw/geographies/` | Raw downloaded TIGER/NHGIS shapefiles | Do not modify | -| `census_raw/population/` | Raw downloaded population CSV tables | Do not modify | -| `census_geographies/` | Merged population+shapefile per year/level (`_.shp`) | Generated; heavy | -| `study_areas/definitions/` | Study area boundary shapefiles | Generated | -| `study_areas//` | Clipped census shapefiles + dual graph JSONs per study area | Generated; heavy | -| `outputs//` | Metric CSVs (`white_black.csv`, `white_poc.csv`), `coverage_stats.csv`, `metric_failures.csv`, figures, `run.log` | Generated outputs | - -`_orig.json` = raw dual graph; `_connected.json` = fully connected, zero-pop nodes contracted (used for metrics). +| Directory | Contents | +|---|---| +| `data/raw/study_area_sources/` | CBSA delineation `.xls` files (required input) | +| `data/raw/geographies/` | Raw downloaded TIGER/NHGIS shapefiles | +| `data/raw/population/` | Raw downloaded population CSV tables | +| `data/processed/census_geographies/` | Population-attributed shapefiles per state/year/level | +| `data/processed/study_area_definitions/` | Study area boundary `.gpkg` + metadata `.json` | +| `data/processed/clipped_geographies/` | Census units clipped to each study area | +| `data/processed/dual_graphs/` | Adjacency graph JSONs (`*_connected.json`) | +| `outputs//` | Metric CSVs, `metric_failures.csv`, figures, `run.log` | ## Important commands -- **Setup:** `bash scripts/setup.sh` -- **Full reproduction:** `bash scripts/reproduce.sh` (takes many hours on full dataset) -- **Generate dual graphs:** `python pipeline/gen_duals.py ` -- **Calculate metrics:** `python pipeline/calculate_metrics.py ` -- **Generate figures:** `python pipeline/generate_figures.py --filename --prefix ` -- **Check overlaps:** `bash scripts/overlaps.sh` -- **Parse output:** `python pipeline/parse_output.py` (consumes CSV format) +All scripts use Typer and accept `--help`. Run from the repo root. + +```bash +# Full reproduction +bash scripts/setup.sh +bash scripts/reproduce.sh + +# Download +poetry run python pipeline/download/download_population_tables.py --level tracts --years "2020 2010 2000" +poetry run python pipeline/download/download_geographies.py --level tracts --years "2020 2010 2000" + +# Preprocessing +poetry run python pipeline/preprocessing/census_geographies.py --level tracts --years "2020 2010 2000" +poetry run python pipeline/preprocessing/study_areas.py --filename data/raw/study_area_sources/list1_march_2020.xls --study-area-type cbsa +poetry run python pipeline/preprocessing/overlaps.py \ + "data/processed/study_area_definitions/cbsa_*_march_2020.gpkg" \ + data/processed/clipped_geographies \ + --census-geography-type tracts \ + --census-geography-years "2020 2010 2000" \ + --definition-vintage march_2020 + +# Graphs and metrics +poetry run python pipeline/graphs.py \ + "data/processed/clipped_geographies/*/tracts_in_cbsa_*_march_2020_vintage.gpkg" +poetry run python pipeline/metrics.py \ + "data/processed/dual_graphs/*/tracts_in_cbsa_*_march_2020_vintage_connected.json" \ + BLACK WHITE TOTPOP outputs/tracts_in_cbsa/white_black.csv + +# Figures +poetry run python pipeline/visualization/generate_figures.py \ + --filename outputs/tracts_in_cbsa/white_black.csv \ + --prefix white_black_cbsa_tracts \ + --geography-type tracts \ + --study-area-type cbsa +``` ## Testing / Verification -- Run `pytest` to execute `tests/test_gen_duals.py` and `tests/test_pipeline_config.py` -- Run `python pipeline/gen_duals.py ...` on a small fixture for a quick sanity check +- Run `pytest` to execute the test suite under `pipeline/tests/` - Run `bash scripts/reproduce.sh` only when full reproduction is needed ## Outputs - Metric CSVs: `outputs//white_black.csv`, `outputs//white_poc.csv` -- Coverage stats: `outputs//coverage_stats.csv` -- Figures: under `outputs//` -- Output formats must not change. +- Run log: `outputs//run.log` +- Figures: `outputs//figures/` +- Output formats must not change. \ No newline at end of file diff --git a/README.md b/README.md index 3fe604e..6450fd7 100644 --- a/README.md +++ b/README.md @@ -10,81 +10,156 @@ capy-bara/ │ ├── raw/ # downloaded source files (gitignored) │ │ ├── geographies/ # TIGER/Line and NHGIS shapefiles │ │ ├── population/ # Census API / NHGIS population tables -│ │ └── ipums_extracts/ # IPUMS extracts (1980, 1990) -│ ├── interim/ # processed intermediates (gitignored) -│ │ ├── census_geographies/ # population-attributed shapefiles per year/level -│ │ ├── cbsas/ # CBSA definitions by decade -│ │ ├── study_areas/ # clipped shapefiles + dual graph JSONs per study area │ │ └── study_area_sources/ # CBSA delineation .xls files -│ └── outputs/ # pipeline run outputs -│ ├── tracts_in_cbsa/ -│ ├── block_groups_in_cbsa/ -│ ├── blocks_in_cbsa/ -│ └── cross_level_comparisons/ # figures comparing results across runs +│ └── processed/ # pipeline intermediates (gitignored) +│ ├── census_geographies/ # population-attributed shapefiles per year/level +│ ├── study_area_definitions/ # study area boundary .gpkg + metadata .json +│ ├── clipped_geographies/ # census units clipped to each study area +│ ├── dual_graphs/ # adjacency graph JSONs per study area +│ └── dropped_nodes/ # zero-population nodes removed from graphs +│ +├── outputs/ # pipeline run outputs (gitignored) +│ ├── tracts_in_cbsa/ # metrics CSVs + figures for this configuration +│ ├── block_groups_in_cbsa/ +│ └── cross_level_comparisons/ # figures comparing results across runs │ ├── pipeline/ # core pipeline modules +│ ├── config.py # config loader; prints shell exports when run directly +│ ├── config.yaml # pipeline configuration +│ ├── graphs.py # dual adjacency graph construction +│ ├── metrics.py # segregation metric calculations +│ ├── process_results.py # enriches metrics CSV with study area metadata │ ├── download/ # download_geographies.py, download_population_tables.py -│ ├── build/ # build_census_geographies.py, build_study_areas.py, -│ │ # overlaps.py, gen_duals.py, filter_cbsas.py -│ ├── metrics/ # calculate_metrics.py, parse_output.py -│ ├── viz/ # generate_figures.py -│ ├── utils/ # definitions.py -│ └── tests/ +│ ├── preprocessing/ # census_geographies.py, study_areas.py, overlaps.py +│ ├── visualization/ # generate_figures.py +│ └── utils/ # definitions.py, pipeline_log.py │ ├── experiments/ # hypothesis-testing experiments -│ ├── notebooks/ # scratch notebooks before a hypothesis becomes a script -│ ├── comparisons/ # cross-experiment analyses and figures -│ └── exp_/ # one folder per experiment: run.py + figures/ +│ └── / # one folder per experiment │ -├── working_paper_reproduction/ # materials for reproducing paper results -│ ├── notebooks/ # reproduction notebooks -│ ├── misc_analysis/ -│ └── figures/ +├── scripts/ # shell scripts +│ ├── reproduce.sh # full pipeline orchestration +│ └── setup.sh # scaffolds directory tree │ -├── scripts/ # shell scripts for running the pipeline └── archive/ # inactive code and old outputs ``` ## Pipeline overview -The full pipeline is driven by `scripts/reproduce.sh`, which sources `scripts/pipeline_config.sh` for all configuration. Steps run in order: +The full pipeline is driven by `scripts/reproduce.sh`. Configuration lives in `pipeline/config.yaml` and is loaded by `pipeline/config.py`. Steps run in order: 1. **`scripts/setup.sh`** — scaffolds the directory tree 2. **`pipeline/download/download_population_tables.py`** — downloads decennial census race/ethnicity counts (TOTPOP, WHITE, BLACK, POC, etc.) via Census API; uses IPUMS/NHGIS extracts for 1980 and 1990 3. **`pipeline/download/download_geographies.py`** — downloads TIGER/Line shapefiles (2000–2020 via Census API; 1980/1990 via IPUMS NHGIS) -4. **`pipeline/build/build_census_geographies.py`** — joins population tables to shapefiles, producing one attributed shapefile per year/level in `data/interim/census_geographies/` -5. **`scripts/build_study_areas.sh`** → **`pipeline/build/build_study_areas.py`** — builds study area boundary polygons (e.g. CBSA outlines from county-component `.xls` files) into `data/interim/study_areas/definitions/` -6. **`scripts/overlaps.sh`** → **`pipeline/build/overlaps.py`** — clips census geography shapefiles to each study area boundary (parallelized over years); outputs clipped shapefiles to `data/interim/study_areas//` and coverage stats to `data/outputs//coverage_stats.csv` -7. **`pipeline/build/gen_duals.py`** — builds the dual adjacency graph from each clipped shapefile; contracts zero-population nodes and ensures full connectivity; outputs `*_orig.json` and `*_connected.json` alongside each shapefile -8. **`pipeline/metrics/calculate_metrics.py`** — computes ~80 segregation metrics per study area / year from each connected graph JSON; outputs one CSV row per area; errors logged to `data/outputs//metric_failures.csv` -9. **`pipeline/viz/generate_figures.py`** — reads aggregated metric CSVs and produces publication figures - -`_orig.json` = raw dual graph; `_connected.json` = fully connected, zero-pop nodes contracted (used for metrics). +4. **`pipeline/preprocessing/census_geographies.py`** — joins population tables to shapefiles, producing one attributed shapefile per state/year/level in `data/processed/census_geographies/` +5. **`pipeline/preprocessing/study_areas.py`** — builds study area boundary polygons (e.g. CBSA outlines from county-component `.xls` files) into `data/processed/study_area_definitions/` +6. **`pipeline/preprocessing/overlaps.py`** — clips census geography shapefiles to each study area boundary; outputs clipped shapefiles to `data/processed/clipped_geographies/` +7. **`pipeline/graphs.py`** — builds the dual adjacency graph from each clipped shapefile; drops zero-population nodes and ensures full connectivity; outputs `*_connected.json` files to `data/processed/dual_graphs/` +8. **`pipeline/metrics.py`** — computes ~80 segregation metrics per study area / year from each connected graph JSON; outputs one CSV row per area; errors logged to `outputs//metric_failures.csv` +9. **`pipeline/process_results.py`** — enriches the metrics CSV with study area metadata (title, population) from the definition JSON files +10. **`pipeline/visualization/generate_figures.py`** — reads the metrics CSV and produces publication figures ## Configuration -All pipeline behavior is controlled by environment variables (with defaults in `scripts/pipeline_config.sh`): +All pipeline behavior is controlled by `pipeline/config.yaml`: -| Variable | Default | Options | +| Key | Example | Options | |---|---|---| -| `STUDY_AREA_TYPE` | `cbsa` | `cbsa`, `county` | -| `CENSUS_GEOGRAPHY_TYPE` | `tracts` | `tracts`, `block_groups`, `blocks`, `counties` | -| `CENSUS_GEOGRAPHY_YEARS` | `2020 2010 2000 1990 1980` | space-separated year list | -| `STUDY_AREA_VINTAGE` | `2020` | year | -| `RUN_NAME` | `_in_` | string | -| `RUN_OUTPUT_DIR` | `data/outputs/` | path | - -For `STUDY_AREA_TYPE=cbsa`, a delineation file matching `list1_*_.xls` must exist in `data/interim/study_area_sources/`. A Census API key and IPUMS API key are required for downloads. - -## Important commands - -- **Setup:** `bash scripts/setup.sh` -- **Full reproduction:** `bash scripts/reproduce.sh` (takes many hours on full dataset) -- **Generate dual graphs:** `python pipeline/build/gen_duals.py ` -- **Calculate metrics:** `python pipeline/metrics/calculate_metrics.py ` -- **Generate figures:** `python pipeline/viz/generate_figures.py --filename --prefix ` -- **Check overlaps:** `bash scripts/overlaps.sh` -- **Parse output:** `python pipeline/metrics/parse_output.py` +| `study_area_type` | `cbsa` | `cbsa`, `max_city`, `max_county`, `county` | +| `census_geography_type` | `tracts` | `tracts`, `block_groups`, `blocks`, `counties` | +| `census_geography_years` | `[2020, 2010, 2000]` | list of years | +| `study_area_vintage` | `2020` | year | + +For `study_area_type: cbsa`, a delineation file matching `list1_*.xls` must exist in `data/raw/study_area_sources/`. A Census API key and IPUMS API key are required for downloads. + +## Running the pipeline + +```bash +# scaffold directories, then run everything +bash scripts/setup.sh +bash scripts/reproduce.sh +``` + +The full run can take from a few minutes to many hours, depending on what level of geography you choose. Each step can also be run standalone, see below. + +## Pipeline scripts + +Run from the repo root with `poetry run python`. + +### `pipeline/config.py` +Prints shell export statements derived from `pipeline/config.yaml`. Used internally by `reproduce.sh`; useful for inspecting resolved config values. +```bash +poetry run python pipeline/config.py +``` + +### `pipeline/download/download_population_tables.py` +Downloads decennial census population tables (race, ethnicity, total) for a given geography level and set of years. +```bash +poetry run python pipeline/download/download_population_tables.py \ + --level tracts \ + --years "2020 2010 2000" +``` + +### `pipeline/download/download_geographies.py` +Downloads TIGER/Line shapefiles (2000–2020) or IPUMS/NHGIS shapefiles (1980–1990) for a given geography level. +```bash +poetry run python pipeline/download/download_geographies.py \ + --level tracts \ + --years "2020 2010 2000" +``` + +### `pipeline/preprocessing/census_geographies.py` +Joins downloaded population tables to shapefiles, writing one `.gpkg` per state/year into `data/processed/census_geographies/`. +```bash +poetry run python pipeline/preprocessing/census_geographies.py \ + --level tracts \ + --years "2020 2010 2000" +``` + +### `pipeline/preprocessing/study_areas.py` +Builds study area boundary files (`.gpkg` + `.json`) from the CBSA definition Excel file. One file pair per study area in `data/processed/study_area_definitions/`. +```bash +poetry run python pipeline/preprocessing/study_areas.py \ + --filename data/raw/study_area_sources/list1_march_2020.xls \ + --study-area-type cbsa +``` + +### `pipeline/preprocessing/overlaps.py` +Clips census geography units to each study area boundary. Writes one `.gpkg` per study area and year to the output directory. +```bash +poetry run python pipeline/preprocessing/overlaps.py \ + "data/processed/study_area_definitions/cbsa_*_march_2020.gpkg" \ + data/processed/clipped_geographies \ + --census-geography-type tracts \ + --census-geography-years "2020 2010 2000" \ + --definition-vintage march_2020 +``` + +### `pipeline/graphs.py` +Builds dual adjacency graphs from clipped shapefiles. Drops zero-population nodes and adds edges between any disconnected components. Writes `*_connected.json` files to `data/processed/dual_graphs/`. +```bash +poetry run python pipeline/graphs.py \ + "data/processed/clipped_geographies/*/tracts_in_cbsa_*_march_2020_vintage.gpkg" +``` + +### `pipeline/metrics.py` +Computes segregation metrics for each study area from connected graph JSONs. Arguments are the glob pattern, group columns, and output CSV path. +```bash +poetry run python pipeline/metrics.py \ + "data/processed/dual_graphs/*/tracts_in_cbsa_*_march_2020_vintage_connected.json" \ + BLACK WHITE TOTPOP \ + outputs/tracts_in_cbsa/white_black.csv +``` + +### `pipeline/visualization/generate_figures.py` +Reads a metrics CSV and writes figures to `outputs//figures/`. +```bash +poetry run python pipeline/visualization/generate_figures.py \ + --filename outputs/tracts_in_cbsa/white_black.csv \ + --prefix white_black_cbsa_tracts \ + --geography-type tracts \ + --study-area-type cbsa +``` ## Dependencies diff --git a/pipeline/config.py b/pipeline/config.py new file mode 100644 index 0000000..e0c47a1 --- /dev/null +++ b/pipeline/config.py @@ -0,0 +1,93 @@ +""" +Pipeline configuration loader. Reads pipeline/config.yaml and returns a plain dict of resolved config values. + +Usage from Python: + from pipeline.config import load_config + cfg = load_config() + +When run directly, prints shell export statements for use in shell scripts: + eval "$(poetry run python pipeline/config.py)" +""" +from __future__ import annotations +import glob +import sys +from pathlib import Path +import yaml + +CONFIG_FILE = Path(__file__).with_name("config.yaml") +REPO_ROOT = Path(__file__).parent.parent + +STUDY_AREA_TYPE_ALIASES = {"counties": "county", + "max_counties": "max_county", + "max_cities": "max_city"} +CENSUS_GEOGRAPHY_TYPE_ALIASES = {"tract": "tracts", + "block_group": "block_groups", + "block": "blocks", + "county": "counties"} + + +def load_config() -> dict: + """Load and validate pipeline/config.yaml. + + Returns a dict with snake_case keys and Python-native values + (list[str] for census_geography_years, Path for run_output_dir). + """ + with open(CONFIG_FILE) as f: + raw = yaml.safe_load(f) + + # study_area_type + study_area_type = STUDY_AREA_TYPE_ALIASES.get(str(raw["study_area_type"]), str(raw["study_area_type"])) + + # census_geography_type + census_geography_type = CENSUS_GEOGRAPHY_TYPE_ALIASES.get( + str(raw["census_geography_type"]), str(raw["census_geography_type"])) + + # census_geography_years + years = [str(y) for y in raw["census_geography_years"]] + if census_geography_type in ("block_groups", "blocks") and "1980" in years: + print(f"Warning: Skipping 1980 for census_geography_type={census_geography_type}. NHGIS does not publish 1980 block group or block boundary shapefiles.", + file=sys.stderr) + years = [y for y in years if y != "1980"] + + study_area_vintage = str(raw.get("study_area_vintage", "2020")) + + # other variables + study_area_definition_geography_type = "places" if study_area_type == "max_city" else "counties" + study_area_definition_geography_year = study_area_vintage + study_area_source_file = None + + if study_area_type in ("cbsa", "max_city", "max_county"): + source_pattern = f"list1_*{study_area_vintage}.xls" + matches = sorted(glob.glob(str(REPO_ROOT / "data" / "raw" / "study_area_sources" / source_pattern))) + if not matches: + raise FileNotFoundError(f"No study area source file found for study_area_type={study_area_type!r}, study_area_vintage={study_area_vintage!r}.") + study_area_source_file = matches[-1] + study_area_definition_vintage = Path(study_area_source_file).stem.removeprefix("list1_") + else: + study_area_definition_vintage = study_area_vintage + + study_area_definition_geographies = (f"data/processed/census_geographies/{study_area_definition_geography_type}/{study_area_definition_geography_year}_{study_area_definition_geography_type}_*.gpkg") + + run_output_dir = REPO_ROOT / "outputs" / f"{census_geography_type}_in_{study_area_type}" + + return {"study_area_type": study_area_type, + "census_geography_type": census_geography_type, + "census_geography_years": years, + "study_area_vintage": study_area_vintage, + "study_area_definition_geography_type": study_area_definition_geography_type, + "study_area_definition_geography_year": study_area_definition_geography_year, + "study_area_source_file": study_area_source_file, + "study_area_definition_vintage": study_area_definition_vintage, + "study_area_definition_geographies": study_area_definition_geographies, + "output_suffix": f"{study_area_type}_{census_geography_type}_{study_area_definition_vintage}", + "run_output_dir": run_output_dir} + + +if __name__ == "__main__": + import shlex + + cfg = load_config() + exports = {**{k.upper(): str(v) for k, v in cfg.items() if v is not None}, + "CENSUS_GEOGRAPHY_YEARS": " ".join(cfg["census_geography_years"])} + for key, value in exports.items(): + print(f"export {key}={shlex.quote(value)}") diff --git a/pipeline/config.yaml b/pipeline/config.yaml index 2f4f45e..6c379a6 100644 --- a/pipeline/config.yaml +++ b/pipeline/config.yaml @@ -3,11 +3,11 @@ # Any value can also be overridden at runtime by setting the corresponding # environment variable (same name, uppercased) before calling reproduce.sh. -# Valid values: cbsa, county +# Valid values: cbsa, county, max_city, max_county study_area_type: max_city # Valid values: tracts, block_groups, blocks, counties -census_geography_type: blocks +census_geography_type: tracts census_geography_years: - 2020 @@ -16,7 +16,4 @@ census_geography_years: - 1990 - 1980 -study_area_vintage: 2020 - -# For cbsa study areas, use counties. -study_area_definition_geography_type: places +study_area_vintage: 2020 \ No newline at end of file diff --git a/pipeline/download/download_geographies.py b/pipeline/download/download_geographies.py index 2934ebd..4f8007f 100644 --- a/pipeline/download/download_geographies.py +++ b/pipeline/download/download_geographies.py @@ -21,15 +21,19 @@ "User-Agent": "capy-bara geography downloader", } +# See https://www.census.gov/library/reference/code-lists/ansi/ansi-codes-for-states.html +# 50 states + DC + Puerto Rico (72). Territories 60/66/69/78 are omitted, they have no CBSA definitions. STATE_FIPS = [ "01", "02", "04", "05", "06", "08", "09", "10", "11", "12", "13", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "44", "45", "46", "47", "48", - "49", "50", "51", "53", "54", "55", "56", - "60", "66", "69", "72", "78", + "49", "50", "51", "53", "54", "55", "56", "72", ] +# "census" values match the Census Bureau Data API's geography `for=` parameter: https://api.census.gov/data/2020/dec/pl/geography.json +# "nhgis_levels" values match the `geographicLevel` field in the IPUMS NHGIS shapefiles metadata API: https://developer.ipums.org/docs/v2/apiprogram/apis/nhgis/ +# Multiple strings per level reflect label changes across decades (e.g. tracts were called Block Numbering Areas before 1990 standardization). LEVELS = { "county": {"label": "counties", "nhgis_levels": ["county"], "census": "county"}, "counties": {"label": "counties", "nhgis_levels": ["county"], "census": "county"}, @@ -331,18 +335,13 @@ def select_shapefile_names(matches: List[dict], year: int, level: str) -> List[s selected = {} for shapefile in sorted( matches, - key=lambda item: (basis_year(item), int(item.get("sequence") or 0)), - ): + key=lambda item: (basis_year(item), int(item.get("sequence") or 0))): selected[shapefile_key(shapefile)] = shapefile return [shapefile["name"] for shapefile in selected.values()] -def nhgis_shapefiles( - client: IpumsApiClient, - year: int, - level: str, -) -> List[str]: +def nhgis_shapefiles(client: IpumsApiClient, year: int, level: str) -> List[str]: config = LEVELS[level] matches = [] @@ -363,11 +362,7 @@ def nhgis_shapefiles( return names -def fetch_nhgis( - year: int, - level: str, - work_dir: Path, -) -> Path: +def fetch_nhgis(year: int, level: str, work_dir: Path) -> Path: api_key = require_env("IPUMS_API_KEY") client = IpumsApiClient(api_key) config = LEVELS[level] @@ -376,8 +371,7 @@ def fetch_nhgis( extract = AggregateDataExtract( collection="nhgis", description=f"{year} {config['label']} shapefiles", - shapefiles=[Shapefile(name) for name in shapefiles], - ) + shapefiles=[Shapefile(name) for name in shapefiles]) submitted = client.submit_extract(extract) client.wait_for_extract(submitted, timeout=10800) @@ -388,8 +382,7 @@ def fetch_nhgis( zip_files = sorted( set(work_dir.glob("*.zip")) - existing_zips, - key=lambda path: path.stat().st_mtime, - ) + key=lambda path: path.stat().st_mtime) if not zip_files: raise FileNotFoundError(f"No new NHGIS zip downloaded to {work_dir}") @@ -410,27 +403,20 @@ def fetch_census(year: int, level: str, output_dir: Path) -> Path: if "{state}" in template: urls = [template.format(state=state) for state in STATE_FIPS] + n_existing = 0 for url in urls: if census_output_exists(url, output_path): - print(f"Skipping existing geography for {url}", flush=True) + n_existing += 1 continue zip_path = download(url, zip_dir / url.split("/")[-1]) extract_zip(zip_path, output_path) + if n_existing: + print(f"Skipped {n_existing} already downloaded {config['label']} shapefiles for {year}", flush=True) return output_path -def main( - level: str = typer.Option( - "tracts", - help="tracts, block_groups, blocks, or counties", - ), - years: Optional[str] = typer.Option(None, "--years", help="Space- or comma-separated years."), - year_values: Optional[List[int]] = typer.Option(None, "--year", "-y"), - output_dir: Path = typer.Option(OUTPUT_DIR), - work_dir: Path = typer.Option(Path("data/raw/geographies/ipums_geography_extracts")), - env_file: Path = typer.Option(Path(".env")), -) -> None: +def main(level: str = typer.Option("tracts", help="tracts, block_groups, blocks, or counties"), years: Optional[str] = typer.Option(None, "--years", help="Space- or comma-separated years."), year_values: Optional[List[int]] = typer.Option(None, "--year", "-y"), output_dir: Path = typer.Option(OUTPUT_DIR), work_dir: Path = typer.Option(Path("data/raw/geographies/ipums_geography_extracts")), env_file: Path = typer.Option(Path(".env"))) -> None: load_dotenv(env_file) if level not in LEVELS: @@ -444,29 +430,19 @@ def main( print(f"Skipping {year} places: only 2020 is used in the pipeline.") continue if year == 1980 and LEVELS[level]["label"] == ("block_groups", "blocks"): - print( - f"Skipping 1980 {level_label}: NHGIS does not publish 1980 block group " - "or block boundary shapefiles. These were not standardized as nationwide " - "geographic units until 1990.", - flush=True, - ) + print(f"Skipping 1980 {level_label}: NHGIS does not publish 1980 block group or block boundary shapefiles. These were not standardized as nationwide geographic units until 1990.", flush=True) continue if year in (1980, 1990): year_work_dir = work_dir / str(year) / level_label path = existing_nhgis_zip(year_work_dir) if path is None: - path = fetch_nhgis( - year, - level, - year_work_dir, - ) + path = fetch_nhgis(year, level, year_work_dir) else: - print(f"Skipping existing NHGIS geography extract {path}", flush=True) + print(f"Skipping existing NHGIS geography {path}", flush=True) elif year in CENSUS_TIGER_URLS: path = fetch_census(year, level, output_dir) else: raise ValueError(f"No geography source is configured for {year}") - print(path) if __name__ == "__main__": diff --git a/pipeline/download/download_population_tables.py b/pipeline/download/download_population_tables.py index 8b4f031..9c26392 100644 --- a/pipeline/download/download_population_tables.py +++ b/pipeline/download/download_population_tables.py @@ -15,13 +15,14 @@ NHGIS_EXTRACTS_DIR = Path("data/raw/population/ipums_population_extracts") DEFAULT_YEARS = [1980, 1990, 2000, 2010, 2020] +# See https://www.census.gov/library/reference/code-lists/ansi/ansi-codes-for-states.html +# 50 states + DC + Puerto Rico (72). Territories 60/66/69/78 are omitted, they have no CBSA definitions. STATES = [ "01", "02", "04", "05", "06", "08", "09", "10", "11", "12", "13", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "44", "45", "46", "47", "48", - "49", "50", "51", "53", "54", "55", "56", "72", -] + "49", "50", "51", "53", "54", "55", "56", "72"] LEVELS = { "tract": { @@ -291,15 +292,9 @@ def fetch_census_state_rows( try: return table.get(variables, geo=geo) except CensusException as exc: - if progress: - print( - " state-wide block group query failed; " - f"falling back to county-by-county: {exc}", - flush=True, - ) + print(f"state {state}: state-wide block group query failed; falling back to county-by-county: {exc}", flush=True) return fetch_county_scoped_census_rows( - table, variables, config, state, progress=progress - ) + table, variables, config, state, progress=progress) return table.get(variables, geo=geo) @@ -320,13 +315,10 @@ def fetch_census(year: int, level: str, states: List[str], output_dir: Path) -> variables = ["NAME"] + list(columns) + print(f"Fetching {year} {config['label']} population ({len(states)} states)", flush=True) rows = [] for state in states: - print( - f"Fetching {year} {config['label']} population for state {state}", - flush=True, - ) - rows.extend(fetch_census_state_rows(table, variables, config, state)) + rows.extend(fetch_census_state_rows(table, variables, config, state, progress=False)) if not rows: raise ValueError(f"Census returned no rows for {year} {config['label']}.") @@ -446,14 +438,12 @@ def main( output_path = output_dir / f"nhgis_{year}_{config['label']}.csv" if output_path.exists(): print(f"Skipping existing {output_path}", flush=True) - print(output_path) continue path = fetch_nhgis(year, level, output_dir, work_dir / str(year)) elif year in CENSUS_COLUMNS: output_path = output_dir / f"census_{year}_{config['label']}.csv" if output_path.exists(): print(f"Skipping existing {output_path}", flush=True) - print(output_path) continue path = fetch_census(year, level, state_fips, output_dir) else: diff --git a/pipeline/graphs.py b/pipeline/graphs.py index d509ad6..8fe8e10 100644 --- a/pipeline/graphs.py +++ b/pipeline/graphs.py @@ -1,75 +1,109 @@ +""" +Create dual graphs of given geometries and save: (1) the original dual graph (2) an edited graph, dropping 0-population nodes and connecting disconnected components. +""" + +import glob import geopandas as gpd import pandas as pd import typer import warnings import gerrychain import networkx as nx +from concurrent.futures import ProcessPoolExecutor +from functools import partial from shapely.strtree import STRtree -from pathlib import Path +from pathlib import Path CONTRACTION_POP_COLS = ("WHITE", "BLACK") -POPULATION_SUM_COLS = ("WHITE", "BLACK", "TOTPOP", "POC") - -def main(filename: str, output_orig: str, output_connected: str, attr: str = "GISJOIN", pop_col: str = "TOTPOP"): - shp = gpd.read_file(filename) - shp = shp.to_crs("esri:102003") #so distances are in meters - # ignore warning about NA values in population columns. These are fine. - warnings.filterwarnings("ignore", message=".*NA values found in column.*") - - right = filename.split("_in_")[1] +def main(input_glob: str, output_base_dir: str = "data/processed/dual_graphs", workers: int = 6, attr: str = "GISJOIN"): + gpkg_files = sorted(glob.glob(input_glob)) + if not gpkg_files: + raise + + worker = partial(_process_file, output_base_dir=output_base_dir, attr=attr) + with ProcessPoolExecutor(max_workers=workers) as pool: + results = list(pool.map(worker, gpkg_files)) + + # infer run output dir from the first matched file's stem: + # e.g. "tracts_in_max_city_35620_2020_vintage" to outputs/tracts_in_max_city/ + stem = Path(gpkg_files[0]).stem + census_geography_type = stem.split("_in_", 1)[0] + right_parts = stem.split("_in_", 1)[1].split("_") + study_area_type = f"max_{right_parts[1]}" if right_parts[0] == "max" else right_parts[0] + dropped_nodes_dir = Path("outputs") / f"{census_geography_type}_in_{study_area_type}" / "dropped_nodes" + + # aggregate dropped nodes by year and write one gpkg per year + by_year = {} + for year, dropped_gdf in results: + if dropped_gdf is not None: + by_year.setdefault(year, []).append(dropped_gdf) + for year, gdfs in by_year.items(): + combined = gpd.GeoDataFrame(pd.concat(gdfs, ignore_index=True), crs=gdfs[0].crs) + dropped_nodes_dir.mkdir(parents=True, exist_ok=True) + combined.to_file(dropped_nodes_dir / f"dropped_nodes_{year}.gpkg", driver="GPKG") + + +def _process_file(gpkg: str, output_base_dir: str, attr: str = "GISJOIN"): + # derive output paths from the filename: + year = Path(gpkg).parent.name + stem = Path(gpkg).stem + out_dir = Path(output_base_dir) / year + out_dir.mkdir(parents=True, exist_ok=True) + + # read and reproject + geofile = gpd.read_file(gpkg) + geofile = geofile.to_crs("esri:102003") # so distances are in meters + warnings.filterwarnings("ignore", message=".*NA values found in column.*") # some fields were introduced in 2000, so they're NA in earlier years. It's expected. + warnings.filterwarnings("ignore", message=".*Found islands.*") # degree-0 nodes are handled explicitly by connect_components. + + # extract area code from the filename + right = gpkg.split("_in_")[1] parts = right.split("_") - area_code = parts[2] if parts[0] == "max" else parts[1] #when the geography is - #cbsa, parts[1] is the code, - #otherwise parts[2] is for - #"max_city" and "max_county" - + area_code = parts[2] if parts[0] == "max" else parts[1] + if geofile.crs is None: + raise ValueError(f"{gpkg} has no CRS defined.") - if shp.crs is None: - raise ValueError(f"{filename} has no CRS defined. Please define a CRS before proceeding.") - - # Compute centroids in esri:102003 so distances are in meters. - centroids = shp.geometry.centroid - + # build a dual graph try: - graph = gerrychain.Graph.from_geodataframe(shp) + graph = gerrychain.Graph.from_geodataframe(geofile) except: - shp["geometry"] = shp["geometry"].buffer(0) - graph = gerrychain.Graph.from_geodataframe(shp) + geofile["geometry"] = geofile["geometry"].buffer(0) + graph = gerrychain.Graph.from_geodataframe(geofile) - for idx in shp.index: + # attach centroid coordinates to each node + centroids = geofile.geometry.centroid + for idx in geofile.index: graph.nodes[idx]["centroid_x"] = centroids.loc[idx].x graph.nodes[idx]["centroid_y"] = centroids.loc[idx].y - graph.to_json(output_orig) - - connected_graph = connect_components(shp, graph, attr) + graph.to_json(str(out_dir / f"{stem}_orig.json")) - zero_nodes = [] + # create an edited version of the graph: + # if the graph has disconnected components, add an edge across the nearest pair of geometries + connected_graph, n_edges_added = connect_components(geofile, graph, attr) + # remove 0-population nodes and their edges + dropped_indices = [] while len(connected_graph.nodes()) != 0 and has_zero_nodes(connected_graph): node_count = len(connected_graph.nodes()) - connected_graph, dropped_nodes = contract_zero_nodes(connected_graph) + connected_graph, dropped = drop_zero_nodes(connected_graph) if len(connected_graph.nodes()) == node_count: - print("No more zero nodes to contract, but graph still has zero nodes. Remaining nodes:", connected_graph.nodes()) break - for _, gisjoin in dropped_nodes: - zero_nodes.append((area_code, gisjoin)) + dropped_indices.extend(n for n, _ in dropped) - connected_graph.to_json(output_connected) + if n_edges_added > 0 or len(dropped_indices) > 0: + print(f"{stem}: +{n_edges_added} edges, {len(dropped_indices)} zero-pop nodes dropped", flush=True) + + connected_graph.to_json(str(out_dir / f"{stem}_connected.json")) - year = Path(filename).parent.name - stem = Path(filename).stem - dropped_dir = Path("data/processed/dropped_nodes") / year - dropped_dir.mkdir(parents=True, exist_ok=True) - - if len(zero_nodes) != 0: - df_zero_nodes = pd.DataFrame(zero_nodes, columns=["area_code", "id"]) - df_zero_nodes.to_csv( - dropped_dir / f"{stem}.csv", index=False - ) + if dropped_indices: + dropped_gdf = geofile.loc[dropped_indices].copy() + dropped_gdf["area_code"] = area_code + return year, dropped_gdf + return year, None def int_attr(attrs, col: str) -> int: @@ -85,19 +119,13 @@ def node_contraction_population(graph: gerrychain.Graph, node) -> int: def has_zero_nodes(graph: gerrychain.Graph): for node in graph.nodes(): - if node_contraction_population(graph, node) == 0: + node_contraction_population = sum(int_attr(graph.nodes[node], col) for col in CONTRACTION_POP_COLS) + if node_contraction_population == 0: return True return False -def add_population_attrs(graph: gerrychain.Graph, target, source): - for col in POPULATION_SUM_COLS: - graph.nodes[target][col] = int_attr(graph.nodes[target], col) + int_attr( - graph.nodes[source], col - ) - - -def contract_zero_nodes(graph: gerrychain.Graph): +def drop_zero_nodes(graph: gerrychain.Graph): zero_nodes = [n for n in graph.nodes() if node_contraction_population(graph, n) == 0] dropped_nodes = [(n, graph.nodes[n].get("GISJOIN", n)) for n in zero_nodes] @@ -106,24 +134,10 @@ def contract_zero_nodes(graph: gerrychain.Graph): return (graph, dropped_nodes) -def select_geom(shp: gpd.GeoDataFrame, geoid: str, attr: str = "GISJOIN"): - filtered_geoms = shp[shp[attr] == geoid] - return filtered_geoms.iloc[0]["geometry"] - - -def distance(shp: gpd.GeoDataFrame, geoid_1: str, geoid_2: str, attr: str = "GISJOIN"): - geom_1 = select_geom(shp, geoid_1, attr) - geom_2 = select_geom(shp, geoid_2, attr) - return geom_1.distance(geom_2) - - -def connect_components(shp: gpd.GeoDataFrame, graph: gerrychain.Graph, attr: str = "GISJOIN"): - geom_by_geoid = dict(zip(shp[attr], shp.geometry)) +def connect_components(geofile: gpd.GeoDataFrame, graph: gerrychain.Graph, attr: str = "GISJOIN"): + geom_by_geoid = dict(zip(geofile[attr], geofile.geometry)) + n_added = 0 while nx.algorithms.components.number_connected_components(graph) != 1: - print( - "Connected components:", - nx.algorithms.components.number_connected_components(graph), - ) cc = list(nx.connected_components(graph))[:2] assert len(cc) == 2 cc_geoids = [] @@ -145,11 +159,8 @@ def connect_components(shp: gpd.GeoDataFrame, graph: gerrychain.Graph, attr: str component_geoms = [geom_by_geoid[geoid] for geoid in cc_geoids[0]] island_geoms = [geom_by_geoid[geoid] for geoid in cc_geoids[1]] tree = STRtree(component_geoms) - pairs, distances = tree.query_nearest( - island_geoms, - return_distance=True, - all_matches=False, - ) + pairs, distances = tree.query_nearest(island_geoms, return_distance=True, + all_matches=False) assert len(distances) > 0 best_index = min(range(len(distances)), key=lambda index: distances[index]) island_index = pairs[0][best_index] @@ -158,9 +169,9 @@ def connect_components(shp: gpd.GeoDataFrame, graph: gerrychain.Graph, attr: str assert min_pair is not None graph.add_edge(geoid_node_mapping[min_pair[0]], geoid_node_mapping[min_pair[1]]) - print("Edge added:", min_pair) + n_added += 1 - return graph + return graph, n_added if __name__ == "__main__": diff --git a/pipeline/metrics.py b/pipeline/metrics.py index 2aa6330..bf0bff5 100644 --- a/pipeline/metrics.py +++ b/pipeline/metrics.py @@ -1,7 +1,9 @@ import geopandas as gpd import typer import os +import csv import glob +import warnings import gerrychain import networkx as nx import matplotlib.pyplot as plt @@ -12,6 +14,27 @@ import scipy.sparse import numpy as np import traceback +from concurrent.futures import ProcessPoolExecutor +from functools import partial +from pathlib import Path + + +def main(input_glob: str, x_col: str, y_col: str, tot_col: str, output: Path, workers: int = 6): + files = sorted(glob.glob(input_glob)) + worker = partial(_process_file, x_col=x_col, y_col=y_col, tot_col=tot_col) + n_ok = 0 + n_failed = 0 + with open(output, "w") as f: + f.write(build_headers(x_col, y_col, tot_col) + "\n") + with ProcessPoolExecutor(max_workers=workers) as pool: + for row in pool.map(worker, files): + if row is not None: + f.write(row + "\n") + n_ok += 1 + else: + n_failed += 1 + print(f"Metrics calculations: {n_ok} processes succeeded, {n_failed} failed. Output: {output}", flush=True) + def study_area_code_from_filename(filename: str) -> str: output_stem = os.path.basename(filename) @@ -52,37 +75,42 @@ def build_headers(x_col: str, y_col: str, tot_col: str) -> str: return ",".join(keys) -def main( - filename: str, x_col: str, y_col: str, tot_col: str, headers_only: bool = False -): - if headers_only: - print(build_headers(x_col, y_col, tot_col)) - return - try: - run_metrics(filename, x_col, y_col, tot_col) - except ZeroDivisionError as e: - metric_failures_file = os.environ.get( - "METRIC_FAILURES_FILE", "outputs/metric_failures.csv" - ) - metric_failures_dir = os.path.dirname(metric_failures_file) - if metric_failures_dir: - os.makedirs(metric_failures_dir, exist_ok=True) - with open(metric_failures_file, "a+") as f: - f.seek(0, os.SEEK_END) - if f.tell() == 0: - print("filename,cbsa_code,error", file=f) - print(f"{filename},{study_area_code_from_filename(filename)},{e}", file=f) - print(filename, e, file=sys.stderr) +FAILURE_FIELDNAMES = ["filename", "study_area_code", "x_col", "y_col", "tot_col", "error_message"] + + +def write_failure(filename: str, x_col: str, y_col: str, tot_col: str, exc: Exception) -> None: + metric_failures_file = os.environ.get("METRIC_FAILURES_FILE", "outputs/metric_failures.csv") + failures_dir = os.path.dirname(metric_failures_file) + os.makedirs(failures_dir, exist_ok=True) + row = {"filename": filename, "study_area_code": study_area_code_from_filename(filename), + "x_col": x_col, "y_col": y_col, "tot_col": tot_col, "error_message": str(exc)} + + # write header once + write_header = not os.path.exists(metric_failures_file) or os.path.getsize(metric_failures_file) == 0 + with open(metric_failures_file, "a", newline="") as f: + writer = csv.DictWriter(f, fieldnames=FAILURE_FIELDNAMES) + if write_header: + writer.writeheader() + writer.writerow(row) + + +def _process_file(filename: str, x_col: str, y_col: str, tot_col: str): + try: + return run_metrics(filename, x_col, y_col, tot_col) + except Exception as e: + write_failure(filename, x_col, y_col, tot_col, e) + print(f"FAILED {filename}: {type(e).__name__}: {e}", file=sys.stderr) + return None def run_metrics(filename: str, x_col: str, y_col: str, tot_col: str): + warnings.filterwarnings("ignore", message=".*Found islands.*") # degree-0 nodes are handled by connect_components in graphs.py. graph = gerrychain.Graph.from_json(filename) for node in graph.nodes(): graph.nodes[node]["white_plus_black"] = ( - int(graph.nodes[node][x_col]) + int(graph.nodes[node][y_col]) - ) + int(graph.nodes[node][x_col]) + int(graph.nodes[node][y_col])) capy_metrics = {} capy_metrics["filename"] = filename @@ -153,7 +181,7 @@ def run_metrics(filename: str, x_col: str, y_col: str, tot_col: str): capy_metrics["total_nodes"] = len(graph.nodes()) capy_metrics["total_edges"] = len(graph.edges()) - print(",".join(map(str, list(capy_metrics.values())))) + return ",".join(map(str, list(capy_metrics.values()))) def angle_1(graph: gerrychain.Graph, x_col: str, y_col: str, lam: float = 1) -> float: diff --git a/pipeline/preprocessing/census_geographies.py b/pipeline/preprocessing/census_geographies.py index fc1b1d7..d173fde 100644 --- a/pipeline/preprocessing/census_geographies.py +++ b/pipeline/preprocessing/census_geographies.py @@ -1,3 +1,7 @@ +""" +This script takes previously downloaded Census/IPUMS geography and population files and creates a .gpkg file with both population and geography information. +""" + import re import tempfile import zipfile @@ -20,8 +24,7 @@ "tract": 6, "block_group": 1, "block": 4, - "place": 5, -} + "place": 5} PART_COLUMNS = { "state": ["STATEFP", "STATEFP20", "STATEFP10", "STATEFP00"], @@ -38,8 +41,7 @@ "tract": "tract", "block_group": "block group", "block": "block", - "place": "place", -} + "place": "place"} LEVELS = { "county": {"label": "counties", "parts": ("state", "county"), "width": 5}, @@ -142,13 +144,20 @@ def is_nhgis_shapefile_for_level(path: Path, year: str, level_label: str) -> boo if level_label == "tracts": return is_original_tract_family_shapefile(path, year) - if is_conflated_nhgis_path(path) or is_county_sidecar_nhgis_path(path): + if is_county_sidecar_nhgis_path(path): return False text = clean_filename(path.stem) compact = text.replace(" ", "") if level_label == "counties": + # County shapefiles may only exist in conflated form for older years + # (e.g. 1980 NHGIS only provides US_county_1980_conflated.shp). + # Accept them — the conflated filter is only meaningful for tracts. return "county" in text and "tract" not in text + + if is_conflated_nhgis_path(path): + return False + if level_label == "block_groups": return is_block_group_name(path) if level_label == "blocks": @@ -245,59 +254,40 @@ def read_nhgis_1990_population(df: pd.DataFrame, path: Path) -> pd.DataFrame: require_columns(df, ["GISJOIN", "STATEA", "COUNTYA"] + race_cols, path) gisjoin = df["GISJOIN"].astype(str) return pd.DataFrame( - { - "JOIN_KEY": gisjoin, - "GISJOIN": gisjoin, - "STATEFP": df["STATEA"].str.zfill(2), - "COUNTYFP": df["COUNTYA"].str.zfill(3), - "WHITE": to_int(df["ET2001"]), - "BLACK": to_int(df["ET2002"]), - "TOTPOP": sum(to_int(df[col]) for col in race_cols), - } - ) + {"JOIN_KEY": gisjoin, + "GISJOIN": gisjoin, + "STATEFP": df["STATEA"].str.zfill(2), + "COUNTYFP": df["COUNTYA"].str.zfill(3), + "WHITE": to_int(df["ET2001"]), + "BLACK": to_int(df["ET2002"]), + "TOTPOP": sum(to_int(df[col]) for col in race_cols)}) -def read_census_population( - df: pd.DataFrame, - path: Path, - year: int, - level_label: str, -) -> pd.DataFrame: +def read_census_population(df: pd.DataFrame, path: Path, year: int, level_label: str) -> pd.DataFrame: config = LEVELS[level_label] part_columns = [POPULATION_PART_COLUMNS[part] for part in config["parts"]] require_columns( df, part_columns + ["TOTPOP", "NH_WHITE", "NH_BLACK"], - path, - ) + path) - parts = [ - normalize_part(df[POPULATION_PART_COLUMNS[part]], part, year) - for part in config["parts"] - ] + parts = [normalize_part(df[POPULATION_PART_COLUMNS[part]], part, year) for part in config["parts"]] join_key = parts[0] for part in parts[1:]: join_key = join_key + part - pop = pd.DataFrame( - { + pop = pd.DataFrame({ "JOIN_KEY": join_key, "GISJOIN": "G" + join_key, "STATEFP": parts[0], "COUNTYFP": parts[1] if len(parts) > 1 else "", "WHITE": to_int(df["NH_WHITE"]), "BLACK": to_int(df["NH_BLACK"]), - "TOTPOP": to_int(df["TOTPOP"]), - } - ) + "TOTPOP": to_int(df["TOTPOP"])}) return pop -def read_population( - year: int, - population_dir: Path, - level_label: str, -) -> pd.DataFrame: +def read_population(year: int, population_dir: Path, level_label: str) -> pd.DataFrame: if year in (1980, 1990): path = population_dir / f"nhgis_{year}_{level_label}.csv" else: @@ -330,10 +320,7 @@ def geography_part(gdf: gpd.GeoDataFrame, part: str) -> pd.Series: return gdf[col].astype(str).str.zfill(PART_WIDTHS[part]) -def standardize_census_geography( - gdf: gpd.GeoDataFrame, - level_label: str, -) -> gpd.GeoDataFrame: +def standardize_census_geography(gdf: gpd.GeoDataFrame, level_label: str) -> gpd.GeoDataFrame: config = LEVELS[level_label] parts = [geography_part(gdf, part) for part in config["parts"]] join_key = parts[0] @@ -348,17 +335,20 @@ def standardize_census_geography( return gdf -def read_census_geography(year: int, geographies_dir: Path, level_label: str) -> gpd.GeoDataFrame: +def read_census_geography(year, geographies_dir, level_label): shape_dir = geographies_dir / f"census_{year}_{level_label}" paths = sorted(path for path in shape_dir.glob("*.shp") if path.is_file()) if not paths: raise FileNotFoundError(f"No Census {level_label} shapefiles found in {shape_dir}") - - frames = [gpd.read_file(path) for path in paths] - crs = frames[0].crs - gdf = gpd.GeoDataFrame(pd.concat(frames, ignore_index=True), crs=crs) - del frames # free the list of frames now that concat is done. crucial for blocks to free the memory - return standardize_census_geography(gdf, level_label) + for path in paths: + gdf = gpd.read_file(path) + gdf = standardize_census_geography(gdf, level_label) + yield gdf.to_crs(TARGET_CRS) + # frames = [gpd.read_file(path) for path in paths] + # crs = frames[0].crs + # gdf = gpd.GeoDataFrame(pd.concat(frames, ignore_index=True), crs=crs) + # del frames # free the list of frames now that concat is done. crucial for blocks to free the memory + # return standardize_census_geography(gdf, level_label) def nested_shapefile_paths(outer_zip: Path, tmp_dir: Path) -> List[Path]: @@ -382,20 +372,13 @@ def nested_shapefile_paths(outer_zip: Path, tmp_dir: Path) -> List[Path]: return shp_paths -def read_nested_nhgis_shapefile( - outer_zip: Path, - year: int, - level_label: str, -) -> Optional[gpd.GeoDataFrame]: +def read_nested_nhgis_shapefile(outer_zip: Path, year: int, level_label: str) -> Optional[gpd.GeoDataFrame]: year_label = str(year) with tempfile.TemporaryDirectory() as tmp_name: tmp_dir = Path(tmp_name) shp_paths = nested_shapefile_paths(outer_zip, tmp_dir) - matches = [ - path - for path in shp_paths - if is_nhgis_shapefile_for_level(path, year_label, level_label) - ] + matches = [path for path in shp_paths + if is_nhgis_shapefile_for_level(path, year_label, level_label)] if not matches: return None @@ -415,10 +398,7 @@ def read_nested_nhgis_shapefile( return gdf -def first_existing_series( - gdf: gpd.GeoDataFrame, - candidates: List[str], -) -> Optional[pd.Series]: +def first_existing_series(gdf: gpd.GeoDataFrame, candidates: List[str]) -> Optional[pd.Series]: for col in candidates: if col in gdf.columns: return gdf[col].astype("string") @@ -439,11 +419,7 @@ def state_county_series(gdf: gpd.GeoDataFrame) -> tuple[pd.Series, pd.Series]: raise ValueError("Missing state/county identifier columns.") -def nhgis_extract_dirs( - geographies_dir: Path, - year: int, - level_label: str, -) -> List[Path]: +def nhgis_extract_dirs(geographies_dir: Path, year: int, level_label: str) -> List[Path]: base_dir = geographies_dir / "ipums_geography_extracts" / str(year) dirs = [base_dir / level_label] if base_dir not in dirs: @@ -451,11 +427,7 @@ def nhgis_extract_dirs( return dirs -def read_nhgis_geography( - year: int, - geographies_dir: Path, - level_label: str, -) -> gpd.GeoDataFrame: +def read_nhgis_geography(year: int, geographies_dir: Path, level_label: str) -> gpd.GeoDataFrame: extract_dirs = nhgis_extract_dirs(geographies_dir, year, level_label) paths = [] for extract_dir in extract_dirs: @@ -481,22 +453,16 @@ def read_nhgis_geography( gdf["COUNTYFP"] = county return gdf - raise ValueError( - f"No {level_label} NHGIS shapefile found in " + raise ValueError(f"No {level_label} NHGIS shapefile found in " f"{', '.join(str(path) for path in extract_dirs)}. " "Rerun download_geographies.py for this year and level after updating " - "the NHGIS selection." - ) + "the NHGIS selection.") TARGET_CRS = "esri:102003" # USA Contiguous Albers Equal Area Conic; meters -def read_geography( - year: int, - geographies_dir: Path, - level_label: str, -) -> gpd.GeoDataFrame: +def read_geography(year: int, geographies_dir: Path, level_label: str) -> gpd.GeoDataFrame: if year in (1980, 1990): gdf = read_nhgis_geography(year, geographies_dir, level_label) else: @@ -504,80 +470,117 @@ def read_geography( return gdf.to_crs(TARGET_CRS) -def join_population( - gdf: gpd.GeoDataFrame, - pop: pd.DataFrame, - year: int, - level_label: str, -) -> gpd.GeoDataFrame: +def join_population(gdf: gpd.GeoDataFrame, pop: pd.DataFrame, year: int, level_label: str) -> gpd.GeoDataFrame: state_fips = set(pop["STATEFP"]) + geo_states = set(gdf["STATEFP"].unique()) gdf = gdf[gdf["STATEFP"].isin(state_fips)].copy() merged = gdf.merge( pop[["JOIN_KEY", "WHITE", "BLACK", "TOTPOP", "POC"]], - on="JOIN_KEY", - how="left", - validate="one_to_one", - ) - + on="JOIN_KEY", how="left", validate="one_to_one") + + if len(merged) == 0: + unmatched = geo_states - state_fips + print( + f" {year} {level_label}: join produced 0 rows — " + f"geography STATEFP {sorted(geo_states)} not found in population data" + + (f" (unmatched: {sorted(unmatched)})" if unmatched != geo_states else ""), + flush=True, + ) + return merged missing = merged["TOTPOP"].isna().sum() - if missing: - # Blocks include many water-only and unpopulated geographic areas that have - # no row in the population CSV; a high miss rate is expected and not an error. - threshold = 0.50 if level_label == "blocks" else 0.01 - if missing / len(merged) > threshold: - raise ValueError( - f"{year}: {missing} {level_label} geometries did not match " - "population rows." - ) - merged = merged[merged["TOTPOP"].notna()].copy() - + if missing / len(merged) > 0.3: + # For 1990 blocks, the population extract omits Block Numbering Area (BNA) + # blocks — rural non-tracted areas identified by tract codes starting with + # "95" in the GISJOIN. Dropping them is expected behaviour. + if level_label == "blocks" and year == 1990: + # ~31% miss rate is expected: zero-pop blocks, alpha-suffix splits, + # and BNA tracts are all absent from the NHGIS extract. + print(f" {year}, {state_fips}: dropping {missing:,} unmatched blocks (zero-population, BNA, or split-block suffixes — expected for 1990).", flush=True) + else: + raise ValueError(f"{year}: {missing} {level_label} geometries did not match population rows.") + + # if level_label == "blocks" and year == 1990: + # unmatched_gisjoin = merged.loc[merged["TOTPOP"].isna(), "GISJOIN"] + # bna_count = unmatched_gisjoin.str[8:10].eq("95").sum() + # if bna_count / missing > 0.8: + # print(f"{year}: dropping {missing:,} unmatched blocks ({bna_count:,} / {missing:,} are BNA blocks absent from the population extract).", flush=True) + # else: + # raise ValueError(f"{year}: {missing} blocks geometries did not match population rows.") + # else: + # raise ValueError(f"{year}: {missing} {level_label} geometries did not match population rows.") + # if missing: + # # Blocks include many water-only and unpopulated geographic areas that have + # # no row in the population CSV; a high miss rate is expected and not an error. + # threshold = 0.50 if level_label == "blocks" else 0.01 + # if missing / len(merged) > threshold: + # raise ValueError( + # f"{year}: {missing} {level_label} geometries did not match " + # "population rows.") + + merged = merged[merged["TOTPOP"].notna()].copy() for col in ["WHITE", "BLACK", "TOTPOP", "POC"]: merged[col] = merged[col].astype("int64") return merged -def write_processed( - gdf: gpd.GeoDataFrame, - year: int, - output_dir: Path, - level_label: str, -) -> Path: - output_dir.mkdir(parents=True, exist_ok=True) - output_path = output_dir / f"{year}_{level_label}.gpkg" +def write_processed(gdf: gpd.GeoDataFrame, year: int, output_dir: Path, level_label: str, statefp: str) -> Path: + output_dir_by_level = output_dir / level_label + output_dir_by_level.mkdir(parents=True, exist_ok=True) + output_path = output_dir_by_level / f"{year}_{level_label}_{statefp}.gpkg" gdf.to_file(output_path, driver="GPKG") return output_path -def main( - level: str = typer.Option( - "tracts", - help="tracts, block_groups, blocks, or counties", - ), +def main(level: str = typer.Option("tracts", + help="tracts, block_groups, blocks, places, or counties"), years: Optional[str] = typer.Option(None, "--years", help="Space- or comma-separated years."), year_values: Optional[List[int]] = typer.Option(None, "--year", "-y"), population_dir: Path = typer.Option(POPULATION_DIR), geographies_dir: Path = typer.Option(GEOGRAPHIES_DIR), - output_dir: Path = typer.Option(OUTPUT_DIR), -) -> None: + output_dir: Path = typer.Option(OUTPUT_DIR)) -> None: level_label = validate_level(level) run_years = parse_years(years, year_values) + states_by_year: dict = {} for year in run_years: if year == 1980 and level_label in ("block_groups", "blocks"): - print( - f"Skipping 1980 {level_label}: NHGIS does not publish 1980 " - "block group or block boundary shapefiles.", - flush=True, - ) + print(f"Skipping 1980 {level_label}: NHGIS does not publish 1980 " + "block group or block boundary shapefiles.", flush=True) continue - print(f"Processing {year} {level_label}") + print(f"Creating {year} {level_label} geography geopackage files", flush=True) pop = read_population(year, population_dir, level_label) - gdf = read_geography(year, geographies_dir, level_label) - merged = join_population(gdf, pop, year, level_label) - path = write_processed(merged, year, output_dir, level_label) - print(path) + if year in (1980, 1990): + gdf = read_nhgis_geography(year, geographies_dir, level_label) + # state_iter = gdf.to_crs(TARGET_CRS).groupby("STATEFP") + state_iter = (group for _, group in gdf.to_crs(TARGET_CRS).groupby("STATEFP")) + else: + state_iter = read_census_geography(year, geographies_dir, level_label) + + n_written = 0 + states_written: set = set() + for state_gdf in state_iter: + statefp = state_gdf["STATEFP"].iloc[0] + merged = join_population(state_gdf, pop, year, level_label) + if len(merged) == 0: + continue + write_processed(merged, year, output_dir, level_label, statefp) + n_written += 1 + states_written.add(statefp) + states_by_year[year] = states_written + print(f" {n_written} {level_label} geopackages written for {year}", flush=True) + + if len(states_by_year) > 1: + all_states = set().union(*states_by_year.values()) + for statefp in sorted(all_states): + missing_in = [y for y, states in sorted(states_by_year.items()) if statefp not in states] + if missing_in: + print(f" Note: FIPS {statefp} has no {level_label} geopackage for {missing_in} " + "(absent from source data for those years)", flush=True) + + # gdf = read_geography(year, geographies_dir, level_label) + # for statefp, state_gdf in merged.groupby("STATEFP"): if __name__ == "__main__": diff --git a/pipeline/preprocessing/overlaps.py b/pipeline/preprocessing/overlaps.py index cad9af0..0a0014b 100644 --- a/pipeline/preprocessing/overlaps.py +++ b/pipeline/preprocessing/overlaps.py @@ -1,82 +1,90 @@ +""" +This script finds which node units (e.g. tracts, blocks) fall within which study areas (e.g. cities, CBSAs). Process: +1. Find a bounding box of each state-level node units file and store. +2. For each study area, find which state-level bounding box oberlap. +3. Load only those files. +4. For these, calculate representative points per node unit. Select units whose representative point falls within the study area. +5. Save these geographies. +""" + +import sys +from pathlib import Path +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) + import typer import tqdm import glob +from pipeline.utils.pipeline_log import tqdm_file import geopandas as gpd -import sys -from pathlib import Path - - -# for flexibility with different versions of geopandas -def union_geometry(gdf: gpd.GeoDataFrame): - geometry = gdf.geometry - if hasattr(geometry, "union_all"): - return geometry.union_all() - return geometry.unary_union +import fiona +import pandas as pd -def output_stem( - study_area_file: str, - prefix: str, - census_geography_type: str, - census_geography_year: str, - definition_vintage: str, -) -> str: +def output_stem(study_area_file: str, prefix: str, census_geography_type: str, census_geography_year: str, definition_vintage: str) -> str: study_area_stem = Path(study_area_file).stem if census_geography_type and census_geography_year and definition_vintage: vintage_suffix = f"_{definition_vintage}" if not study_area_stem.endswith(vintage_suffix): - raise ValueError( - f"{study_area_file} does not end with vintage {definition_vintage}" - ) + raise ValueError(f"{study_area_file} does not end with vintage {definition_vintage}") study_area_identity = study_area_stem.removesuffix(vintage_suffix) - return ( - f"{prefix}{census_geography_type}_in_{study_area_identity}_" - f"{census_geography_year}_{definition_vintage}_vintage" - ) + return (f"{prefix}{census_geography_type}_in_{study_area_identity}_" + f"{census_geography_year}_{definition_vintage}_vintage") return f"{prefix}{study_area_stem}_geographies" -def main( - census_geographies_file: str, - study_area_glob: str, - output_dir: str, - prefix: str = "", - census_geography_type: str = "", - census_geography_year: str = "", - definition_vintage: str = "", -): - """ - Writes census geographies whose representative points fall within each study area. - """ - census_geographies = gpd.read_file(census_geographies_file) - geography_points = census_geographies.geometry.representative_point() +def _run_year(study_area_glob: str, output_dir: str, prefix: str, census_geography_type: str, census_geography_year: str, definition_vintage: str, census_geographies_dir: str) -> None: + """Run overlap clipping for a single census geography year.""" + # get 4 bounds of each state block collection; keep state_files and state_bounds in sync + state_files = [] + state_bounds = [] + for f in sorted((Path(census_geographies_dir) / census_geography_type).glob(f"{census_geography_year}_{census_geography_type}_*.gpkg")): + with fiona.open(f) as src: + if len(src) == 0: + print("Skipping empty file:", f) + else: + state_files.append(f) + state_bounds.append(src.bounds) # (minx, miny, maxx, maxy) Path(output_dir).mkdir(parents=True, exist_ok=True) - for study_area_file in tqdm.tqdm(sorted(glob.glob(study_area_glob))): - study_area_gdf = gpd.read_file(study_area_file).to_crs(census_geographies.crs) - study_area_boundary = union_geometry(study_area_gdf) + n_written = 0 + n_skipped = 0 + for study_area_file in tqdm.tqdm(sorted(glob.glob(study_area_glob)), desc=census_geography_year, file=tqdm_file): + study_area_gdf = gpd.read_file(study_area_file).to_crs("esri:102003") + study_area_boundary = study_area_gdf.union_all() + minx, miny, maxx, maxy = study_area_boundary.bounds + + # select block files if their bounds are within study_area_boundary + needed = [f for f, b in zip(state_files, state_bounds) + if b[0] <= maxx and b[2] >= minx and b[1] <= maxy and b[3] >= miny] + if not needed: + print("No state files intersect", study_area_file, file=sys.stderr) + n_skipped += 1 + continue + + frames = [gpd.read_file(f) for f in needed] + census_geographies = gpd.GeoDataFrame(pd.concat(frames, ignore_index=True), crs=frames[0].crs) - geography_indices = geography_points.sindex.query( - study_area_boundary, predicate="covers") + geography_points = census_geographies.geometry.representative_point() + geography_indices = geography_points.sindex.query(study_area_boundary, predicate="covers") selected_geographies = census_geographies.iloc[sorted(geography_indices)] if len(selected_geographies) != 0: - selected_geographies_stem = output_stem( - study_area_file, - prefix, - census_geography_type, - census_geography_year, - definition_vintage, - ) - selected_geographies.to_file( - f"{output_dir}/{selected_geographies_stem}.gpkg", driver="GPKG") + selected_geographies_stem = output_stem(study_area_file, prefix, census_geography_type, census_geography_year, definition_vintage) + selected_geographies.to_file(f"{output_dir}/{selected_geographies_stem}.gpkg", driver="GPKG") + n_written += 1 else: - print( - "empty overlaps computed:", - census_geographies_file, - study_area_file, - output_dir, - file=sys.stderr) + print("Empty overlaps between study area file:", study_area_file, "and state files:", [str(file_link) for file_link in needed], file=sys.stderr) + n_skipped += 1 + print(f"Overlaps between node units and study area in {census_geography_year}: {n_written} written, {n_skipped} skipped (empty)", flush=True) + + +def main(study_area_glob: str, output_base_dir: str, prefix: str = "", census_geography_type: str = "", census_geography_years: str = "", definition_vintage: str = "2020", census_geographies_dir: str = "data/processed/census_geographies"): + """ + Writes census geographies whose representative points fall within each study area, + for each year in census_geography_years (space-separated string). + """ + for year in census_geography_years.split(): + _run_year(study_area_glob, f"{output_base_dir}/{year}", prefix, census_geography_type, year, definition_vintage, census_geographies_dir) if __name__ == "__main__": diff --git a/pipeline/preprocessing/study_areas.py b/pipeline/preprocessing/study_areas.py index 1a004a5..c7d574d 100644 --- a/pipeline/preprocessing/study_areas.py +++ b/pipeline/preprocessing/study_areas.py @@ -1,7 +1,15 @@ +""" +The script builds study area definition files (.gpkg + .json), one per study area. It writes two files per study area into "data/processed/study_area_definitions": +{type}_{code}_{vintage}.gpkg with the boundary geometry +{type}_{code}_{vintage}.json with metadata (CBSA code, title, component counties, total population) +The .gpkg files are what overlaps.py reads as study_area_glob. +""" + import sys from pathlib import Path sys.path.insert(0, str(Path(__file__).resolve().parents[2])) -from pipeline.utils.definitions import CBSA +from pipeline.utils.definitions import StudyArea +from pipeline.utils.pipeline_log import tqdm_file import tqdm import pandas as pd @@ -11,110 +19,80 @@ from pathlib import Path -def cbsa_to_dict(cbsa: CBSA) -> dict: - if hasattr(cbsa, "model_dump"): - return cbsa.model_dump(exclude={"geometry"}) - return json.loads(cbsa.json(exclude={"geometry"})) - - -def main( - filename: str = "", - definition_geographies: str = "data/processed/census_geographies/2020_tracts.gpkg", - output_dir: str = "data/processed/study_area_definitions", - study_area_type: str = "cbsa", - definition_vintage: str = "", - cbsa_geographies: str = None -): +def main(filename: str = "data/raw/study_area_sources/list1_march_2020.xls", definition_geographies: str = None, output_dir: str = "data/processed/study_area_definitions", study_area_type: str = "cbsa", definition_vintage: str = "march_2020", cbsa_geographies: str = None): if study_area_type == "counties": study_area_type = "county" if study_area_type not in {"cbsa", "max_county", "county", "max_city"}: - raise ValueError( - f"Unsupported study area type {study_area_type!r}. " - "Use 'cbsa', 'max_county', 'max_city', or 'county'." - ) + raise ValueError(f"Unsupported study area type {study_area_type}. Use 'cbsa', 'max_county', 'max_city', or 'county'.") + + if definition_geographies is None: + if study_area_type == "max_city": + definition_geographies = "data/processed/census_geographies/places/2020_places_*.gpkg" + else: + definition_geographies = "data/processed/census_geographies/counties/2020_counties_*.gpkg" + if cbsa_geographies is None and study_area_type == "max_city": + cbsa_geographies = "data/processed/census_geographies/counties/2020_counties_*.gpkg" if study_area_type == "county": - build_county_definitions( - definition_geographies, - output_dir, - definition_vintage or Path(definition_geographies).stem.split("_", 1)[0], - ) + build_county_definitions(definition_geographies, output_dir, definition_vintage or Path(definition_geographies).stem.split("_", 1)[0]) return - if study_area_type == "max_county": - build_max_county_definitions( - filename, - definition_geographies, - output_dir, - definition_vintage or Path(definition_geographies).stem.split("_", 1)[0], - ) + build_max_county_definitions(filename, definition_geographies, output_dir, definition_vintage or Path(definition_geographies).stem.split("_", 1)[0]) return if study_area_type == "max_city": if not filename: raise ValueError("max_city study areas require --filename.") - build_max_city_definitions( - filename, - definition_geographies, - output_dir, - definition_vintage or Path(definition_geographies).stem.split("_", 1)[0], - cbsa_geographies=cbsa_geographies - ) + build_max_city_definitions(filename, definition_geographies, output_dir, definition_vintage or Path(definition_geographies).stem.split("_", 1)[0], + cbsa_geographies=cbsa_geographies) return - if not filename: raise ValueError("CBSA study areas require --filename.") - source_stem = Path(filename).stem - if not definition_vintage: - if source_stem.startswith("list1_"): - definition_vintage = source_stem.removeprefix("list1_") - elif source_stem.startswith(f"{study_area_type}_"): - definition_vintage = source_stem.removeprefix(f"{study_area_type}_") - else: - definition_vintage = source_stem - metro_mappings = create_metro_mappings(fetch_metro_areas(filename)) - country = gpd.read_file(definition_geographies) + country = load_census_geography(definition_geographies) country["STATEFP"] = country["STATEFP"].astype(str).str.zfill(2) country["COUNTYFP"] = country["COUNTYFP"].astype(str).str.zfill(3) country["STCNTYFP"] = country["STATEFP"] + country["COUNTYFP"] Path(output_dir).mkdir(parents=True, exist_ok=True) - for cbsa_code, cbsa in tqdm.tqdm(metro_mappings.items()): + for cbsa_code, cbsa in tqdm.tqdm(metro_mappings.items(), file=tqdm_file): cbsa = add_cbsa_pop_and_geometry(country, cbsa) output_stem = f"{study_area_type}_{cbsa_code}_{definition_vintage}" with open(f"{output_dir}/{output_stem}.json", "w") as w: - json.dump(cbsa_to_dict(cbsa), w) + json.dump(area_to_dict(cbsa), w) cbsa.geometry.to_file(f"{output_dir}/{output_stem}.gpkg", driver="GPKG") def fetch_metro_areas(filename) -> pd.DataFrame: + """ + Reads the Census Bureau's CBSA Excel file (the delineation file), filters to Metropolitan Statistical Areas only, returns a df of CBSA-codes and FIPS rows + """ cbsa_counties = pd.read_excel(filename, skiprows=2) cbsa_counties = cbsa_counties[~cbsa_counties["FIPS County Code"].isna()] cbsa_counties["FIPS County Code"] = ( cbsa_counties["FIPS County Code"] .astype(int) .astype(str) - .str.zfill(3) - ) + .str.zfill(3)) cbsa_counties["FIPS State Code"] = ( cbsa_counties["FIPS State Code"] .astype(int) .astype(str) - .str.zfill(2) - ) + .str.zfill(2)) metro_areas = cbsa_counties[ cbsa_counties["Metropolitan/Micropolitan Statistical Area"] - == "Metropolitan Statistical Area" - ] + == "Metropolitan Statistical Area"] return metro_areas -def create_metro_mappings(metro_areas: pd.DataFrame) -> dict[str, CBSA]: +def create_metro_mappings(metro_areas: pd.DataFrame) -> dict[str, StudyArea]: + """ + Groups metro area rows into a dict of cbsa_code mapped to CBSA objects, each holding a list of component county FIPS codes + """ metro_mappings = {} for _, row in metro_areas.iterrows(): cbsa_code = row["CBSA Code"] @@ -123,22 +101,17 @@ def create_metro_mappings(metro_areas: pd.DataFrame) -> dict[str, CBSA]: if cbsa_code in metro_mappings: metro_mappings[cbsa_code].component_counties_fips.append(fips_code) else: - metro_mappings[cbsa_code] = CBSA( - area_code=cbsa_code, - cbsa_title=cbsa_title, - component_counties_fips=[fips_code], - total_population=None, - ) - + metro_mappings[cbsa_code] = StudyArea(area_code=cbsa_code, area_title=cbsa_title, component_counties_fips=[fips_code], total_population=None) return metro_mappings -def add_cbsa_pop_and_geometry(country: gpd.GeoDataFrame, cbsa: CBSA) -> CBSA: +def add_cbsa_pop_and_geometry(country: gpd.GeoDataFrame, cbsa: StudyArea) -> StudyArea: + """ + Filters the national counties GDF to this CBSA's component counties, dissolves them into one polygon, sums population, then writes .gpkg and .json + """ assert cbsa.total_population is None - cbsa_components = country[ - country["STCNTYFP"].apply(lambda x: x in cbsa.component_counties_fips) - ] + cbsa_components = country[country["STCNTYFP"].apply(lambda x: x in cbsa.component_counties_fips)] cbsa.geometry = cbsa_components.dissolve() cbsa.total_population = int(cbsa_components["TOTPOP"].sum()) @@ -162,65 +135,52 @@ def county_title(row: pd.Series) -> str: return f"County {row['STATEFP']}{row['COUNTYFP']}" -def build_county_definitions( - definition_geographies: str, - output_dir: str, - definition_vintage: str, -) -> None: - counties = gpd.read_file(definition_geographies) +def build_county_definitions(definition_geographies: str, output_dir: str, definition_vintage: str) -> None: + counties = load_census_geography(definition_geographies) state_col = first_existing_column( counties, - ["STATEFP", "STATEFP20", "STATEFP10", "STATEFP00"], - ) + ["STATEFP", "STATEFP20", "STATEFP10", "STATEFP00"]) county_col = first_existing_column( counties, - ["COUNTYFP", "COUNTYFP20", "COUNTYFP10", "COUNTYFP00"], - ) + ["COUNTYFP", "COUNTYFP20", "COUNTYFP10", "COUNTYFP00"]) counties["STATEFP"] = counties[state_col].astype(str).str.zfill(2) counties["COUNTYFP"] = counties[county_col].astype(str).str.zfill(3) counties["STCNTYFP"] = counties["STATEFP"] + counties["COUNTYFP"] Path(output_dir).mkdir(parents=True, exist_ok=True) - for _, county in tqdm.tqdm(counties.iterrows(), total=len(counties)): + for _, county in tqdm.tqdm(counties.iterrows(), total=len(counties), file=tqdm_file): county_fips = county["STCNTYFP"] output_stem = f"county_{county_fips}_{definition_vintage}" county_gdf = gpd.GeoDataFrame( [county], columns=counties.columns, - crs=counties.crs, - ) - study_area = CBSA( + crs=counties.crs) + study_area = StudyArea( area_code=county_fips, - cbsa_title=county_title(county), + area_title=county_title(county), component_counties_fips=[county_fips], total_population=( int(county["TOTPOP"]) if "TOTPOP" in county and pd.notna(county["TOTPOP"]) else None ), - geometry=county_gdf, - ) + geometry=county_gdf) with open(f"{output_dir}/{output_stem}.json", "w") as w: - json.dump(cbsa_to_dict(study_area), w) + json.dump(area_to_dict(study_area), w) county_gdf.to_file(f"{output_dir}/{output_stem}.gpkg", driver="GPKG") -def build_max_county_definitions( - filename: str, - definition_geographies: str, - output_dir: str, - definition_vintage: str, -) -> None: +def build_max_county_definitions(filename: str, definition_geographies: str, output_dir: str, definition_vintage: str) -> None: metro_mappings = create_metro_mappings(fetch_metro_areas(filename)) - counties = gpd.read_file(definition_geographies) - + counties = load_census_geography(definition_geographies) + counties["STATEFP"] = counties["STATEFP"].astype(str).str.zfill(2) counties["COUNTYFP"] = counties["COUNTYFP"].astype(str).str.zfill(3) counties["STCNTYFP"] = counties["STATEFP"] + counties["COUNTYFP"] Path(output_dir).mkdir(parents=True, exist_ok=True) - for cbsa_code, cbsa in tqdm.tqdm(metro_mappings.items()): + for cbsa_code, cbsa in tqdm.tqdm(metro_mappings.items(), file=tqdm_file): components = counties[counties["STCNTYFP"].isin(cbsa.component_counties_fips)] try: ##guard agains the extremely unlikely possibility of nonexistent counties in a cbsa. max_county = components.loc[components["TOTPOP"].idxmax()] @@ -235,32 +195,24 @@ def build_max_county_definitions( columns=counties.columns, crs=counties.crs, ) - study_area = CBSA( + study_area = StudyArea( area_code=county_fips, - cbsa_title=county_title(max_county), + area_title=county_title(max_county), component_counties_fips=[county_fips], total_population=( - int(max_county["TOTPOP"]) if "TOTPOP" in max_county and pd.notna(max_county["TOTPOP"]) else None - ), - geometry=county_gdf, - ) + int(max_county["TOTPOP"]) if "TOTPOP" in max_county and pd.notna(max_county["TOTPOP"]) else None), + geometry=county_gdf) with open(f"{output_dir}/{output_stem}.json", "w") as w: - json.dump(cbsa_to_dict(study_area), w) + json.dump(area_to_dict(study_area), w) county_gdf.to_file(f"{output_dir}/{output_stem}.gpkg", driver= "GPKG") -def build_max_city_definitions( - filename: str, - definition_geographies: str, - output_dir: str, - definition_vintage: str, - cbsa_geographies: str = None -) -> None: +def build_max_city_definitions(filename: str, definition_geographies: str, output_dir: str, definition_vintage: str, cbsa_geographies: str = None) -> None: metro_mappings = create_metro_mappings(fetch_metro_areas(filename)) - places = gpd.read_file(definition_geographies).to_crs("esri:102003") - counties = gpd.read_file(cbsa_geographies).to_crs("esri:102003") if cbsa_geographies else gpd.read_file(definition_geographies).to_crs("esri:102003") + places = load_census_geography(definition_geographies).to_crs("esri:102003") + counties = load_census_geography(cbsa_geographies).to_crs("esri:102003") counties["STATEFP"] = counties["STATEFP"].astype(str).str.zfill(2) counties["COUNTYFP"] = counties["COUNTYFP"].astype(str).str.zfill(3) @@ -268,7 +220,7 @@ def build_max_city_definitions( Path(output_dir).mkdir(parents=True, exist_ok=True) - for cbsa_code, cbsa in tqdm.tqdm(metro_mappings.items()): + for cbsa_code, cbsa in tqdm.tqdm(metro_mappings.items(), file=tqdm_file): components = counties[counties["STCNTYFP"].isin(cbsa.component_counties_fips)] cbsa_boundary = components.dissolve() @@ -281,20 +233,36 @@ def build_max_city_definitions( output_stem = f"max_city_{max_place['GEOID'].iloc[0]}_{definition_vintage}" - study_area = CBSA( + study_area = StudyArea( area_code=max_place["GEOID"].iloc[0], - cbsa_title=str(max_place["NAMELSAD"].iloc[0]), + area_title=str(max_place["NAMELSAD"].iloc[0]), component_counties_fips=[max_place["GEOID"].iloc[0]], - total_population=( - int(max_place["TOTPOP"].iloc[0]) if "TOTPOP" in max_place.columns and pd.notna(max_place["TOTPOP"].iloc[0]) else None - ), - geometry=max_place, - ) + total_population=(int(max_place["TOTPOP"].iloc[0]) if "TOTPOP" in max_place.columns and pd.notna(max_place["TOTPOP"].iloc[0]) else None), + geometry=max_place) with open(f"{output_dir}/{output_stem}.json", "w") as w: - json.dump(cbsa_to_dict(study_area), w) + json.dump(area_to_dict(study_area), w) max_place.to_file(f"{output_dir}/{output_stem}.gpkg", driver = "GPKG") +def area_to_dict(cbsa: StudyArea) -> dict: + if hasattr(cbsa, "model_dump"): + return cbsa.model_dump(exclude={"geometry"}) + return json.loads(cbsa.json(exclude={"geometry"})) + + +def load_census_geography(path_or_glob: str) -> gpd.GeoDataFrame: + """ + Helper that accepts either a file path or a glob pattern, and loads + concatenates accordingly + """ + p = Path(path_or_glob) + if p.is_file(): + return gpd.read_file(p) + files = sorted(p.parent.glob(p.name)) + if not files: + raise FileNotFoundError(f"No files matching {path_or_glob}") + frames = [gpd.read_file(f) for f in files] + return gpd.GeoDataFrame(pd.concat(frames, ignore_index=True), crs=frames[0].crs) + if __name__ == "__main__": typer.run(main) diff --git a/pipeline/process_results.py b/pipeline/process_results.py index 36c8a19..39e6646 100644 --- a/pipeline/process_results.py +++ b/pipeline/process_results.py @@ -7,7 +7,7 @@ from pipeline.utils import definitions -def parse_cbsa(config_loc: str) -> definitions.CBSA: +def parse_cbsa(config_loc: str) -> definitions.StudyArea: with open(config_loc) as f: data = json.load(f) @@ -17,9 +17,9 @@ def parse_cbsa(config_loc: str) -> definitions.CBSA: data.setdefault("geometry", None) data.setdefault("total_population", None) - if hasattr(definitions.CBSA, "model_validate"): - return definitions.CBSA.model_validate(data) - return definitions.CBSA.parse_obj(data) + if hasattr(definitions.StudyArea, "model_validate"): + return definitions.StudyArea.model_validate(data) + return definitions.StudyArea.parse_obj(data) def _strip_graph_suffix(filename: str) -> str: @@ -74,8 +74,8 @@ def enrich_metrics(df: pd.DataFrame) -> pd.DataFrame: lambda x: output_name_parts(x)[2] ) df["year"] = df["filename"].apply(lambda x: int(output_name_parts(x)[1])) - df["cbsa_title"] = cbsa_infos.apply(lambda x: x.cbsa_title) - df["cbsa_code"] = cbsa_infos.apply(lambda x: x.area_code) + df["area_title"] = cbsa_infos.apply(lambda x: x.area_title) + df["area_code"] = cbsa_infos.apply(lambda x: x.area_code) df["total_population_2020"] = cbsa_infos.apply(lambda x: x.total_population) return df diff --git a/pipeline/utils/definitions.py b/pipeline/utils/definitions.py index 52a3eef..e2e3aa0 100644 --- a/pipeline/utils/definitions.py +++ b/pipeline/utils/definitions.py @@ -4,20 +4,20 @@ from typing import List, Optional, Dict -class CBSA(pydantic.BaseModel): +class StudyArea(pydantic.BaseModel): class Config: arbitrary_types_allowed = True area_code: str - cbsa_title: str + area_title: str component_counties_fips: List[str] total_population: Optional[int] = None geometry: Optional[gpd.GeoDataFrame] = None if hasattr(pydantic, "RootModel"): - class CBSADict(pydantic.RootModel[Dict[str, CBSA]]): + class StudyAreaDict(pydantic.RootModel[Dict[str, StudyArea]]): pass else: - class CBSADict(pydantic.BaseModel): - __root__: Dict[str, CBSA] + class StudyAreaDict(pydantic.BaseModel): + __root__: Dict[str, StudyArea] diff --git a/pipeline/utils/pipeline_log.py b/pipeline/utils/pipeline_log.py new file mode 100644 index 0000000..d306ad7 --- /dev/null +++ b/pipeline/utils/pipeline_log.py @@ -0,0 +1,32 @@ +""" +Wrapper around print() for pipeline scripts. +When running a python script directly (i.e. when PIPELINE_LOG_FILE not set), `log()` behaves like `print()`. +""" + +import os +import sys + + +def log(msg: str, *, file=None) -> None: + """Print *msg* with immediate flush. + Parameters + ---------- + msg: + The message to emit. + file: + Passed to ``print()`` (default: stdout). + """ + print(msg, file=file, flush=True) + + +# In pipeline mode (PIPELINE_LOG_FILE set by reproduce.sh), tqdm should write +# directly to the terminal device so progress bars appear on screen but are not +# captured by the log file's tee redirect. Pass this as `file=` to every tqdm +# call. Falls back to None (tqdm default → stderr) when /dev/tty is +# unavailable (e.g. headless CI). +tqdm_file = None +if os.environ.get("PIPELINE_LOG_FILE"): + try: + tqdm_file = open("/dev/tty", "w") + except OSError: + pass # no controlling terminal — tqdm will use stderr diff --git a/pipeline/utils/visualization_settings.py b/pipeline/utils/visualization_settings.py index 54f5888..4eb51cb 100644 --- a/pipeline/utils/visualization_settings.py +++ b/pipeline/utils/visualization_settings.py @@ -88,8 +88,8 @@ def _build_metric_labels() -> dict: def _short_name(cbsa_title: str) -> str: - city_part, _, state = cbsa_title.rpartition(", ") - return f"{city_part.split('-')[0]}, {state}" + city_part, sep, state = cbsa_title.rpartition(", ") + return f"{city_part.split('-')[0]}, {state}" if sep else cbsa_title _NARROW_Y_THRESHOLD = 0.1 diff --git a/pipeline/visualization/generate_figures.py b/pipeline/visualization/generate_figures.py index fd188bb..0df8166 100644 --- a/pipeline/visualization/generate_figures.py +++ b/pipeline/visualization/generate_figures.py @@ -9,41 +9,24 @@ import typer from pipeline.process_results import enrich_metrics -from pipeline.utils.visualization_settings import ( - METRIC_LABELS, - METRICS, - PALETTE, - _apply_panel_style, - _shorten_prefix, - _short_name, -) +from pipeline.utils.visualization_settings import (METRIC_LABELS, METRICS, PALETTE, _apply_panel_style, _shorten_prefix, _short_name) from pipeline.visualization.plot_family_grids import plot_family_grids from pipeline.visualization.plot_single_metric import plot_single_metric from pipeline.visualization.plot_grid_top10 import plot_grid_top10 -def ensure_metadata(df: pd.DataFrame) -> pd.DataFrame: - required = {"definition_month_year", "year", "cbsa_title", "cbsa_code", "total_population_2020"} - if required.issubset(df.columns): - return df - return enrich_metrics(df) - - -def main( - filename: str = "outputs/tracts_in_cbsa/white_poc.csv", - n: int = 10, - prefix: str = "white_poc", - geography_type: Optional[str] = None, - fixed_y: bool = False, - study_area_type: Optional[str] = None -): +def main(filename: str = "", n: int = 10, prefix: str = "white_poc", geography_type: Optional[str] = None, fixed_y: bool = False, study_area_type: Optional[str] = None): if study_area_type == "max_county": - area_label = "Most Populous Counties within CBSAs" + area_label = "most populous counties within CBSAs" elif study_area_type == "max_city": - area_label = "Most Populous Cities within CBSAs" + area_label = "most populous cities within CBSAs" else: area_label = "CBSAs" - + + if not filename: + run_name = f"tracts_in_{study_area_type or 'cbsa'}" + filename = f"outputs/{run_name}/white_poc.csv" + if geography_type is None: for geo in ("block_groups", "blocks", "tracts", "counties"): if geo in prefix: @@ -67,11 +50,9 @@ def main( for month_year in set(df["definition_month_year"]): month_year_df = df[df["definition_month_year"] == month_year] - top_n_metros = list(month_year_df["cbsa_code"].drop_duplicates()[:n]) - code_to_title = month_year_df.drop_duplicates("cbsa_code").set_index("cbsa_code")["cbsa_title"] - top_n_df = month_year_df[ - month_year_df["cbsa_code"].isin(top_n_metros) - ].sort_values(["cbsa_code", "year"]) + top_n_metros = list(month_year_df["area_code"].drop_duplicates()[:n]) + code_to_title = month_year_df.drop_duplicates("area_code").set_index("area_code")["area_title"] + top_n_df = month_year_df[month_year_df["area_code"].isin(top_n_metros)].sort_values(["area_code", "year"]) color_map = {cbsa: PALETTE[i % len(PALETTE)] for i, cbsa in enumerate(top_n_metros)} years = sorted(top_n_df["year"].unique()) @@ -85,11 +66,10 @@ def main( _apply_panel_style(ax, years, None, y_range=y_range) for cbsa in top_n_metros: - cbsa_df = top_n_df[top_n_df["cbsa_code"] == cbsa] + cbsa_df = top_n_df[top_n_df["area_code"] == cbsa] ax.plot( cbsa_df["year"], cbsa_df[metric], - color=color_map[cbsa], linewidth=1.8, marker="o", markersize=4, zorder=2, - ) + color=color_map[cbsa], linewidth=1.8, marker="o", markersize=4, zorder=2) title, subtitle = METRIC_LABELS.get(metric, (metric.replace("_", " ").title(), "")) ax.set_title(title, fontsize=13, fontweight="bold", @@ -97,36 +77,17 @@ def main( if subtitle: ax.text(0.5, 1.04, subtitle, transform=ax.transAxes, ha="center", va="bottom", fontsize=9, color="#777777") - fig.suptitle( - f"Segregation over time: {pair_label}", - fontsize=14, fontweight="bold", color="#111111", y=1.06, - ) - fig.text( - 0.5, 1, - f"Top {n} U.S. metros by 2020 population · Census {geography_label} in {area_label}", - ha="center", fontsize=9, color="#555555", - ) - - handles = [ - plt.Line2D([0], [0], color=color_map[c], linewidth=2.5, label=_short_name(code_to_title[c])) - for c in top_n_metros - ] - fig.legend( - handles=handles, - loc="lower center", - ncol=min(5, len(top_n_metros)), - bbox_to_anchor=(0.5, -0.1), - frameon=False, - fontsize=8, - handlelength=1.5, - columnspacing=1.0, - labelcolor="#333333", - ) - - fig.savefig( - output_dir / "lineplots" / f"{prefix}_{metric}.png", - dpi=150, bbox_inches="tight", facecolor=BG, - ) + fig.suptitle(f"Segregation over time: {pair_label}", + fontsize=14, fontweight="bold", color="#111111", y=1.06) + fig.text(0.5, 1, f"Top {n} U.S. metros by 2020 population. Census {geography_label} in {area_label}", + ha="center", fontsize=9, color="#555555") + + handles = [plt.Line2D([0], [0], color=color_map[c], linewidth=2.5, label=_short_name(code_to_title[c])) + for c in top_n_metros] + fig.legend(handles=handles, loc="lower center", ncol=min(5, len(top_n_metros)), bbox_to_anchor=(0.5, -0.1), frameon=False, fontsize=8, handlelength=1.5, columnspacing=1.0, labelcolor="#333333") + + fig.savefig(output_dir / "lineplots" / f"{prefix}_{metric}.png", + dpi=150, bbox_inches="tight", facecolor=BG) plt.close(fig) plot_grid_top10(df, prefix, month_year, output_dir, n, geography_label=geography_label, area_label= area_label, fixed_y=fixed_y) @@ -134,5 +95,12 @@ def main( plot_family_grids(df, prefix, month_year, output_dir, n, geography_label=geography_label, area_label= area_label, fixed_y=fixed_y) +def ensure_metadata(df: pd.DataFrame) -> pd.DataFrame: + required = {"definition_month_year", "year", "area_title", "area_code", "total_population_2020"} + if required.issubset(df.columns): + return df + return enrich_metrics(df) + + if __name__ == "__main__": typer.run(main) diff --git a/pipeline/visualization/plot_family_grids.py b/pipeline/visualization/plot_family_grids.py index 93baded..85de308 100644 --- a/pipeline/visualization/plot_family_grids.py +++ b/pipeline/visualization/plot_family_grids.py @@ -3,34 +3,15 @@ import matplotlib.pyplot as plt import pandas as pd -from pipeline.utils.visualization_settings import ( - METRIC_LABELS, - METRICS, - PALETTE, - _apply_panel_style, - _short_name, -) +from pipeline.utils.visualization_settings import (METRIC_LABELS, METRICS, PALETTE, _apply_panel_style, _short_name) -def plot_family_grids( - df: pd.DataFrame, - prefix: str, - month_year: str, - output_dir: Path, - n: int = 10, - n_cols: int = 6, - geography_label: str = "tracts", - area_label: str = "CBSA", - fixed_y: bool = False, -) -> None: +def plot_family_grids(df: pd.DataFrame, prefix: str, month_year: str, output_dir: Path, n: int = 10, n_cols: int = 6, geography_label: str = "tracts", area_label: str = "CBSA", fixed_y: bool = False) -> None: BG = "#fafafa" month_year_df = df[df["definition_month_year"] == month_year] - top_n_metros = list(month_year_df["cbsa_code"].drop_duplicates()[:n]) - code_to_title = month_year_df.drop_duplicates("cbsa_code").set_index("cbsa_code")["cbsa_title"] - plot_df = ( - month_year_df[month_year_df["cbsa_code"].isin(top_n_metros)] - .sort_values(["cbsa_code", "year"]) - ) + top_n_metros = list(month_year_df["area_code"].drop_duplicates()[:n]) + code_to_title = month_year_df.drop_duplicates("area_code").set_index("area_code")["area_title"] + plot_df = (month_year_df[month_year_df["area_code"].isin(top_n_metros)].sort_values(["area_code", "year"])) color_map = {cbsa: PALETTE[i % len(PALETTE)] for i, cbsa in enumerate(top_n_metros)} years = sorted(plot_df["year"].unique()) @@ -46,10 +27,8 @@ def plot_family_grids( family_dir = output_dir / "metric_family_grids" family_dir.mkdir(parents=True, exist_ok=True) - handles = [ - plt.Line2D([0], [0], color=color_map[c], linewidth=2.5, label=_short_name(code_to_title[c])) - for c in top_n_metros - ] + handles = [plt.Line2D([0], [0], color=color_map[c], linewidth=2.5, label=_short_name(code_to_title[c])) + for c in top_n_metros] for family_title, members in families.items(): n_metrics = len(members) @@ -63,8 +42,7 @@ def plot_family_grids( figsize=(5 * cols, 5 * rows), facecolor=BG, sharey=False, - squeeze=False, - ) + squeeze=False) for idx, (metric, subtitle) in enumerate(members): ax = axes[idx // cols][idx % cols] @@ -72,48 +50,25 @@ def plot_family_grids( _apply_panel_style(ax, years, ylim, y_range=y_range) ax.set_title( subtitle if subtitle else family_title, - fontsize=10, fontweight="bold", pad=8, color="#111111", - ) + fontsize=10, fontweight="bold", pad=8, color="#111111") for cbsa in top_n_metros: - cbsa_df = plot_df[plot_df["cbsa_code"] == cbsa] + cbsa_df = plot_df[plot_df["area_code"] == cbsa] ax.plot( cbsa_df["year"], cbsa_df[metric], - color=color_map[cbsa], linewidth=1.8, marker="o", markersize=4, zorder=2, - ) + color=color_map[cbsa], linewidth=1.8, marker="o", markersize=4, zorder=2) for idx in range(n_metrics, rows * cols): axes[idx // cols][idx % cols].set_visible(False) SUPTITLE_Y = 1.02 - fig.suptitle( - f"{family_title} · Segregation over time: {pair_label}", - fontsize=14, fontweight="bold", color="#111111", y=SUPTITLE_Y, - ) - fig.legend( - handles=handles, - loc="lower center", - ncol=min(5, len(top_n_metros)), - bbox_to_anchor=(0.5, -0.03), - frameon=False, - fontsize=8, - handlelength=1.5, - columnspacing=1.0, - labelcolor="#333333", - ) + fig.suptitle(f"{family_title}. Segregation over time: {pair_label}", fontsize=14, fontweight="bold", color="#111111", y=SUPTITLE_Y) + fig.legend(handles=handles, loc="lower center", ncol=min(5, len(top_n_metros)), bbox_to_anchor=(0.5, -0.03), frameon=False, fontsize=8, handlelength=1.5, columnspacing=1.0, labelcolor="#333333") subtitle_y = SUPTITLE_Y - 20 / (72 * fig.get_figheight()) - fig.text( - 0.5, subtitle_y, - f"Top {n} U.S. metros by 2020 population · Census {geography_label} in {area_label}", - ha="center", va="top", fontsize=9, color="#555555", - ) + fig.text(0.5, subtitle_y, + f"Top {n} U.S. metros by 2020 population. Census {geography_label} in {area_label}", + ha="center", va="top", fontsize=9, color="#555555") - safe_name = ( - family_title.lower() - .replace("'", "").replace("(", "").replace(")", "") - .replace(" ", "_") - ) - fig.savefig( - family_dir / f"{prefix}_{safe_name}.png", - dpi=150, bbox_inches="tight", facecolor=BG, - ) + safe_name = (family_title.lower().replace("'", "").replace("(", "").replace(")", "").replace(" ", "_")) + fig.savefig(family_dir / f"{prefix}_{safe_name}.png", + dpi=150, bbox_inches="tight", facecolor=BG) plt.close(fig) diff --git a/pipeline/visualization/plot_grid_top10.py b/pipeline/visualization/plot_grid_top10.py index 3231749..bb113a7 100644 --- a/pipeline/visualization/plot_grid_top10.py +++ b/pipeline/visualization/plot_grid_top10.py @@ -3,32 +3,17 @@ import matplotlib.pyplot as plt import pandas as pd -from pipeline.utils.visualization_settings import ( - GRID_METRICS, - PALETTE, - _apply_panel_style, - _short_name, -) +from pipeline.utils.visualization_settings import (GRID_METRICS, PALETTE, _apply_panel_style, _short_name) -def plot_grid_top10( - df: pd.DataFrame, - prefix: str, - month_year: str, - output_dir: Path, - n: int = 10, - geography_label: str = "tracts", - area_label: str = "CBSA", - fixed_y: bool = False, -) -> None: +def plot_grid_top10(df: pd.DataFrame, prefix: str, month_year: str, output_dir: Path, n: int = 10, geography_label: str = "tracts", area_label: str = "CBSA", fixed_y: bool = False) -> None: BG = "#fafafa" month_year_df = df[df["definition_month_year"] == month_year] - top_n_metros = list(month_year_df["cbsa_code"].drop_duplicates()[:n]) - code_to_title = month_year_df.drop_duplicates("cbsa_code").set_index("cbsa_code")["cbsa_title"] + top_n_metros = list(month_year_df["area_code"].drop_duplicates()[:n]) + code_to_title = month_year_df.drop_duplicates("area_code").set_index("area_code")["area_title"] plot_df = ( - month_year_df[month_year_df["cbsa_code"].isin(top_n_metros)] - .sort_values(["cbsa_code", "year"]) - ) + month_year_df[month_year_df["area_code"].isin(top_n_metros)] + .sort_values(["area_code", "year"])) available = [m for m in GRID_METRICS if m in plot_df.columns] if not available: @@ -48,49 +33,29 @@ def plot_grid_top10( _apply_panel_style(ax, years, ylim, y_range=y_range) ax.set_title(GRID_METRICS[metric], fontsize=11, fontweight="bold", pad=8, color="#111111") for cbsa in top_n_metros: - cbsa_df = plot_df[plot_df["cbsa_code"] == cbsa] + cbsa_df = plot_df[plot_df["area_code"] == cbsa] ax.plot( cbsa_df["year"], cbsa_df[metric], - color=color_map[cbsa], linewidth=1.8, marker="o", markersize=4, zorder=2, - ) + color=color_map[cbsa], linewidth=1.8, marker="o", markersize=4, zorder=2) pair_label = "White–Black" if prefix.startswith("wb") else "White–POC" fig.suptitle( f"Segregation over time: {pair_label}", - fontsize=14, fontweight="bold", color="#111111", y=1.04, - ) - fig.text( - 0.5, 0.97, - f"Top {n} U.S. metros by 2020 population · Census {geography_label} within {area_label}", - ha="center", fontsize=9, color="#555555", - ) - - handles = [ - plt.Line2D([0], [0], color=color_map[c], linewidth=2.5, label=_short_name(code_to_title[c])) - for c in top_n_metros - ] - fig.legend( - handles=handles, - loc="lower center", - ncol=min(5, len(top_n_metros)), - bbox_to_anchor=(0.5, -0.13), - frameon=False, - fontsize=8, - handlelength=1.5, - columnspacing=1.0, - labelcolor="#333333", - ) - fig.text( - 0.5, -0.22, - "Notes: Moran's I uses weights matrix P. Half Edge uses λ=1.\n" + fontsize=14, fontweight="bold", color="#111111", y=1.04) + fig.text(0.5, 0.97, + f"Segregation metrics in top {n} {area_label} by 2020 population.", ha="center", fontsize=9, color="#555555") + + handles = [plt.Line2D([0], [0], color=color_map[c], linewidth=2.5, label=_short_name(code_to_title[c])) + for c in top_n_metros] + fig.legend(handles=handles, loc="lower center", ncol=min(5, len(top_n_metros)), bbox_to_anchor=(0.5, -0.13), frameon=False, fontsize=8, handlelength=1.5, columnspacing=1.0, labelcolor="#333333") + fig.text(0.5, -0.22, + f"Notes: Calculated using Census {geography_label} in {area_label}.\n" "Sources: Decennial census and TIGER/Line shapefiles via Census API (2000-2020) and NHGIS (before 2000).", - ha="center", fontsize=7, color="#383838", linespacing=1.6, - ) + ha="center", fontsize=7, color="#383838", linespacing=1.6) grid_dir = output_dir / "grid_lineplots" grid_dir.mkdir(parents=True, exist_ok=True) fig.savefig( grid_dir / f"{prefix}_top10.png", - dpi=150, bbox_inches="tight", facecolor=BG, - ) + dpi=150, bbox_inches="tight", facecolor=BG) plt.close(fig) diff --git a/pipeline/visualization/plot_single_metric.py b/pipeline/visualization/plot_single_metric.py index 7e34205..823ad47 100644 --- a/pipeline/visualization/plot_single_metric.py +++ b/pipeline/visualization/plot_single_metric.py @@ -6,16 +6,10 @@ from pipeline.utils.visualization_settings import GRID_METRICS, _apply_panel_style -def plot_single_metric( - df: pd.DataFrame, - prefix: str, - month_year: str, - output_dir: Path, - geography_label: str = "tracts", - area_label: str = "CBSA", - fixed_y: bool = False, -) -> None: +def plot_single_metric(df: pd.DataFrame, prefix: str, month_year: str, output_dir: Path, geography_label: str = "tracts", area_label: str = "CBSA", fixed_y: bool = False) -> None: MIN_POPULATION = 100_000 + if "Cities" in area_label: + MIN_POPULATION = 0 BG = "#fafafa" month_year_df = df[df["definition_month_year"] == month_year] @@ -26,12 +20,12 @@ def plot_single_metric( years = sorted(month_year_df["year"].unique()) - cbsa_year_counts = month_year_df.groupby("cbsa_code")["year"].nunique() + cbsa_year_counts = month_year_df.groupby("area_code")["year"].nunique() complete_cbsas = cbsa_year_counts[cbsa_year_counts == len(years)].index - cbsa_pop = month_year_df.drop_duplicates("cbsa_code").set_index("cbsa_code")["total_population_2020"] + cbsa_pop = month_year_df.drop_duplicates("area_code").set_index("area_code")["total_population_2020"] eligible_cbsas = complete_cbsas[cbsa_pop.reindex(complete_cbsas).fillna(0) >= MIN_POPULATION] - month_year_df = month_year_df[month_year_df["cbsa_code"].isin(eligible_cbsas)] + month_year_df = month_year_df[month_year_df["area_code"].isin(eligible_cbsas)] all_cbsas = eligible_cbsas ylim = (month_year_df[available].min().min(), month_year_df[available].max().max()) if fixed_y else None @@ -49,53 +43,30 @@ def plot_single_metric( ax.set_title(GRID_METRICS[metric], fontsize=11, fontweight="bold", pad=8, color="#111111") for cbsa in all_cbsas: - cbsa_df = month_year_df[month_year_df["cbsa_code"] == cbsa].sort_values("year") + cbsa_df = month_year_df[month_year_df["area_code"] == cbsa].sort_values("year") ax.plot( cbsa_df["year"], cbsa_df[metric], - color="#aaaaaa", linewidth=0.7, alpha=0.4, zorder=1, - ) - ax.plot( - yearly_mean.index, yearly_mean[metric], - color="#0072b2", linewidth=2.4, marker="o", markersize=5, zorder=3, - ) + color="#aaaaaa", linewidth=0.7, alpha=0.4, zorder=1) + ax.plot(yearly_mean.index, yearly_mean[metric], + color="#0072b2", linewidth=2.4, marker="o", markersize=5, zorder=3) pair_label = "White–Black" if prefix.startswith("wb") else "White–POC" - fig.suptitle( - f"Segregation over time: {pair_label}", - fontsize=14, fontweight="bold", color="#111111", y=1.04, - ) - fig.text( - 0.5, 0.95, - f"U.S. CBSAs ≥100k pop., present in all years · Census {geography_label} in {area_label} · Mean in blue", - ha="center", fontsize=9, color="#555555", - ) - - handles = [ - plt.Line2D([0], [0], color="#aaaaaa", linewidth=1.5, alpha=0.6, label="Individual CBSA"), - plt.Line2D([0], [0], color="#0072b2", linewidth=2.4, marker="o", markersize=5, label="Mean across all CBSAs"), - ] - fig.legend( - handles=handles, - loc="lower center", - ncol=2, - bbox_to_anchor=(0.5, -0.08), - frameon=False, - fontsize=8, - handlelength=1.5, - columnspacing=1.5, - labelcolor="#333333", - ) - fig.text( - 0.5, -0.16, - "Notes: Moran's I uses weights matrix P. Half Edge uses λ=1.\n" - "Sources: Decennial census and TIGER/Line shapefiles via Census API (2000-2020) and NHGIS (before 2000).", - ha="center", fontsize=7, color="#383838", linespacing=1.6, - ) + fig.suptitle(f"Segregation over time: {pair_label}", + fontsize=14, fontweight="bold", color="#111111", y=1.04) + if "cities" in area_label: + fig.text(0.5, 0.95, f"Segregation metrics in {area_label}, present in all years ({len(eligible_cbsas)}). Mean in blue.", ha="center", fontsize=9, color="#555555") + else: + fig.text(0.5, 0.95, f"{len(eligible_cbsas)} {area_label} ≥100k pop., present in all years. Census {geography_label} in {area_label}. Mean in blue.", ha="center", fontsize=9, color="#555555") + + handles = [plt.Line2D([0], [0], color="#aaaaaa", linewidth=1.5, alpha=0.6, label="Individual area"), + plt.Line2D([0], [0], color="#0072b2", linewidth=2.4, marker="o", markersize=5, label="Mean across all areas")] + fig.legend(handles=handles, loc="lower center", ncol=2, bbox_to_anchor=(0.5, -0.08), frameon=False, fontsize=8, handlelength=1.5, columnspacing=1.5, labelcolor="#333333") + fig.text(0.5, -0.16, + f"Notes: Calculated using Census {geography_label} in {area_label}.\n" + # Moran's I uses weights matrix P. Half Edge uses λ=1.\n" + "Sources: Decennial census and TIGER/Line shapefiles via Census API (2000-2020) and NHGIS (before 2000).", ha="center", fontsize=7, color="#383838", linespacing=1.6) grid_dir = output_dir / "grid_lineplots" grid_dir.mkdir(parents=True, exist_ok=True) - fig.savefig( - grid_dir / f"{prefix}_all_cbsa.png", - dpi=150, bbox_inches="tight", facecolor=BG, - ) + fig.savefig(grid_dir / f"{prefix}_all_cbsa.png", dpi=150, bbox_inches="tight", facecolor=BG) plt.close(fig) diff --git a/scripts/build_study_areas.sh b/scripts/build_study_areas.sh deleted file mode 100644 index 13aa81a..0000000 --- a/scripts/build_study_areas.sh +++ /dev/null @@ -1,29 +0,0 @@ -#!/usr/bin/env bash - - -SCRIPT_DIR="$( - cd -- "$(dirname -- "${BASH_SOURCE[0]}")" >/dev/null 2>&1 - pwd -P -)" -TOP_DIR="$(cd -- "${SCRIPT_DIR}/.." >/dev/null 2>&1 && pwd -P)" -cd "${TOP_DIR}" - -_config="$(poetry run python scripts/resolve_config.py)" || exit 1 -eval "${_config}" - -if [ "${STUDY_AREA_TYPE}" = "cbsa" ] || [ "${STUDY_AREA_TYPE}" = "max_county" ] || [ "${STUDY_AREA_TYPE}" = "max_city" ]; then - poetry run python pipeline/preprocessing/study_areas.py \ - --filename "${STUDY_AREA_SOURCE_FILE}" \ - --definition-geographies "${STUDY_AREA_DEFINITION_GEOGRAPHIES}" \ - --output-dir "data/processed/study_area_definitions" \ - --study-area-type "${STUDY_AREA_TYPE}" \ - --definition-vintage "${STUDY_AREA_DEFINITION_VINTAGE}" \ - --cbsa-geographies "data/processed/census_geographies/${STUDY_AREA_DEFINITION_GEOGRAPHY_YEAR}_counties.gpkg" - -else - poetry run python pipeline/preprocessing/study_areas.py \ - --definition-geographies "${STUDY_AREA_DEFINITION_GEOGRAPHIES}" \ - --output-dir "data/processed/study_area_definitions" \ - --study-area-type "${STUDY_AREA_TYPE}" \ - --definition-vintage "${STUDY_AREA_DEFINITION_VINTAGE}" -fi diff --git a/scripts/overlaps.sh b/scripts/overlaps.sh deleted file mode 100644 index 188f3e7..0000000 --- a/scripts/overlaps.sh +++ /dev/null @@ -1,20 +0,0 @@ -#!/usr/bin/env bash - -SCRIPT_DIR="$( - cd -- "$(dirname -- "${BASH_SOURCE[0]}")" >/dev/null 2>&1 - pwd -P -)" -TOP_DIR="$(cd -- "${SCRIPT_DIR}/.." >/dev/null 2>&1 && pwd -P)" -cd "${TOP_DIR}" - -_config="$(poetry run python scripts/resolve_config.py)" || exit 1 -eval "${_config}" -IFS=" " read -r -a census_geography_years <<< "${CENSUS_GEOGRAPHY_YEARS}" - -parallel --bar \ - bash "${SCRIPT_DIR}/run_overlap.sh" \ - {} \ - "${CENSUS_GEOGRAPHY_TYPE}" \ - "${STUDY_AREA_TYPE}" \ - "${STUDY_AREA_DEFINITION_VINTAGE}" \ - ::: "${census_geography_years[@]}" diff --git a/scripts/reproduce.sh b/scripts/reproduce.sh index cf13188..2a131c0 100644 --- a/scripts/reproduce.sh +++ b/scripts/reproduce.sh @@ -1,114 +1,77 @@ #!/usr/bin/env bash +cd "$(dirname "${BASH_SOURCE[0]}")/.." +config="$(poetry run python pipeline/config.py)" || exit 1 +eval "${config}" -SCRIPT_DIR="$( - cd -- "$(dirname -- "${BASH_SOURCE[0]}")" >/dev/null 2>&1 - pwd -P -)" -TOP_DIR="$(cd -- "${SCRIPT_DIR}/.." >/dev/null 2>&1 && pwd -P)" -cd "${TOP_DIR}" - -_config="$(poetry run python scripts/resolve_config.py)" || exit 1 -eval "${_config}" -IFS=" " read -r -a census_geography_years <<< "${CENSUS_GEOGRAPHY_YEARS}" - -build_geography_inputs() { - local geography_type="$1" - local geography_years="$2" - - poetry run python pipeline/download/download_population_tables.py \ - --level "${geography_type}" \ - --years "${geography_years}" - poetry run python pipeline/download/download_geographies.py \ - --level "${geography_type}" \ - --years "${geography_years}" - poetry run python pipeline/preprocessing/census_geographies.py \ - --level "${geography_type}" \ - --years "${geography_years}" -} - -year_list_contains() { - local needle="$1" - shift - local year - for year in "$@"; do - if [ "${year}" = "${needle}" ]; then - return 0 - fi - done - return 1 -} - -calculate_csv() { - local output_file="$1"; shift - poetry run python pipeline/metrics.py "" "$@" --headers-only > "${output_file}" - for year in "${census_geography_years[@]}"; do - find "data/processed/dual_graphs/${year}" -type f \ - -name "${CENSUS_GEOGRAPHY_TYPE}_in_${STUDY_AREA_TYPE}_*_${year}_${STUDY_AREA_DEFINITION_VINTAGE}_vintage_connected.json" | - parallel --bar -j 6 poetry run python pipeline/metrics.py {} "$@" >> "${output_file}" - done -} - -# Set up folder structure. -bash "${SCRIPT_DIR}/setup.sh" +# Set up folder structure +bash scripts/setup.sh # Save a log of the run configuration RUN_STARTED_AT="$(date -u +"%Y-%m-%dT%H:%M:%SZ")" +export PIPELINE_LOG_FILE="${RUN_OUTPUT_DIR}/run.log" export METRIC_FAILURES_FILE="${RUN_OUTPUT_DIR}/metric_failures.csv" { - echo "run_name=${RUN_NAME}" echo "start_timestamp=${RUN_STARTED_AT}" - echo "" echo "graph_area_type=${STUDY_AREA_TYPE}" echo "nodes_area_type=${CENSUS_GEOGRAPHY_TYPE}" - echo "" echo "study_area_source_file=${STUDY_AREA_SOURCE_FILE:-}" - echo "" echo "census_geography_years=${CENSUS_GEOGRAPHY_YEARS}" + echo "" +} > "${PIPELINE_LOG_FILE}" + +# Write to the log file +exec > >(tee -a "${PIPELINE_LOG_FILE}") 2>&1 -} > "${RUN_OUTPUT_DIR}/run.log" +echo "=== 1. Download ===" +poetry run python pipeline/download/download_population_tables.py --level "${CENSUS_GEOGRAPHY_TYPE}" +poetry run python pipeline/download/download_geographies.py --level "${CENSUS_GEOGRAPHY_TYPE}" +poetry run python pipeline/preprocessing/census_geographies.py --level "${CENSUS_GEOGRAPHY_TYPE}" -# Download and join population values to census geography shapefiles. -build_geography_inputs "${CENSUS_GEOGRAPHY_TYPE}" "${CENSUS_GEOGRAPHY_YEARS}" -# Only download study-area-definition geographies separately when they differ from -# the census geography type/year already fetched above (avoids a redundant download). +# Only download study-area-definition geographies separately when they differ from the census geography type/year already fetched above (avoids a redundant download). if [ "${STUDY_AREA_DEFINITION_GEOGRAPHY_TYPE}" != "${CENSUS_GEOGRAPHY_TYPE}" ] || - ! year_list_contains "${STUDY_AREA_DEFINITION_GEOGRAPHY_YEAR}" "${census_geography_years[@]}"; then - build_geography_inputs \ - "${STUDY_AREA_DEFINITION_GEOGRAPHY_TYPE}" \ - "${STUDY_AREA_DEFINITION_GEOGRAPHY_YEAR}" + [[ " ${CENSUS_GEOGRAPHY_YEARS} " != *" ${STUDY_AREA_DEFINITION_GEOGRAPHY_YEAR} "* ]]; then + poetry run python pipeline/download/download_population_tables.py \ + --level "${STUDY_AREA_DEFINITION_GEOGRAPHY_TYPE}" --years "${STUDY_AREA_DEFINITION_GEOGRAPHY_YEAR}" + poetry run python pipeline/download/download_geographies.py \ + --level "${STUDY_AREA_DEFINITION_GEOGRAPHY_TYPE}" --years "${STUDY_AREA_DEFINITION_GEOGRAPHY_YEAR}" + poetry run python pipeline/preprocessing/census_geographies.py \ + --level "${STUDY_AREA_DEFINITION_GEOGRAPHY_TYPE}" --years "${STUDY_AREA_DEFINITION_GEOGRAPHY_YEAR}" fi -# Generate study area definition shapefiles. -bash "${SCRIPT_DIR}/build_study_areas.sh" - -# Select census geographies that overlap with study area definition shapefiles. -bash "${SCRIPT_DIR}/overlaps.sh" - -# Generate dual graphs. -_build_graph() { - local shp="$1" - local stem - stem="$(basename "${shp%.gpkg}")" - local dual_dir - dual_dir="data/processed/dual_graphs/$(basename "$(dirname "$shp")")" - poetry run python pipeline/graphs.py "$shp" "${dual_dir}/${stem}_orig.json" "${dual_dir}/${stem}_connected.json" -} -export -f _build_graph - -for year in "${census_geography_years[@]}"; do - find "data/processed/clipped_geographies/${year}" \ - -type f \ - -name "${CENSUS_GEOGRAPHY_TYPE}_in_${STUDY_AREA_TYPE}_*_${year}_${STUDY_AREA_DEFINITION_VINTAGE}_vintage.gpkg" -done | - parallel --bar -j 6 _build_graph {} - -# Calculate metrics. -calculate_csv "${RUN_OUTPUT_DIR}/white_black.csv" BLACK WHITE TOTPOP -calculate_csv "${RUN_OUTPUT_DIR}/white_poc.csv" POC WHITE TOTPOP - -# Generate figures under the configured run output. The figure script enriches -# raw metric CSVs with study-area metadata as needed. +echo "" +echo "=== 2. Study areas ===" +poetry run python pipeline/preprocessing/study_areas.py \ + ${STUDY_AREA_SOURCE_FILE:+--filename "${STUDY_AREA_SOURCE_FILE}"} \ + --study-area-type "${STUDY_AREA_TYPE}" +echo "Study areas are saved to `data/processed/study_area_definitions`." + +echo "" +echo "=== 3. Overlaps ===" +poetry run python pipeline/preprocessing/overlaps.py \ + "data/processed/study_area_definitions/${STUDY_AREA_TYPE}_*_${STUDY_AREA_DEFINITION_VINTAGE}.gpkg" \ + "data/processed/clipped_geographies" \ + --census-geography-type "${CENSUS_GEOGRAPHY_TYPE}" \ + --census-geography-years "${CENSUS_GEOGRAPHY_YEARS}" \ + --definition-vintage "${STUDY_AREA_DEFINITION_VINTAGE}" + +echo "" +echo "=== 4. Graphs ===" +poetry run python pipeline/graphs.py \ + "data/processed/clipped_geographies/*/${CENSUS_GEOGRAPHY_TYPE}_in_${STUDY_AREA_TYPE}_*_${STUDY_AREA_DEFINITION_VINTAGE}_vintage.gpkg" + +echo "" +echo "=== 5. Metrics ===" +poetry run python pipeline/metrics.py \ + "data/processed/dual_graphs/*/${CENSUS_GEOGRAPHY_TYPE}_in_${STUDY_AREA_TYPE}_*_${STUDY_AREA_DEFINITION_VINTAGE}_vintage_connected.json" \ + BLACK WHITE TOTPOP "${RUN_OUTPUT_DIR}/white_black.csv" + +poetry run python pipeline/metrics.py \ + "data/processed/dual_graphs/*/${CENSUS_GEOGRAPHY_TYPE}_in_${STUDY_AREA_TYPE}_*_${STUDY_AREA_DEFINITION_VINTAGE}_vintage_connected.json" \ + POC WHITE TOTPOP "${RUN_OUTPUT_DIR}/white_poc.csv" + +echo "" +echo "=== 6. Figures ===" for metric in white_black white_poc; do poetry run python pipeline/visualization/generate_figures.py \ --filename "${RUN_OUTPUT_DIR}/${metric}.csv" \ @@ -116,3 +79,4 @@ for metric in white_black white_poc; do --geography-type "${CENSUS_GEOGRAPHY_TYPE}" \ --study-area-type "${STUDY_AREA_TYPE}" done +echo "Saved to ${RUN_OUTPUT_DIR}/figures" diff --git a/scripts/resolve_config.py b/scripts/resolve_config.py deleted file mode 100644 index bb7e05c..0000000 --- a/scripts/resolve_config.py +++ /dev/null @@ -1,140 +0,0 @@ -#!/usr/bin/env python3 -""" -Read pipeline/config.yaml, apply normalization/validation/file-finding logic, -and print shell export statements to stdout. Intended to be consumed via: - - eval "$(poetry run python scripts/resolve_config.py)" - -Must be run from the repository root (all caller scripts guarantee this -via cd "${TOP_DIR}" before invoking). -""" -import glob -import os -import shlex -import sys - -import yaml - - -def _normalize(value: str, aliases: dict[str, str], valid: set[str], name: str) -> str: - value = aliases.get(value, value) - if value not in valid: - print(f"Unsupported {name}={value!r}. Valid values: {', '.join(sorted(valid))}.", file=sys.stderr) - sys.exit(1) - return value - - -def main() -> None: - with open("pipeline/config.yaml") as f: - cfg = yaml.safe_load(f) - - study_area_type = os.environ.get("STUDY_AREA_TYPE", str(cfg["study_area_type"])) - study_area_type = _normalize( - study_area_type, - {"counties": "county", - "max_counties": "max_county", - "max_cities": "max_city"}, - {"cbsa", "county", "max_city", "max_county"}, - "STUDY_AREA_TYPE", - ) - - census_geography_type = os.environ.get("CENSUS_GEOGRAPHY_TYPE", str(cfg["census_geography_type"])) - census_geography_type = _normalize( - census_geography_type, - {"tract": "tracts", "block_group": "block_groups", "block": "blocks", "county": "counties"}, - {"tracts", "block_groups", "blocks", "counties"}, - "CENSUS_GEOGRAPHY_TYPE", - ) - - years_raw = os.environ.get("CENSUS_GEOGRAPHY_YEARS", "") - years = years_raw.split() if years_raw else [str(y) for y in cfg["census_geography_years"]] - if census_geography_type in ("block_groups", "blocks") and "1980" in years: - print( - f"Warning: Skipping 1980 for CENSUS_GEOGRAPHY_TYPE={census_geography_type} — " - "NHGIS does not publish 1980 block group or block boundary shapefiles.", - file=sys.stderr, - ) - years = [y for y in years if y != "1980"] - census_geography_years = " ".join(years) - - study_area_vintage = os.environ.get("STUDY_AREA_VINTAGE", str(cfg.get("study_area_vintage", "2020"))) - - study_area_definition_geography_type = os.environ.get( - "STUDY_AREA_DEFINITION_GEOGRAPHY_TYPE", - str(cfg.get("study_area_definition_geography_type", - "places" if study_area_type == "max_city" else "counties")), - ) - - study_area_definition_geography_type = _normalize( - study_area_definition_geography_type, - {"tract": "tracts", "block_group": "block_groups", "block": "blocks", "county": "counties", "place": "places"}, - {"tracts", "block_groups", "blocks", "counties", "places"}, - "STUDY_AREA_DEFINITION_GEOGRAPHY_TYPE", - ) - - if study_area_type == "max_city" and study_area_definition_geography_type != "places": - print(f"Warning: STUDY_AREA_DEFINITION_GEOGRAPHY_TYPE={study_area_definition_geography_type!r} " - "is incompatible with STUDY_AREA_TYPE=max_city — switching to 'places'.", - file=sys.stderr) - - study_area_definition_geography_type = "places" - - study_area_definition_geography_year = os.environ.get( - "STUDY_AREA_DEFINITION_GEOGRAPHY_YEAR", study_area_vintage - ) - - study_area_source_pattern = os.environ.get( - "STUDY_AREA_SOURCE_PATTERN", f"list1_*{study_area_vintage}.xls" - ) - study_area_source_file = os.environ.get("STUDY_AREA_SOURCE_FILE", "") - study_area_definition_vintage = os.environ.get("STUDY_AREA_DEFINITION_VINTAGE", "") - - if study_area_type in ("cbsa", "max_city", "max_county"): - if not study_area_source_file: - matches = sorted(glob.glob(f"data/raw/study_area_sources/{study_area_source_pattern}")) - if not matches: - print( - f"No study area source file found for STUDY_AREA_TYPE={study_area_type}, " - f"STUDY_AREA_VINTAGE={study_area_vintage}, " - f"STUDY_AREA_SOURCE_PATTERN={study_area_source_pattern}", - file=sys.stderr, - ) - sys.exit(1) - study_area_source_file = matches[-1] - - if not study_area_definition_vintage: - stem = os.path.splitext(os.path.basename(study_area_source_file))[0] - study_area_definition_vintage = stem.removeprefix("list1_") - else: - study_area_definition_vintage = study_area_definition_vintage or study_area_vintage - - study_area_definition_geographies = os.environ.get( - "STUDY_AREA_DEFINITION_GEOGRAPHIES", - f"data/processed/census_geographies/{study_area_definition_geography_year}_{study_area_definition_geography_type}.gpkg", - ) - - run_name = os.environ.get("RUN_NAME", f"{census_geography_type}_in_{study_area_type}") - run_output_dir = os.environ.get("RUN_OUTPUT_DIR", f"outputs/{run_name}") - - exports = { - "STUDY_AREA_TYPE": study_area_type, - "CENSUS_GEOGRAPHY_TYPE": census_geography_type, - "CENSUS_GEOGRAPHY_YEARS": census_geography_years, - "STUDY_AREA_VINTAGE": study_area_vintage, - "STUDY_AREA_DEFINITION_GEOGRAPHY_TYPE": study_area_definition_geography_type, - "STUDY_AREA_DEFINITION_GEOGRAPHY_YEAR": study_area_definition_geography_year, - "STUDY_AREA_SOURCE_PATTERN": study_area_source_pattern, - "STUDY_AREA_SOURCE_FILE": study_area_source_file, - "STUDY_AREA_DEFINITION_VINTAGE": study_area_definition_vintage, - "STUDY_AREA_DEFINITION_GEOGRAPHIES": study_area_definition_geographies, - "OUTPUT_SUFFIX": f"{study_area_type}_{census_geography_type}_{study_area_definition_vintage}", - "RUN_NAME": run_name, - "RUN_OUTPUT_DIR": run_output_dir, - } - - for key, value in exports.items(): - print(f"export {key}={shlex.quote(value)}") - - -if __name__ == "__main__": - main() diff --git a/scripts/run_overlap.sh b/scripts/run_overlap.sh deleted file mode 100644 index 65373bd..0000000 --- a/scripts/run_overlap.sh +++ /dev/null @@ -1,32 +0,0 @@ -#!/usr/bin/env bash - -if (( $# != 4 )); then - printf 'Usage: %s YEAR GEOGRAPHY_TYPE STUDY_AREA_TYPE DEFINITION_VINTAGE\n' "$0" >&2 - exit 2 -fi - -year=$1 -census_geography_type=$2 -study_area_type=$3 -definition_vintage=$4 - -SCRIPT_DIR="$( - cd -- "$(dirname -- "${BASH_SOURCE[0]}")" >/dev/null 2>&1 - pwd -P -)" -TOP_DIR="$(cd -- "${SCRIPT_DIR}/.." >/dev/null 2>&1 && pwd -P)" -cd "${TOP_DIR}" - -output_dir="data/processed/clipped_geographies/${year}" -census_file="data/processed/census_geographies/${year}_${census_geography_type}.gpkg" -definition_pattern="data/processed/study_area_definitions/${study_area_type}_*_${definition_vintage}.gpkg" - -mkdir -p "${output_dir}" - -poetry run python pipeline/preprocessing/overlaps.py \ - "${census_file}" \ - "${definition_pattern}" \ - "${output_dir}" \ - --census-geography-type "${census_geography_type}" \ - --census-geography-year "${year}" \ - --definition-vintage "${definition_vintage}" diff --git a/scripts/setup.sh b/scripts/setup.sh index f6c1b47..9024f83 100644 --- a/scripts/setup.sh +++ b/scripts/setup.sh @@ -1,6 +1,5 @@ #!/usr/bin/env bash - SCRIPT_DIR="$( cd -- "$(dirname -- "${BASH_SOURCE[0]}")" >/dev/null 2>&1 pwd -P @@ -8,7 +7,7 @@ SCRIPT_DIR="$( TOP_DIR="$(cd -- "${SCRIPT_DIR}/.." >/dev/null 2>&1 && pwd -P)" cd "${TOP_DIR}" -_config="$(poetry run python scripts/resolve_config.py)" || exit 1 +_config="$(poetry run python pipeline/config.py)" || exit 1 eval "${_config}" IFS=" " read -r -a census_geography_years <<< "${CENSUS_GEOGRAPHY_YEARS}"