Python script for normalizing demographic data across ZIP codes

import os
import logging
import pandas as pd
import numpy as np
import geopandas as gpd
from pathlib import Path
from typing import Optional, Dict, List
import requests
from requests.adapters import HTTPAdapter
from urllib3.util.retry import Retry
from sklearn.impute import SimpleImputer
from shapely.geometry import Point

logging.basicConfig(
    level=os.getenv("LOG_LEVEL", "INFO"),
    format="%(asctime)s [%(levelname)s] %(name)s: %(message)s"
)
logger = logging.getLogger(__name__)


class CensusAPIError(Exception):
    pass


class SpatialAlignmentError(Exception):
    pass


class NormalizationError(Exception):
    pass


class DemographicNormalizer:
    def __init__(self, census_api_key: str, cache_dir: str = "./data_cache"):
        self.api_key = census_api_key
        self.cache_dir = Path(cache_dir)
        self.cache_dir.mkdir(parents=True, exist_ok=True)
        self.session = self._configure_session()
        # ACS 5-year 2023 release; update year annually
        self.base_url = "https://api.census.gov/data/2023/acs/acs5"
        self.zcta_gdf: Optional[gpd.GeoDataFrame] = None

    def _configure_session(self) -> requests.Session:
        # Retry on throttling (429) and transient server errors with exponential backoff.
        session = requests.Session()
        retry_strategy = Retry(
            total=5,
            backoff_factor=2,
            status_forcelist=[429, 500, 502, 503, 504],
            allowed_methods=["GET"]
        )
        session.mount("https://", HTTPAdapter(max_retries=retry_strategy))
        return session

    def fetch_acs_for_zctas(self, variables: List[str], zcta_list: List[str]) -> pd.DataFrame:
        """
        Fetch ACS variables for a list of ZCTAs.
        The Census API supports querying specific ZCTAs via 'for=zip code tabulation area:XXXXX,...'
        Keep batches under ~50 ZCTAs to stay within URL length limits.
        """
        if not zcta_list:
            raise ValueError("Empty ZCTA list provided.")

        chunk_size = 50
        all_data = []

        for i in range(0, len(zcta_list), chunk_size):
            chunk = zcta_list[i:i + chunk_size]
            params = {
                "get": ",".join(variables),
                "for": f"zip code tabulation area:{','.join(chunk)}",
                "key": self.api_key,
            }
            try:
                resp = self.session.get(self.base_url, params=params, timeout=30)
                resp.raise_for_status()
                data = resp.json()
                if len(data) <= 1:  # header row only -> no records
                    logger.warning("No data returned for ZCTA chunk starting at index %d.", i)
                    continue
                df = pd.DataFrame(data[1:], columns=data[0])
                all_data.append(df)
            except requests.exceptions.RequestException as e:
                logger.error("API request failed for chunk %d: %s", i, e)
                raise CensusAPIError(f"Failed to fetch ACS data: {e}")

        if not all_data:
            return pd.DataFrame()
        return pd.concat(all_data, ignore_index=True)

    def load_zcta_boundaries(self, shapefile_path: str) -> gpd.GeoDataFrame:
        # Cache the parsed shapefile to Parquet; reading the cache avoids re-parsing on every run.
        cache_path = self.cache_dir / "zcta_index.parquet"
        if cache_path.exists():
            logger.info("Loading cached ZCTA boundaries.")
            return gpd.read_parquet(cache_path)

        logger.info("Loading ZCTA shapefile from %s", shapefile_path)
        gdf = gpd.read_file(shapefile_path)
        gdf = gdf.to_crs(epsg=4326)  # assert geographic CRS for the intersect join
        gdf.to_parquet(cache_path)
        return gdf

    def zip_to_zcta(self, zip_codes: List[str], zcta_gdf: gpd.GeoDataFrame) -> pd.DataFrame:
        """
        Map postal ZIPs to ZCTAs.
        The preferred production approach is the HUD USPS ZIP-to-ZCTA crosswalk, a
        deterministic table join without spatial computation. The spatial approach
        below is a fallback when centroids for input ZIPs are available.
        """
        # In production: load a ZIP centroid file and do a point-in-polygon join.
        # Here we demonstrate the structure; replace dummy Point(0,0) with real centroids.
        zip_points = gpd.GeoDataFrame(
            {"zip_code": zip_codes},
            geometry=[Point(0, 0) for _ in zip_codes],  # replace with real ZIP centroids
            crs="EPSG:4326"
        )
        logger.warning(
            "Using placeholder ZIP centroids. Replace with real centroid data or "
            "HUD USPS crosswalk for production use."
        )
        try:
            joined = gpd.sjoin(zip_points, zcta_gdf, how="left", predicate="intersects")
            # TIGER/Line ZCTA shapefile uses ZCTA5CE20 for 2020-vintage boundaries
            zcta_col = "ZCTA5CE20" if "ZCTA5CE20" in joined.columns else joined.columns[-1]
            return joined[["zip_code", zcta_col]].drop_duplicates().rename(
                columns={zcta_col: "zcta"}
            )
        except Exception as e:
            raise SpatialAlignmentError(f"Spatial join failed: {e}")

    def impute_nulls(self, df: pd.DataFrame, numeric_cols: List[str]) -> pd.DataFrame:
        """Applies median imputation for sparse rural geographies."""
        imputer = SimpleImputer(strategy="median")
        df = df.copy()
        df[numeric_cols] = imputer.fit_transform(df[numeric_cols])
        return df

    def normalize_and_weight(
        self,
        df: pd.DataFrame,
        numeric_cols: List[str],
        population_col: str,
        audience_weights: Optional[Dict[str, float]] = None
    ) -> pd.DataFrame:
        """Applies population-weighted min-max normalization and custom audience scaling."""
        try:
            weights = np.array(df[population_col].values, dtype=float)
            weights = np.where(weights <= 0, 1e-6, weights)  # guard against divide-by-zero

            norm_df = df.copy()
            for col in numeric_cols:
                col_min = norm_df[col].min()
                col_max = norm_df[col].max()
                denom = col_max - col_min
                if denom == 0:  # constant column -> no spread to scale
                    norm_df[col] = 0.0
                else:
                    norm_df[col] = (norm_df[col] - col_min) / denom

                # Population-weighted score
                norm_df[f"{col}_weighted"] = norm_df[col] * weights

            if audience_weights:
                for col, weight in audience_weights.items():
                    if col in norm_df.columns:
                        norm_df[col] *= weight
                        logger.info("Applied audience multiplier %.2f to %s", weight, col)

            return norm_df
        except Exception as e:
            raise NormalizationError(f"Normalization failed: {e}")

    def run_pipeline(
        self,
        target_zips: List[str],
        acs_variables: List[str],
        population_col: str = "B01003_001E",
        audience_weights: Optional[Dict[str, float]] = None
    ) -> pd.DataFrame:
        logger.info("Starting demographic normalization pipeline.")
        zcta_gdf = self.load_zcta_boundaries(
            str(self.cache_dir / "zcta_boundaries/tl_2022_us_zcta520.shp")
        )

        # 1. Map ZIPs to ZCTAs
        zip_zcta_map = self.zip_to_zcta(target_zips, zcta_gdf)
        target_zctas = zip_zcta_map["zcta"].dropna().unique().tolist()

        # 2. Fetch raw ACS data for the resolved ZCTAs
        raw_df = self.fetch_acs_for_zctas(acs_variables, target_zctas)
        if raw_df.empty:
            logger.critical("Pipeline halted: No ACS data retrieved.")
            return pd.DataFrame()

        # 3. Join ZIP -> ZCTA mapping back to ACS data
        zcta_col = "zip code tabulation area"
        merged = raw_df.merge(zip_zcta_map, left_on=zcta_col, right_on="zcta", how="inner")

        # 4. Imputation & normalization
        numeric_cols = [c for c in acs_variables if c != population_col]
        merged[numeric_cols] = merged[numeric_cols].apply(pd.to_numeric, errors="coerce")
        merged = self.impute_nulls(merged, numeric_cols)
        final_df = self.normalize_and_weight(merged, numeric_cols, population_col, audience_weights)

        logger.info("Pipeline complete. Processed %d geographies.", len(final_df))
        return final_df


if __name__ == "__main__":
    normalizer = DemographicNormalizer(
        census_api_key=os.getenv("CENSUS_API_KEY", ""),
        cache_dir=os.getenv("CACHE_DIR", "./data_cache")
    )

    output = normalizer.run_pipeline(
        target_zips=["90210", "10001", "60601"],
        acs_variables=["B01003_001E", "B19013_001E", "B15001_001E"],
        population_col="B01003_001E",
        audience_weights={"B19013_001E": 1.25}
    )
    print(output.head())