Eastern North Carolina Loss Figures¶

This notebook reproduces the seven Eastern North Carolina figures of the CHEERSafe application section from the package's exposure and loss-estimation modules:

Figure	Content	Package numerics
1	Modal structural archetype `M` per 1 km cell	`calculate_cell_index`
2	Modal resistance configuration `C` per 1 km cell (ten most common + Other)	`calculate_cell_index`
3	Probability-weighted mean wind speed per hazard cell	`calculate_probability_weighted_intensity`
4	Probability-weighted mean flood depth conditional on inundation	`calculate_probability_weighted_intensity`
5-7	Expected annual loss (wind, flood, combined) per LEAPHI-503 zone	`LossEstimation`

The plotting helpers live in this notebook; every number comes from cheersafe.exposure and cheersafe.loss_estimation.

!!! warning "Execution profiles" Smoke mode (default) draws all seven figures at small scale on a synthetic window (2,000 buildings, a 12 x 10 km grid, five events, six zones) without external data or published-value assertions.

**Full mode** (`CHEERSAFE_FIGURES_SMOKE=0`) reads from the `Loss Estimates` folder
(`CHEERSAFE_FIGURES_DATA`): the 97 `Old97` hazard layers, the frozen 07/07/2026 export
(1,154,657 buildings, memory-mapped), the frozen lookup tables (or the 997 MB damage library when the
tables are absent) and the LEAPHI zones. It asserts the published modal-C legend and the three published
EAL colorbar maxima within 0.5%. Allow about 2 GB of RAM and a few minutes. County boundaries are
optional: point `CHEERSAFE_FIGURES_COUNTIES` at a TIGER 2020 `tl_2020_us_county.shp`
(https://www2.census.gov/geo/tiger/TIGER2020/COUNTY/tl_2020_us_county.zip, 130 MB, not shipped);
without it the maps show only the 1 km study-area outline.

The frozen export is the full-mode default because the exposure stage of the source run used an
unseeded configuration draw; re-running the exposure chain from the raw inventory reproduces the
archetype distributions but not the row-by-row `C` values behind figures 2 and 5-7.

Plan¶

Choose the execution profile, the aggregation unit, the value basis and the index convention.
Define the plotting helpers (grid rasters, study-area outline, categorical and continuous maps).
Load (or generate) the exposure export, the hazard catalogue, the zones and the lookup source.
Figures 1-2: modal M and C per 1 km cell.
Figures 3-4: probability-weighted wind and conditional flood per hazard cell.
Figures 5-7: expected annual loss per zone, with the index-convention and value-basis sensitivity table.
Write the numeric checks next to the figures.

In [1]:

Copied!





# Repository setup and imports
import collections
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 matplotlib.cm import ScalarMappable
from matplotlib.collections import LineCollection
from matplotlib.colors import BoundaryNorm, ListedColormap, Normalize
from shapely.geometry import box
from shapely.ops import unary_union

from cheersafe.config import Config
from cheersafe.exposure import BaseExposure, Exposure, SyntheticExposureData
from cheersafe.exposure.hazard_mapping import (
    calculate_cell_index,
    calculate_grid_mask,
    calculate_probability_weighted_intensity,
)
from cheersafe.exposure.synthetic import GRID_CELL_SIZE, GRID_CRS, GRID_ORIGIN
from cheersafe.loss_estimation import LossEstimation
from cheersafe.tools import Utils

print(f"repo_root={repo_root}")
# Repository setup and imports
import collections
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 matplotlib.cm import ScalarMappable
from matplotlib.collections import LineCollection
from matplotlib.colors import BoundaryNorm, ListedColormap, Normalize
from shapely.geometry import box
from shapely.ops import unary_union

from cheersafe.config import Config
from cheersafe.exposure import BaseExposure, Exposure, SyntheticExposureData
from cheersafe.exposure.hazard_mapping import (
    calculate_cell_index,
    calculate_grid_mask,
    calculate_probability_weighted_intensity,
)
from cheersafe.exposure.synthetic import GRID_CELL_SIZE, GRID_CRS, GRID_ORIGIN
from cheersafe.loss_estimation import 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 window by default; full mode reads the frozen Eastern North Carolina export
FULL_DATA_ROOT = Path(os.environ.get("CHEERSAFE_FIGURES_DATA", "../Loss Estimates"))
FULL_HAZARD_DIR = FULL_DATA_ROOT / "Hazards" / "Old97"
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_LIBRARY = FULL_DATA_ROOT / "CHEER_Safe_data" / "ProbDS_cwfd_MaxIter3000_V2.pkl"
FULL_ZONES = FULL_DATA_ROOT / "Zones" / "leaphi_zones_20230920.shp"
PROB_JSON = Path("Inputs/Hurricane_ProbWindFlood.json")
COUNTY_BOUNDARIES = os.environ.get("CHEERSAFE_FIGURES_COUNTIES")  # TIGER 2020 tl_2020_us_county.shp; None draws no county lines

AGG_UNIT = "i_lh503"              # the published EAL maps aggregate to LEAPHI-503 zones
ZONE_ID = "i_lh740"               # zone column of the exposure export
VALUE_FIELD = "structure_value"
HAZARD_INDEX = "floor"
COMBINE_MODE = "additive"
MODAL_C_LEGEND = [1, 3, 5, 7, 25, 49, 97, 145, 169, 175]
PUBLISHED_VMAX = {"wind": 6.487e6, "flood": 68.784e6, "combined": 71.806e6}
VMAX_TOLERANCE = 0.005

SMOKE_N_BUILDINGS = 2000
SMOKE_N_EVENTS = 5
SMOKE_GRID = (12, 10)
SMOKE_SEED = 42

smoke_flag = os.environ.get("CHEERSAFE_FIGURES_SMOKE", "1")
if smoke_flag not in {"0", "1"}:
    raise ValueError("CHEERSAFE_FIGURES_SMOKE must be '0' or '1'")
smoke_mode = smoke_flag == "1"
output_override = os.environ.get("CHEERSAFE_FIGURES_OUTPUT")
output_path = (
    Path(output_override).expanduser().resolve()
    if output_override
    else (repo_root / "Outputs/eastern_nc_loss_figures").resolve()
)
figures_path = output_path / "figures"
checks_path = output_path / "checks"
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'}, unit={AGG_UNIT}, value={VALUE_FIELD}, index={HAZARD_INDEX}, output={output_path}")
# Execution profile: smoke mode on a synthetic window by default; full mode reads the frozen Eastern North Carolina export
FULL_DATA_ROOT = Path(os.environ.get("CHEERSAFE_FIGURES_DATA", "../Loss Estimates"))
FULL_HAZARD_DIR = FULL_DATA_ROOT / "Hazards" / "Old97"
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_LIBRARY = FULL_DATA_ROOT / "CHEER_Safe_data" / "ProbDS_cwfd_MaxIter3000_V2.pkl"
FULL_ZONES = FULL_DATA_ROOT / "Zones" / "leaphi_zones_20230920.shp"
PROB_JSON = Path("Inputs/Hurricane_ProbWindFlood.json")
COUNTY_BOUNDARIES = os.environ.get("CHEERSAFE_FIGURES_COUNTIES")  # TIGER 2020 tl_2020_us_county.shp; None draws no county lines

AGG_UNIT = "i_lh503"              # the published EAL maps aggregate to LEAPHI-503 zones
ZONE_ID = "i_lh740"               # zone column of the exposure export
VALUE_FIELD = "structure_value"
HAZARD_INDEX = "floor"
COMBINE_MODE = "additive"
MODAL_C_LEGEND = [1, 3, 5, 7, 25, 49, 97, 145, 169, 175]
PUBLISHED_VMAX = {"wind": 6.487e6, "flood": 68.784e6, "combined": 71.806e6}
VMAX_TOLERANCE = 0.005

SMOKE_N_BUILDINGS = 2000
SMOKE_N_EVENTS = 5
SMOKE_GRID = (12, 10)
SMOKE_SEED = 42

smoke_flag = os.environ.get("CHEERSAFE_FIGURES_SMOKE", "1")
if smoke_flag not in {"0", "1"}:
    raise ValueError("CHEERSAFE_FIGURES_SMOKE must be '0' or '1'")
smoke_mode = smoke_flag == "1"
output_override = os.environ.get("CHEERSAFE_FIGURES_OUTPUT")
output_path = (
    Path(output_override).expanduser().resolve()
    if output_override
    else (repo_root / "Outputs/eastern_nc_loss_figures").resolve()
)
figures_path = output_path / "figures"
checks_path = output_path / "checks"
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'}, unit={AGG_UNIT}, value={VALUE_FIELD}, index={HAZARD_INDEX}, output={output_path}")

profile=full, unit=i_lh503, value=structure_value, index=floor, output=/home/mammad/Documents/EDU/cope/CHEERSafe/LossModel-exposure/Outputs/eastern_nc_loss_figures

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
Config.AUTO_EXECUTE_METHODS = False
Config.INPUTS_PATH = "Inputs"
figures_path.mkdir(parents=True, exist_ok=True)
checks_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
Config.AUTO_EXECUTE_METHODS = False
Config.INPUTS_PATH = "Inputs"
figures_path.mkdir(parents=True, exist_ok=True)
checks_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}

Plotting helpers¶

Ported from the figure-reproduction scripts: rasters on the 1 km grid, the staircase outline of the union of hazard cells (the study area), optional county lines, categorical maps with tab10 / tab20 palettes, continuous viridis maps and zone choropleths with a dollar-formatted colorbar.

In [4]:

Copied!





LAND_GREY = "#eeeeee"
LINE_GREY = "#696969"
COUNTY_GREY = "#bfbfbf"
FIGSIZE = (9.0, 7.5)
SAVE_DPI = 150
X0, Y0 = GRID_ORIGIN
DX = GRID_CELL_SIZE
HALF = DX / 2.0


def raster_extent(nx, ny):
    """imshow extent of an (ny, nx) raster whose cells are centred on the grid nodes."""
    return (X0 - HALF, X0 + nx * DX - HALF, Y0 - HALF, Y0 + ny * DX - HALF)


def to_raster(values, ix, iy, nx, ny):
    """Scatter per-cell values onto an (ny, nx) raster, masked where no cell exists."""
    raster = np.full((ny, nx), np.nan, dtype="float64")
    raster[iy, ix] = values
    return np.ma.masked_invalid(raster)


def mask_boundary_segments(mask, extent, dx):
    """Edges of the union of the True cells of `mask`, as (n, 2, 2) line segments."""
    x0, _, y0, _ = extent
    segments = []

    def push(ci, ri, horizontal, offset):
        gx = x0 + ci * dx
        gy = y0 + ri * dx
        if horizontal:
            segments.append(np.stack([np.stack([gx, gy + offset], 1), np.stack([gx + dx, gy + offset], 1)], 1))
        else:
            segments.append(np.stack([np.stack([gx + offset, gy], 1), np.stack([gx + offset, gy + dx], 1)], 1))

    left = mask & ~np.pad(mask, ((0, 0), (1, 0)))[:, :-1]
    right = mask & ~np.pad(mask, ((0, 0), (0, 1)))[:, 1:]
    below = mask & ~np.pad(mask, ((1, 0), (0, 0)))[:-1, :]
    above = mask & ~np.pad(mask, ((0, 1), (0, 0)))[1:, :]
    for edge_mask, horizontal, offset in ((left, False, 0.0), (right, False, dx), (below, True, 0.0), (above, True, dx)):
        rows, cols = np.nonzero(edge_mask)
        if rows.size:
            push(cols.astype("float64"), rows.astype("float64"), horizontal, offset)
    return np.concatenate(segments, axis=0)


def study_area(mask, extent, dx):
    """The staircase study area: the union of the hazard cell squares."""
    ys, xs = np.nonzero(mask)
    x0, _, y0, _ = extent
    cells = [box(x0 + x * dx, y0 + y * dx, x0 + (x + 1) * dx, y0 + (y + 1) * dx) for y, x in zip(ys.tolist(), xs.tolist())]
    return unary_union(cells)


def clip_to_study(gdf, area):
    """Clip a layer to the study area, dropping anything fully outside."""
    out = gdf.copy()
    out["geometry"] = out.geometry.intersection(area)
    keep = (~out.geometry.is_empty) & out.geometry.notna()
    return out[keep]


def load_counties(path, extent, crs, area):
    """North Carolina counties (TIGER 2020) clipped to the study area, or None."""
    if not path or not Path(path).exists():
        return None
    counties = gpd.read_file(path, columns=["STATEFP", "COUNTYFP", "GEOID"])
    counties = counties[counties["STATEFP"] == "37"].to_crs(crs)
    counties = counties[counties.intersects(box(extent[0], extent[2], extent[1], extent[3]))]
    return clip_to_study(counties, area)


def frame_axes(ax, extent):
    width = extent[1] - extent[0]
    height = extent[3] - extent[2]
    ax.set_xlim(extent[0] - 0.05 * width, extent[1] + 0.05 * width)
    ax.set_ylim(extent[2] - 0.05 * height, extent[3] + 0.05 * height)
    ax.set_aspect("equal")
    ax.set_axis_off()


def draw_basemap(ax, mask, extent, counties):
    ax.imshow(np.ma.masked_where(~mask, np.zeros(mask.shape)), origin="lower", extent=extent,
              cmap=ListedColormap([LAND_GREY]), vmin=0, vmax=1, interpolation="nearest", zorder=1)
    if counties is not None and len(counties):
        counties.boundary.plot(ax=ax, color=COUNTY_GREY, linewidth=0.6, zorder=3)


def draw_outline(ax, segments):
    ax.add_collection(LineCollection(segments, colors=LINE_GREY, linewidths=0.8, zorder=4, capstyle="projecting"))


def finish(fig, out_png):
    fig.patch.set_facecolor("white")
    fig.savefig(out_png, dpi=SAVE_DPI, bbox_inches="tight", facecolor="white")
    plt.show()
    plt.close(fig)


def draw_raster_map(raster, extent, mask, counties, segments, label, out_png, vmin=None, vmax=None, with_base=True):
    fig, ax = plt.subplots(figsize=FIGSIZE)
    if with_base:
        draw_basemap(ax, mask, extent, counties)
    image = ax.imshow(raster, extent=extent, origin="lower", cmap="viridis", vmin=vmin, vmax=vmax,
                      interpolation="nearest", zorder=2)
    if counties is not None and len(counties):
        counties.boundary.plot(ax=ax, color="white", linewidth=0.5, zorder=3)
    draw_outline(ax, segments)
    frame_axes(ax, extent)
    colorbar = fig.colorbar(image, ax=ax, fraction=0.033, pad=0.02, aspect=26)
    colorbar.set_label(label, fontsize=14)
    finish(fig, out_png)


def categorical_map(codes, occupied, cell_extent, mask, extent, segments, counties, cmap, labels, cbar_label, out_png):
    data = np.ma.masked_where(~occupied, codes)
    n_categories = len(labels)
    norm = BoundaryNorm(np.arange(-0.5, n_categories + 0.5, 1.0), n_categories)
    fig, ax = plt.subplots(figsize=FIGSIZE)
    draw_basemap(ax, mask, extent, counties)
    image = ax.imshow(data, origin="lower", extent=cell_extent, cmap=cmap, norm=norm, interpolation="nearest", zorder=2)
    draw_outline(ax, segments)
    frame_axes(ax, extent)
    colorbar = fig.colorbar(image, ax=ax, fraction=0.040, pad=0.02, ticks=np.arange(n_categories))
    colorbar.ax.set_yticklabels(labels)
    colorbar.ax.tick_params(length=0, labelsize=11)
    colorbar.set_label(cbar_label, rotation=270, labelpad=18, fontsize=14)
    finish(fig, out_png)


def money_fmt(value, _pos=None):
    if value == 0:
        return "$0"
    if abs(value) >= 1e6:
        return f"${value / 1e6:.0f}M"
    return f"${value / 1e3:.0f}k"


def draw_zone_map(polygons, values, counties, segments, extent, out_png):
    """Zone choropleth on a full-range viridis scale; returns the colorbar maximum."""
    layer = polygons.copy()
    layer["val"] = values.reindex(layer["unit"]).fillna(0.0).to_numpy()
    vmax = float(layer["val"].max())
    norm = Normalize(vmin=0.0, vmax=vmax)
    fig, ax = plt.subplots(figsize=FIGSIZE)
    layer.plot(ax=ax, column="val", cmap="viridis", norm=norm, edgecolor="white", linewidth=0.15, zorder=2)
    if counties is not None and len(counties):
        counties.boundary.plot(ax=ax, color="white", linewidth=0.9, zorder=3)
    ax.add_collection(LineCollection(segments, colors="0.45", linewidths=1.0, zorder=4, capstyle="projecting"))
    frame_axes(ax, extent)
    colorbar = fig.colorbar(ScalarMappable(norm=norm, cmap="viridis"), ax=ax, fraction=0.046, pad=0.02)
    colorbar.ax.yaxis.set_major_formatter(plt.FuncFormatter(money_fmt))
    colorbar.outline.set_visible(False)
    finish(fig, out_png)
    return vmax


def modal_grids(x, y, archetypes, configurations):
    """Per-cell modal M and C on the 1 km grid; ties go to the smallest value."""
    ix, iy = calculate_cell_index(x, y, method="nearest")
    assert ix.min() >= 0 and iy.min() >= 0, "a building falls left of or below the grid origin"
    nx = int(ix.max()) + 1
    ny = int(iy.max()) + 1
    flat = iy * nx + ix
    n_cells = nx * ny

    def mode_of(values, n_values):
        counts = np.bincount(flat * (n_values + 1) + values, minlength=n_cells * (n_values + 1)).reshape(n_cells, n_values + 1)
        best = counts.argmax(axis=1)
        top = counts.max(axis=1)
        total = counts.sum(axis=1)
        tied = ((counts == top[:, None]).sum(axis=1) > 1) & (total > 0)
        return best, total, tied

    mode_m, total, tied_m = mode_of(archetypes, 8)
    mode_c, _, tied_c = mode_of(configurations, 192)
    shape = (ny, nx)
    return dict(nx=nx, ny=ny, mode_m=mode_m.reshape(shape), mode_c=mode_c.reshape(shape),
                count=total.reshape(shape), occupied=(total > 0).reshape(shape),
                tied_m=int(tied_m.sum()), tied_c=int(tied_c.sum()))


def make_zone_grid(origin, extent_cells, cell, n_x=3, n_y=2):
    """Synthetic LEAPHI-style zone grid: n_x * n_y fine zones nesting pairwise into coarse zones."""
    x0, y0 = origin
    nx, ny = extent_cells
    xmin, ymin = x0 - cell / 2.0, y0 - cell / 2.0
    width, height = nx * cell, ny * cell
    rows = []
    for j in range(n_y):
        for i in range(n_x):
            k = j * n_x + i
            rows.append({
                "COUNTYFP00": f"{101 + k:03d}", "zone": 1.0, "i_lh503": float(k // 2 + 1),
                "i_lh1509": float(k + 1), "i_lh708": float(k + 1), "i_lh740": float(k + 1),
                "geometry": box(xmin + i * width / n_x, ymin + j * height / n_y,
                                xmin + (i + 1) * width / n_x, ymin + (j + 1) * height / n_y),
            })
    return gpd.GeoDataFrame(rows, crs=GRID_CRS)


def run_exposure(inventory_path, hazard_dir, zones_path, seed):
    """Run the exposure chain with explicit calls; every result is assigned here."""
    exposure = Exposure(inventory_path=str(inventory_path), hazard_dir=str(hazard_dir), zones_path=str(zones_path),
                        average_mode="Nearest", 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)
    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
LAND_GREY = "#eeeeee"
LINE_GREY = "#696969"
COUNTY_GREY = "#bfbfbf"
FIGSIZE = (9.0, 7.5)
SAVE_DPI = 150
X0, Y0 = GRID_ORIGIN
DX = GRID_CELL_SIZE
HALF = DX / 2.0


def raster_extent(nx, ny):
    """imshow extent of an (ny, nx) raster whose cells are centred on the grid nodes."""
    return (X0 - HALF, X0 + nx * DX - HALF, Y0 - HALF, Y0 + ny * DX - HALF)


def to_raster(values, ix, iy, nx, ny):
    """Scatter per-cell values onto an (ny, nx) raster, masked where no cell exists."""
    raster = np.full((ny, nx), np.nan, dtype="float64")
    raster[iy, ix] = values
    return np.ma.masked_invalid(raster)


def mask_boundary_segments(mask, extent, dx):
    """Edges of the union of the True cells of `mask`, as (n, 2, 2) line segments."""
    x0, _, y0, _ = extent
    segments = []

    def push(ci, ri, horizontal, offset):
        gx = x0 + ci * dx
        gy = y0 + ri * dx
        if horizontal:
            segments.append(np.stack([np.stack([gx, gy + offset], 1), np.stack([gx + dx, gy + offset], 1)], 1))
        else:
            segments.append(np.stack([np.stack([gx + offset, gy], 1), np.stack([gx + offset, gy + dx], 1)], 1))

    left = mask & ~np.pad(mask, ((0, 0), (1, 0)))[:, :-1]
    right = mask & ~np.pad(mask, ((0, 0), (0, 1)))[:, 1:]
    below = mask & ~np.pad(mask, ((1, 0), (0, 0)))[:-1, :]
    above = mask & ~np.pad(mask, ((0, 1), (0, 0)))[1:, :]
    for edge_mask, horizontal, offset in ((left, False, 0.0), (right, False, dx), (below, True, 0.0), (above, True, dx)):
        rows, cols = np.nonzero(edge_mask)
        if rows.size:
            push(cols.astype("float64"), rows.astype("float64"), horizontal, offset)
    return np.concatenate(segments, axis=0)


def study_area(mask, extent, dx):
    """The staircase study area: the union of the hazard cell squares."""
    ys, xs = np.nonzero(mask)
    x0, _, y0, _ = extent
    cells = [box(x0 + x * dx, y0 + y * dx, x0 + (x + 1) * dx, y0 + (y + 1) * dx) for y, x in zip(ys.tolist(), xs.tolist())]
    return unary_union(cells)


def clip_to_study(gdf, area):
    """Clip a layer to the study area, dropping anything fully outside."""
    out = gdf.copy()
    out["geometry"] = out.geometry.intersection(area)
    keep = (~out.geometry.is_empty) & out.geometry.notna()
    return out[keep]


def load_counties(path, extent, crs, area):
    """North Carolina counties (TIGER 2020) clipped to the study area, or None."""
    if not path or not Path(path).exists():
        return None
    counties = gpd.read_file(path, columns=["STATEFP", "COUNTYFP", "GEOID"])
    counties = counties[counties["STATEFP"] == "37"].to_crs(crs)
    counties = counties[counties.intersects(box(extent[0], extent[2], extent[1], extent[3]))]
    return clip_to_study(counties, area)


def frame_axes(ax, extent):
    width = extent[1] - extent[0]
    height = extent[3] - extent[2]
    ax.set_xlim(extent[0] - 0.05 * width, extent[1] + 0.05 * width)
    ax.set_ylim(extent[2] - 0.05 * height, extent[3] + 0.05 * height)
    ax.set_aspect("equal")
    ax.set_axis_off()


def draw_basemap(ax, mask, extent, counties):
    ax.imshow(np.ma.masked_where(~mask, np.zeros(mask.shape)), origin="lower", extent=extent,
              cmap=ListedColormap([LAND_GREY]), vmin=0, vmax=1, interpolation="nearest", zorder=1)
    if counties is not None and len(counties):
        counties.boundary.plot(ax=ax, color=COUNTY_GREY, linewidth=0.6, zorder=3)


def draw_outline(ax, segments):
    ax.add_collection(LineCollection(segments, colors=LINE_GREY, linewidths=0.8, zorder=4, capstyle="projecting"))


def finish(fig, out_png):
    fig.patch.set_facecolor("white")
    fig.savefig(out_png, dpi=SAVE_DPI, bbox_inches="tight", facecolor="white")
    plt.show()
    plt.close(fig)


def draw_raster_map(raster, extent, mask, counties, segments, label, out_png, vmin=None, vmax=None, with_base=True):
    fig, ax = plt.subplots(figsize=FIGSIZE)
    if with_base:
        draw_basemap(ax, mask, extent, counties)
    image = ax.imshow(raster, extent=extent, origin="lower", cmap="viridis", vmin=vmin, vmax=vmax,
                      interpolation="nearest", zorder=2)
    if counties is not None and len(counties):
        counties.boundary.plot(ax=ax, color="white", linewidth=0.5, zorder=3)
    draw_outline(ax, segments)
    frame_axes(ax, extent)
    colorbar = fig.colorbar(image, ax=ax, fraction=0.033, pad=0.02, aspect=26)
    colorbar.set_label(label, fontsize=14)
    finish(fig, out_png)


def categorical_map(codes, occupied, cell_extent, mask, extent, segments, counties, cmap, labels, cbar_label, out_png):
    data = np.ma.masked_where(~occupied, codes)
    n_categories = len(labels)
    norm = BoundaryNorm(np.arange(-0.5, n_categories + 0.5, 1.0), n_categories)
    fig, ax = plt.subplots(figsize=FIGSIZE)
    draw_basemap(ax, mask, extent, counties)
    image = ax.imshow(data, origin="lower", extent=cell_extent, cmap=cmap, norm=norm, interpolation="nearest", zorder=2)
    draw_outline(ax, segments)
    frame_axes(ax, extent)
    colorbar = fig.colorbar(image, ax=ax, fraction=0.040, pad=0.02, ticks=np.arange(n_categories))
    colorbar.ax.set_yticklabels(labels)
    colorbar.ax.tick_params(length=0, labelsize=11)
    colorbar.set_label(cbar_label, rotation=270, labelpad=18, fontsize=14)
    finish(fig, out_png)


def money_fmt(value, _pos=None):
    if value == 0:
        return "$0"
    if abs(value) >= 1e6:
        return f"${value / 1e6:.0f}M"
    return f"${value / 1e3:.0f}k"


def draw_zone_map(polygons, values, counties, segments, extent, out_png):
    """Zone choropleth on a full-range viridis scale; returns the colorbar maximum."""
    layer = polygons.copy()
    layer["val"] = values.reindex(layer["unit"]).fillna(0.0).to_numpy()
    vmax = float(layer["val"].max())
    norm = Normalize(vmin=0.0, vmax=vmax)
    fig, ax = plt.subplots(figsize=FIGSIZE)
    layer.plot(ax=ax, column="val", cmap="viridis", norm=norm, edgecolor="white", linewidth=0.15, zorder=2)
    if counties is not None and len(counties):
        counties.boundary.plot(ax=ax, color="white", linewidth=0.9, zorder=3)
    ax.add_collection(LineCollection(segments, colors="0.45", linewidths=1.0, zorder=4, capstyle="projecting"))
    frame_axes(ax, extent)
    colorbar = fig.colorbar(ScalarMappable(norm=norm, cmap="viridis"), ax=ax, fraction=0.046, pad=0.02)
    colorbar.ax.yaxis.set_major_formatter(plt.FuncFormatter(money_fmt))
    colorbar.outline.set_visible(False)
    finish(fig, out_png)
    return vmax


def modal_grids(x, y, archetypes, configurations):
    """Per-cell modal M and C on the 1 km grid; ties go to the smallest value."""
    ix, iy = calculate_cell_index(x, y, method="nearest")
    assert ix.min() >= 0 and iy.min() >= 0, "a building falls left of or below the grid origin"
    nx = int(ix.max()) + 1
    ny = int(iy.max()) + 1
    flat = iy * nx + ix
    n_cells = nx * ny

    def mode_of(values, n_values):
        counts = np.bincount(flat * (n_values + 1) + values, minlength=n_cells * (n_values + 1)).reshape(n_cells, n_values + 1)
        best = counts.argmax(axis=1)
        top = counts.max(axis=1)
        total = counts.sum(axis=1)
        tied = ((counts == top[:, None]).sum(axis=1) > 1) & (total > 0)
        return best, total, tied

    mode_m, total, tied_m = mode_of(archetypes, 8)
    mode_c, _, tied_c = mode_of(configurations, 192)
    shape = (ny, nx)
    return dict(nx=nx, ny=ny, mode_m=mode_m.reshape(shape), mode_c=mode_c.reshape(shape),
                count=total.reshape(shape), occupied=(total > 0).reshape(shape),
                tied_m=int(tied_m.sum()), tied_c=int(tied_c.sum()))


def make_zone_grid(origin, extent_cells, cell, n_x=3, n_y=2):
    """Synthetic LEAPHI-style zone grid: n_x * n_y fine zones nesting pairwise into coarse zones."""
    x0, y0 = origin
    nx, ny = extent_cells
    xmin, ymin = x0 - cell / 2.0, y0 - cell / 2.0
    width, height = nx * cell, ny * cell
    rows = []
    for j in range(n_y):
        for i in range(n_x):
            k = j * n_x + i
            rows.append({
                "COUNTYFP00": f"{101 + k:03d}", "zone": 1.0, "i_lh503": float(k // 2 + 1),
                "i_lh1509": float(k + 1), "i_lh708": float(k + 1), "i_lh740": float(k + 1),
                "geometry": box(xmin + i * width / n_x, ymin + j * height / n_y,
                                xmin + (i + 1) * width / n_x, ymin + (j + 1) * height / n_y),
            })
    return gpd.GeoDataFrame(rows, crs=GRID_CRS)


def run_exposure(inventory_path, hazard_dir, zones_path, seed):
    """Run the exposure chain with explicit calls; every result is assigned here."""
    exposure = Exposure(inventory_path=str(inventory_path), hazard_dir=str(hazard_dir), zones_path=str(zones_path),
                        average_mode="Nearest", 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)
    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

Inputs¶

Smoke mode builds a synthetic window (inventory spread over a 12 x 10 km grid, five hazard events, six zones nesting into three), runs the exposure chain and generates a damage library for the (M, C) pairs that occur. Full mode reads the frozen Eastern North Carolina export and the hazard catalogue.

In [5]:

Copied!





data_available = True
exposure = None
library_path = None
lookup_paths = None
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"),
}
if smoke_mode:
    case_root = output_path / "synthetic_case"
    probabilities = SyntheticExposureData.make_hurricane_probabilities(SMOKE_N_EVENTS, seed=SMOKE_SEED)
    zones_gdf = make_zone_grid(GRID_ORIGIN, SMOKE_GRID, GRID_CELL_SIZE)
    written = SyntheticExposureData.save_inputs(
        str(case_root),
        inventory=SyntheticExposureData.make_inventory(
            SMOKE_N_BUILDINGS, seed=SMOKE_SEED, extent_cells=SMOKE_GRID, include_edge_cases=False, spread_full_window=True
        ),
        hazard_events=SyntheticExposureData.make_hazard_events(
            tuple(f"h{index}" for index in range(1, SMOKE_N_EVENTS + 1)), nx=SMOKE_GRID[0], ny=SMOKE_GRID[1], seed=SMOKE_SEED
        ),
        zones=zones_gdf,
        hurricane_probabilities=probabilities,
        **canonical_tables,
    )
    exposure = run_exposure(written["inventory"], written["hazard_dir"], written["zones"], SMOKE_SEED)
    building_table, wind_ids, flood_ids, event_names = exposure.calculate_starr_export()
    hazard_points_xy, hazard_wind, hazard_flood = exposure.hazard_points_xy, exposure.hazard_wind, exposure.hazard_flood
    pairs = sorted(set(zip(building_table["M"].tolist(), building_table["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 window: {len(building_table)} buildings, {hazard_wind.shape[1]} control points, {len(event_names)} events, {len(pairs)} (M, C) pairs")
else:
    hazard_dir = full_data_root / FULL_HAZARD_DIR.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)
    lookup_dir = full_data_root / FULL_LOOKUP_DIR.relative_to(FULL_DATA_ROOT)
    candidate_library = full_data_root / FULL_LIBRARY.relative_to(FULL_DATA_ROOT)
    missing = [str(p) for p in [hazard_dir, zones_path, *export_paths.values()] if not p.exists()]
    if (lookup_dir / "L_mcw_wd.npy").exists() and (lookup_dir / "L_mcf_fld.npy").exists():
        lookup_paths = (str(lookup_dir / "L_mcw_wd.npy"), str(lookup_dir / "L_mcf_fld.npy"))
    elif candidate_library.exists():
        library_path = candidate_library
    else:
        missing.append(f"{lookup_dir} or {candidate_library}")
    if missing:
        data_available = False
        print("skipped: missing", missing)
    else:
        started = time.perf_counter()
        catalogue = BaseExposure(hazard_dir=str(hazard_dir))
        catalogue.load_hazard_events()
        hazard_points_xy, hazard_wind, hazard_flood = catalogue.hazard_points_xy, catalogue.hazard_wind, catalogue.hazard_flood
        event_names = catalogue.event_names
        building_table = pd.read_csv(
            export_paths["inventory_csv"],
            usecols=["building_ID", "M", "C", "structure_value", "contents_value", "home_value", "x", "y", ZONE_ID],
        )
        wind_ids = np.load(export_paths["wind_ids"], mmap_mode="r")
        flood_ids = np.load(export_paths["flood_ids"], mmap_mode="r")
        zones_gdf = gpd.read_file(zones_path).to_crs(GRID_CRS)
        probabilities = np.asarray(Utils.from_json(str(PROB_JSON), dict_key="Hurricane_Prob"))
        assert wind_ids.shape == (len(building_table), len(probabilities))
        print(f"frozen export: {len(building_table):,} buildings x {len(event_names)} events; "
              f"{hazard_wind.shape[1]:,} control points; loaded in {time.perf_counter() - started:.1f} s")
if data_available:
    print(f"sum of event probabilities = {probabilities[:, 1].sum():.6f} per year")
data_available = True
exposure = None
library_path = None
lookup_paths = None
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"),
}
if smoke_mode:
    case_root = output_path / "synthetic_case"
    probabilities = SyntheticExposureData.make_hurricane_probabilities(SMOKE_N_EVENTS, seed=SMOKE_SEED)
    zones_gdf = make_zone_grid(GRID_ORIGIN, SMOKE_GRID, GRID_CELL_SIZE)
    written = SyntheticExposureData.save_inputs(
        str(case_root),
        inventory=SyntheticExposureData.make_inventory(
            SMOKE_N_BUILDINGS, seed=SMOKE_SEED, extent_cells=SMOKE_GRID, include_edge_cases=False, spread_full_window=True
        ),
        hazard_events=SyntheticExposureData.make_hazard_events(
            tuple(f"h{index}" for index in range(1, SMOKE_N_EVENTS + 1)), nx=SMOKE_GRID[0], ny=SMOKE_GRID[1], seed=SMOKE_SEED
        ),
        zones=zones_gdf,
        hurricane_probabilities=probabilities,
        **canonical_tables,
    )
    exposure = run_exposure(written["inventory"], written["hazard_dir"], written["zones"], SMOKE_SEED)
    building_table, wind_ids, flood_ids, event_names = exposure.calculate_starr_export()
    hazard_points_xy, hazard_wind, hazard_flood = exposure.hazard_points_xy, exposure.hazard_wind, exposure.hazard_flood
    pairs = sorted(set(zip(building_table["M"].tolist(), building_table["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 window: {len(building_table)} buildings, {hazard_wind.shape[1]} control points, {len(event_names)} events, {len(pairs)} (M, C) pairs")
else:
    hazard_dir = full_data_root / FULL_HAZARD_DIR.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)
    lookup_dir = full_data_root / FULL_LOOKUP_DIR.relative_to(FULL_DATA_ROOT)
    candidate_library = full_data_root / FULL_LIBRARY.relative_to(FULL_DATA_ROOT)
    missing = [str(p) for p in [hazard_dir, zones_path, *export_paths.values()] if not p.exists()]
    if (lookup_dir / "L_mcw_wd.npy").exists() and (lookup_dir / "L_mcf_fld.npy").exists():
        lookup_paths = (str(lookup_dir / "L_mcw_wd.npy"), str(lookup_dir / "L_mcf_fld.npy"))
    elif candidate_library.exists():
        library_path = candidate_library
    else:
        missing.append(f"{lookup_dir} or {candidate_library}")
    if missing:
        data_available = False
        print("skipped: missing", missing)
    else:
        started = time.perf_counter()
        catalogue = BaseExposure(hazard_dir=str(hazard_dir))
        catalogue.load_hazard_events()
        hazard_points_xy, hazard_wind, hazard_flood = catalogue.hazard_points_xy, catalogue.hazard_wind, catalogue.hazard_flood
        event_names = catalogue.event_names
        building_table = pd.read_csv(
            export_paths["inventory_csv"],
            usecols=["building_ID", "M", "C", "structure_value", "contents_value", "home_value", "x", "y", ZONE_ID],
        )
        wind_ids = np.load(export_paths["wind_ids"], mmap_mode="r")
        flood_ids = np.load(export_paths["flood_ids"], mmap_mode="r")
        zones_gdf = gpd.read_file(zones_path).to_crs(GRID_CRS)
        probabilities = np.asarray(Utils.from_json(str(PROB_JSON), dict_key="Hurricane_Prob"))
        assert wind_ids.shape == (len(building_table), len(probabilities))
        print(f"frozen export: {len(building_table):,} buildings x {len(event_names)} events; "
              f"{hazard_wind.shape[1]:,} control points; loaded in {time.perf_counter() - started:.1f} s")
if data_available:
    print(f"sum of event probabilities = {probabilities[:, 1].sum():.6f} per year")

frozen export: 1,154,657 buildings x 97 events; 63,094 control points; loaded in 3.3 s
sum of event probabilities = 0.323470 per year

Buildings are binned onto the 1 km grid by their nearest control point (calculate_cell_index, the same rule as nearest-neighbour hazard mapping). Each cell shows its most common M (tab10) and its most common C among the ten most common modal values plus "Other" (tab20); ties go to the smallest value.

In [6]:

Copied!





if data_available:
    mask, extent = calculate_grid_mask(hazard_points_xy[:, 0], hazard_points_xy[:, 1])
    segments = mask_boundary_segments(mask, extent, GRID_CELL_SIZE)
    area = study_area(mask, extent, GRID_CELL_SIZE)
    counties = load_counties(COUNTY_BOUNDARIES, extent, GRID_CRS, area)
    grid = modal_grids(
        building_table["x"].to_numpy("float64"), building_table["y"].to_numpy("float64"),
        building_table["M"].to_numpy("int64"), building_table["C"].to_numpy("int64"),
    )
    occupied = grid["occupied"]
    cell_extent = raster_extent(grid["nx"], grid["ny"])
    counts_m = collections.Counter(grid["mode_m"][occupied].tolist())
    counts_c = collections.Counter(grid["mode_c"][occupied].tolist())
    top10 = sorted(value for value, _ in counts_c.most_common(10))
    print(f"{int(occupied.sum()):,} occupied cells; modal-M cells {dict(sorted(counts_m.items()))}; modal-C top-10 {top10}")
    if not smoke_mode:
        assert top10 == MODAL_C_LEGEND, f"modal-C top-10 {top10} differs from the published legend {MODAL_C_LEGEND}"
        print("modal-C legend matches the published figure")
    cmap_m = plt.colormaps["tab10"].resampled(8)
    categorical_map(grid["mode_m"] - 1, occupied, cell_extent, mask, extent, segments, counties, cmap_m,
                    [str(v) for v in range(1, 9)], "Mode M",
                    figures_path / "eastern_nc_inventory_modal_building_category.png")
    lut = np.full(193, len(top10), dtype="int64")
    for rank, value in enumerate(top10):
        lut[value] = rank
    cmap_c = plt.colormaps["tab20"].resampled(len(top10) + 1)
    categorical_map(lut[grid["mode_c"]], occupied, cell_extent, mask, extent, segments, counties, cmap_c,
                    [str(v) for v in top10] + ["Other"], "Mode C",
                    figures_path / "eastern_nc_inventory_modal_resistance_configuration.png")
if data_available:
    mask, extent = calculate_grid_mask(hazard_points_xy[:, 0], hazard_points_xy[:, 1])
    segments = mask_boundary_segments(mask, extent, GRID_CELL_SIZE)
    area = study_area(mask, extent, GRID_CELL_SIZE)
    counties = load_counties(COUNTY_BOUNDARIES, extent, GRID_CRS, area)
    grid = modal_grids(
        building_table["x"].to_numpy("float64"), building_table["y"].to_numpy("float64"),
        building_table["M"].to_numpy("int64"), building_table["C"].to_numpy("int64"),
    )
    occupied = grid["occupied"]
    cell_extent = raster_extent(grid["nx"], grid["ny"])
    counts_m = collections.Counter(grid["mode_m"][occupied].tolist())
    counts_c = collections.Counter(grid["mode_c"][occupied].tolist())
    top10 = sorted(value for value, _ in counts_c.most_common(10))
    print(f"{int(occupied.sum()):,} occupied cells; modal-M cells {dict(sorted(counts_m.items()))}; modal-C top-10 {top10}")
    if not smoke_mode:
        assert top10 == MODAL_C_LEGEND, f"modal-C top-10 {top10} differs from the published legend {MODAL_C_LEGEND}"
        print("modal-C legend matches the published figure")
    cmap_m = plt.colormaps["tab10"].resampled(8)
    categorical_map(grid["mode_m"] - 1, occupied, cell_extent, mask, extent, segments, counties, cmap_m,
                    [str(v) for v in range(1, 9)], "Mode M",
                    figures_path / "eastern_nc_inventory_modal_building_category.png")
    lut = np.full(193, len(top10), dtype="int64")
    for rank, value in enumerate(top10):
        lut[value] = rank
    cmap_c = plt.colormaps["tab20"].resampled(len(top10) + 1)
    categorical_map(lut[grid["mode_c"]], occupied, cell_extent, mask, extent, segments, counties, cmap_c,
                    [str(v) for v in top10] + ["Other"], "Mode C",
                    figures_path / "eastern_nc_inventory_modal_resistance_configuration.png")

37,315 occupied cells; modal-M cells {1: 162, 2: 1917, 3: 1704, 4: 31006, 5: 39, 6: 348, 7: 272, 8: 1867}; modal-C top-10 [1, 3, 5, 7, 25, 49, 97, 145, 169, 175]
modal-C legend matches the published figure

No description has been provided for this image

Figures 3-4: probability-weighted hazard¶

Per control point, the wind map shows sum_h p_h V_h / sum_h p_h over all events and the flood map sum_h p_h z_h / sum_h p_h over the events that inundate the point (cells never inundated stay grey). Colour limits are the 2nd and 98th percentiles, as in the published figures.

In [7]:

Copied!





if data_available:
    annual = probabilities[:, 1]
    wind_mean = calculate_probability_weighted_intensity(hazard_wind.T, annual)
    flood_mean = calculate_probability_weighted_intensity(hazard_flood.T, annual, conditional_on_positive=True)
    ix, iy = calculate_cell_index(hazard_points_xy[:, 0], hazard_points_xy[:, 1], method="nearest")
    ny_mask, nx_mask = mask.shape
    wet = ~np.isnan(flood_mean)
    wind_limits = (float(np.percentile(wind_mean, 2)), float(np.percentile(wind_mean, 98)))
    flood_limits = (float(np.percentile(flood_mean[wet], 2)), float(np.percentile(flood_mean[wet], 98))) if wet.any() else (None, None)
    print(f"mean wind {wind_mean.min():.2f}-{wind_mean.max():.2f} m/s; conditional flood over {int(wet.sum())} of {len(wet)} cells "
          f"({flood_mean[wet].min() if wet.any() else float('nan'):.3f}-{flood_mean[wet].max() if wet.any() else float('nan'):.3f} m)")
    draw_raster_map(to_raster(wind_mean, ix, iy, nx_mask, ny_mask), extent, mask, counties, segments, "Wind speed (m/s)",
                    figures_path / "eastern_nc_hazard_probability_weighted_wind_speed.png",
                    vmin=wind_limits[0], vmax=wind_limits[1], with_base=False)
    draw_raster_map(to_raster(flood_mean, ix, iy, nx_mask, ny_mask), extent, mask, counties, segments, "Flood depth (m)",
                    figures_path / "eastern_nc_hazard_probability_weighted_flood_depth.png",
                    vmin=flood_limits[0], vmax=flood_limits[1], with_base=True)
if data_available:
    annual = probabilities[:, 1]
    wind_mean = calculate_probability_weighted_intensity(hazard_wind.T, annual)
    flood_mean = calculate_probability_weighted_intensity(hazard_flood.T, annual, conditional_on_positive=True)
    ix, iy = calculate_cell_index(hazard_points_xy[:, 0], hazard_points_xy[:, 1], method="nearest")
    ny_mask, nx_mask = mask.shape
    wet = ~np.isnan(flood_mean)
    wind_limits = (float(np.percentile(wind_mean, 2)), float(np.percentile(wind_mean, 98)))
    flood_limits = (float(np.percentile(flood_mean[wet], 2)), float(np.percentile(flood_mean[wet], 98))) if wet.any() else (None, None)
    print(f"mean wind {wind_mean.min():.2f}-{wind_mean.max():.2f} m/s; conditional flood over {int(wet.sum())} of {len(wet)} cells "
          f"({flood_mean[wet].min() if wet.any() else float('nan'):.3f}-{flood_mean[wet].max() if wet.any() else float('nan'):.3f} m)")
    draw_raster_map(to_raster(wind_mean, ix, iy, nx_mask, ny_mask), extent, mask, counties, segments, "Wind speed (m/s)",
                    figures_path / "eastern_nc_hazard_probability_weighted_wind_speed.png",
                    vmin=wind_limits[0], vmax=wind_limits[1], with_base=False)
    draw_raster_map(to_raster(flood_mean, ix, iy, nx_mask, ny_mask), extent, mask, counties, segments, "Flood depth (m)",
                    figures_path / "eastern_nc_hazard_probability_weighted_flood_depth.png",
                    vmin=flood_limits[0], vmax=flood_limits[1], with_base=True)

mean wind 15.54-28.31 m/s; conditional flood over 10270 of 63094 cells (0.004-6.858 m)

Figures 5-7: expected annual loss per zone¶

LossEstimation maps the hazard IDs onto the lookup tables (floor), adds wind and flood, weights by the annual event probabilities and multiplies by the structure value; the building EAL is summed within LEAPHI-503 zones (i_lh740 -> i_lh503 through the zone layer). Each panel uses its own full-range colour scale; full mode checks the colorbar maxima against the published figures.

In [8]:

Copied!





if data_available:
    try:
        started = time.perf_counter()
        estimator = LossEstimation(
            prob_ds_file_path=str(library_path) if library_path else None,
            hurricane_probabilities=probabilities,
            zones_path=str(zones_path),
            zone_id_col=ZONE_ID,
            aggregation_col=AGG_UNIT,
            index_mode=HAZARD_INDEX,
            combine_mode=COMBINE_MODE,
            value_basis=VALUE_FIELD,
        )
        estimator.load_hurricane_probabilities()
        estimator.load_zones()
        if lookup_paths is not None:
            estimator.load_lookup_tables(*lookup_paths)
            print("lookup tables: frozen STARR_lookup_tables")
        else:
            estimator.load_damage_library()
            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("lookup tables: built from the damage library")
        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()
        annual = estimator.hurricane_prob[:, 1]
        estimator.building_eal = estimator.calculate_building_eal(estimator.building_table, estimator.wind_ids, estimator.flood_ids, annual)
        estimator.aggregation_key = estimator.calculate_aggregation_key(estimator.building_table, estimator.zones, ZONE_ID, AGG_UNIT)
        estimator.zone_eal = estimator.calculate_zone_eal(estimator.building_eal, estimator.aggregation_key)
        assigned = estimator.zone_eal.drop(index=-1, errors="ignore")
        print(f"EAL computed for {len(estimator.building_eal):,} buildings in {len(assigned)} {AGG_UNIT} units ({time.perf_counter() - started:.1f} s)")

        polygons = estimator.zones.to_crs(GRID_CRS).dissolve(by=AGG_UNIT).reset_index().rename(columns={AGG_UNIT: "unit"})
        polygons["unit"] = polygons["unit"].astype("int64")
        polygons = clip_to_study(polygons[["unit", "geometry"]], area)
        eal_figures = {
            "wind": ("eal_wind", "eastern_nc_expected_annual_loss_wind.png"),
            "flood": ("eal_flood", "eastern_nc_expected_annual_loss_flood.png"),
            "combined": ("eal_combined", "eastern_nc_expected_annual_loss_combined.png"),
        }
        reproduced_vmax = {}
        for panel, (column, file_name) in eal_figures.items():
            reproduced_vmax[panel] = draw_zone_map(polygons, assigned[column], counties, segments, extent,
                                                   figures_path / file_name)
        comparison = pd.DataFrame({"reproduced": reproduced_vmax, "published": PUBLISHED_VMAX})
        comparison["relative_difference"] = comparison["reproduced"] / comparison["published"] - 1.0
        display(comparison)
        if not smoke_mode:
            worst = comparison["relative_difference"].abs().max()
            assert worst <= VMAX_TOLERANCE, f"colorbar maxima differ from the published figures by up to {worst:.2%}"
            print(f"all three colorbar maxima within {worst:.2%} of the published figures")

        # Sensitivity of the colorbar maxima to the index convention and the value basis
        sensitivity = {}
        for basis in ("structure_value", "home_value"):
            for mode in ("floor", "ceil", "round", "linear"):
                estimator.value_basis, estimator.index_mode = basis, mode
                eal_variant = estimator.calculate_building_eal(estimator.building_table, estimator.wind_ids, estimator.flood_ids, annual)
                zone_variant = estimator.calculate_zone_eal(eal_variant, estimator.aggregation_key).drop(index=-1, errors="ignore")
                sensitivity[f"{basis}|{mode}"] = {panel: zone_variant[column].max() for panel, column in
                                                  (("wind", "eal_wind"), ("flood", "eal_flood"), ("combined", "eal_combined"))}
        estimator.value_basis, estimator.index_mode = VALUE_FIELD, HAZARD_INDEX
        sensitivity_table = pd.DataFrame(sensitivity).T
        reference = sensitivity_table.loc[f"{VALUE_FIELD}|{HAZARD_INDEX}"]
        display((sensitivity_table / reference - 1.0).map("{:+.1%}".format))
    except BaseException:
        restore_state()
        raise
if data_available:
    try:
        started = time.perf_counter()
        estimator = LossEstimation(
            prob_ds_file_path=str(library_path) if library_path else None,
            hurricane_probabilities=probabilities,
            zones_path=str(zones_path),
            zone_id_col=ZONE_ID,
            aggregation_col=AGG_UNIT,
            index_mode=HAZARD_INDEX,
            combine_mode=COMBINE_MODE,
            value_basis=VALUE_FIELD,
        )
        estimator.load_hurricane_probabilities()
        estimator.load_zones()
        if lookup_paths is not None:
            estimator.load_lookup_tables(*lookup_paths)
            print("lookup tables: frozen STARR_lookup_tables")
        else:
            estimator.load_damage_library()
            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("lookup tables: built from the damage library")
        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()
        annual = estimator.hurricane_prob[:, 1]
        estimator.building_eal = estimator.calculate_building_eal(estimator.building_table, estimator.wind_ids, estimator.flood_ids, annual)
        estimator.aggregation_key = estimator.calculate_aggregation_key(estimator.building_table, estimator.zones, ZONE_ID, AGG_UNIT)
        estimator.zone_eal = estimator.calculate_zone_eal(estimator.building_eal, estimator.aggregation_key)
        assigned = estimator.zone_eal.drop(index=-1, errors="ignore")
        print(f"EAL computed for {len(estimator.building_eal):,} buildings in {len(assigned)} {AGG_UNIT} units ({time.perf_counter() - started:.1f} s)")

        polygons = estimator.zones.to_crs(GRID_CRS).dissolve(by=AGG_UNIT).reset_index().rename(columns={AGG_UNIT: "unit"})
        polygons["unit"] = polygons["unit"].astype("int64")
        polygons = clip_to_study(polygons[["unit", "geometry"]], area)
        eal_figures = {
            "wind": ("eal_wind", "eastern_nc_expected_annual_loss_wind.png"),
            "flood": ("eal_flood", "eastern_nc_expected_annual_loss_flood.png"),
            "combined": ("eal_combined", "eastern_nc_expected_annual_loss_combined.png"),
        }
        reproduced_vmax = {}
        for panel, (column, file_name) in eal_figures.items():
            reproduced_vmax[panel] = draw_zone_map(polygons, assigned[column], counties, segments, extent,
                                                   figures_path / file_name)
        comparison = pd.DataFrame({"reproduced": reproduced_vmax, "published": PUBLISHED_VMAX})
        comparison["relative_difference"] = comparison["reproduced"] / comparison["published"] - 1.0
        display(comparison)
        if not smoke_mode:
            worst = comparison["relative_difference"].abs().max()
            assert worst <= VMAX_TOLERANCE, f"colorbar maxima differ from the published figures by up to {worst:.2%}"
            print(f"all three colorbar maxima within {worst:.2%} of the published figures")

        # Sensitivity of the colorbar maxima to the index convention and the value basis
        sensitivity = {}
        for basis in ("structure_value", "home_value"):
            for mode in ("floor", "ceil", "round", "linear"):
                estimator.value_basis, estimator.index_mode = basis, mode
                eal_variant = estimator.calculate_building_eal(estimator.building_table, estimator.wind_ids, estimator.flood_ids, annual)
                zone_variant = estimator.calculate_zone_eal(eal_variant, estimator.aggregation_key).drop(index=-1, errors="ignore")
                sensitivity[f"{basis}|{mode}"] = {panel: zone_variant[column].max() for panel, column in
                                                  (("wind", "eal_wind"), ("flood", "eal_flood"), ("combined", "eal_combined"))}
        estimator.value_basis, estimator.index_mode = VALUE_FIELD, HAZARD_INDEX
        sensitivity_table = pd.DataFrame(sensitivity).T
        reference = sensitivity_table.loc[f"{VALUE_FIELD}|{HAZARD_INDEX}"]
        display((sensitivity_table / reference - 1.0).map("{:+.1%}".format))
    except BaseException:
        restore_state()
        raise

lookup tables: frozen STARR_lookup_tables

EAL computed for 1,154,657 buildings in 502 i_lh503 units (3.5 s)

	reproduced	published	relative_difference
wind	6.492232e+06	6487000.0	0.000806
flood	6.888614e+07	68784000.0	0.001485
combined	7.190775e+07	71806000.0	0.001417

all three colorbar maxima within 0.15% of the published figures

	wind	flood	combined
structure_value\|floor	+0.0%	+0.0%	+0.0%
structure_value\|ceil	+40.6%	+2.6%	+4.5%
structure_value\|round	+21.1%	+1.4%	+2.7%
structure_value\|linear	+20.1%	+1.4%	+2.5%
home_value\|floor	+50.0%	+50.0%	+50.0%
home_value\|ceil	+110.9%	+53.9%	+56.8%
home_value\|round	+81.6%	+52.2%	+54.1%
home_value\|linear	+80.2%	+52.1%	+53.8%

Numeric checks¶

The quantities behind the figures are written next to them so a rerun can be compared without opening the PNGs.

In [9]:

Copied!





if data_available:
    checks = {
        "profile": "smoke" if smoke_mode else "full",
        "inventory": {
            "n_buildings": int(len(building_table)),
            "n_occupied_cells": int(occupied.sum()),
            "modal_M_cell_counts": {str(k): int(v) for k, v in sorted(counts_m.items())},
            "modal_C_top10": [int(v) for v in top10],
            "published_legend": MODAL_C_LEGEND,
            "tied_cells_modal_M": grid["tied_m"],
            "tied_cells_modal_C": grid["tied_c"],
        },
        "hazard": {
            "n_cells": int(hazard_points_xy.shape[0]),
            "n_events": int(len(probabilities)),
            "sum_p_h": float(probabilities[:, 1].sum()),
            "wind_mean_min_max": [float(wind_mean.min()), float(wind_mean.max())],
            "n_ever_inundated": int(wet.sum()),
            "wind_colour_limits_p2_p98": list(wind_limits),
            "flood_colour_limits_p2_p98": [None if v is None else float(v) for v in flood_limits],
        },
        "eal": {
            "aggregation_unit": AGG_UNIT,
            "value_field": VALUE_FIELD,
            "index_mode": HAZARD_INDEX,
            "combine_mode": COMBINE_MODE,
            "n_units_with_buildings": int((assigned["n_buildings"] > 0).sum()),
            "reproduced_vmax": {k: float(v) for k, v in reproduced_vmax.items()},
            "published_vmax": PUBLISHED_VMAX,
            "argmax_units": {panel: int(assigned[column].idxmax()) for panel, column in
                             (("wind", "eal_wind"), ("flood", "eal_flood"), ("combined", "eal_combined"))},
            "totals": {column: float(assigned[column].sum()) for column in ("eal_wind", "eal_flood", "eal_combined")},
            "sensitivity_vmax": {k: {p: float(v) for p, v in row.items()} for k, row in sensitivity.items()},
        },
    }
    for name in ("inventory", "hazard", "eal"):
        with open(checks_path / f"{name}.json", "w", encoding="utf-8") as handle:
            json.dump(checks[name], handle, indent=2)
    print(f"checks written to {checks_path}")
    print(json.dumps(checks["eal"]["reproduced_vmax"], indent=2))
if data_available:
    checks = {
        "profile": "smoke" if smoke_mode else "full",
        "inventory": {
            "n_buildings": int(len(building_table)),
            "n_occupied_cells": int(occupied.sum()),
            "modal_M_cell_counts": {str(k): int(v) for k, v in sorted(counts_m.items())},
            "modal_C_top10": [int(v) for v in top10],
            "published_legend": MODAL_C_LEGEND,
            "tied_cells_modal_M": grid["tied_m"],
            "tied_cells_modal_C": grid["tied_c"],
        },
        "hazard": {
            "n_cells": int(hazard_points_xy.shape[0]),
            "n_events": int(len(probabilities)),
            "sum_p_h": float(probabilities[:, 1].sum()),
            "wind_mean_min_max": [float(wind_mean.min()), float(wind_mean.max())],
            "n_ever_inundated": int(wet.sum()),
            "wind_colour_limits_p2_p98": list(wind_limits),
            "flood_colour_limits_p2_p98": [None if v is None else float(v) for v in flood_limits],
        },
        "eal": {
            "aggregation_unit": AGG_UNIT,
            "value_field": VALUE_FIELD,
            "index_mode": HAZARD_INDEX,
            "combine_mode": COMBINE_MODE,
            "n_units_with_buildings": int((assigned["n_buildings"] > 0).sum()),
            "reproduced_vmax": {k: float(v) for k, v in reproduced_vmax.items()},
            "published_vmax": PUBLISHED_VMAX,
            "argmax_units": {panel: int(assigned[column].idxmax()) for panel, column in
                             (("wind", "eal_wind"), ("flood", "eal_flood"), ("combined", "eal_combined"))},
            "totals": {column: float(assigned[column].sum()) for column in ("eal_wind", "eal_flood", "eal_combined")},
            "sensitivity_vmax": {k: {p: float(v) for p, v in row.items()} for k, row in sensitivity.items()},
        },
    }
    for name in ("inventory", "hazard", "eal"):
        with open(checks_path / f"{name}.json", "w", encoding="utf-8") as handle:
            json.dump(checks[name], handle, indent=2)
    print(f"checks written to {checks_path}")
    print(json.dumps(checks["eal"]["reproduced_vmax"], indent=2))

checks written to /home/mammad/Documents/EDU/cope/CHEERSafe/LossModel-exposure/Outputs/eastern_nc_loss_figures/checks
{
  "wind": 6492231.506356327,
  "flood": 68886135.98809117,
  "combined": 71907750.37074493
}

In [10]:

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¶

Figures are written to Outputs/eastern_nc_loss_figures/figures/ and the numbers behind them to checks/*.json. In full mode the modal-C legend and the three EAL colorbar maxima match the published figures within 0.5%.
The EAL maps aggregate to LEAPHI-503 zones with structure_value and floor indexing; the paper caption says "census tract", but the published polygons are the 503-zone system. The sensitivity table shows how far the other index conventions and home_value move the colorbar maxima.
Full mode starts from the frozen 07/07/2026 export because the source's configuration draw was unseeded; a seeded rerun of the exposure chain reproduces the archetype distributions, not the row-by-row C values. Smoke mode runs the whole chain on a synthetic window instead.
County boundaries need the TIGER 2020 county shapefile (CHEERSAFE_FIGURES_COUNTIES); without it the maps show the 1 km study-area outline only.