CHEERSafe Loss Estimation Workflow¶

This notebook runs the cheersafe.loss_estimation chain: the damage library becomes wind-only and flood-only expected loss-ratio lookup tables, the exposure export's continuous hazard IDs select table bins under a chosen index convention, wind and flood are combined, and the annual event probabilities turn the per-event loss ratios into expected annual loss (EAL) per building and per zone.

Every calculation is called explicitly and its return value assigned by the notebook. See the loss-estimation manual for the methodology.

!!! warning "Execution profiles" Smoke mode (default) builds a synthetic exposure case and a synthetic damage library for the (M, C) pairs that occur in it (about 20 seconds, no external data).

**Full mode** (`CHEERSAFE_LOSS_SMOKE=0`) reads the 997 MB Eastern North Carolina damage library, the
frozen 07/07/2026 export (1,154,657 buildings x 97 events, memory-mapped) and the LEAPHI zones from the
`Loss Estimates` folder (`CHEERSAFE_LOSS_DATA`). Allow about 2 GB of RAM and a few minutes. Missing
files are reported and the dependent cells are skipped.

Plan¶

Choose the execution profile, the index convention, the combination mode and the value basis.
Load (or generate) the exposure export, the damage library and the event probabilities.
Query the damage library through DamageProbabilityProvider.
Build the lookup tables with LookupTableBuilder (full mode checks them against the frozen tables).
Compare the four index conventions.
Compute building EAL for each value basis.
Aggregate EAL to zones and compare with the orchestrator result.

In [1]:

Copied!





# Repository setup and imports
import json
import os
import pickle
import sys
import time
from pathlib import Path

starting_directory = Path.cwd().resolve()
repo_root = starting_directory
while not (repo_root / "pyproject.toml").exists() and repo_root != repo_root.parent:
    repo_root = repo_root.parent
os.chdir(repo_root)
sys.path.insert(0, str(repo_root / "src"))

import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from cheersafe.config import Config
from cheersafe.exposure import Exposure, SyntheticExposureData
from cheersafe.loss_estimation import DamageProbabilityProvider, LossEstimation
from cheersafe.tools import Utils

print(f"repo_root={repo_root}")
# Repository setup and imports
import json
import os
import pickle
import sys
import time
from pathlib import Path

starting_directory = Path.cwd().resolve()
repo_root = starting_directory
while not (repo_root / "pyproject.toml").exists() and repo_root != repo_root.parent:
    repo_root = repo_root.parent
os.chdir(repo_root)
sys.path.insert(0, str(repo_root / "src"))

import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from cheersafe.config import Config
from cheersafe.exposure import Exposure, SyntheticExposureData
from cheersafe.loss_estimation import DamageProbabilityProvider, LossEstimation
from cheersafe.tools import Utils

print(f"repo_root={repo_root}")

/home/mammad/.local/lib/python3.12/site-packages/matplotlib/projections/__init__.py:63: UserWarning: Unable to import Axes3D. This may be due to multiple versions of Matplotlib being installed (e.g. as a system package and as a pip package). As a result, the 3D projection is not available.
  warnings.warn("Unable to import Axes3D. This may be due to multiple versions of "

repo_root=/home/mammad/Documents/EDU/cope/CHEERSafe/LossModel-exposure

In [2]:

parameters

Copied!





# Execution profile: smoke mode on a synthetic case by default; full mode reads the external Eastern North Carolina data
FULL_DATA_ROOT = Path(os.environ.get("CHEERSAFE_LOSS_DATA", "../Loss Estimates"))
FULL_LIBRARY = FULL_DATA_ROOT / "CHEER_Safe_data" / "ProbDS_cwfd_MaxIter3000_V2.pkl"
FULL_EXPORT_DIR = FULL_DATA_ROOT / "User_Output_07072026" / "Old97"
FULL_EXPORT_SUFFIX = "_valid_zone_order_07072026"
FULL_LOOKUP_DIR = FULL_DATA_ROOT / "User_Output_07072026" / "STARR_lookup_tables"
FULL_ZONES = FULL_DATA_ROOT / "Zones" / "leaphi_zones_20230920.shp"
PROB_JSON = Path("Inputs/Hurricane_ProbWindFlood.json")

SMOKE_N_BUILDINGS = 2000
SMOKE_N_EVENTS = 5
SMOKE_SEED = 42

INDEX_MODE = "floor"              # floor (published figures) | ceil (STARR export) | round | linear
COMBINE_MODE = "additive"         # additive (wind + flood) | joint (full wind-flood surface)
VALUE_BASIS = "structure_value"   # structure_value | contents_value | home_value
ZONE_ID = "i_lh740"
TARGET_ZONE_ID = "i_lh503"
AVERAGE_MODE = "Nearest"

smoke_flag = os.environ.get("CHEERSAFE_LOSS_SMOKE", "1")
if smoke_flag not in {"0", "1"}:
    raise ValueError("CHEERSAFE_LOSS_SMOKE must be '0' or '1'")
smoke_mode = smoke_flag == "1"
output_override = os.environ.get("CHEERSAFE_LOSS_OUTPUT")
output_path = (
    Path(output_override).expanduser().resolve()
    if output_override
    else (repo_root / "Outputs/loss_estimation_workflow").resolve()
)
full_data_root = (repo_root / FULL_DATA_ROOT).resolve() if not FULL_DATA_ROOT.is_absolute() else FULL_DATA_ROOT
print(f"profile={'smoke' if smoke_mode else 'full'}, index_mode={INDEX_MODE}, combine_mode={COMBINE_MODE}, "
      f"value_basis={VALUE_BASIS}, output={output_path}")
# Execution profile: smoke mode on a synthetic case by default; full mode reads the external Eastern North Carolina data
FULL_DATA_ROOT = Path(os.environ.get("CHEERSAFE_LOSS_DATA", "../Loss Estimates"))
FULL_LIBRARY = FULL_DATA_ROOT / "CHEER_Safe_data" / "ProbDS_cwfd_MaxIter3000_V2.pkl"
FULL_EXPORT_DIR = FULL_DATA_ROOT / "User_Output_07072026" / "Old97"
FULL_EXPORT_SUFFIX = "_valid_zone_order_07072026"
FULL_LOOKUP_DIR = FULL_DATA_ROOT / "User_Output_07072026" / "STARR_lookup_tables"
FULL_ZONES = FULL_DATA_ROOT / "Zones" / "leaphi_zones_20230920.shp"
PROB_JSON = Path("Inputs/Hurricane_ProbWindFlood.json")

SMOKE_N_BUILDINGS = 2000
SMOKE_N_EVENTS = 5
SMOKE_SEED = 42

INDEX_MODE = "floor"              # floor (published figures) | ceil (STARR export) | round | linear
COMBINE_MODE = "additive"         # additive (wind + flood) | joint (full wind-flood surface)
VALUE_BASIS = "structure_value"   # structure_value | contents_value | home_value
ZONE_ID = "i_lh740"
TARGET_ZONE_ID = "i_lh503"
AVERAGE_MODE = "Nearest"

smoke_flag = os.environ.get("CHEERSAFE_LOSS_SMOKE", "1")
if smoke_flag not in {"0", "1"}:
    raise ValueError("CHEERSAFE_LOSS_SMOKE must be '0' or '1'")
smoke_mode = smoke_flag == "1"
output_override = os.environ.get("CHEERSAFE_LOSS_OUTPUT")
output_path = (
    Path(output_override).expanduser().resolve()
    if output_override
    else (repo_root / "Outputs/loss_estimation_workflow").resolve()
)
full_data_root = (repo_root / FULL_DATA_ROOT).resolve() if not FULL_DATA_ROOT.is_absolute() else FULL_DATA_ROOT
print(f"profile={'smoke' if smoke_mode else 'full'}, index_mode={INDEX_MODE}, combine_mode={COMBINE_MODE}, "
      f"value_basis={VALUE_BASIS}, output={output_path}")

profile=smoke, index_mode=floor, combine_mode=additive, value_basis=structure_value, output=/home/mammad/Documents/EDU/cope/CHEERSafe/LossModel-exposure/Outputs/loss_estimation_workflow

In [3]:

Copied!





# Snapshot the global configuration touched below, then switch to manual control
config_snapshot = {
    name: getattr(Config, name)
    for name in ("AUTO_LOAD_DATASETS", "AUTO_EXECUTE_METHODS", "INPUTS_PATH", "OUTPUTS_PATH", "RANDOM_SEED")
}


def restore_state() -> None:
    """Restore the configuration snapshot and the starting directory."""
    for name, value in config_snapshot.items():
        setattr(Config, name, value)
    os.chdir(starting_directory)


Config.AUTO_LOAD_DATASETS = False      # every load_* call below is explicit
Config.AUTO_EXECUTE_METHODS = False    # every calculate_* result is assigned by this notebook
Config.INPUTS_PATH = "Inputs"          # the package's JSON tables
output_path.mkdir(parents=True, exist_ok=True)
print(config_snapshot)
# Snapshot the global configuration touched below, then switch to manual control
config_snapshot = {
    name: getattr(Config, name)
    for name in ("AUTO_LOAD_DATASETS", "AUTO_EXECUTE_METHODS", "INPUTS_PATH", "OUTPUTS_PATH", "RANDOM_SEED")
}


def restore_state() -> None:
    """Restore the configuration snapshot and the starting directory."""
    for name, value in config_snapshot.items():
        setattr(Config, name, value)
    os.chdir(starting_directory)


Config.AUTO_LOAD_DATASETS = False      # every load_* call below is explicit
Config.AUTO_EXECUTE_METHODS = False    # every calculate_* result is assigned by this notebook
Config.INPUTS_PATH = "Inputs"          # the package's JSON tables
output_path.mkdir(parents=True, exist_ok=True)
print(config_snapshot)

{'AUTO_LOAD_DATASETS': True, 'AUTO_EXECUTE_METHODS': True, 'INPUTS_PATH': 'Inputs', 'OUTPUTS_PATH': 'Outputs', 'RANDOM_SEED': 42}

Inputs¶

Smoke mode runs the exposure chain on a synthetic case and generates a damage library for the (M, C) pairs that occur, so no lookup row is NaN. Full mode uses the frozen Eastern North Carolina export, the real damage library and the 97-event probabilities bundled with the package.

In [4]:

Copied!





canonical_tables = {
    "wind_speeds": Utils.from_json("Inputs/WindSpeed.json"),
    "flood_depths": Utils.from_json("Inputs/FloodDepth.json"),
    "damage_states": Utils.from_json("Inputs/DamageState.json"),
    "component_combinations": Utils.from_json("Inputs/ComponentCombination_All.json"),
}

def run_exposure(inventory_path, hazard_dir, zones_path, seed):
    """Run the exposure chain with explicit calls; every result is assigned by the caller."""
    exposure = Exposure(
        inventory_path=str(inventory_path),
        hazard_dir=str(hazard_dir),
        zones_path=str(zones_path) if zones_path else None,
        average_mode=AVERAGE_MODE,
        zone_id_col=ZONE_ID,
        seed=seed,
    )
    exposure.load_hazard_tables()
    exposure.load_component_combinations()
    exposure.load_inventory()
    exposure.load_hazard_events()
    exposure.load_zones()
    exposure.archetype_inventory = exposure.calculate_archetypes(exposure.inventory)
    exposure.clipped_inventory = exposure.calculate_clipped_inventory(exposure.archetype_inventory)
    centroids = exposure.clipped_inventory.geometry.centroid
    exposure.neighbor_index, exposure.neighbor_distance = exposure.calculate_neighbor_match(
        np.column_stack([centroids.x, centroids.y]), exposure.hazard_points_xy, 1 if AVERAGE_MODE == "Nearest" else 4
    )
    exposure.wind_ids, exposure.flood_ids = exposure.calculate_hazard_ids(
        exposure.neighbor_index, exposure.neighbor_distance, exposure.hazard_wind, exposure.hazard_flood
    )
    exposure.zone_ids = exposure.calculate_zone_assignment(exposure.clipped_inventory, exposure.zones, ZONE_ID)
    exposure.starr_inventory = exposure.calculate_starr_inventory(exposure.clipped_inventory, exposure.zone_ids)
    return exposure

data_available = True
exposure = None
if smoke_mode:
    case_root = output_path / "synthetic_case"
    probabilities = SyntheticExposureData.make_hurricane_probabilities(SMOKE_N_EVENTS, seed=SMOKE_SEED)
    written = SyntheticExposureData.save_inputs(
        str(case_root),
        inventory=SyntheticExposureData.make_inventory(SMOKE_N_BUILDINGS, seed=SMOKE_SEED, include_edge_cases=False),
        hazard_events=SyntheticExposureData.make_hazard_events(
            tuple(f"h{index}" for index in range(1, SMOKE_N_EVENTS + 1)), nx=4, ny=4, seed=SMOKE_SEED
        ),
        zones=SyntheticExposureData.make_zones(extent_cells=(4, 4)),
        hurricane_probabilities=probabilities,
        **canonical_tables,
    )
    exposure = run_exposure(written["inventory"], written["hazard_dir"], written["zones"], SMOKE_SEED)
    pairs = sorted(set(zip(exposure.starr_inventory["M"].tolist(), exposure.starr_inventory["C"].tolist())))
    library_path = case_root / "Inputs" / "damage_library.pkl"
    with open(library_path, "wb") as handle:
        pickle.dump(SyntheticExposureData.make_damage_library(pairs, seed=SMOKE_SEED, unnormalized_rows=0), handle)
    zones_path = Path(written["zones"])
    export_paths = None
    print(f"synthetic case: {len(exposure.starr_inventory)} buildings, {len(exposure.event_names)} events, {len(pairs)} (M, C) pairs")
else:
    library_path = full_data_root / FULL_LIBRARY.relative_to(FULL_DATA_ROOT)
    export_dir = full_data_root / FULL_EXPORT_DIR.relative_to(FULL_DATA_ROOT)
    export_paths = {
        "inventory_csv": export_dir / f"Inventory{FULL_EXPORT_SUFFIX}.csv",
        "wind_ids": export_dir / f"Wind_Tr{FULL_EXPORT_SUFFIX}.npy",
        "flood_ids": export_dir / f"Flood_Tr{FULL_EXPORT_SUFFIX}.npy",
    }
    zones_path = full_data_root / FULL_ZONES.relative_to(FULL_DATA_ROOT)
    probabilities = np.asarray(Utils.from_json(str(PROB_JSON), dict_key="Hurricane_Prob"))
    missing = [str(p) for p in [library_path, zones_path, *export_paths.values()] if not p.exists()]
    if missing:
        data_available = False
        print("skipped: missing", missing)
print(f"library={library_path}\nzones={zones_path}\nevents={len(probabilities)} (sum p = {probabilities[:, 1].sum():.4f}/yr)\navailable={data_available}")
canonical_tables = {
    "wind_speeds": Utils.from_json("Inputs/WindSpeed.json"),
    "flood_depths": Utils.from_json("Inputs/FloodDepth.json"),
    "damage_states": Utils.from_json("Inputs/DamageState.json"),
    "component_combinations": Utils.from_json("Inputs/ComponentCombination_All.json"),
}

def run_exposure(inventory_path, hazard_dir, zones_path, seed):
    """Run the exposure chain with explicit calls; every result is assigned by the caller."""
    exposure = Exposure(
        inventory_path=str(inventory_path),
        hazard_dir=str(hazard_dir),
        zones_path=str(zones_path) if zones_path else None,
        average_mode=AVERAGE_MODE,
        zone_id_col=ZONE_ID,
        seed=seed,
    )
    exposure.load_hazard_tables()
    exposure.load_component_combinations()
    exposure.load_inventory()
    exposure.load_hazard_events()
    exposure.load_zones()
    exposure.archetype_inventory = exposure.calculate_archetypes(exposure.inventory)
    exposure.clipped_inventory = exposure.calculate_clipped_inventory(exposure.archetype_inventory)
    centroids = exposure.clipped_inventory.geometry.centroid
    exposure.neighbor_index, exposure.neighbor_distance = exposure.calculate_neighbor_match(
        np.column_stack([centroids.x, centroids.y]), exposure.hazard_points_xy, 1 if AVERAGE_MODE == "Nearest" else 4
    )
    exposure.wind_ids, exposure.flood_ids = exposure.calculate_hazard_ids(
        exposure.neighbor_index, exposure.neighbor_distance, exposure.hazard_wind, exposure.hazard_flood
    )
    exposure.zone_ids = exposure.calculate_zone_assignment(exposure.clipped_inventory, exposure.zones, ZONE_ID)
    exposure.starr_inventory = exposure.calculate_starr_inventory(exposure.clipped_inventory, exposure.zone_ids)
    return exposure

data_available = True
exposure = None
if smoke_mode:
    case_root = output_path / "synthetic_case"
    probabilities = SyntheticExposureData.make_hurricane_probabilities(SMOKE_N_EVENTS, seed=SMOKE_SEED)
    written = SyntheticExposureData.save_inputs(
        str(case_root),
        inventory=SyntheticExposureData.make_inventory(SMOKE_N_BUILDINGS, seed=SMOKE_SEED, include_edge_cases=False),
        hazard_events=SyntheticExposureData.make_hazard_events(
            tuple(f"h{index}" for index in range(1, SMOKE_N_EVENTS + 1)), nx=4, ny=4, seed=SMOKE_SEED
        ),
        zones=SyntheticExposureData.make_zones(extent_cells=(4, 4)),
        hurricane_probabilities=probabilities,
        **canonical_tables,
    )
    exposure = run_exposure(written["inventory"], written["hazard_dir"], written["zones"], SMOKE_SEED)
    pairs = sorted(set(zip(exposure.starr_inventory["M"].tolist(), exposure.starr_inventory["C"].tolist())))
    library_path = case_root / "Inputs" / "damage_library.pkl"
    with open(library_path, "wb") as handle:
        pickle.dump(SyntheticExposureData.make_damage_library(pairs, seed=SMOKE_SEED, unnormalized_rows=0), handle)
    zones_path = Path(written["zones"])
    export_paths = None
    print(f"synthetic case: {len(exposure.starr_inventory)} buildings, {len(exposure.event_names)} events, {len(pairs)} (M, C) pairs")
else:
    library_path = full_data_root / FULL_LIBRARY.relative_to(FULL_DATA_ROOT)
    export_dir = full_data_root / FULL_EXPORT_DIR.relative_to(FULL_DATA_ROOT)
    export_paths = {
        "inventory_csv": export_dir / f"Inventory{FULL_EXPORT_SUFFIX}.csv",
        "wind_ids": export_dir / f"Wind_Tr{FULL_EXPORT_SUFFIX}.npy",
        "flood_ids": export_dir / f"Flood_Tr{FULL_EXPORT_SUFFIX}.npy",
    }
    zones_path = full_data_root / FULL_ZONES.relative_to(FULL_DATA_ROOT)
    probabilities = np.asarray(Utils.from_json(str(PROB_JSON), dict_key="Hurricane_Prob"))
    missing = [str(p) for p in [library_path, zones_path, *export_paths.values()] if not p.exists()]
    if missing:
        data_available = False
        print("skipped: missing", missing)
print(f"library={library_path}\nzones={zones_path}\nevents={len(probabilities)} (sum p = {probabilities[:, 1].sum():.4f}/yr)\navailable={data_available}")

synthetic case: 2000 buildings, 5 events, 370 (M, C) pairs
library=/home/mammad/Documents/EDU/cope/CHEERSafe/LossModel-exposure/Outputs/loss_estimation_workflow/synthetic_case/Inputs/damage_library.pkl
zones=/home/mammad/Documents/EDU/cope/CHEERSafe/LossModel-exposure/Outputs/loss_estimation_workflow/synthetic_case/Inputs/zones.gpkg
events=5 (sum p = 0.3200/yr)
available=True

Damage library queries¶

DamageProbabilityProvider.calculate_damage_ratio returns the expected damage ratio for one (M, C) at a wind speed (m/s) and flood depth (m), selecting the first bin at or above each value.

In [5]:

Copied!





if data_available:
    provider = DamageProbabilityProvider(
        damage_state_path="Inputs/DamageState.json",
        wind_speed_path="Inputs/WindSpeed.json",
        flood_depth_path="Inputs/FloodDepth.json",
        prob_ds_file_path=str(library_path),
    )
    sample_m, sample_c = sorted(provider.prob_damage_states)[0]
    scenarios = pd.DataFrame({"wind_mps": [0.0, 22.4, 42.5, 62.6, 96.1], "flood_m": [0.0, 0.0, 0.3, 1.0, 3.0]})
    scenarios["damage_ratio"] = [
        provider.calculate_damage_ratio(row.wind_mps, row.flood_m, m=sample_m, c=sample_c)
        for row in scenarios.itertuples(index=False)
    ]
    print(f"{len(provider.prob_damage_states)} (M, C) pairs in the library; showing M={sample_m}, C={sample_c}")
    display(scenarios)
    fig, ax = plt.subplots(figsize=(6, 4))
    wind_axis = provider.wind_speeds[:, 2]
    ax.plot(wind_axis, [provider.calculate_damage_ratio(w, 0.0, sample_m, sample_c) for w in wind_axis], marker="o")
    ax.set_xlabel("Wind speed (m/s)"); ax.set_ylabel("Expected damage ratio (no flood)")
    ax.set_title(f"Provider lookup, M={sample_m}, C={sample_c}"); ax.grid(True, alpha=0.3)
    plt.show()
if data_available:
    provider = DamageProbabilityProvider(
        damage_state_path="Inputs/DamageState.json",
        wind_speed_path="Inputs/WindSpeed.json",
        flood_depth_path="Inputs/FloodDepth.json",
        prob_ds_file_path=str(library_path),
    )
    sample_m, sample_c = sorted(provider.prob_damage_states)[0]
    scenarios = pd.DataFrame({"wind_mps": [0.0, 22.4, 42.5, 62.6, 96.1], "flood_m": [0.0, 0.0, 0.3, 1.0, 3.0]})
    scenarios["damage_ratio"] = [
        provider.calculate_damage_ratio(row.wind_mps, row.flood_m, m=sample_m, c=sample_c)
        for row in scenarios.itertuples(index=False)
    ]
    print(f"{len(provider.prob_damage_states)} (M, C) pairs in the library; showing M={sample_m}, C={sample_c}")
    display(scenarios)
    fig, ax = plt.subplots(figsize=(6, 4))
    wind_axis = provider.wind_speeds[:, 2]
    ax.plot(wind_axis, [provider.calculate_damage_ratio(w, 0.0, sample_m, sample_c) for w in wind_axis], marker="o")
    ax.set_xlabel("Wind speed (m/s)"); ax.set_ylabel("Expected damage ratio (no flood)")
    ax.set_title(f"Provider lookup, M={sample_m}, C={sample_c}"); ax.grid(True, alpha=0.3)
    plt.show()

370 (M, C) pairs in the library; showing M=1, C=1

	wind_mps	flood_m	damage_ratio
0	0.0	0.0	0.000000
1	22.4	0.0	0.007125
2	42.5	0.3	0.046568
3	62.6	1.0	0.477198
4	96.1	3.0	1.332117

No description has been provided for this image

Lookup tables¶

L_mcw_wd[m-1, c-1, w-1] is the wind-only expected loss ratio (flood at the zero-hazard ID) and L_mcf_fld[m-1, c-1, f-1] the flood-only one. Pairs absent from the library are NaN. In full mode the tables must equal the frozen STARR_lookup_tables byte for byte.

In [6]:

Copied!





if data_available:
    estimator = LossEstimation(
        prob_ds_file_path=str(library_path),
        hurricane_probabilities=probabilities,
        zones_path=str(zones_path),
        zone_id_col=ZONE_ID,
        aggregation_col=TARGET_ZONE_ID,
        index_mode=INDEX_MODE,
        combine_mode=COMBINE_MODE,
        value_basis=VALUE_BASIS,
        output_dir=str(output_path / "loss"),
    )
    estimator.load_damage_library()
    estimator.load_hurricane_probabilities()
    estimator.load_zones()
    estimator.wind_lookup, estimator.flood_lookup = estimator.calculate_marginal_lookups(
        estimator.prob_damage_states, estimator.damage_states[:, 1],
        estimator.n_building_types, estimator.n_component_configurations,
    )
    print("wind_lookup", estimator.wind_lookup.shape, "flood_lookup", estimator.flood_lookup.shape,
          f"NaN pairs: {int(np.isnan(estimator.wind_lookup[:, :, 0]).sum())}")
    if not smoke_mode:
        lookup_dir = full_data_root / FULL_LOOKUP_DIR.relative_to(FULL_DATA_ROOT)
        frozen_wind = np.load(lookup_dir / "L_mcw_wd.npy")
        frozen_flood = np.load(lookup_dir / "L_mcf_fld.npy")
        print("byte-identical to the frozen tables:",
              np.array_equal(estimator.wind_lookup, frozen_wind) and np.array_equal(estimator.flood_lookup, frozen_flood))
    lookup_paths = estimator.save_lookup_tables(str(output_path / "loss"))
if data_available:
    estimator = LossEstimation(
        prob_ds_file_path=str(library_path),
        hurricane_probabilities=probabilities,
        zones_path=str(zones_path),
        zone_id_col=ZONE_ID,
        aggregation_col=TARGET_ZONE_ID,
        index_mode=INDEX_MODE,
        combine_mode=COMBINE_MODE,
        value_basis=VALUE_BASIS,
        output_dir=str(output_path / "loss"),
    )
    estimator.load_damage_library()
    estimator.load_hurricane_probabilities()
    estimator.load_zones()
    estimator.wind_lookup, estimator.flood_lookup = estimator.calculate_marginal_lookups(
        estimator.prob_damage_states, estimator.damage_states[:, 1],
        estimator.n_building_types, estimator.n_component_configurations,
    )
    print("wind_lookup", estimator.wind_lookup.shape, "flood_lookup", estimator.flood_lookup.shape,
          f"NaN pairs: {int(np.isnan(estimator.wind_lookup[:, :, 0]).sum())}")
    if not smoke_mode:
        lookup_dir = full_data_root / FULL_LOOKUP_DIR.relative_to(FULL_DATA_ROOT)
        frozen_wind = np.load(lookup_dir / "L_mcw_wd.npy")
        frozen_flood = np.load(lookup_dir / "L_mcf_fld.npy")
        print("byte-identical to the frozen tables:",
              np.array_equal(estimator.wind_lookup, frozen_wind) and np.array_equal(estimator.flood_lookup, frozen_flood))
    lookup_paths = estimator.save_lookup_tables(str(output_path / "loss"))

1166 of 1536 (m, c) pairs are missing from the damage library; their lookup rows are NaN.

wind_lookup (8, 192, 25) flood_lookup (8, 192, 136) NaN pairs: 1166

In [7]:

Copied!





if data_available:
    legend_c = [1, 3, 5, 7, 25, 49, 97, 145, 169, 175]
    present = sorted({(m, c) for (m, c) in estimator.prob_damage_states if c in legend_c}) or sorted(estimator.prob_damage_states)[:6]
    fig, axes = plt.subplots(1, 2, figsize=(12, 4))
    for m, c in present[:8]:
        axes[0].plot(estimator.wind_speeds[:, 1], estimator.wind_lookup[m - 1, c - 1], label=f"M={m}, C={c}")
        axes[1].plot(estimator.flood_depths[:, 1], estimator.flood_lookup[m - 1, c - 1], label=f"M={m}, C={c}")
    axes[0].set_xlabel("Wind speed (mph)"); axes[0].set_ylabel("Expected loss ratio"); axes[0].set_title("L_mcw_wd (wind only)")
    axes[1].set_xlabel("Flood depth (inches)"); axes[1].set_title("L_mcf_fld (flood only)")
    for ax in axes:
        ax.grid(True, alpha=0.3); ax.legend(fontsize=7)
    plt.tight_layout()
    plt.show()
if data_available:
    legend_c = [1, 3, 5, 7, 25, 49, 97, 145, 169, 175]
    present = sorted({(m, c) for (m, c) in estimator.prob_damage_states if c in legend_c}) or sorted(estimator.prob_damage_states)[:6]
    fig, axes = plt.subplots(1, 2, figsize=(12, 4))
    for m, c in present[:8]:
        axes[0].plot(estimator.wind_speeds[:, 1], estimator.wind_lookup[m - 1, c - 1], label=f"M={m}, C={c}")
        axes[1].plot(estimator.flood_depths[:, 1], estimator.flood_lookup[m - 1, c - 1], label=f"M={m}, C={c}")
    axes[0].set_xlabel("Wind speed (mph)"); axes[0].set_ylabel("Expected loss ratio"); axes[0].set_title("L_mcw_wd (wind only)")
    axes[1].set_xlabel("Flood depth (inches)"); axes[1].set_title("L_mcf_fld (flood only)")
    for ax in axes:
        ax.grid(True, alpha=0.3); ax.legend(fontsize=7)
    plt.tight_layout()
    plt.show()

Index conventions¶

Wind_Tr and Flood_Tr are continuous 1-based IDs. floor (int(id) - 1) reproduces the published Eastern North Carolina maps; ceil is the verified STARR export's searchsorted mode; round and linear are alternatives. On the frozen run the wind EAL maximum under ceil is 40.7% above floor, so the convention must accompany every result.

In [8]:

Copied!





if data_available:
    try:
        if smoke_mode:
            estimator.load_from_exposure(exposure)
        else:
            estimator.inventory_csv_path = str(export_paths["inventory_csv"])
            estimator.wind_ids_path = str(export_paths["wind_ids"])
            estimator.flood_ids_path = str(export_paths["flood_ids"])
            estimator.load_exposure_inputs()
        print(f"{len(estimator.building_table)} buildings x {len(estimator.event_names)} events loaded")
        annual_probabilities = estimator.hurricane_prob[:, 1]
        mode_totals = {}
        for mode in ("floor", "ceil", "round", "linear"):
            estimator.index_mode = mode
            eal_mode = estimator.calculate_building_eal(
                estimator.building_table, estimator.wind_ids, estimator.flood_ids, annual_probabilities
            )
            mode_totals[mode] = {"eal_wind": eal_mode["eal_wind"].sum(), "eal_flood": eal_mode["eal_flood"].sum(),
                                 "eal_combined": eal_mode["eal_combined"].sum()}
        estimator.index_mode = INDEX_MODE
        mode_table = pd.DataFrame(mode_totals).T
        mode_table["combined_vs_floor"] = mode_table["eal_combined"] / mode_table.loc["floor", "eal_combined"] - 1.0
        shown = mode_table.copy()
        for column in ("eal_wind", "eal_flood", "eal_combined"):
            shown[column] = shown[column].map("${:,.0f}".format)
        shown["combined_vs_floor"] = shown["combined_vs_floor"].map("{:+.1%}".format)
        display(shown)
    except BaseException:
        restore_state()
        raise
if data_available:
    try:
        if smoke_mode:
            estimator.load_from_exposure(exposure)
        else:
            estimator.inventory_csv_path = str(export_paths["inventory_csv"])
            estimator.wind_ids_path = str(export_paths["wind_ids"])
            estimator.flood_ids_path = str(export_paths["flood_ids"])
            estimator.load_exposure_inputs()
        print(f"{len(estimator.building_table)} buildings x {len(estimator.event_names)} events loaded")
        annual_probabilities = estimator.hurricane_prob[:, 1]
        mode_totals = {}
        for mode in ("floor", "ceil", "round", "linear"):
            estimator.index_mode = mode
            eal_mode = estimator.calculate_building_eal(
                estimator.building_table, estimator.wind_ids, estimator.flood_ids, annual_probabilities
            )
            mode_totals[mode] = {"eal_wind": eal_mode["eal_wind"].sum(), "eal_flood": eal_mode["eal_flood"].sum(),
                                 "eal_combined": eal_mode["eal_combined"].sum()}
        estimator.index_mode = INDEX_MODE
        mode_table = pd.DataFrame(mode_totals).T
        mode_table["combined_vs_floor"] = mode_table["eal_combined"] / mode_table.loc["floor", "eal_combined"] - 1.0
        shown = mode_table.copy()
        for column in ("eal_wind", "eal_flood", "eal_combined"):
            shown[column] = shown[column].map("${:,.0f}".format)
        shown["combined_vs_floor"] = shown["combined_vs_floor"].map("{:+.1%}".format)
        display(shown)
    except BaseException:
        restore_state()
        raise

2000 buildings x 5 events loaded

	eal_wind	eal_flood	eal_combined	combined_vs_floor
floor	$1,811,293	$9,624,792	$11,436,084	+0.0%
ceil	$2,677,880	$10,012,771	$12,690,651	+11.0%
round	$2,339,602	$9,694,365	$12,033,968	+5.2%
linear	$2,294,029	$9,759,553	$12,053,582	+5.4%

Expected annual loss per building¶

EAL_b = sum_h p_h * lr_bh * value_b with p_h the annual event probability (not normalized). The value basis changes only the monetary scale; home_value = structure_value + contents_value.

In [9]:

Copied!





if data_available:
    basis_totals = {}
    for basis in ("structure_value", "contents_value", "home_value"):
        estimator.value_basis = basis
        eal_basis = estimator.calculate_building_eal(
            estimator.building_table, estimator.wind_ids, estimator.flood_ids, annual_probabilities
        )
        basis_totals[basis] = eal_basis[["eal_wind", "eal_flood", "eal_combined"]].sum()
    estimator.value_basis = VALUE_BASIS
    estimator.building_eal = estimator.calculate_building_eal(
        estimator.building_table, estimator.wind_ids, estimator.flood_ids, annual_probabilities
    )
    display(pd.DataFrame(basis_totals).T.map("${:,.0f}".format))
    fig, ax = plt.subplots(figsize=(7, 4))
    positive = estimator.building_eal["eal_combined"][estimator.building_eal["eal_combined"] > 0]
    ax.hist(np.log10(positive), bins=40, color="#55a868")
    ax.set_xlabel(f"log10 combined EAL per building ($/yr, {VALUE_BASIS}, {INDEX_MODE})"); ax.set_ylabel("Buildings")
    ax.set_title("Building expected annual loss"); ax.grid(True, alpha=0.3)
    plt.show()
if data_available:
    basis_totals = {}
    for basis in ("structure_value", "contents_value", "home_value"):
        estimator.value_basis = basis
        eal_basis = estimator.calculate_building_eal(
            estimator.building_table, estimator.wind_ids, estimator.flood_ids, annual_probabilities
        )
        basis_totals[basis] = eal_basis[["eal_wind", "eal_flood", "eal_combined"]].sum()
    estimator.value_basis = VALUE_BASIS
    estimator.building_eal = estimator.calculate_building_eal(
        estimator.building_table, estimator.wind_ids, estimator.flood_ids, annual_probabilities
    )
    display(pd.DataFrame(basis_totals).T.map("${:,.0f}".format))
    fig, ax = plt.subplots(figsize=(7, 4))
    positive = estimator.building_eal["eal_combined"][estimator.building_eal["eal_combined"] > 0]
    ax.hist(np.log10(positive), bins=40, color="#55a868")
    ax.set_xlabel(f"log10 combined EAL per building ($/yr, {VALUE_BASIS}, {INDEX_MODE})"); ax.set_ylabel("Buildings")
    ax.set_title("Building expected annual loss"); ax.grid(True, alpha=0.3)
    plt.show()

	eal_wind	eal_flood	eal_combined
structure_value	$1,811,293	$9,624,792	$11,436,084
contents_value	$905,646	$4,812,396	$5,718,042
home_value	$2,716,939	$14,437,187	$17,154,126

Zone aggregation¶

Building EAL is summed within zones. aggregation_col remaps the inventory's zone column onto a coarser nested system through the zone layer (i_lh740 -> i_lh503 for the Eastern North Carolina zones); buildings outside every zone are reported as unit -1.

In [10]:

Copied!





if data_available:
    estimator.aggregation_key = estimator.calculate_aggregation_key(
        estimator.building_table, estimator.zones, ZONE_ID, TARGET_ZONE_ID
    )
    estimator.zone_eal = estimator.calculate_zone_eal(estimator.building_eal, estimator.aggregation_key)
    estimator.event_losses = estimator.calculate_event_losses(
        estimator.building_table, estimator.wind_ids, estimator.flood_ids, annual_probabilities
    )
    top_zones = estimator.zone_eal.drop(index=-1, errors="ignore").sort_values("eal_combined", ascending=False).head(15)
    shown_zones = top_zones.copy()
    for column in ("eal_wind", "eal_flood", "eal_combined", "value"):
        shown_zones[column] = shown_zones[column].map("${:,.0f}".format)
    display(shown_zones)
    fig, ax = plt.subplots(figsize=(9, 4))
    top_zones[["eal_wind", "eal_flood"]].plot.bar(stacked=True, ax=ax, color=["#4c72b0", "#c44e52"])
    ax.set_xlabel(f"Zone ({TARGET_ZONE_ID})"); ax.set_ylabel("Expected annual loss ($/yr)")
    ax.set_title("Zones with the largest combined EAL"); ax.grid(True, axis="y", alpha=0.3)
    plt.tight_layout()
    plt.show()
    annual_from_events = (estimator.event_losses["loss_combined"] * estimator.event_losses["probability"]).sum()
    print(f"sum_h p_h * regional loss_h = ${annual_from_events:,.0f}/yr; sum of building EAL = ${estimator.building_eal['eal_combined'].sum():,.0f}/yr")
if data_available:
    estimator.aggregation_key = estimator.calculate_aggregation_key(
        estimator.building_table, estimator.zones, ZONE_ID, TARGET_ZONE_ID
    )
    estimator.zone_eal = estimator.calculate_zone_eal(estimator.building_eal, estimator.aggregation_key)
    estimator.event_losses = estimator.calculate_event_losses(
        estimator.building_table, estimator.wind_ids, estimator.flood_ids, annual_probabilities
    )
    top_zones = estimator.zone_eal.drop(index=-1, errors="ignore").sort_values("eal_combined", ascending=False).head(15)
    shown_zones = top_zones.copy()
    for column in ("eal_wind", "eal_flood", "eal_combined", "value"):
        shown_zones[column] = shown_zones[column].map("${:,.0f}".format)
    display(shown_zones)
    fig, ax = plt.subplots(figsize=(9, 4))
    top_zones[["eal_wind", "eal_flood"]].plot.bar(stacked=True, ax=ax, color=["#4c72b0", "#c44e52"])
    ax.set_xlabel(f"Zone ({TARGET_ZONE_ID})"); ax.set_ylabel("Expected annual loss ($/yr)")
    ax.set_title("Zones with the largest combined EAL"); ax.grid(True, axis="y", alpha=0.3)
    plt.tight_layout()
    plt.show()
    annual_from_events = (estimator.event_losses["loss_combined"] * estimator.event_losses["probability"]).sum()
    print(f"sum_h p_h * regional loss_h = ${annual_from_events:,.0f}/yr; sum of building EAL = ${estimator.building_eal['eal_combined'].sum():,.0f}/yr")

	eal_wind	eal_flood	eal_combined	n_buildings	value
unit
1	$1,811,293	$9,624,792	$11,436,084	2000	$600,842,304

sum_h p_h * regional loss_h = $11,436,084/yr; sum of building EAL = $11,436,084/yr

Orchestrator¶

LossEstimation.estimate_all_losses() recomputes every result from the current inputs and returns them in one dictionary; it must agree with the stepwise results above.

In [11]:

Copied!





if data_available:
    results = estimator.estimate_all_losses()
    pd.testing.assert_frame_equal(results["building_eal"], estimator.building_eal)
    pd.testing.assert_frame_equal(results["zone_eal"], estimator.zone_eal)
    pd.testing.assert_frame_equal(results["event_losses"], estimator.event_losses)
    saved = estimator.save_all_outputs(str(output_path / "loss"))
    print("orchestrator results match the stepwise results")
    print(json.dumps(saved, indent=2))
if data_available:
    results = estimator.estimate_all_losses()
    pd.testing.assert_frame_equal(results["building_eal"], estimator.building_eal)
    pd.testing.assert_frame_equal(results["zone_eal"], estimator.zone_eal)
    pd.testing.assert_frame_equal(results["event_losses"], estimator.event_losses)
    saved = estimator.save_all_outputs(str(output_path / "loss"))
    print("orchestrator results match the stepwise results")
    print(json.dumps(saved, indent=2))

orchestrator results match the stepwise results
{
  "wind_lookup": "/home/mammad/Documents/EDU/cope/CHEERSafe/LossModel-exposure/Outputs/loss_estimation_workflow/loss/L_mcw_wd.npy",
  "flood_lookup": "/home/mammad/Documents/EDU/cope/CHEERSafe/LossModel-exposure/Outputs/loss_estimation_workflow/loss/L_mcf_fld.npy",
  "metadata": "/home/mammad/Documents/EDU/cope/CHEERSafe/LossModel-exposure/Outputs/loss_estimation_workflow/loss/STARR_lookup_metadata.txt",
  "building_eal": "/home/mammad/Documents/EDU/cope/CHEERSafe/LossModel-exposure/Outputs/loss_estimation_workflow/loss/building_eal.csv",
  "zone_eal": "/home/mammad/Documents/EDU/cope/CHEERSafe/LossModel-exposure/Outputs/loss_estimation_workflow/loss/zone_eal.csv",
  "event_losses": "/home/mammad/Documents/EDU/cope/CHEERSafe/LossModel-exposure/Outputs/loss_estimation_workflow/loss/event_losses.csv"
}

In [12]:

Copied!

# Restore the configuration changed above and return to the starting directory
restore_state()
print("configuration restored; working directory:", Path.cwd())
# Restore the configuration changed above and return to the starting directory
restore_state()
print("configuration restored; working directory:", Path.cwd())

configuration restored; working directory: /home/mammad/Documents/EDU/cope/CHEERSafe/LossModel-exposure/examples

Results, pitfalls, and next steps¶

Outputs are under Outputs/loss_estimation_workflow/loss/: the lookup tables, building_eal.csv, zone_eal.csv and event_losses.csv.
Always report the index convention with a result; floor reproduces the published maps, ceil the verified STARR export. The additive combination discards the wind-flood interaction of the library; use COMBINE_MODE = "joint" to keep it.
home_value overshoots the published Eastern North Carolina panels by 50% to 111% because contents are half the structure value in that inventory.
The Eastern North Carolina figures notebook draws the published hazard, inventory and EAL maps from the same chain.