diff --git a/.env.example b/.env.example index 9018d4bd..e6136e4f 100644 --- a/.env.example +++ b/.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= diff --git a/context/code/data/resstock_sb_release_pipeline_main_py.md b/context/code/data/resstock_sb_release_pipeline_main_py.md index 7e7622ec..2dc016b9 100644 --- a/context/code/data/resstock_sb_release_pipeline_main_py.md +++ b/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 | diff --git a/data/resstock/load_curve/adjust_gas_usage.py b/data/resstock/load_curve/adjust_gas_usage.py new file mode 100644 index 00000000..8d20277f --- /dev/null +++ b/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 diff --git a/data/resstock/main.py b/data/resstock/main.py index a9fcfa5b..b22eb867 100644 --- a/data/resstock/main.py +++ b/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 diff --git a/data/resstock/manifest.py b/data/resstock/manifest.py index d3a44fbf..c3bc42f1 100644 --- a/data/resstock/manifest.py +++ b/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"}), diff --git a/data/resstock/utility/assign_utility_ct.py b/data/resstock/utility/assign_utility_ct.py index 1ee64353..2925d5ea 100644 --- a/data/resstock/utility/assign_utility_ct.py +++ b/data/resstock/utility/assign_utility_ct.py @@ -25,6 +25,16 @@ per-PUMA utility probability distributions, which are then used to sample a utility assignment for each ResStock building. +Electric probabilities are each utility's share of the PUMA's land area. Gas +probabilities are each utility's share of the PUMA's homes heated with utility +gas (ACS 2016-2020 table B25040 by block group, spread to 2020 census blocks +by housing units): every block is placed, by its internal point, in a PUMA and +a gas territory. The same gas probabilities are used for every building with a +gas connection, including ones that do not heat with gas, since gas-heated +homes mark where the mains are. Land area and total housing units both +over-assign Yankee Gas, whose territory covers much of the state but has a low +share of homes on gas. + CT electric utilities (IOUs + municipals): - Eversource Energy CT (std_name ``ct_eversource``; HIFLD still ``CONNECTICUT LIGHT & POWER CO``) @@ -51,23 +61,27 @@ import geopandas as gpd import polars as pl -from data.resstock.utils import ( - load_state_configs, - select_puma_and_heating_fuel_metadata, -) from data.resstock.utility.utils import ( GIS_CACHE_DIR, + add_block_gas_heated_homes, calculate_prior_distributions, + calculate_puma_utility_housing_overlap, calculate_puma_utility_overlap, calculate_utility_probabilities, + fetch_acs_gas_heated_homes_by_block_group, fill_missing_puma_probabilities, filter_hifld_for_state, + load_census_blocks, load_national_hifld, load_pumas, print_comparison_summary, sample_utility_per_building, zero_excluded_gas_utilities_and_renormalize, ) +from data.resstock.utils import ( + load_state_configs, + select_puma_and_heating_fuel_metadata, +) # ── CT-specific constants ───────────────────────────────────────────────────── @@ -77,6 +91,7 @@ CT_STATE_CRS: int = _CT_CFG["state_crs"] CT_PUMA_YEAR: int = _CT_CFG["puma_year"] +CT_STATE_FIPS: str = _STATE_CONFIGS[_STATE]["state_fips"] # ── Pipeline entry point ────────────────────────────────────────────────────── @@ -131,12 +146,19 @@ def assign_utility( ) pumas = pumas.to_crs(epsg=state_crs) + print(" Loading CT 2020 Census blocks and ACS gas-heated homes ...", flush=True) + blocks = add_block_gas_heated_homes( + load_census_blocks(_STATE), + fetch_acs_gas_heated_homes_by_block_group(CT_STATE_FIPS), + ) + return assign_utility_ct( input_metadata=metadata, electric_polygons=elec_ct, gas_polygons=gas_ct, pumas=pumas, state_crs=state_crs, + gas_blocks=blocks, excluded_gas_utilities=frozenset(excluded_gas_utilities) if excluded_gas_utilities is not None else frozenset(), @@ -153,11 +175,13 @@ def assign_utility_ct( pumas: gpd.GeoDataFrame, state_crs: int, excluded_gas_utilities: frozenset[str] = frozenset(), + gas_blocks: gpd.GeoDataFrame | None = None, ) -> pl.LazyFrame: """Assign electric and gas utilities to ResStock buildings in CT. - Both electric and gas utilities are assigned via PUMA-polygon area - overlap on HIFLD service territory shapes. + Electric utilities are assigned via PUMA-polygon area overlap on HIFLD + service territory shapes. Gas utilities use gas-heated-home overlap when + ``gas_blocks`` is given, and area overlap otherwise. Args: input_metadata: ResStock metadata LazyFrame. @@ -171,6 +195,9 @@ def assign_utility_ct( state_crs: EPSG code for the CT projected CRS. excluded_gas_utilities: Gas utility names whose PUMA probabilities are zeroed before sampling (default: empty). + gas_blocks: Census block internal points with ``gas_heated_homes`` + (from :func:`load_census_blocks` and + :func:`add_block_gas_heated_homes`). Returns: LazyFrame with all original metadata columns plus @@ -190,7 +217,14 @@ def assign_utility_ct( puma_elec_overlap = calculate_puma_utility_overlap( pumas, electric_polygons, state_crs ) - puma_gas_overlap = calculate_puma_utility_overlap(pumas, gas_polygons, state_crs) + if gas_blocks is not None: + puma_gas_overlap = calculate_puma_utility_housing_overlap( + pumas, gas_polygons, gas_blocks, state_crs, weight_col="gas_heated_homes" + ) + else: + puma_gas_overlap = calculate_puma_utility_overlap( + pumas, gas_polygons, state_crs + ) puma_elec_probs = calculate_utility_probabilities( puma_elec_overlap, diff --git a/data/resstock/utility/utils.py b/data/resstock/utility/utils.py index c7f98933..3d2252ee 100644 --- a/data/resstock/utility/utils.py +++ b/data/resstock/utility/utils.py @@ -35,16 +35,16 @@ import polars as pl import yaml from cloudpathlib import S3Path +from dotenv import load_dotenv +from pygris import blocks as pygris_blocks +from pygris import counties as pygris_counties from pygris import pumas as pygris_pumas from ruamel.yaml import YAML -from pygris import counties as pygris_counties - from data.resstock.constants import CONFIG_PATH, STATE_CONFIGS_PATH -from utils import get_aws_region +from utils import get_aws_region, get_project_root from utils.utility_codes import get_hifld_to_std_name - # --------------------------------------------------------------------------- # GIS / I/O helpers # --------------------------------------------------------------------------- @@ -222,6 +222,164 @@ def calculate_puma_utility_overlap( return pl.from_pandas(puma_overlap).lazy() +def load_census_blocks(state: str, year: int = 2020) -> gpd.GeoDataFrame: + """Fetch Census tabulation blocks with full-count housing units via pygris. + + TIGER/Line ``tabblock20`` carries ``HOUSING20`` (2020 decennial housing + units) and the Census internal point (``INTPTLAT20`` / ``INTPTLON20``), + so no Census API key is needed. + + Returns: + GeoDataFrame with one point geometry per block (its internal point) + and a ``housing_units`` column, in the TIGER CRS. + """ + if year != 2020: + raise ValueError(f"Only 2020 blocks (HOUSING20) are supported, got {year}.") + raw = cast(gpd.GeoDataFrame, pygris_blocks(state=state, year=year, cache=True)) + points = gpd.points_from_xy( + raw["INTPTLON20"].astype(float), raw["INTPTLAT20"].astype(float) + ) + return gpd.GeoDataFrame( + { + "block_geoid": raw["GEOID20"].astype(str), + "housing_units": raw["HOUSING20"].astype(int), + }, + geometry=points, + crs=raw.crs, + ) + + +ACS5_BLOCK_GROUP_YEAR = 2020 +ACS_HEATING_FUEL_TOTAL = "B25040_001E" +ACS_HEATING_FUEL_UTILITY_GAS = "B25040_002E" + + +def fetch_acs_gas_heated_homes_by_block_group( + state_fips: str, year: int = ACS5_BLOCK_GROUP_YEAR +) -> pl.DataFrame: + """ACS 5-year occupied homes heated with utility gas, by block group. + + Table B25040 (House Heating Fuel). The 2016-2020 release is tabulated on + 2020 block groups, so 2020 census blocks nest inside them exactly. + Requires ``CENSUS_API_KEY`` in the environment or the project ``.env``. + + Returns: + DataFrame with ``block_group_geoid`` (12 chars) and + ``gas_heated_homes``. + """ + load_dotenv(dotenv_path=get_project_root() / ".env") + key = os.getenv("CENSUS_API_KEY") + if not key: + raise ValueError( + "CENSUS_API_KEY not set. Add it to the project .env file " + "(free key: https://api.census.gov/data/key_signup.html)." + ) + resp = get_with_retry( + f"https://api.census.gov/data/{year}/acs/acs5", + { + "get": f"{ACS_HEATING_FUEL_TOTAL},{ACS_HEATING_FUEL_UTILITY_GAS}", + "for": "block group:*", + "in": [f"state:{state_fips}", "county:*", "tract:*"], + "key": key, + }, + ) + rows = resp.json() + return pl.DataFrame(rows[1:], schema=rows[0], orient="row").select( + pl.concat_str("state", "county", "tract", "block group").alias( + "block_group_geoid" + ), + pl.col(ACS_HEATING_FUEL_UTILITY_GAS) + .cast(pl.Int64) + .clip(lower_bound=0) + .alias("gas_heated_homes"), + ) + + +def add_block_gas_heated_homes( + blocks: gpd.GeoDataFrame, block_groups: pl.DataFrame +) -> gpd.GeoDataFrame: + """Spread each block group's gas-heated homes over its blocks. + + Each block gets the block group's count times the block's share of the + block group's housing units. A block group with gas-heated homes but no + housing units in the 2020 blocks spreads them equally over its blocks. + """ + block_df = pl.DataFrame( + { + "block_geoid": blocks["block_geoid"].tolist(), + "housing_units": blocks["housing_units"].tolist(), + } + ).with_columns(pl.col("block_geoid").str.slice(0, 12).alias("block_group_geoid")) + bg_hu = pl.col("housing_units").sum().over("block_group_geoid") + allocated = block_df.join( + block_groups, on="block_group_geoid", how="left" + ).with_columns( + ( + pl.col("gas_heated_homes").fill_null(0) + * pl.when(bg_hu > 0) + .then(pl.col("housing_units") / bg_hu) + .otherwise(1 / pl.len().over("block_group_geoid")) + ).alias("gas_heated_homes") + ) + out = blocks.copy() + out["gas_heated_homes"] = allocated["gas_heated_homes"].to_numpy() + return out + + +def calculate_puma_utility_housing_overlap( + pumas: gpd.GeoDataFrame, + utility_gdf: gpd.GeoDataFrame, + blocks: gpd.GeoDataFrame, + state_crs: int, + weight_col: str = "housing_units", +) -> pl.LazyFrame: + """PUMA-utility overlap weighted by a block count instead of land area. + + Each block is placed whole, by its internal point, in one PUMA and in + every utility polygon that contains that point. A block inside more than + one utility polygon splits its ``weight_col`` equally among them. + + Returns the same columns as :func:`calculate_puma_utility_overlap` + (``puma_id``, ``utility``, ``pct_overlap``), where ``pct_overlap`` is the + utility's share of the PUMA's ``weight_col`` total, in percent. + """ + pumas_proj = pumas[["PUMACE10", "geometry"]].to_crs(epsg=state_crs) + utilities_proj = ( + utility_gdf[["utility", "geometry"]].to_crs(epsg=state_crs).dissolve("utility") + ).reset_index() + blocks_proj = blocks[["block_geoid", weight_col, "geometry"]].to_crs(epsg=state_crs) + + blocks_in_puma = gpd.sjoin( + blocks_proj, pumas_proj, how="inner", predicate="within" + ).drop(columns="index_right") + puma_totals = ( + pl.from_pandas(blocks_in_puma[["PUMACE10", weight_col]]) + .group_by("PUMACE10") + .agg(pl.col(weight_col).sum().alias("puma_total")) + ) + + in_utility = gpd.sjoin( + blocks_in_puma, utilities_proj, how="inner", predicate="within" + ) + block_utilities = pl.from_pandas( + in_utility[["block_geoid", "PUMACE10", "utility", weight_col]] + ).with_columns((pl.col(weight_col) / pl.len().over("block_geoid")).alias("share")) + + return ( + block_utilities.group_by("PUMACE10", "utility") + .agg(pl.col("share").sum()) + .join(puma_totals, on="PUMACE10", how="left") + .filter(pl.col("share") > 0) + .select( + pl.col("PUMACE10").alias("puma_id"), + pl.col("utility"), + (pl.col("share") / pl.col("puma_total") * 100).alias("pct_overlap"), + ) + .sort("puma_id", "utility") + .lazy() + ) + + def calculate_puma_county_utility_overlap( pumas: gpd.GeoDataFrame, county_gdf: gpd.GeoDataFrame, diff --git a/data/resstock/validations.py b/data/resstock/validations.py index dca3960c..7f76f73c 100644 --- a/data/resstock/validations.py +++ b/data/resstock/validations.py @@ -178,6 +178,42 @@ def validate_utility_assignment_args( ) +def validate_gas_usage_inputs( + states: list[str], + upgrade_ids: list[str], + file_types: list[str], + adjust_gas_usage: bool, + path_raw: Path, +) -> None: + """Raise RuntimeError if the gas scale factor's upgrade-00 annual file is unavailable. + + The factor is computed from the raw upgrade-00 ``load_curve_annual``. It is + available if this run fetches it (upgrade 00 in ``upgrade_ids`` and + ``load_curve_annual`` in ``file_types``) or if it is already on disk for + every requested state. Skipped when the step would not run anyway + (``load_curve_hourly`` not requested). + """ + if not adjust_gas_usage or "load_curve_hourly" not in file_types: + return + if "00" in [u.zfill(2) for u in upgrade_ids] and "load_curve_annual" in file_types: + return + + missing: list[str] = [] + for s in states: + lca_dir = path_raw / "load_curve_annual" / f"state={s}" / "upgrade=00" + if not any(lca_dir.glob("*.parquet")): + missing.append(str(lca_dir)) + if missing: + raise RuntimeError( + "--adjust-gas-usage needs the raw upgrade-00 load_curve_annual to " + "compute the gas scale factor, but it is not on disk and this run " + "will not fetch it:\n" + + "\n".join(f" {m}" for m in missing) + + "\nAdd 0 to --upgrade-ids and 'load_curve_annual' to --file-types, " + "or disable the step with --adjust-gas-usage False." + ) + + def validate_no_stale_aggregate_loads( state: list[str], upgrade_ids: list[str], diff --git a/data/resstock/warnings.py b/data/resstock/warnings.py index 93d498b2..66eebdba 100644 --- a/data/resstock/warnings.py +++ b/data/resstock/warnings.py @@ -22,6 +22,7 @@ def collect_run_warnings( assign_utility: bool, add_monthly_loads: bool, add_annual_loads: bool = True, + adjust_gas_usage: bool = False, ) -> list[str]: """Run all pre-run argument/file-type mismatch checks. @@ -48,6 +49,8 @@ def collect_run_warnings( Whether the monthly load aggregation step is enabled. add_annual_loads: Whether the annual load aggregation step is enabled. + adjust_gas_usage: + Whether the natural gas scaling to EIA-176 step is enabled. """ warnings: list[str] = [] @@ -96,6 +99,23 @@ def _warn(msg: str) -> None: f"step will be skipped. Check --upgrade-ids or mf_adj_upgrade_ids in config.yaml." ) + if adjust_gas_usage and "load_curve_hourly" not in file_types: + _warn( + "--adjust-gas-usage is enabled but 'load_curve_hourly' is not in " + "--file-types. The gas usage adjustment step will be skipped." + ) + + if adjust_gas_usage and ( + "load_curve_annual" not in file_types + or "00" not in [u.zfill(2) for u in upgrade_ids] + ): + _warn( + "--adjust-gas-usage is enabled but this run does not fetch the raw " + "upgrade-00 load_curve_annual (needs 0 in --upgrade-ids and " + "'load_curve_annual' in --file-types). The gas scale factor will use " + "the copy already on disk; pre-flight fails if there is none." + ) + if assign_utility and "metadata" not in file_types: _warn( "--assign-utility is enabled but 'metadata' is not in --file-types. " diff --git a/tests/test_adjust_gas_usage.py b/tests/test_adjust_gas_usage.py new file mode 100644 index 00000000..b784fe9a --- /dev/null +++ b/tests/test_adjust_gas_usage.py @@ -0,0 +1,171 @@ +"""Tests for data/resstock/load_curve/adjust_gas_usage.py.""" + +from __future__ import annotations + +from pathlib import Path + +import polars as pl +import pytest + +from data.resstock.load_curve.adjust_gas_usage import ( + ANNUAL_GAS_TOTAL_COL, + HOURLY_GAS_COLS, + KWH_PER_MMBTU, + adjust_gas_usage_hourly_dir, + load_eia176_residential_kwh, + resstock_weighted_gas_kwh, + scale_gas_hourly, +) +from data.resstock.validations import validate_gas_usage_inputs + +ELEC_COL = "out.electricity.total.energy_consumption" + + +def _hourly_frame() -> pl.DataFrame: + data: dict[str, list[float]] = {c: [1.0, 2.0] for c in HOURLY_GAS_COLS} + data[ELEC_COL] = [5.0, 6.0] + return pl.DataFrame(data) + + +def test_scale_gas_hourly_scales_all_gas_columns_only() -> None: + out = scale_gas_hourly(_hourly_frame().lazy(), 1.5).collect() + for col in HOURLY_GAS_COLS: + assert out[col].to_list() == [1.5, 3.0] + assert out[ELEC_COL].to_list() == [5.0, 6.0] + + +def test_hourly_gas_cols_cover_consumption_and_intensity() -> None: + assert len(HOURLY_GAS_COLS) == 22 + assert "out.natural_gas.total.energy_consumption" in HOURLY_GAS_COLS + assert "out.natural_gas.heating.energy_consumption_intensity" in HOURLY_GAS_COLS + + +def test_scale_gas_hourly_raises_on_missing_gas_column() -> None: + frame = _hourly_frame().drop("out.natural_gas.grill.energy_consumption") + with pytest.raises(ValueError, match="grill"): + scale_gas_hourly(frame.lazy(), 1.1) + + +def test_resstock_weighted_gas_kwh_treats_null_as_zero() -> None: + annual = pl.LazyFrame( + { + "bldg_id": [1, 2, 3], + ANNUAL_GAS_TOTAL_COL: [100.0, None, 50.0], + "weight": [2.0, 3.0, 4.0], + } + ) + assert resstock_weighted_gas_kwh(annual) == pytest.approx(400.0) + + +def test_resstock_weighted_gas_kwh_requires_weight() -> None: + annual = pl.LazyFrame({"bldg_id": [1], ANNUAL_GAS_TOTAL_COL: [1.0]}) + with pytest.raises(ValueError, match="weight"): + resstock_weighted_gas_kwh(annual) + + +def _write_eia(tmp_path: Path) -> tuple[str, str]: + by_consumer = tmp_path / "by_consumer.parquet" + disposition = tmp_path / "disposition.parquet" + pl.DataFrame( + { + "report_year": [2018, 2018, 2018, 2018, 2019, 2018], + "operator_id_eia": ["A", "A", "B", "B", "A", "C"], + "operating_state": ["CT", "CT", "CT", "CT", "CT", "NY"], + "customer_class": [ + "residential", + "residential", + "residential", + "commercial", + "residential", + "residential", + ], + "revenue_class": ["sales", "transport", "sales", "sales", "sales", "sales"], + "volume_mcf": [100.0, 10.0, 200.0, 999.0, 999.0, 999.0], + } + ).write_parquet(by_consumer) + pl.DataFrame( + { + "report_year": [2018, 2018, 2018], + "operator_id_eia": ["A", "B", "C"], + "operating_state": ["CT", "CT", "NY"], + "delivered_gas_heat_content_mmbtu_per_mcf": [1.0, 1.05, 1.0], + } + ).write_parquet(disposition) + return str(by_consumer), str(disposition) + + +def test_eia176_residential_includes_transport_and_excludes_other_rows( + tmp_path: Path, +) -> None: + by_consumer, disposition = _write_eia(tmp_path) + kwh = load_eia176_residential_kwh( + "CT", 2018, path_by_consumer=by_consumer, path_disposition=disposition + ) + assert kwh == pytest.approx(((100.0 + 10.0) * 1.0 + 200.0 * 1.05) * KWH_PER_MMBTU) + + +def test_eia176_residential_raises_when_no_rows(tmp_path: Path) -> None: + by_consumer, disposition = _write_eia(tmp_path) + with pytest.raises(ValueError, match="no residential deliveries"): + load_eia176_residential_kwh( + "RI", 2018, path_by_consumer=by_consumer, path_disposition=disposition + ) + + +def _validate(tmp_path: Path, upgrade_ids: list[str], file_types: list[str]) -> None: + validate_gas_usage_inputs( + states=["CT"], + upgrade_ids=upgrade_ids, + file_types=file_types, + adjust_gas_usage=True, + path_raw=tmp_path, + ) + + +def test_validate_gas_usage_passes_when_run_fetches_u00_annual( + tmp_path: Path, +) -> None: + _validate(tmp_path, ["0", "2"], ["load_curve_hourly", "load_curve_annual"]) + + +@pytest.mark.parametrize( + ("upgrade_ids", "file_types"), + [ + (["2"], ["load_curve_hourly", "load_curve_annual"]), + (["0", "2"], ["load_curve_hourly"]), + ], +) +def test_validate_gas_usage_raises_when_u00_annual_missing( + tmp_path: Path, upgrade_ids: list[str], file_types: list[str] +) -> None: + with pytest.raises(RuntimeError, match="upgrade-00 load_curve_annual"): + _validate(tmp_path, upgrade_ids, file_types) + + +def test_validate_gas_usage_passes_when_u00_annual_on_disk(tmp_path: Path) -> None: + lca_dir = tmp_path / "load_curve_annual" / "state=CT" / "upgrade=00" + lca_dir.mkdir(parents=True) + pl.DataFrame({"bldg_id": [1]}).write_parquet(lca_dir / "annual.parquet") + _validate(tmp_path, ["2"], ["load_curve_hourly"]) + + +def test_validate_gas_usage_skips_when_disabled_or_no_hourly(tmp_path: Path) -> None: + validate_gas_usage_inputs( + states=["CT"], + upgrade_ids=["2"], + file_types=["load_curve_hourly"], + adjust_gas_usage=False, + path_raw=tmp_path, + ) + _validate(tmp_path, ["2"], ["metadata"]) + + +def test_adjust_gas_usage_hourly_dir_rewrites_every_file(tmp_path: Path) -> None: + for bldg_id in (1, 2): + _hourly_frame().write_parquet(tmp_path / f"{bldg_id}-0.parquet") + n = adjust_gas_usage_hourly_dir(tmp_path, 2.0, max_workers=2) + assert n == 2 + for bldg_id in (1, 2): + out = pl.read_parquet(tmp_path / f"{bldg_id}-0.parquet") + assert out["out.natural_gas.total.energy_consumption"].to_list() == [2.0, 4.0] + assert out[ELEC_COL].to_list() == [5.0, 6.0] diff --git a/tests/test_compare_resstock_eia176_gas.py b/tests/test_compare_resstock_eia176_gas.py new file mode 100644 index 00000000..6fd4d0da --- /dev/null +++ b/tests/test_compare_resstock_eia176_gas.py @@ -0,0 +1,122 @@ +"""Tests for the ResStock vs EIA-176 gas comparison.""" + +from __future__ import annotations + +from pathlib import Path + +import polars as pl +import pytest + +from utils.post.compare_resstock_eia176_gas import ( + ANNUAL_GAS_COL, + EIA_RESIDENTIAL_CUSTOMERS, + EIA_RESIDENTIAL_KWH, + KWH_PER_MMBTU, + RESSTOCK_CUSTOMERS, + RESSTOCK_TOTAL_KWH, + compare_resstock_eia176_gas, + eia_year_for_resstock_release, +) + + +def test_eia_year_comes_from_amy_in_the_release_name() -> None: + assert eia_year_for_resstock_release("res_2024_amy2018_2_sb") == 2018 + assert eia_year_for_resstock_release("res_2024_AMY2012") == 2012 + with pytest.raises(ValueError, match="amyYYYY"): + eia_year_for_resstock_release("res_2024_tmy3") + + +def test_compare_resstock_eia176_gas(tmp_path: Path) -> None: + annual = pl.DataFrame( + { + "bldg_id": [1, 2, 3], + ANNUAL_GAS_COL: [1000.0, 2000.0, 500.0], + "weight": [10.0, 5.0, 2.0], + } + ) + path_annual = tmp_path / "annual.parquet" + annual.write_parquet(path_annual) + path_assignment = tmp_path / "utility_assignment.parquet" + pl.DataFrame( + { + "bldg_id": [1, 2, 3, 4], + "sb.gas_utility": ["yankee_gas", "yankee_gas", None, "norwich_muni"], + } + ).write_parquet(path_assignment) + + # 10 Mcf * 1.0 MMBtu/Mcf * KWH_PER_MMBTU, sales + transport. + by_consumer = pl.DataFrame( + { + "report_year": [2018, 2018, 2018, 2018, 2018], + "operator_id_eia": [ + "17602792CT", + "17602792CT", + "17610359CT", + "17619018CT", + "17619865CT", + ], + "operating_state": ["CT", "CT", "CT", "CT", "CT"], + "customer_class": ["residential"] * 5, + "revenue_class": ["sales", "transport", "sales", "sales", "sales"], + "consumers": [100, 4, 8, 1, 1], + "revenue": [1.0, 1.0, 1.0, 1.0, 1.0], + "volume_mcf": [8.0, 2.0, 1.0, 1.0, 1.0], + } + ) + path_consumer = tmp_path / "by_consumer.parquet" + by_consumer.write_parquet(path_consumer) + path_disposition = tmp_path / "disposition.parquet" + pl.DataFrame( + { + "operator_id_eia": [ + "17602792CT", + "17619018CT", + "17619865CT", + "17610359CT", + ], + "report_year": [2018, 2018, 2018, 2018], + "operating_state": ["CT", "CT", "CT", "CT"], + "delivered_gas_heat_content_mmbtu_per_mcf": [1.0, 1.0, 1.0, 2.0], + } + ).write_parquet(path_disposition) + + result = compare_resstock_eia176_gas( + path_annual=str(path_annual), + path_utility_assignment=str(path_assignment), + path_by_consumer=str(path_consumer), + path_disposition=str(path_disposition), + state="CT", + year=2018, + ) + + yankee = result.filter(pl.col("utility_code") == "yankee_gas").to_dicts()[0] + assert yankee[RESSTOCK_TOTAL_KWH] == pytest.approx(1000 * 10 + 2000 * 5) + assert yankee[EIA_RESIDENTIAL_KWH] == pytest.approx(10.0 * KWH_PER_MMBTU) + assert yankee[EIA_RESIDENTIAL_CUSTOMERS] == 104 + assert yankee[RESSTOCK_CUSTOMERS] == pytest.approx(15.0) + assert yankee["kwh_ratio"] == pytest.approx( + yankee[RESSTOCK_TOTAL_KWH] / yankee[EIA_RESIDENTIAL_KWH] + ) + + norwich = result.filter(pl.col("utility_code") == "norwich_muni") + assert norwich.height == 0 + + +def test_unknown_state_raises(tmp_path: Path) -> None: + path_annual = tmp_path / "annual.parquet" + pl.DataFrame( + {"bldg_id": [1], ANNUAL_GAS_COL: [1.0], "weight": [1.0]} + ).write_parquet(path_annual) + path_assignment = tmp_path / "utility_assignment.parquet" + pl.DataFrame({"bldg_id": [1], "sb.gas_utility": ["yankee_gas"]}).write_parquet( + path_assignment + ) + with pytest.raises(ValueError, match="No EIA-176"): + compare_resstock_eia176_gas( + path_annual=str(path_annual), + path_utility_assignment=str(path_assignment), + path_by_consumer=str(path_annual), + path_disposition=str(path_annual), + state="NY", + year=2018, + ) diff --git a/tests/test_puma_utility_housing_overlap.py b/tests/test_puma_utility_housing_overlap.py new file mode 100644 index 00000000..88694ab0 --- /dev/null +++ b/tests/test_puma_utility_housing_overlap.py @@ -0,0 +1,108 @@ +"""Tests for housing-unit-weighted PUMA-utility overlap.""" + +from __future__ import annotations + +import geopandas as gpd +import polars as pl +import pytest +from shapely.geometry import Point, box + +from data.resstock.utility.utils import ( + add_block_gas_heated_homes, + calculate_puma_utility_housing_overlap, + calculate_utility_probabilities, +) + +CRS = 2234 + + +def test_housing_units_set_probabilities_not_area() -> None: + pumas = gpd.GeoDataFrame( + {"PUMACE10": ["00100"]}, geometry=[box(0, 0, 100, 100)], crs=CRS + ) + # Yankee covers 90% of the land; the small city polygon holds most homes. + utilities = gpd.GeoDataFrame( + {"utility": ["yankee_gas", "norwich_muni"]}, + geometry=[box(10, 0, 100, 100), box(0, 0, 10, 100)], + crs=CRS, + ) + blocks = gpd.GeoDataFrame( + { + "block_geoid": ["a", "b", "c", "d"], + "housing_units": [300, 500, 200, 1000], + }, + geometry=[Point(5, 50), Point(6, 60), Point(50, 50), Point(200, 200)], + crs=CRS, + ) + + overlap = calculate_puma_utility_housing_overlap(pumas, utilities, blocks, CRS) + probs = calculate_utility_probabilities( + overlap, + pl.DataFrame( + { + "state_name": pl.Series([], dtype=pl.Utf8), + "std_name": pl.Series([], dtype=pl.Utf8), + } + ).lazy(), + handle_municipal=False, + ).collect() + + row = probs.to_dicts()[0] + assert row["norwich_muni"] == pytest.approx(0.8) + assert row["yankee_gas"] == pytest.approx(0.2) + + +def test_block_in_two_polygons_splits_its_housing_units() -> None: + pumas = gpd.GeoDataFrame( + {"PUMACE10": ["00100"]}, geometry=[box(0, 0, 100, 100)], crs=CRS + ) + utilities = gpd.GeoDataFrame( + {"utility": ["yankee_gas", "ct_natural_gas"]}, + geometry=[box(0, 0, 60, 100), box(40, 0, 100, 100)], + crs=CRS, + ) + blocks = gpd.GeoDataFrame( + {"block_geoid": ["shared"], "housing_units": [100]}, + geometry=[Point(50, 50)], + crs=CRS, + ) + + overlap = calculate_puma_utility_housing_overlap( + pumas, utilities, blocks, CRS + ).collect() + + assert overlap["pct_overlap"].to_list() == pytest.approx([50.0, 50.0]) + + +def test_gas_heated_homes_spread_by_housing_units_and_weight_overlap() -> None: + bg = "090010001001" + blocks = gpd.GeoDataFrame( + { + "block_geoid": [f"{bg}001", f"{bg}002", "090010001002001"], + "housing_units": [30, 10, 50], + }, + geometry=[Point(5, 50), Point(50, 50), Point(60, 60)], + crs=CRS, + ) + block_groups = pl.DataFrame( + {"block_group_geoid": [bg, "090010001002"], "gas_heated_homes": [20, 0]} + ) + + with_gas = add_block_gas_heated_homes(blocks, block_groups) + assert with_gas["gas_heated_homes"].tolist() == pytest.approx([15.0, 5.0, 0.0]) + + pumas = gpd.GeoDataFrame( + {"PUMACE10": ["00100"]}, geometry=[box(0, 0, 100, 100)], crs=CRS + ) + utilities = gpd.GeoDataFrame( + {"utility": ["norwich_muni", "yankee_gas"]}, + geometry=[box(0, 0, 10, 100), box(10, 0, 100, 100)], + crs=CRS, + ) + overlap = calculate_puma_utility_housing_overlap( + pumas, utilities, with_gas, CRS, weight_col="gas_heated_homes" + ).collect() + + shares = dict(zip(overlap["utility"], overlap["pct_overlap"], strict=True)) + assert shares["norwich_muni"] == pytest.approx(75.0) + assert shares["yankee_gas"] == pytest.approx(25.0) diff --git a/tests/test_resstock_warnings.py b/tests/test_resstock_warnings.py index 51e3a5cf..474b850b 100644 --- a/tests/test_resstock_warnings.py +++ b/tests/test_resstock_warnings.py @@ -33,3 +33,36 @@ def test_approx_does_not_warn_when_one_configured_upgrade_is_requested() -> None add_annual_loads=False, ) assert not any("approx_upgrade_ids" in w for w in warnings) + + +def test_adjust_gas_usage_warns_without_hourly_or_annual() -> None: + warnings = collect_run_warnings( + file_types=["metadata"], + upgrade_ids=["0"], + approximate_non_hp_load=False, + approx_upgrades=["01", "02"], + adjust_mf_electricity=False, + mf_adj_upgrades=["00", "02"], + assign_utility=False, + add_monthly_loads=False, + add_annual_loads=False, + adjust_gas_usage=True, + ) + gas_warnings = [w for w in warnings if "--adjust-gas-usage" in w] + assert len(gas_warnings) == 2 + + +def test_adjust_gas_usage_warns_when_upgrade_00_not_requested() -> None: + warnings = collect_run_warnings( + file_types=["load_curve_hourly", "load_curve_annual"], + upgrade_ids=["2"], + approximate_non_hp_load=False, + approx_upgrades=["01", "02"], + adjust_mf_electricity=False, + mf_adj_upgrades=["00", "02"], + assign_utility=False, + add_monthly_loads=False, + add_annual_loads=False, + adjust_gas_usage=True, + ) + assert any("upgrade-00 load_curve_annual" in w for w in warnings) diff --git a/utils/post/compare_resstock_eia176_gas.py b/utils/post/compare_resstock_eia176_gas.py new file mode 100644 index 00000000..6a84c94f --- /dev/null +++ b/utils/post/compare_resstock_eia176_gas.py @@ -0,0 +1,349 @@ +"""Compare ResStock natural gas use per gas utility to EIA-176 residential deliveries. + +ResStock ``res_2024_amy2018_2`` is simulated on Actual Meteorological Year 2018, +so the default EIA year is 2018 (parsed from the release name). Each building's +annual gas kWh (``out.natural_gas.total.energy_consumption.kwh``) is multiplied +by its sample weight and summed by ``sb.gas_utility``. + +EIA-176 (via the PUDL release pinned in ``data.eia.constants``) reports company- +and state-level deliveries. The comparison uses residential sales plus +residential transportation: ResStock is a residential building sample, and in a +retail-choice state the transportation lines are gas the utility delivers for a +marketer. Volumes are Mcf, converted to kWh with that company's delivered-gas +heat content. + +Connecticut operator ids are the EIA-176 company ids from the 2018 company list: +Yankee Gas, Connecticut Natural Gas, Southern Connecticut Gas, and Norwich +Public Utilities. +""" + +from __future__ import annotations + +import argparse +import re +import sys +from pathlib import Path + +import polars as pl + +from data.eia.constants import PUDL_STABLE_VERSION +from utils import get_aws_region + +BLDG_ID_COL = "bldg_id" +GAS_UTILITY_COL = "sb.gas_utility" +WEIGHT_COL = "weight" +ANNUAL_GAS_COL = "out.natural_gas.total.energy_consumption.kwh" +RESSTOCK_TOTAL_KWH = "resstock_total_kwh" +RESSTOCK_CUSTOMERS = "resstock_customers" +EIA_RESIDENTIAL_KWH = "eia_residential_kwh" +EIA_RESIDENTIAL_CUSTOMERS = "eia_residential_customers" + +# 1 MMBtu = 1e6 Btu; 1 kWh = 3412.141633 Btu. +KWH_PER_MMBTU = 1_000_000 / 3412.141633 + +DEFAULT_RESSTOCK_RELEASE = "res_2024_amy2018_2_sb" +S3_BASE_RESSTOCK = "s3://data.sb/nrel/resstock" +_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" + +# EIA-176 operator_id_eia -> platform gas utility, by state. +# Connecticut names are from eia176_2018_company_list.csv. +GAS_OPERATORS_BY_STATE: dict[str, dict[str, str]] = { + "CT": { + "17602792CT": "yankee_gas", # YANKEE GAS SVC CO + "17619018CT": "ct_natural_gas", # CONNECTICUT NAT GAS CORP + "17619865CT": "southern_ct_gas", # SOUTHERN CONNECTICUT GAS COMPANY + "17610359CT": "norwich_muni", # NORWICH PUBLIC UTILITIES + }, +} + +_AMY_YEAR = re.compile(r"amy(\d{4})", re.IGNORECASE) + + +def eia_year_for_resstock_release(release: str) -> int: + """Return the AMY weather year embedded in a ResStock release name.""" + match = _AMY_YEAR.search(release) + if match is None: + raise ValueError( + f"Cannot infer the EIA year from ResStock release {release!r}. " + "Expected a name containing amyYYYY, such as res_2024_amy2018_2_sb. " + "Pass --eia-year." + ) + return int(match.group(1)) + + +def _storage_options() -> dict[str, str]: + return {"aws_region": get_aws_region()} + + +def _scan(path: str, storage_options: dict[str, str] | None) -> pl.LazyFrame: + opts = storage_options if path.startswith("s3://") else None + return pl.scan_parquet(path, storage_options=opts) + + +def _default_annual_path(release: str, state: str, upgrade: str) -> str: + base = ( + f"{S3_BASE_RESSTOCK}/{release}/load_curve_annual/" + f"state={state}/upgrade={upgrade}" + ) + return f"{base}/{state}_upgrade{upgrade}_metadata_and_annual_results.parquet" + + +def _default_utility_assignment_path(release: str, state: str) -> str: + return ( + f"{S3_BASE_RESSTOCK}/{release}/metadata_utility/" + f"state={state}/utility_assignment.parquet" + ) + + +def load_resstock_gas_kwh_per_utility( + path_annual: str, + path_utility_assignment: str, + storage_options: dict[str, str] | None, + gas_column: str = ANNUAL_GAS_COL, +) -> pl.DataFrame: + """Sum weighted annual gas kWh and sample weights by gas utility.""" + annual_schema = _scan(path_annual, storage_options).collect_schema().names() + for column in (gas_column, BLDG_ID_COL, WEIGHT_COL): + if column not in annual_schema: + raise ValueError( + f"Annual parquet at {path_annual!r} is missing column {column!r}." + ) + annual = ( + _scan(path_annual, storage_options) + .select( + pl.col(BLDG_ID_COL), + pl.col(gas_column).alias("annual_kwh"), + pl.col(WEIGHT_COL), + ) + .collect() + ) + assignment_schema = ( + _scan(path_utility_assignment, storage_options).collect_schema().names() + ) + if GAS_UTILITY_COL not in assignment_schema: + raise ValueError( + f"Utility assignment at {path_utility_assignment!r} is missing " + f"{GAS_UTILITY_COL!r}." + ) + assignment = ( + _scan(path_utility_assignment, storage_options) + .select(BLDG_ID_COL, GAS_UTILITY_COL) + .collect() + ) + joined = annual.join(assignment, on=BLDG_ID_COL, how="inner").filter( + pl.col(GAS_UTILITY_COL).is_not_null() + ) + return joined.group_by(GAS_UTILITY_COL).agg( + (pl.col("annual_kwh").fill_null(0) * pl.col(WEIGHT_COL)) + .sum() + .alias(RESSTOCK_TOTAL_KWH), + pl.col(WEIGHT_COL).sum().alias(RESSTOCK_CUSTOMERS), + ) + + +def load_eia176_residential( + path_by_consumer: str, + path_disposition: str, + state: str, + year: int, + operator_to_utility: dict[str, str], +) -> pl.DataFrame: + """Residential deliveries (sales + transport) in kWh and customer counts.""" + operator_ids = list(operator_to_utility) + deliveries = ( + _scan(path_by_consumer, None) + .filter( + (pl.col("operating_state") == state) + & (pl.col("report_year") == year) + & (pl.col("customer_class") == "residential") + & pl.col("operator_id_eia").is_in(operator_ids) + ) + .group_by("operator_id_eia") + .agg( + pl.col("volume_mcf").sum().alias("volume_mcf"), + pl.col("consumers").sum().alias(EIA_RESIDENTIAL_CUSTOMERS), + ) + .collect() + ) + if not isinstance(deliveries, pl.DataFrame): + raise TypeError("Expected DataFrame from EIA-176 deliveries collect()") + missing = sorted(set(operator_ids) - set(deliveries["operator_id_eia"].to_list())) + if missing: + raise ValueError( + f"EIA-176 has no residential deliveries for {state} {year} " + f"operator(s) {missing}." + ) + heat = ( + _scan(path_disposition, None) + .filter( + (pl.col("operating_state") == state) + & (pl.col("report_year") == year) + & pl.col("operator_id_eia").is_in(operator_ids) + ) + .select("operator_id_eia", "delivered_gas_heat_content_mmbtu_per_mcf") + .collect() + ) + if not isinstance(heat, pl.DataFrame): + raise TypeError("Expected DataFrame from EIA-176 heat content collect()") + mapping = pl.DataFrame( + { + "operator_id_eia": operator_ids, + "utility_code": [operator_to_utility[i] for i in operator_ids], + } + ) + out = deliveries.join(heat, on="operator_id_eia", how="left").join( + mapping, on="operator_id_eia", how="left" + ) + if out["delivered_gas_heat_content_mmbtu_per_mcf"].null_count() > 0: + raise ValueError( + f"EIA-176 is missing delivered-gas heat content for {state} {year}." + ) + return out.select( + "utility_code", + ( + pl.col("volume_mcf") + * pl.col("delivered_gas_heat_content_mmbtu_per_mcf") + * KWH_PER_MMBTU + ).alias(EIA_RESIDENTIAL_KWH), + pl.col(EIA_RESIDENTIAL_CUSTOMERS), + ) + + +def compare_resstock_eia176_gas( + path_annual: str, + path_utility_assignment: str, + path_by_consumer: str, + path_disposition: str, + state: str, + year: int, + storage_options: dict[str, str] | None = None, + gas_column: str = ANNUAL_GAS_COL, +) -> pl.DataFrame: + """Build the per-utility ResStock vs EIA-176 comparison table.""" + operators = GAS_OPERATORS_BY_STATE.get(state) + if operators is None: + supported = ", ".join(sorted(GAS_OPERATORS_BY_STATE)) + raise ValueError( + f"No EIA-176 gas-utility map for state {state}. Supported: {supported}." + ) + resstock = load_resstock_gas_kwh_per_utility( + path_annual, + path_utility_assignment, + storage_options, + gas_column=gas_column, + ).rename({GAS_UTILITY_COL: "utility_code"}) + eia = load_eia176_residential( + path_by_consumer, path_disposition, state, year, operators + ) + return ( + resstock.join(eia, on="utility_code", how="left") + .with_columns( + (pl.col(RESSTOCK_TOTAL_KWH) / pl.col(EIA_RESIDENTIAL_KWH)).alias( + "kwh_ratio" + ), + ( + (pl.col(RESSTOCK_TOTAL_KWH) - pl.col(EIA_RESIDENTIAL_KWH)) + / pl.col(EIA_RESIDENTIAL_KWH) + * 100 + ).alias("kwh_pct_diff"), + (pl.col(RESSTOCK_CUSTOMERS) / pl.col(EIA_RESIDENTIAL_CUSTOMERS)).alias( + "customers_ratio" + ), + ( + (pl.col(RESSTOCK_CUSTOMERS) - pl.col(EIA_RESIDENTIAL_CUSTOMERS)) + / pl.col(EIA_RESIDENTIAL_CUSTOMERS) + * 100 + ).alias("customers_pct_diff"), + ) + .sort("utility_code") + ) + + +def main() -> None: + parser = argparse.ArgumentParser( + description=( + "Compare weighted ResStock natural gas use by gas utility to " + "EIA-176 residential deliveries for the ResStock weather year." + ) + ) + parser.add_argument("--state", required=True, help="State abbreviation, e.g. CT.") + parser.add_argument( + "--resstock-release", + default=DEFAULT_RESSTOCK_RELEASE, + help="ResStock release (default: res_2024_amy2018_2_sb).", + ) + parser.add_argument("--upgrade", default="00", help="Zero-padded upgrade id.") + parser.add_argument("--path-annual", default=None) + parser.add_argument("--path-utility-assignment", default=None) + parser.add_argument( + "--eia-year", + type=int, + default=None, + help="EIA-176 report year. Default: AMY year parsed from the release name.", + ) + parser.add_argument("--path-eia176-by-consumer", default=EIA176_BY_CONSUMER) + parser.add_argument("--path-eia176-disposition", default=EIA176_DISPOSITION) + parser.add_argument( + "--output", + type=Path, + default=None, + help="Write comparison CSV here. Prints CSV to stdout when omitted.", + ) + parser.add_argument( + "--load-column", + default=ANNUAL_GAS_COL, + help=f"Annual gas kWh column (default: {ANNUAL_GAS_COL}).", + ) + args = parser.parse_args() + + state = args.state.strip().upper() + year = ( + args.eia_year + if args.eia_year is not None + else eia_year_for_resstock_release(args.resstock_release) + ) + path_annual = args.path_annual or _default_annual_path( + args.resstock_release, state, args.upgrade + ) + path_utility_assignment = ( + args.path_utility_assignment + or _default_utility_assignment_path(args.resstock_release, state) + ) + needs_aws = path_annual.startswith("s3://") or path_utility_assignment.startswith( + "s3://" + ) + comparison = compare_resstock_eia176_gas( + path_annual=path_annual, + path_utility_assignment=path_utility_assignment, + path_by_consumer=args.path_eia176_by_consumer, + path_disposition=args.path_eia176_disposition, + state=state, + year=year, + storage_options=_storage_options() if needs_aws else None, + gas_column=args.load_column, + ) + print( + f"ResStock {args.resstock_release} upgrade {args.upgrade} vs EIA-176 {year} " + f"residential deliveries (sales + transport) in {state}.", + file=sys.stderr, + ) + if args.output is not None: + args.output.parent.mkdir(parents=True, exist_ok=True) + comparison.write_csv(args.output) + print(f"Wrote {args.output}", file=sys.stderr) + else: + print(comparison.write_csv(), end="") + + +if __name__ == "__main__": + try: + main() + except (FileNotFoundError, ValueError) as exc: + print(f"Error: {exc}", file=sys.stderr) + sys.exit(1)