Validate ResStock gas consumption against EIA and add adjustment factor

.env.example

@@ -10,6 +10,9 @@
 # EIA data (e.g. data/eia/hourly_loads/, data/eia/861/)
 # EIA_API_KEY=
 
+# Census data API (ACS block-group heating fuel for CT gas utility assignment)
+# CENSUS_API_KEY=
+
 # OpenEI / NREL (e.g. URDB tariff data)
 # OPENEI_API_KEY=
 

context/code/data/resstock_sb_release_pipeline_main_py.md

@@ -47,6 +47,8 @@ Release-level defaults are loaded from `data/resstock/config.yaml`. State-specif
 | `--add-vulnerability-columns`  | per-state from `state_configs.yaml`            | Add LMI columns; defaults True for NY, False for RI. Pass True/False to override.   |
 | `--approximate-non-hp-load`    | `True`                                         | Run k-nearest-neighbor HVAC substitution for each upgrade in `approx_upgrade_ids`   |
 | `--adjust-mf-electricity`      | `True`                                         | Apply MF non-HVAC electricity adjustment (00 and 02)                                |
+| `--adjust-gas-usage`           | `True`                                         | Step 2c-iii: scale hourly gas of all requested upgrades by one statewide factor     |
+| `--gas-usage-eia-year`         | `2018`                                         | EIA-176 year for the gas factor (res. sales + transport ÷ weighted u00 ResStock)    |
 | `--assign-utility`             | `True`                                         | Assign electric/gas utilities (NY, RI only)                                         |
 | `--electric-poly-filename`     | from `state_configs.yaml`                      | Electric utility polygon CSV; overrides config default                              |
 | `--gas-poly-filename`          | from `state_configs.yaml`                      | Gas utility polygon CSV; overrides config default                                   |

data/resstock/load_curve/adjust_gas_usage.py

@@ -0,0 +1,152 @@
+"""Scale ResStock hourly natural gas load curves to EIA-176 residential deliveries.
+
+The scale factor is statewide: EIA-176 residential sales plus transportation
+(all operators in the state, Mcf converted to kWh with each operator's
+delivered-gas heat content) divided by weighted upgrade-00 ResStock gas use. The same factor is applied to
+every requested upgrade so that retrofit gas savings are scaled consistently
+with the baseline.
+"""
+
+from __future__ import annotations
+
+from concurrent.futures import ThreadPoolExecutor
+from pathlib import Path
+
+import polars as pl
+
+from data.eia.constants import PUDL_STABLE_VERSION
+
+BLDG_ID_COL = "bldg_id"
+WEIGHT_COL = "weight"
+ANNUAL_GAS_TOTAL_COL = "out.natural_gas.total.energy_consumption.kwh"
+
+GAS_END_USES = (
+    "clothes_dryer",
+    "fireplace",
+    "grill",
+    "heating",
+    "heating_hp_bkup",
+    "hot_water",
+    "lighting",
+    "permanent_spa_heat",
+    "pool_heater",
+    "range_oven",
+    "total",
+)
+HOURLY_GAS_COLS = tuple(
+    f"out.natural_gas.{end_use}.{suffix}"
+    for end_use in GAS_END_USES
+    for suffix in ("energy_consumption", "energy_consumption_intensity")
+)
+
+# 1 MMBtu = 1e6 Btu; 1 kWh = 3412.141633 Btu.
+KWH_PER_MMBTU = 1_000_000 / 3412.141633
+
+_PUDL_PARQUET = (
+    f"https://s3.us-west-2.amazonaws.com/pudl.catalyst.coop/{PUDL_STABLE_VERSION}"
+)
+EIA176_BY_CONSUMER = (
+    f"{_PUDL_PARQUET}/core_eia176__yearly_gas_disposition_by_consumer.parquet"
+)
+EIA176_DISPOSITION = f"{_PUDL_PARQUET}/core_eia176__yearly_gas_disposition.parquet"
+
+
+def load_eia176_residential_kwh(
+    state: str,
+    year: int,
+    path_by_consumer: str = EIA176_BY_CONSUMER,
+    path_disposition: str = EIA176_DISPOSITION,
+) -> float:
+    """Statewide EIA-176 residential deliveries (sales + transportation) in kWh.
+
+    Transportation volumes are gas the utility delivers to residential customers
+    who buy the commodity from a marketer; ResStock models their use too.
+    """
+    sales = (
+        pl.scan_parquet(path_by_consumer)
+        .filter(
+            (pl.col("operating_state") == state)
+            & (pl.col("report_year") == year)
+            & (pl.col("customer_class") == "residential")
+            & pl.col("revenue_class").cast(pl.String).is_in(["sales", "transport"])
+        )
+        .group_by("operator_id_eia")
+        .agg(pl.col("volume_mcf").cast(pl.Float64).sum())
+    )
+    heat = (
+        pl.scan_parquet(path_disposition)
+        .filter((pl.col("operating_state") == state) & (pl.col("report_year") == year))
+        .select(
+            "operator_id_eia",
+            pl.col("delivered_gas_heat_content_mmbtu_per_mcf").cast(pl.Float64),
+        )
+    )
+    joined = sales.join(heat, on="operator_id_eia", how="left").collect()
+    if joined.height == 0:
+        raise ValueError(f"EIA-176 has no residential deliveries for {state} {year}.")
+    missing_heat = joined.filter(
+        pl.col("delivered_gas_heat_content_mmbtu_per_mcf").is_null()
+    )["operator_id_eia"].to_list()
+    if missing_heat:
+        raise ValueError(
+            f"EIA-176 is missing delivered-gas heat content for {state} {year} "
+            f"operator(s) {missing_heat}."
+        )
+    return float(
+        joined.select(
+            (
+                pl.col("volume_mcf")
+                * pl.col("delivered_gas_heat_content_mmbtu_per_mcf")
+                * KWH_PER_MMBTU
+            ).sum()
+        ).item()
+    )
+
+
+def resstock_weighted_gas_kwh(load_curve_annual: pl.LazyFrame) -> float:
+    """Weighted sum of annual total natural gas kWh across all buildings."""
+    schema = load_curve_annual.collect_schema().names()
+    for col in (ANNUAL_GAS_TOTAL_COL, WEIGHT_COL):
+        if col not in schema:
+            raise ValueError(f"Load curve annual is missing column {col!r}.")
+    return float(
+        load_curve_annual.select(
+            (pl.col(ANNUAL_GAS_TOTAL_COL).fill_null(0) * pl.col(WEIGHT_COL)).sum()
+        )
+        .collect()
+        .item()
+    )
+
+
+def scale_gas_hourly(load_curve_hourly: pl.LazyFrame, factor: float) -> pl.LazyFrame:
+    """Multiply every natural gas consumption and intensity column by ``factor``."""
+    schema = load_curve_hourly.collect_schema().names()
+    missing = [c for c in HOURLY_GAS_COLS if c not in schema]
+    if missing:
+        raise ValueError(f"Hourly load curve is missing gas columns {missing}.")
+    return load_curve_hourly.with_columns(
+        (pl.col(c) * factor).alias(c) for c in HOURLY_GAS_COLS
+    )
+
+
+def adjust_gas_usage_hourly_dir(
+    load_curve_hourly_dir: Path, factor: float, max_workers: int = 128
+) -> int:
+    """Scale gas columns in every hourly parquet in the directory, in place.
+
+    Returns the number of files rewritten.
+    """
+    paths = sorted(load_curve_hourly_dir.glob("*.parquet"))
+
+    def _scale_one(path: Path) -> None:
+        # Collect before writing: source and sink are the same file.
+        scale_gas_hourly(pl.scan_parquet(str(path)), factor).collect().write_parquet(
+            str(path)
+        )
+
+    n_files = len(paths)
+    with ThreadPoolExecutor(max_workers=max_workers) as executor:
+        for i, _ in enumerate(executor.map(_scale_one, paths), 1):
+            if i % 500 == 0 or i == n_files:
+                print(f"    {i} out of {n_files} files scaled", flush=True)
+    return n_files

data/resstock/main.py

@@ -51,6 +51,11 @@
     BUILDING_TYPE_RECS_COL,
     adjust_mf_electricity_parquet,
 )
+from data.resstock.load_curve.adjust_gas_usage import (
+    adjust_gas_usage_hourly_dir,
+    load_eia176_residential_kwh,
+    resstock_weighted_gas_kwh,
+)
 from data.resstock.load_curve.approximate_non_hp_load import (
     _find_nearest_neighbors,
     _identify_non_hp_mf,
@@ -89,6 +94,7 @@
     validate_metadata_columns,
     validate_metadata_output,
     validate_metadata_readable,
+    validate_gas_usage_inputs,
     validate_no_stale_aggregate_loads,
     validate_s3_objects,
     validate_utility_assignment_args,
@@ -413,6 +419,70 @@ def _adjust_mf_electricity(
     return processed_upgrades
 
 
+def _adjust_gas_usage(
+    *,
+    states: list[str],
+    path_raw: Path,
+    path_sb: Path,
+    upgrade_ids: list[str],
+    eia_year: int,
+    sample: int,
+) -> dict[str, float]:
+    """Scale _sb hourly natural gas columns so weighted ResStock gas matches EIA-176.
+
+    For each state, the factor is EIA-176 residential sales plus transportation
+    in ``eia_year``
+    divided by weighted upgrade-00 gas from the raw ``load_curve_annual``. That
+    one factor is applied to the hourly gas columns of every requested upgrade.
+
+    Returns the factor applied per state.
+    """
+    factors: dict[str, float] = {}
+    for s in states:
+        lca_dir = path_raw / "load_curve_annual" / f"state={s}" / "upgrade=00"
+        if not any(lca_dir.glob("*.parquet")):
+            raise RuntimeError(
+                f"Upgrade 00 load_curve_annual not found at {lca_dir}; cannot "
+                f"compute the gas scale factor for state={s}."
+            )
+        if sample > 0:
+            print(
+                "  WARNING: --sample active — weighted ResStock gas covers only the "
+                "sampled buildings, so the gas scale factor will be too large. "
+                "Run without --sample for production use.",
+                flush=True,
+            )
+
+        resstock_kwh = resstock_weighted_gas_kwh(pl.scan_parquet(str(lca_dir)))
+        eia_kwh = load_eia176_residential_kwh(state=s, year=eia_year)
+        factor = eia_kwh / resstock_kwh
+        print(
+            f"  state={s}: EIA-176 {eia_year} residential sales + transport "
+            f"{eia_kwh:,.0f} kWh / "
+            f"ResStock upgrade 00 {resstock_kwh:,.0f} kWh = factor {factor:.4f}",
+            flush=True,
+        )
+
+        for uid in upgrade_ids:
+            uid_padded = uid.zfill(2)
+            lc_hourly_dir = (
+                path_sb / "load_curve_hourly" / f"state={s}" / f"upgrade={uid_padded}"
+            )
+            if not lc_hourly_dir.exists():
+                print(
+                    f"  WARNING: {lc_hourly_dir} not found, skipping "
+                    f"state={s} upgrade={uid_padded}.",
+                    flush=True,
+                )
+                continue
+            print(f"  Scaling gas for state={s} upgrade={uid_padded}...", flush=True)
+            adjust_gas_usage_hourly_dir(lc_hourly_dir, factor)
+        factors[s] = factor
+        gc.collect()
+
+    return factors
+
+
 def _assign_utility(
     *,
     states: list[str],
@@ -852,6 +922,27 @@ def _parse_args(argv: list[str] | None = None) -> argparse.Namespace:
             "load_curve_hourly and load_curve_annual in --file-types."
         ),
     )
+    parser.add_argument(
+        "--adjust-gas-usage",
+        type=parse_bool,
+        default=True,
+        metavar="BOOL",
+        help=(
+            "Scale hourly natural gas columns of every requested upgrade by one "
+            "statewide factor: EIA-176 residential sales + transport in "
+            "--gas-usage-eia-year "
+            "divided by weighted upgrade-00 ResStock gas (default: True). Runs "
+            "after adjust-mf-electricity. Requires load_curve_hourly and "
+            "load_curve_annual in --file-types."
+        ),
+    )
+    parser.add_argument(
+        "--gas-usage-eia-year",
+        type=int,
+        default=2018,
+        metavar="YEAR",
+        help="EIA-176 report year to scale ResStock gas to (default: 2018).",
+    )
     parser.add_argument(
         "--assign-utility",
         type=parse_bool,
@@ -974,6 +1065,8 @@ def main(argv: list[str] | None = None) -> None:
             },
             "approximate_non_hp_load": args.approximate_non_hp_load,
             "adjust_mf_electricity": args.adjust_mf_electricity,
+            "adjust_gas_usage": args.adjust_gas_usage,
+            "gas_usage_eia_year": args.gas_usage_eia_year,
             "assign_utility": args.assign_utility,
             "add_monthly_loads": args.add_monthly_loads,
             "add_annual_loads": args.add_annual_loads,
@@ -1005,6 +1098,7 @@ def main(argv: list[str] | None = None) -> None:
             approx_upgrades=_APPROX_UPGRADES,
             adjust_mf_electricity=args.adjust_mf_electricity,
             mf_adj_upgrades=_MF_ADJ_UPGRADES,
+            adjust_gas_usage=args.adjust_gas_usage,
             assign_utility=args.assign_utility,
             add_monthly_loads=args.add_monthly_loads,
             add_annual_loads=args.add_annual_loads,
@@ -1026,6 +1120,13 @@ def main(argv: list[str] | None = None) -> None:
             utility_assign_upgrade=_UTILITY_ASSIGN_UPGRADE,
             supported_states=SUPPORTED_UTILITY_STATES,
         )
+        validate_gas_usage_inputs(
+            states=args.state,
+            upgrade_ids=args.upgrade_ids,
+            file_types=args.file_types,
+            adjust_gas_usage=args.adjust_gas_usage,
+            path_raw=path_raw,
+        )
 
         # ── 1. Fetch ──────────────────────────────────────────────────────────
         # ── 1a. Fetch raw ResStock data ────────────────────────────────────────
@@ -1216,6 +1317,27 @@ def main(argv: list[str] | None = None) -> None:
             upsert_run(path_sb, run)
             gc.collect()
 
+        # ── 2c-iii. Scale natural gas to EIA-176 residential deliveries ───────
+        if args.adjust_gas_usage and "load_curve_hourly" in args.file_types:
+            print("Adjusting natural gas usage to EIA-176...", flush=True)
+            gas_factors = _adjust_gas_usage(
+                states=args.state,
+                path_raw=path_raw,
+                path_sb=path_sb,
+                upgrade_ids=args.upgrade_ids,
+                eia_year=args.gas_usage_eia_year,
+                sample=args.sample,
+            )
+            record_step(
+                run,
+                "adjust_gas_usage",
+                eia_year=args.gas_usage_eia_year,
+                scale_factors=gas_factors,
+                upgrades=[u.zfill(2) for u in args.upgrade_ids],
+            )
+            upsert_run(path_sb, run)
+            gc.collect()
+
         # ── 2d. Aggregate load curves (monthly / annual) ────────────────────
         _wants_aggregation = (
             args.add_monthly_loads or args.add_annual_loads

data/resstock/manifest.py

@@ -499,6 +499,7 @@ def print_status(
     "assign_utility": frozenset({"metadata_utility"}),
     "approximate_non_hp_load": frozenset({"load_curve_hourly"}),
     "adjust_mf_electricity": frozenset({"load_curve_hourly"}),
+    "adjust_gas_usage": frozenset({"load_curve_hourly"}),
     "add_monthly_loads": frozenset({"load_curve_monthly"}),
     "add_annual_loads": frozenset({"load_curve_annual"}),
     "add_aggregate_loads": frozenset({"load_curve_monthly", "load_curve_annual"}),