diff --git a/src/hazus/flood/aal.py b/src/hazus/flood/aal.py index cb0a1f9..167cbc3 100644 --- a/src/hazus/flood/aal.py +++ b/src/hazus/flood/aal.py @@ -14,10 +14,12 @@ from __future__ import annotations +import math from decimal import ROUND_HALF_UP, Decimal, getcontext from fractions import Fraction from typing import Iterable, Mapping +import numpy as np import pandas as pd getcontext().prec = 38 @@ -48,34 +50,49 @@ def aal_coefficients_decimal(return_periods: Iterable[int]) -> dict[int, Decimal return coefs -def _round_half_away(x: Fraction) -> int: - n = abs(x) - r = int(n + Fraction(1, 2)) - return r if x >= 0 else -r - - def aal_by_soccup(ecloss_soccup: pd.DataFrame, coefficients: Mapping[int, Decimal]) -> pd.DataFrame: - """``flFRGBSEcLossBySOccup`` rows with ``ReturnPeriodId = 'Annual'``.""" - df = ecloss_soccup.copy() - rp = pd.to_numeric(df["ReturnPeriodId"].astype(str).str.strip(), errors="coerce") - df = df[rp.isin(list(coefficients))].assign(ReturnPeriod=rp.astype("Int64")) + """``flFRGBSEcLossBySOccup`` rows with ``ReturnPeriodId = 'Annual'``. + + The sum is exact: with coefficients ``C / L`` (``L`` their common denominator), ``2 L × AAL`` is the + integer ``Σ 2 C v`` (longest period) ``+ Σ C (v + v_longer)``, rounded half away from zero by integer + division. Python integers take over if int64 could overflow. + """ + rp = pd.to_numeric(ecloss_soccup["ReturnPeriodId"].astype(str).str.strip(), errors="coerce") + df = ecloss_soccup[rp.isin(list(coefficients))].assign(ReturnPeriod=rp.astype("Int64")) part_keys = ["StudyCaseId", "CensusBlock", "ControlHazard", "SOccup", "Pre_PostFirm", "AnalysisOptId"] + columns = part_keys + ["ReturnPeriodId"] + ECON_COLS df = df.groupby(part_keys + ["ReturnPeriod"], as_index=False)[ECON_COLS].sum() - coef = {int(k): Fraction(v) for k, v in coefficients.items()} - - rows = [] - for key, g in df.sort_values("ReturnPeriod", ascending=False).groupby(part_keys, sort=False): - totals = dict.fromkeys(ECON_COLS, Fraction(0)) - prev = None - for r in g.itertuples(index=False): - c = coef[int(r.ReturnPeriod)] - for col in ECON_COLS: - v = int(getattr(r, col)) - totals[col] += c * v if prev is None else c * (v + int(getattr(prev, col))) / 2 - prev = r - rows.append({**dict(zip(part_keys, key)), "ReturnPeriodId": "Annual", - **{col: _round_half_away(totals[col]) for col in ECON_COLS}}) - return pd.DataFrame(rows, columns=part_keys + ["ReturnPeriodId"] + ECON_COLS) + if df.empty: + return pd.DataFrame(columns=columns) + + # groups in order of first appearance after sorting by return period, longest period first within each + df = df.sort_values("ReturnPeriod", ascending=False) + gid = df.groupby(part_keys, sort=False).ngroup().to_numpy() + periods = df["ReturnPeriod"].to_numpy("int64") + order = np.lexsort((-periods, gid)) + gid, periods = gid[order], periods[order] + values = df[ECON_COLS].to_numpy("int64")[order] + + fractions = {int(k): Fraction(v) for k, v in coefficients.items()} + denominator = math.lcm(*(f.denominator for f in fractions.values())) + scaled = {k: f.numerator * (denominator // f.denominator) for k, f in fractions.items()} + bound = 2 * int(np.abs(values).max(initial=0)) * sum(abs(c) for c in scaled.values()) + denominator + dtype = "int64" if bound < 2 ** 62 else object + values = values.astype(dtype) + coef = np.array([scaled[int(p)] for p in periods], dtype=dtype)[:, None] + + first = np.r_[True, gid[1:] != gid[:-1]] + longer = np.roll(values, 1, axis=0) + longer[first] = values[first] # the longest period counts twice: 2 C v + starts = np.flatnonzero(first) + totals = np.add.reduceat(coef * (values + longer), starts, axis=0) + magnitude = (np.abs(totals) + denominator) // (2 * denominator) + rounded = np.where(totals >= 0, magnitude, -magnitude).astype("int64") + + out = df.iloc[order[starts]][part_keys].reset_index(drop=True).assign(ReturnPeriodId="Annual") + for j, col in enumerate(ECON_COLS): + out[col] = rounded[:, j] + return out[columns] def aal_gbs(ecloss_soccup: pd.DataFrame, coefficients: Mapping[int, Decimal], region: Mapping[str, pd.DataFrame], diff --git a/src/hazus/flood/debris.py b/src/hazus/flood/debris.py index c9ddb8e..2bd4f33 100644 --- a/src/hazus/flood/debris.py +++ b/src/hazus/flood/debris.py @@ -71,7 +71,7 @@ def calculate_debris(region: Mapping[str, pd.DataFrame], area_weighted: pd.DataF scheme["CensusBlock"] = _strip(scheme["CensusBlock"]) demo = region["hzDemographicsB"].copy() demo["CensusBlock"] = _strip(demo["CensusBlock"]) - sc_blocks = aw_t.loc[aw_t["StudyCaseID"] == study_case_id, "CensusBlock"] + sc_blocks = aw_t.loc[aw_t["StudyCaseID"] == study_case_id, "CensusBlock"].unique() ppp = _pre_post_pct(_pre_post_bldg(demo[demo["CensusBlock"].isin(sc_blocks)], scheme), scheme) gbs = aw[(aw["StudyCaseID"] == study_case_id) & (aw["AnalysisOptId"] == analysis_opt_id)].merge(ppp, on="CensusBlock") diff --git a/src/hazus/flood/exposure.py b/src/hazus/flood/exposure.py index a53d702..d23925b 100644 --- a/src/hazus/flood/exposure.py +++ b/src/hazus/flood/exposure.py @@ -94,25 +94,33 @@ def snapped_grid(bounds: tuple[float, float, float, float], snap_transform, cell return Grid(left=ll_x, top=ll_y + height * cell, cell=cell, width=width, height=height) -def _block_counts(blocks, grid: Grid, block_index: np.ndarray) -> np.ndarray: - """Total cell count of each block (the ``BLKGRID`` value table), rasterizing each block on its own window.""" +def _block_counts(blocks, grid: Grid, block_index: np.ndarray, cache: dict | None = None) -> np.ndarray: + """Total cell count of each block (the ``BLKGRID`` value table), rasterizing each block on its own window. + + ``cache`` (block index → count) is reused by later calls on the same blocks and grid. + """ + import rasterio from rasterio.features import rasterize - counts = np.zeros(len(block_index), dtype="int64") + cache = {} if cache is None else cache geoms = blocks.geometry.values - for k, i in enumerate(block_index): - geom = geoms[i] - r0, r1, c0, c1 = grid.window_for_bounds(*geom.bounds) - if r1 <= r0 or c1 <= c0: - continue - burned = rasterize([(geom, 1)], out_shape=(r1 - r0, c1 - c0), transform=grid.sub_transform(r0, c0), - fill=0, dtype="uint8") - counts[k] = int(burned.sum()) - return counts + with rasterio.Env(): # one GDAL environment for all the small rasterize calls + for i in block_index: + if i in cache: + continue + geom = geoms[i] + r0, r1, c0, c1 = grid.window_for_bounds(*geom.bounds) + if r1 <= r0 or c1 <= c0: + cache[i] = 0 + continue + burned = rasterize([(geom, 1)], out_shape=(r1 - r0, c1 - c0), transform=grid.sub_transform(r0, c0), + fill=0, dtype="uint8") + cache[i] = int(burned.sum()) + return np.array([cache[i] for i in block_index], dtype="int64") def area_weighted_depths(blocks, flood: FloodRaster, *, block_id: str = "CensusBlock", boundary=None, - cell: float = CELL_SIZE, tile: int = 2048) -> pd.DataFrame: + cell: float = CELL_SIZE, tile: int = 2048, block_counts: dict | None = None) -> pd.DataFrame: """Percent of each census block's cells in each (zone, half-foot depth bin). Parameters @@ -124,6 +132,9 @@ def area_weighted_depths(blocks, flood: FloodRaster, *, block_id: str = "CensusB Depth grid in feet, EPSG:4326 (see :func:`hazus.flood.hazard.prepare_depth_grid`). boundary: Optional study-region boundary geometry (``hzBoundary``). Defaults to the union of the blocks. + block_counts: + Optional dict to share block cell counts between calls on the same blocks (e.g. the grids of a + scenario). Counts are kept per block grid, so grids snapped differently do not mix. Returns ------- @@ -131,6 +142,7 @@ def area_weighted_depths(blocks, flood: FloodRaster, *, block_id: str = "CensusB (the long form of ``flAnGBSCombineAWTemp``). """ import rasterio + import shapely from rasterio.features import rasterize blocks = blocks.reset_index(drop=True) @@ -138,8 +150,12 @@ def area_weighted_depths(blocks, flood: FloodRaster, *, block_id: str = "CensusB blocks = blocks.to_crs(4326) hazard = HazardType(flood.hazard_type) sindex = blocks.sindex + geoms = blocks.geometry.values + local = np.zeros(len(blocks) + 1, dtype="int64") # block value (index + 1) -> position in the tile's hits + if boundary is not None: + shapely.prepare(boundary) - with rasterio.open(flood.path) as ds: + with rasterio.open(flood.path) as ds, rasterio.Env(): grid = snapped_grid(tuple(blocks.total_bounds), ds.transform, cell) dl, dt = ds.transform.c, ds.transform.f drx, dry = abs(ds.transform.a), abs(ds.transform.e) @@ -155,50 +171,59 @@ def area_weighted_depths(blocks, flood: FloodRaster, *, block_id: str = "CensusB tx0 = grid.left + tc * grid.cell ty0 = grid.top - tr * grid.cell tbox = (tx0, ty0 - h * grid.cell, tx0 + w * grid.cell, ty0) - hits = sindex.query(_box(*tbox)) + tile_box = _box(*tbox) + if boundary is not None and not boundary.intersects(tile_box): + continue # no cell centre of this tile is inside the study region + hits = sindex.query(tile_box) if len(hits) == 0: continue transform = grid.sub_transform(tr, tc) - blk = rasterize(((blocks.geometry.values[i], int(i) + 1) for i in hits), out_shape=(h, w), + blk = rasterize(((geoms[i], int(i) + 1) for i in hits), out_shape=(h, w), transform=transform, fill=0, dtype="int32") if not blk.any(): continue - if boundary is not None: - inside = rasterize([(boundary, 1)], out_shape=(h, w), transform=transform, fill=0, - dtype="uint8").astype(bool) - else: - inside = blk > 0 - # Nearest-neighbour resample: depth cell containing each block-grid cell centre. + # Nearest-neighbour resample: depth cell containing each block-grid cell centre. The rows and + # columns that fall on the depth grid are contiguous (ra:rb, ca:cb); every other cell is NoData. xs = tx0 + (np.arange(w) + 0.5) * grid.cell ys = ty0 - (np.arange(h) + 0.5) * grid.cell dcols = np.floor((xs - dl) / drx).astype("int64") drows = np.floor((dt - ys) / dry).astype("int64") - col_ok = (dcols >= 0) & (dcols < ds.width) - row_ok = (drows >= 0) & (drows < ds.height) - if not col_ok.any() or not row_ok.any(): + ci = np.flatnonzero((dcols >= 0) & (dcols < ds.width)) + ri = np.flatnonzero((drows >= 0) & (drows < ds.height)) + if len(ci) == 0 or len(ri) == 0: continue - wc0, wc1 = dcols[col_ok].min(), dcols[col_ok].max() + 1 - wr0, wr1 = drows[row_ok].min(), drows[row_ok].max() + 1 + ca, cb, ra, rb = ci[0], ci[-1] + 1, ri[0], ri[-1] + 1 + wc0, wc1 = dcols[ca], dcols[cb - 1] + 1 + wr0, wr1 = drows[ra], drows[rb - 1] + 1 band = ds.read(1, window=((wr0, wr1), (wc0, wc1))).astype("float64") if nodata is not None and not np.isnan(nodata): band[band == nodata] = np.nan - depth = np.full((h, w), np.nan) - rr = np.clip(drows - wr0, 0, wr1 - wr0 - 1) - cc = np.clip(dcols - wc0, 0, wc1 - wc0 - 1) - depth_full = band[np.ix_(rr, cc)] - ok2d = np.outer(row_ok, col_ok) - depth[ok2d] = depth_full[ok2d] - + # zone << 8 | bin, or -1 where not counted. Bin whichever is smaller: the depth cells (then + # resample the codes) or the resampled cells. + take = np.ix_(drows[ra:rb] - wr0, dcols[ca:cb] - wc0) + depth = band if band.size < (rb - ra) * (cb - ca) else band[take] bins = depth_bin(depth) - valid = (blk > 0) & inside & (bins >= 0) + code = (zone_code(depth, hazard).astype("int64") << 8) | (bins // 50) + code[bins < 0] = -1 + if depth is band: + code = code[take] + + b = blk[ra:rb, ca:cb] + valid = (b > 0) & (code >= 0) + if boundary is not None and not boundary.contains(tile_box): + inside = rasterize([(boundary, 1)], out_shape=(h, w), transform=transform, fill=0, + dtype="uint8") + valid &= inside[ra:rb, ca:cb].astype(bool) if not valid.any(): continue - d = depth[valid] - zones = zone_code(d, hazard) - key = (blk[valid].astype("int64") << 16) | (zones.astype("int64") << 8) | (bins[valid] // 50) - uniq, counts = np.unique(key, return_counts=True) - parts.append(pd.DataFrame({"key": uniq, "Count": counts})) + # count per (block, zone, bin) with a bincount over the tile's own blocks + local[hits + 1] = np.arange(len(hits)) + width = 5 << 8 # zone codes 1-4, bins 0-48 + counts = np.bincount(local[b[valid]] * width + code[valid], minlength=len(hits) * width) + nz = np.flatnonzero(counts) + key = ((hits[nz // width].astype("int64") + 1) << 16) | (nz % width) + parts.append(pd.DataFrame({"key": key, "Count": counts[nz]})) if not parts: return pd.DataFrame(columns=["CensusBlock", "CoastalZone", "ReClass", "Count", "TotalCount", "PctAffected"]) @@ -210,7 +235,8 @@ def area_weighted_depths(blocks, flood: FloodRaster, *, block_id: str = "CensusB reclass = (keys & 0xFF) * 50 flooded = np.unique(idx) - totals = pd.Series(_block_counts(blocks, grid, flooded), index=flooded) + cache = None if block_counts is None else block_counts.setdefault(grid, {}) + totals = pd.Series(_block_counts(blocks, grid, flooded, cache), index=flooded) out = pd.DataFrame({ "CensusBlock": blocks[block_id].to_numpy()[idx], diff --git a/src/hazus/flood/facility_analysis.py b/src/hazus/flood/facility_analysis.py index b7d6e47..1693742 100644 --- a/src/hazus/flood/facility_analysis.py +++ b/src/hazus/flood/facility_analysis.py @@ -114,8 +114,22 @@ def get(name: str) -> pd.DataFrame: "cont": (get("flEssntFltyContDmgXRef"), "ContDmgFnId"), "util": (get("flUtilFltyDmgXRef"), "UtilFltyDmgFnId"), } + # default ids per (kind, hazard flag, occupancy), in table order + self._defaults: dict[tuple[str, str], dict[str, list]] = {} + for kind, (xref, id_col) in self._xrefs.items(): + x = xref[xref["Selected"].astype(bool) | xref["DefaultFn"].astype(bool)] + occ = x["SOccupId"].map(_strip) + for flag in ("HazardR", "HazardCA", "HazardCV"): + ids: dict[str, list] = {} + hit = x[flag].astype(bool).to_numpy() + for o, i in zip(occ[hit], x.loc[hit, id_col]): + ids.setdefault(o, []).append(i) + self._defaults[kind, flag] = ids rs = get("flRsFnEssntFlty") - self._restoration = rs.assign(EfClass=rs["EfClass"].map(_strip)) + self._restoration: dict[str, list[tuple]] = {} # rows per EF class, in table order + for row in zip(rs["EfClass"].map(_strip), rs["Min_Depth"].to_numpy(), rs["Max_Depth"].to_numpy(), + rs["Max_Restor_Days"].to_numpy(), rs["FunctionalDepth"].to_numpy()): + self._restoration.setdefault(row[0], []).append(row[1:]) details = get("flUtilFltyDmgUnionDetails") self._functional_depth = {_strip(o): d for o, d in zip(details["Occupancy"], details["Fnctdepth"])} ratios = get("hzPctContentOfStructureValue") @@ -123,13 +137,11 @@ def get(name: str) -> pd.DataFrame: def default_id(self, kind: str, hazard: str, depth: float, occupancy: str) -> int | None: """The ``flGetDamageXRefId*`` procedure: the xref default for an occupancy, hazard and depth.""" - xref, id_col = self._xrefs[kind] flag = "HazardR" if hazard == "R" else "HazardCV" if hazard == "C" and depth > 6 else "HazardCA" - rows = xref[(xref["SOccupId"].map(_strip) == _strip(occupancy)) & xref[flag].astype(bool) - & (xref["Selected"].astype(bool) | xref["DefaultFn"].astype(bool))] - if len(rows) > 1: # Dapper's QuerySingleOrDefault throws here - raise LookupError(f"{len(rows)} default {kind} damage functions for {occupancy!r} ({flag})") - return int(rows[id_col].iloc[0]) if len(rows) else None + ids = self._defaults[kind, flag].get(_strip(occupancy), []) + if len(ids) > 1: # Dapper's QuerySingleOrDefault throws here + raise LookupError(f"{len(ids)} default {kind} damage functions for {occupancy!r} ({flag})") + return int(ids[0]) if ids else None def function(self, kind: str, fn_id, hazard: str, depth: float, occupancy: str) -> dict[str, float] | None: """The function with id ``fn_id``, else the xref default (``kind`` is "struct", "cont" or "util").""" @@ -141,16 +153,11 @@ def function(self, kind: str, fn_id, hazard: str, depth: float, occupancy: str) def restoration(self, ef_class: str, flooded_depth: float) -> dict | None: """``flRsFnEssntFlty`` row with ``Min_Depth <= flooded_depth < Max_Depth`` as ``MaxDays``/``FunctionalDepth``.""" - rs = self._restoration - rows = rs[(rs["EfClass"] == _strip(ef_class)) & (rs["Min_Depth"] <= flooded_depth) - & (rs["Max_Depth"] > flooded_depth)] - if rows.empty: - return None - r = rows.iloc[0] - days = r["Max_Restor_Days"] - fd = r["FunctionalDepth"] - return {"MaxDays": None if pd.isna(days) else int(days), - "FunctionalDepth": None if pd.isna(fd) else float(np.float32(fd))} + for lo, hi, days, fd in self._restoration.get(_strip(ef_class), ()): + if lo <= flooded_depth and hi > flooded_depth: + return {"MaxDays": None if pd.isna(days) else int(days), + "FunctionalDepth": None if pd.isna(fd) else float(np.float32(fd))} + return None def functional_depth(self, utility_class: str) -> float | None: return self._functional_depth.get(_strip(utility_class)) @@ -324,18 +331,24 @@ def run_facilities(region: Mapping[str, pd.DataFrame], rasters, *, study_case_id if utilities: kinds += [(k, utility_results) for k in UTILITIES.values()] + kinds = [(kind, fn) for kind, fn in kinds if kind.base in region] + # sample every grid once for the facilities of all kinds + points = [facility_points(region[kind.base]) for kind, _ in kinds] + xs = np.concatenate([p[0] for p in points]) if points else np.empty(0) + ys = np.concatenate([p[1] for p in points]) if points else np.empty(0) + bounds = np.cumsum([0] + [len(p[0]) for p in points]) + sampled = {period: max_depth_at_points([r for r in rasters if int(r.return_period) == period], xs, ys) + for period in sorted(ids)} + out: dict[str, pd.DataFrame] = {} - for kind, fn in kinds: - if kind.base not in region: - continue + for k, (kind, fn) in enumerate(kinds): base = region[kind.base] flood = region[kind.flood] if kind.flood in region else flood_defaults(kind, base) - xs, ys = facility_points(base) + part = slice(bounds[k], bounds[k + 1]) parts = [fn(kind, base.iloc[:0], flood, [], [], study_case_id=study_case_id, functions=functions)] for period in sorted(ids): - grids = [r for r in rasters if int(r.return_period) == period] - depth, hazard = max_depth_at_points(grids, xs, ys) - parts.append(fn(kind, base, flood, depth, hazard, return_period=ids[period], + depth, hazard = sampled[period] + parts.append(fn(kind, base, flood, depth[part], hazard[part], return_period=ids[period], study_case_id=study_case_id, functions=functions)) out[kind.result] = pd.concat(parts, ignore_index=True) if len(parts) > 1 else parts[0] return out diff --git a/src/hazus/flood/gbs.py b/src/hazus/flood/gbs.py index 8d2294a..a531da3 100644 --- a/src/hazus/flood/gbs.py +++ b/src/hazus/flood/gbs.py @@ -80,13 +80,40 @@ def _occupancies() -> pd.DataFrame: return cls[["SOccup", "SortOrder"]] -def _occ_column(region: Region, name: str, suffix: str, blocks: pd.Series, occ: pd.Series) -> np.ndarray: - """Per-row value of the block's ``{SOccup}{suffix}`` column in a ``*OccupB`` table (NaN if missing).""" - t = region[name].copy() - t["CensusBlock"] = _strip(t["CensusBlock"]) - long = t.melt(id_vars="CensusBlock", var_name="col", value_name="v") - key = pd.DataFrame({"CensusBlock": blocks.to_numpy(), "col": (occ + suffix).to_numpy()}) - return key.merge(long, on=["CensusBlock", "col"], how="left")["v"].to_numpy(dtype="float64") +def _occ_positions(occid: np.ndarray) -> tuple[list[str], np.ndarray]: + """Specific occupancies in ``ClSOccupancy`` order, and each row's position among them (-1 if unknown).""" + occ = _occupancies() + lut = np.full(int(max(occ["SortOrder"].max(), occid.max(initial=0))) + 1, -1) + lut[occ["SortOrder"].to_numpy("int64")] = np.arange(len(occ)) + return list(occ["SOccup"]), lut[occid] + + +def _occ_column(region: Region, name: str, suffix: str, occs: list[str], uniques, codes: np.ndarray, + occ_pos: np.ndarray) -> np.ndarray: + """Per-row value of the block's ``{SOccup}{suffix}`` column in a ``*OccupB`` table (NaN if missing). + + Rows are given as ``codes`` into the block ids ``uniques`` (``pd.factorize``) and positions in ``occs``. + """ + t = region[name] + index = pd.Index(_strip(t["CensusBlock"])) + values = np.full((len(t), len(occs)), np.nan) + for j, o in enumerate(occs): + if f"{o}{suffix}" in t.columns: + values[:, j] = t[f"{o}{suffix}"].to_numpy("float64", na_value=np.nan) + rows = index.get_indexer(uniques)[codes] + ok = (rows >= 0) & (occ_pos >= 0) + out = np.full(len(codes), np.nan) + out[ok] = values[rows[ok], occ_pos[ok]] + return out + + +def _blocks_in_all(region: Region, blocks: np.ndarray) -> np.ndarray: + """Which of ``blocks`` appear in every exposure table (the ``flGBSXResults*`` inner joins).""" + keep = np.ones(len(blocks), dtype=bool) + for name in ("hzExposureContentOccupB", "hzExposureOccupB", "hzBldgCountOccupB", "hzSqFootageOccupB", + "flDepExposureOccupB", "flDepExposureContentOccupB"): + keep &= pd.Index(blocks).isin(_strip(region[name]["CensusBlock"]).unique()) + return keep # ------------------------------------------------------------------------------------------- base prep @@ -215,21 +242,21 @@ def _merge_calcs(nrds: pd.DataFrame, ppaw: pd.DataFrame, region: Region, with_de g["DecadeID"] = 0 # Blocks must appear in every exposure table (flGBSXResults* inner joins). - blocks = set(g["CensusBlock"]) - for name in ("hzExposureContentOccupB", "hzExposureOccupB", "hzBldgCountOccupB", "hzSqFootageOccupB", - "flDepExposureOccupB", "flDepExposureContentOccupB"): - blocks &= set(_strip(region[name]["CensusBlock"])) - g = g[g["CensusBlock"].isin(blocks)].reset_index(drop=True) - - occ = _occupancies().set_index("SortOrder")["SOccup"] - soccup = g["OCCID"].map(occ) - blk = g["CensusBlock"] - sqft = _occ_column(region, "hzSqFootageOccupB", "F", blk, soccup).astype(F32).astype("float64") # REAL column - eob = _occ_column(region, "hzExposureOccupB", "I", blk, soccup) - ecob = _occ_column(region, "hzExposureContentOccupB", "I", blk, soccup) - bcob = _occ_column(region, "hzBldgCountOccupB", "I", blk, soccup) - deob = _occ_column(region, "flDepExposureOccupB", "I", blk, soccup) - decob = _occ_column(region, "flDepExposureContentOccupB", "I", blk, soccup) + codes, uniques = pd.factorize(g["CensusBlock"]) + g = g[(codes >= 0) & _blocks_in_all(region, uniques)[codes]].reset_index(drop=True) + codes, uniques = pd.factorize(g["CensusBlock"]) + occid = g["OCCID"].to_numpy("int64") + occs, occ_pos = _occ_positions(occid) + + def occ_column(name, suffix): + return _occ_column(region, name, suffix, occs, uniques, codes, occ_pos) + + sqft = occ_column("hzSqFootageOccupB", "F").astype(F32).astype("float64") # REAL column + eob = occ_column("hzExposureOccupB", "I") + ecob = occ_column("hzExposureContentOccupB", "I") + bcob = occ_column("hzBldgCountOccupB", "I") + deob = occ_column("flDepExposureOccupB", "I") + decob = occ_column("flDepExposureContentOccupB", "I") pre = g["PrePct"].to_numpy("float64") post = g["PostPct"].to_numpy("float64") @@ -263,10 +290,13 @@ def amount(expo, share, factor): if inventory: psi = table("flBldgEconParamSalesAndInv") psi["Occupancy"] = _strip(psi["Occupancy"]) - p = pd.DataFrame({"SOccup": soccup}).merge(psi, left_on="SOccup", right_on="Occupancy", how="left") - sales = p["AnnualSalesPerSqFt"].fillna(0).to_numpy("float64") - invpct = p["BusinessInvPctofSales"].fillna(0).to_numpy("float64") / 100.0 - is_inv = g["OCCID"].isin(INV_OCCIDS).to_numpy() + p = pd.DataFrame({"SOccup": occs}).merge(psi, left_on="SOccup", right_on="Occupancy", how="left") + if len(p) != len(occs): + raise ValueError("flBldgEconParamSalesAndInv has duplicate occupancies") + # per occupancy, plus a trailing 0 for rows without one (position -1) + sales = np.append(p["AnnualSalesPerSqFt"].fillna(0).to_numpy("float64"), 0.0)[occ_pos] + invpct = np.append(p["BusinessInvPctofSales"].fillna(0).to_numpy("float64"), 0.0)[occ_pos] / 100.0 + is_inv = np.isin(occid, list(INV_OCCIDS)) def inv_amount(share): # SqFt*(share/100)*(PctAffected/100)*PctDmgInvC*Sales*(InvPct/100), strictly left to right @@ -305,7 +335,8 @@ def base_prep(region: Region, area_weighted: pd.DataFrame, dmg_final: pd.DataFra demo = region["hzDemographicsB"].copy() demo["CensusBlock"] = _strip(demo["CensusBlock"]) - used_schemes = set(scheme.loc[scheme["CensusBlock"].isin(aw_t["CensusBlock"]), "SchemeId"]) + aw_blocks = aw_t["CensusBlock"].unique() + used_schemes = set(scheme.loc[scheme["CensusBlock"].isin(aw_blocks), "SchemeId"]) used_haz = set(aw["HazTypeID"]) dec = table("flDecadeDescription") myb = demo[["CensusBlock", "MedianYearBuilt"]].dropna() @@ -323,7 +354,7 @@ def base_prep(region: Region, area_weighted: pd.DataFrame, dmg_final: pd.DataFra rds = _dist_damage(dmg_final, dist_r, used_schemes, used_haz, with_decade=False) # pre/post-FIRM split per block (only blocks in the study case, demographics and scheme mapping) - ppb = _pre_post_bldg(demo[demo["CensusBlock"].isin(aw_t["CensusBlock"])], scheme) + ppb = _pre_post_bldg(demo[demo["CensusBlock"].isin(aw_blocks)], scheme) ppp = _pre_post_pct(ppb, scheme) typeofhaz = aw["HazTypeID"].map({1: 1, 2: 2, 3: 2, 4: 2, 5: 5, 6: 6, 7: 7}) @@ -422,18 +453,23 @@ def dmg_counts(final: pd.DataFrame, region: Region) -> dict[str, pd.DataFrame]: def base_tab(final: pd.DataFrame, region: Region) -> pd.DataFrame: """``absf_GBSBaseTabT``: exposure, content exposure, count, depreciated values and sq ft per block × occupancy.""" - blocks = pd.Series(sorted(set(final["CensusBlock"])), name="CensusBlock") - for name in ("hzExposureOccupB", "hzExposureContentOccupB", "hzBldgCountOccupB", "flDepExposureOccupB", - "flDepExposureContentOccupB", "hzSqFootageOccupB"): - blocks = blocks[blocks.isin(set(_strip(region[name]["CensusBlock"])))] + blocks = np.sort(final["CensusBlock"].unique().astype(object)) + blocks = pd.Series(blocks[_blocks_in_all(region, blocks)].tolist(), name="CensusBlock") occ = _occupancies() + occs = list(occ["SOccup"]) x = pd.DataFrame({"CensusBlock": blocks}).merge(occ[["SOccup"]], how="cross") - x["Exposure"] = _occ_column(region, "hzExposureOccupB", "I", x["CensusBlock"], x["SOccup"]) - x["ExpCont"] = _occ_column(region, "hzExposureContentOccupB", "I", x["CensusBlock"], x["SOccup"]) - x["BldgCnt"] = _occ_column(region, "hzBldgCountOccupB", "I", x["CensusBlock"], x["SOccup"]) - x["DepExp"] = _occ_column(region, "flDepExposureOccupB", "I", x["CensusBlock"], x["SOccup"]) - x["DepExpCont"] = _occ_column(region, "flDepExposureContentOccupB", "I", x["CensusBlock"], x["SOccup"]) - x["SqFt"] = _occ_column(region, "hzSqFootageOccupB", "F", x["CensusBlock"], x["SOccup"]).astype(F32) + codes, uniques = pd.factorize(x["CensusBlock"]) + occ_pos = pd.Index(occs).get_indexer(x["SOccup"]) + + def occ_column(name, suffix): + return _occ_column(region, name, suffix, occs, uniques, codes, occ_pos) + + x["Exposure"] = occ_column("hzExposureOccupB", "I") + x["ExpCont"] = occ_column("hzExposureContentOccupB", "I") + x["BldgCnt"] = occ_column("hzBldgCountOccupB", "I") + x["DepExp"] = occ_column("flDepExposureOccupB", "I") + x["DepExpCont"] = occ_column("flDepExposureContentOccupB", "I") + x["SqFt"] = occ_column("hzSqFootageOccupB", "F").astype(F32) return x diff --git a/src/hazus/flood/hazard.py b/src/hazus/flood/hazard.py index f40521a..d57f557 100644 --- a/src/hazus/flood/hazard.py +++ b/src/hazus/flood/hazard.py @@ -11,6 +11,7 @@ from __future__ import annotations +import os from dataclasses import dataclass from enum import Enum from pathlib import Path @@ -20,6 +21,7 @@ FEET_PER_METER = 5000.0 / 1524.0 WGS84 = "EPSG:4326" +_POINT_TILE = 512 # cells: point reads are grouped into windows of at most this size class HazardType(str, Enum): @@ -75,16 +77,17 @@ def to_feet(depth: np.ndarray, units: str) -> np.ndarray: raise ValueError(f"unknown depth units {units!r}") -def prepare_depth_grid(src: str | Path, dst: str | Path, units: str = "feet", *, block: int = 1024) -> Path: +def prepare_depth_grid(src: str | Path, dst: str | Path, units: str = "feet", *, block: int = 1024, + threads: int | None = None) -> Path: """Write ``src`` as a float32 GeoTIFF in EPSG:4326 with depths in feet. Hazus uses ArcGIS ``ProjectRaster`` with nearest-neighbour resampling. Here the output grid size comes from GDAL's default transform, so cell boundaries can differ slightly from ArcGIS. Each output cell takes the source cell under its centre, computed exactly (no approximate transformer). The result - therefore does not depend on how the grid is split. + therefore does not depend on how the grid is split, nor on ``threads``. The grid is processed in ``block`` × ``block`` output windows, so memory use stays small for large - (e.g. 1 m lidar-based) depth grids. + (e.g. 1 m lidar-based) depth grids. GDAL warps and compresses with ``threads`` threads (default: all CPUs). """ import rasterio from rasterio.warp import Resampling, calculate_default_transform, reproject, transform_bounds @@ -100,7 +103,7 @@ def prepare_depth_grid(src: str | Path, dst: str | Path, units: str = "feet", *, transform, width, height = calculate_default_transform(ds.crs, WGS84, ds.width, ds.height, *ds.bounds) profile = dict(driver="GTiff", height=height, width=width, count=1, dtype="float32", crs=WGS84, transform=transform, nodata=np.float32(np.nan), compress="deflate", tiled=True, - blockxsize=256, blockysize=256, BIGTIFF="IF_SAFER") + blockxsize=256, blockysize=256, BIGTIFF="IF_SAFER", NUM_THREADS=threads or "ALL_CPUS") whole = Window(0, 0, ds.width, ds.height) with rasterio.open(dst, "w", **profile) as out_ds: for row in range(0, height, block): @@ -123,7 +126,8 @@ def prepare_depth_grid(src: str | Path, dst: str | Path, units: str = "feet", *, data = ds.read(1, window=swin, masked=True).astype("float64").filled(np.nan) reproject(source=data, destination=out, src_transform=ds.window_transform(swin), src_crs=ds.crs, src_nodata=np.nan, dst_transform=window_transform(win, transform), - dst_crs=WGS84, dst_nodata=np.nan, resampling=Resampling.nearest, tolerance=0) + dst_crs=WGS84, dst_nodata=np.nan, resampling=Resampling.nearest, tolerance=0, + num_threads=threads or os.cpu_count() or 1) out_ds.write(to_feet(out, units).astype("float32"), 1, window=win) return dst @@ -145,7 +149,8 @@ def merge_max(paths: Sequence[str | Path], dst: str | Path) -> Path: s.close() dst = Path(dst) profile = dict(driver="GTiff", height=mosaic.shape[1], width=mosaic.shape[2], count=1, dtype="float32", - crs=WGS84, transform=transform, nodata=np.float32(np.nan), compress="deflate", tiled=True) + crs=WGS84, transform=transform, nodata=np.float32(np.nan), compress="deflate", tiled=True, + NUM_THREADS="ALL_CPUS") with rasterio.open(dst, "w", **profile) as out: out.write(mosaic[0], 1) return dst @@ -229,9 +234,16 @@ def max_depth_at_points(rasters: Iterable[FloodRaster], xs: Sequence[float], cols = np.asarray(cols).astype("int64") inside = (rows >= 0) & (cols >= 0) & (rows < ds.height) & (cols < ds.width) values = np.full(len(xs), np.nan) - for k in np.flatnonzero(inside): # one cell per point; the grid is never read whole - cell = ds.read(1, window=((rows[k], rows[k] + 1), (cols[k], cols[k] + 1)), masked=True) - values[k] = np.ma.filled(cell.astype("float64"), np.nan)[0, 0] + # one read per tile that holds points (the window spanning that tile's points); never the whole grid + k = np.flatnonzero(inside) + tile = (rows[k] // _POINT_TILE) * (ds.width // _POINT_TILE + 1) + cols[k] // _POINT_TILE + order = np.argsort(tile, kind="stable") + k, tile = k[order], tile[order] + for group in np.split(k, np.flatnonzero(np.diff(tile)) + 1) if len(k) else []: + r, c = rows[group], cols[group] + r0, c0 = r.min(), c.min() + cells = ds.read(1, window=((r0, r.max() + 1), (c0, c.max() + 1)), masked=True) + values[group] = np.ma.filled(cells.astype("float64"), np.nan)[r - r0, c - c0] higher = np.nan_to_num(values, nan=-np.inf) > depth depth[higher] = values[higher] hazard[higher] = fr.hazard_type.controlling_code diff --git a/src/hazus/flood/scenario.py b/src/hazus/flood/scenario.py index 0efc0ba..ea9c618 100644 --- a/src/hazus/flood/scenario.py +++ b/src/hazus/flood/scenario.py @@ -50,11 +50,12 @@ def area_weighting(blocks, rasters: Sequence[FloodRaster], *, boundary=None, stu ids = return_period_ids(rasters) aw: dict[tuple[str, str], pd.DataFrame] = {} sr: dict[tuple[str, str], pd.DataFrame] = {} + block_counts: dict = {} # shared by grids that snap to the same block grid for period in sorted(ids): for raster in (r for r in rasters if int(r.return_period) == period): rp = ids[period] hazard = HazardType(raster.hazard_type).controlling_code - long = area_weighted_depths(blocks, raster, boundary=boundary) + long = area_weighted_depths(blocks, raster, boundary=boundary, block_counts=block_counts) # absp_GBSInsertAW replaces earlier rows of the same return period and hazard aw[rp, hazard] = area_weighted_rows(long, study_case_id=study_case_id, return_period=rp, controlling_hazard=hazard, analysis_opt_id=analysis_opt_id) diff --git a/tests/test_vectorized_paths.py b/tests/test_vectorized_paths.py new file mode 100644 index 0000000..c458a12 --- /dev/null +++ b/tests/test_vectorized_paths.py @@ -0,0 +1,421 @@ +"""Vectorized code paths against the per-row code they replaced. + +Each reference below is the earlier implementation, kept here as an oracle. The fast paths must give the same +result exactly (values, dtypes and row order) on inputs chosen to hit the edge cases: duplicate rows, missing +return periods, half-way rounding, int64 overflow, NoData, misaligned and partial grids, study boundaries that +cut through tiles, points outside every grid, ambiguous or missing damage-function defaults. +""" + +import sys +from fractions import Fraction +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + +from hazus import data +from hazus.flood.aal import ECON_COLS, aal_by_soccup, aal_coefficients, aal_coefficients_decimal +from hazus.flood.facility_analysis import ( + ESSENTIAL_FACILITIES, + UTILITIES, + FacilityDamageFunctions, + _strip, + run_facilities, +) +from hazus.flood.gbs import _occ_column, _occ_positions, _occupancies + +sys.path.insert(0, str(Path(__file__).parent)) +import synth # noqa: E402 + +PART_KEYS = ["StudyCaseId", "CensusBlock", "ControlHazard", "SOccup", "Pre_PostFirm", "AnalysisOptId"] + + +# --- AAL (absp_AALDoIt) --- + +def aal_reference(ecloss: pd.DataFrame, coefficients) -> pd.DataFrame: + def round_half_away(x: Fraction) -> int: + r = int(abs(x) + Fraction(1, 2)) + return r if x >= 0 else -r + + df = ecloss.copy() + rp = pd.to_numeric(df["ReturnPeriodId"].astype(str).str.strip(), errors="coerce") + df = df[rp.isin(list(coefficients))].assign(ReturnPeriod=rp.astype("Int64")) + df = df.groupby(PART_KEYS + ["ReturnPeriod"], as_index=False)[ECON_COLS].sum() + coef = {int(k): Fraction(v) for k, v in coefficients.items()} + rows = [] + for key, g in df.sort_values("ReturnPeriod", ascending=False).groupby(PART_KEYS, sort=False): + totals = dict.fromkeys(ECON_COLS, Fraction(0)) + prev = None + for r in g.itertuples(index=False): + c = coef[int(r.ReturnPeriod)] + for col in ECON_COLS: + v = int(getattr(r, col)) + totals[col] += c * v if prev is None else c * (v + int(getattr(prev, col))) / 2 + prev = r + rows.append({**dict(zip(PART_KEYS, key)), "ReturnPeriodId": "Annual", + **{col: round_half_away(totals[col]) for col in ECON_COLS}}) + return pd.DataFrame(rows, columns=PART_KEYS + ["ReturnPeriodId"] + ECON_COLS) + + +def ecloss_rows(seed: int, periods, n: int = 2000, scale: int = 1) -> pd.DataFrame: + """Random flFRGBSEcLossBySOccup rows: duplicates, unknown and padded return periods, negative values.""" + rng = np.random.default_rng(seed) + labels = [str(p) for p in periods] + [f" {periods[0]} ", "Annual", "999"] + df = pd.DataFrame({ + "StudyCaseId": 1, "ReturnPeriodId": rng.choice(labels, n), + "CensusBlock": [f"240010001001{i:03d}" for i in rng.integers(0, 80, n)], + "SOccup": rng.choice(["RES1", "COM1", "IND2"], n), "ControlHazard": rng.choice(["R", "C"], n), + "Pre_PostFirm": rng.integers(0, 2, n), "AnalysisOptId": 0}) + for col in ECON_COLS: + df[col] = rng.integers(-500, 20000, n) * scale + df.loc[:9, "BuildingLoss"] = 250 # 500-year coefficient 0.002 x 250 = 0.5: exactly half way + return df + + +@pytest.mark.parametrize("periods, coefficients, scale", [ + ((10, 25, 50, 100, 500), "decimal", 1), + ((7, 13, 333), "decimal", 1), + ((2, 5, 10, 25, 50, 100, 200, 500, 1000), "decimal", 1), + ((10, 100, 500), "float", 1), # binary fractions: Python integers, not int64 + ((10, 100, 500), "decimal", 10 ** 12), # int64 would overflow +]) +def test_aal_matches_fraction_loop(periods, coefficients, scale): + coef = aal_coefficients_decimal(periods) if coefficients == "decimal" else aal_coefficients(periods) + ecloss = ecloss_rows(len(periods), periods, scale=scale) + expected = aal_reference(ecloss, coef) + pd.testing.assert_frame_equal(aal_by_soccup(ecloss, coef), expected, check_exact=True) + assert len(expected) > 100 + + +def test_aal_without_matching_return_periods(): + ecloss = ecloss_rows(0, (10, 100, 500)) + coef = aal_coefficients_decimal((20, 200, 2000)) + pd.testing.assert_frame_equal(aal_by_soccup(ecloss, coef), aal_reference(ecloss, coef), check_exact=True) + + +# --- GBS per-block x occupancy lookups --- + +def occ_column_reference(region, name, suffix, blocks: pd.Series, occ: pd.Series) -> np.ndarray: + t = region[name].copy() + t["CensusBlock"] = t["CensusBlock"].astype(str).str.rstrip() + long = t.melt(id_vars="CensusBlock", var_name="col", value_name="v") + key = pd.DataFrame({"CensusBlock": blocks.to_numpy(), "col": (occ + suffix).to_numpy()}) + return key.merge(long, on=["CensusBlock", "col"], how="left")["v"].to_numpy(dtype="float64") + + +@pytest.mark.parametrize("name, suffix", [("hzExposureOccupB", "I"), ("hzSqFootageOccupB", "F")]) +def test_occ_column_matches_melt_merge(name, suffix): + rng = np.random.default_rng(3) + region = synth.make_region(seed=3, n_blocks=12) + t = region[name].drop(columns=f"RES2{suffix}").astype({f"COM1{suffix}": "float64"}) + t.loc[2, f"COM1{suffix}"] = np.nan + t["CensusBlock"] = t["CensusBlock"] + " " # padded, as char(15) can be + region[name] = t + blocks = pd.Series(list(rng.choice(list(t["CensusBlock"].str.strip()) + ["UNKNOWN"], 500))) + occid = rng.choice(list(_occupancies()["SortOrder"]) + [999], 500) + occs, occ_pos = _occ_positions(occid) + codes, uniques = pd.factorize(blocks) + got = _occ_column(region, name, suffix, occs, uniques, codes, occ_pos) + soccup = pd.Series(occid).map(_occupancies().set_index("SortOrder")["SOccup"]) + expected = occ_column_reference(region, name, suffix, blocks, soccup) + np.testing.assert_array_equal(got, expected) + assert np.isnan(expected).any() and not np.isnan(expected).all() + + +# --- area weighting (Analysis Area Weighted) --- + +def area_weighted_reference(blocks, flood, *, block_id="CensusBlock", boundary=None, tile=2048): + import rasterio + from rasterio.features import rasterize + + from hazus.flood.exposure import CELL_SIZE, ZONE_LETTERS, _box, depth_bin, snapped_grid, zone_code + from hazus.flood.hazard import HazardType + + blocks = blocks.reset_index(drop=True) + hazard = HazardType(flood.hazard_type) + sindex = blocks.sindex + with rasterio.open(flood.path) as ds: + grid = snapped_grid(tuple(blocks.total_bounds), ds.transform, CELL_SIZE) + dl, dt = ds.transform.c, ds.transform.f + drx, dry = abs(ds.transform.a), abs(ds.transform.e) + nodata = ds.nodata + r0, r1, c0, c1 = grid.window_for_bounds(ds.bounds.left, ds.bounds.bottom, ds.bounds.right, ds.bounds.top) + parts = [] + for tr in range(r0, r1, tile): + for tc in range(c0, c1, tile): + h, w = min(tile, r1 - tr), min(tile, c1 - tc) + tx0, ty0 = grid.left + tc * grid.cell, grid.top - tr * grid.cell + hits = sindex.query(_box(tx0, ty0 - h * grid.cell, tx0 + w * grid.cell, ty0)) + if len(hits) == 0: + continue + transform = grid.sub_transform(tr, tc) + blk = rasterize(((blocks.geometry.values[i], int(i) + 1) for i in hits), out_shape=(h, w), + transform=transform, fill=0, dtype="int32") + if not blk.any(): + continue + if boundary is not None: + inside = rasterize([(boundary, 1)], out_shape=(h, w), transform=transform, fill=0, + dtype="uint8").astype(bool) + else: + inside = blk > 0 + xs = tx0 + (np.arange(w) + 0.5) * grid.cell + ys = ty0 - (np.arange(h) + 0.5) * grid.cell + dcols = np.floor((xs - dl) / drx).astype("int64") + drows = np.floor((dt - ys) / dry).astype("int64") + col_ok = (dcols >= 0) & (dcols < ds.width) + row_ok = (drows >= 0) & (drows < ds.height) + if not col_ok.any() or not row_ok.any(): + continue + wc0, wc1 = dcols[col_ok].min(), dcols[col_ok].max() + 1 + wr0, wr1 = drows[row_ok].min(), drows[row_ok].max() + 1 + band = ds.read(1, window=((wr0, wr1), (wc0, wc1))).astype("float64") + if nodata is not None and not np.isnan(nodata): + band[band == nodata] = np.nan + depth = np.full((h, w), np.nan) + full = band[np.ix_(np.clip(drows - wr0, 0, wr1 - wr0 - 1), np.clip(dcols - wc0, 0, wc1 - wc0 - 1))] + ok2d = np.outer(row_ok, col_ok) + depth[ok2d] = full[ok2d] + bins = depth_bin(depth) + valid = (blk > 0) & inside & (bins >= 0) + if not valid.any(): + continue + zones = zone_code(depth[valid], hazard) + key = (blk[valid].astype("int64") << 16) | (zones.astype("int64") << 8) | (bins[valid] // 50) + uniq, counts = np.unique(key, return_counts=True) + parts.append(pd.DataFrame({"key": uniq, "Count": counts})) + if not parts: + return pd.DataFrame(columns=["CensusBlock", "CoastalZone", "ReClass", "Count", "TotalCount", "PctAffected"]) + agg = pd.concat(parts).groupby("key", as_index=False)["Count"].sum() + keys = agg["key"].to_numpy(dtype="int64") + idx = (keys >> 16) - 1 + flooded = np.unique(idx) + totals = np.zeros(len(flooded), dtype="int64") + for k, i in enumerate(flooded): + geom = blocks.geometry.values[i] + b0, b1, a0, a1 = grid.window_for_bounds(*geom.bounds) + if b1 > b0 and a1 > a0: + totals[k] = int(rasterize([(geom, 1)], out_shape=(b1 - b0, a1 - a0), transform=grid.sub_transform(b0, a0), + fill=0, dtype="uint8").sum()) + out = pd.DataFrame({ + "CensusBlock": blocks[block_id].to_numpy()[idx], + "CoastalZone": [ZONE_LETTERS.get(int(z), "N") for z in (keys >> 8) & 0xFF], + "ReClass": ((keys & 0xFF) * 50).astype("int32"), + "Count": agg["Count"].to_numpy(), + "TotalCount": pd.Series(totals, index=flooded).loc[idx].to_numpy(), + }) + out["PctAffected"] = out["Count"] / out["TotalCount"] * 100.0 + return out.sort_values(["CensusBlock", "CoastalZone", "ReClass"], ignore_index=True) + + +def jittered_blocks(n_side: int = 10, side: float = 0.001, seed: int = 3): + """Census blocks that tile the region: quadrilaterals on a jittered lattice.""" + gpd = pytest.importorskip("geopandas") + from shapely.geometry import Polygon + + rng = np.random.default_rng(seed) + gx = -77.0 + np.arange(n_side + 1)[None, :] * side + rng.uniform(-0.3, 0.3, (n_side + 1, n_side + 1)) * side + gy = 38.8 + np.arange(n_side + 1)[:, None] * side + rng.uniform(-0.3, 0.3, (n_side + 1, n_side + 1)) * side + polys = [Polygon([(gx[r, c], gy[r, c]), (gx[r, c + 1], gy[r, c + 1]), (gx[r + 1, c + 1], gy[r + 1, c + 1]), + (gx[r + 1, c], gy[r + 1, c])]) for r in range(n_side) for c in range(n_side)] + return gpd.GeoDataFrame({"CensusBlock": [f"B{i:03d}" for i in range(len(polys))]}, geometry=polys, crs="EPSG:4326") + + +def depth_grid(path, bounds, res, seed, *, nodata=np.nan, scale=1.0, offset=(0.0, 0.0), shape=None): + """Depth (ft) falling from west to east with noise; <= 0 and 2 % of cells at random are NoData.""" + rasterio = pytest.importorskip("rasterio") + from rasterio.transform import from_origin + + left, top = bounds[0] + offset[0], bounds[3] + offset[1] + h, w = shape or (int(np.ceil((bounds[3] - bounds[1]) / res)), int(np.ceil((bounds[2] - bounds[0]) / res))) + rng = np.random.default_rng(seed) + d = (np.linspace(12, -3, w)[None, :] + rng.normal(0, 1.5, (h, w))) * scale + d[(d <= 0) | (rng.random((h, w)) < 0.02)] = nodata + with rasterio.open(path, "w", driver="GTiff", height=h, width=w, count=1, dtype="float32", crs="EPSG:4326", + transform=from_origin(left, top, res, res), nodata=nodata, tiled=True, blockxsize=256, + blockysize=256) as ds: + ds.write(d.astype("float32"), 1) + return path + + +AW_CASES = { # name: (grid options, hazard, boundary, tile) + "coarse": (dict(res=0.0001, seed=1), "riverine", None, 257), + "finer_than_block_grid_misaligned": (dict(res=0.0000073, seed=2, offset=(0.0003, -0.0002)), "coastal", None, 2048), + "nodata_value_whole_boundary": (dict(res=0.000031, seed=4, nodata=-9999.0, scale=1.7), "coastal", "whole", 300), + "partial_boundary": (dict(res=0.0000413, seed=5, offset=(-0.001, 0.0011)), "riverine", "partial", 300), + "grid_inside_region": (dict(res=0.00005, seed=6, offset=(0.004, -0.004), shape=(60, 70)), "coastal", None, 517), +} + + +@pytest.mark.parametrize("case", list(AW_CASES)) +def test_area_weighting_matches_reference(case, tmp_path): + from shapely.geometry import box + + from hazus.flood.exposure import area_weighted_depths + from hazus.flood.hazard import FloodRaster + + opts, hazard, boundary, tile = AW_CASES[case] + blocks = jittered_blocks() + bounds = tuple(blocks.total_bounds) + union = blocks.geometry.union_all() + boundary = {None: None, "whole": union.buffer(0.00001), + "partial": union.intersection(box(bounds[0], bounds[1], (bounds[0] + bounds[2]) / 2 + 0.0007, + bounds[3])).buffer(0.00001)}[boundary] + flood = FloodRaster(depth_grid(tmp_path / "d.tif", bounds, **opts), hazard, 100) + expected = area_weighted_reference(blocks, flood, boundary=boundary, tile=tile) + got = area_weighted_depths(blocks, flood, boundary=boundary, tile=tile) + pd.testing.assert_frame_equal(got, expected, check_exact=True) + assert len(expected) > 20 + + +def test_area_weighting_shares_block_counts(tmp_path): + from hazus.flood.exposure import area_weighted_depths + from hazus.flood.hazard import FloodRaster + + blocks = jittered_blocks() + bounds = tuple(blocks.total_bounds) + counts: dict = {} + for k, res in enumerate([0.0001, 0.0001, 0.00003]): # the last grid snaps the block grid differently + flood = FloodRaster(depth_grid(tmp_path / f"d{k}.tif", bounds, res=res, seed=k), "riverine", 100) + shared = area_weighted_depths(blocks, flood, block_counts=counts) + pd.testing.assert_frame_equal(shared, area_weighted_depths(blocks, flood), check_exact=True) + assert len(counts) == 2 + + +# --- facilities --- + +def default_id_reference(xref, id_col, hazard, depth, occupancy): + flag = "HazardR" if hazard == "R" else "HazardCV" if hazard == "C" and depth > 6 else "HazardCA" + rows = xref[(xref["SOccupId"].map(_strip) == _strip(occupancy)) & xref[flag].astype(bool) + & (xref["Selected"].astype(bool) | xref["DefaultFn"].astype(bool))] + if len(rows) > 1: + raise LookupError + return int(rows[id_col].iloc[0]) if len(rows) else None + + +def outcome(fn, *args): + try: + return fn(*args) + except LookupError: + return LookupError + + +XREFS = {"struct": ("flEssntFltyStructDmgXRef", "BldgDmgFnId"), "cont": ("flEssntFltyContDmgXRef", "ContDmgFnId"), + "util": ("flUtilFltyDmgXRef", "UtilFltyDmgFnId")} + + +def test_default_ids_match_table_filter(): + custom = pd.DataFrame({"SOccupId": ["PWTL", "PWTL ", "PWTL", "WWTL"], "UtilFltyDmgFnId": [7, 10, 11, 12], + "DefaultFn": [1, 0, 0, 0], "Selected": [0, 1, 0, 1], "HazardR": [1, 0, 0, 1], + "HazardCA": [0, 1, np.nan, 0], "HazardCV": [0, 0, 1, 1]}) + for tables in ({}, {"flUtilFltyDmgXRef": custom}): + f = FacilityDamageFunctions(tables) + for kind, (name, id_col) in XREFS.items(): + xref = tables.get(name, data.table(name)) + for occ in [*set(xref["SOccupId"]), " PWTL", "NOPE"]: + for hazard, depth in [("R", 1.0), ("C", 6.0), ("C", 6.5)]: + assert outcome(f.default_id, kind, hazard, depth, occ) == \ + outcome(default_id_reference, xref, id_col, hazard, depth, occ), (kind, occ, hazard, depth) + + +def test_restoration_matches_table_filter(): + rs = data.table("flRsFnEssntFlty").assign(EfClass=lambda t: t["EfClass"].map(_strip)) + f = FacilityDamageFunctions() + + def reference(ef_class, d): + rows = rs[(rs["EfClass"] == _strip(ef_class)) & (rs["Min_Depth"] <= d) & (rs["Max_Depth"] > d)] + if rows.empty: + return None + r = rows.iloc[0] + return {"MaxDays": None if pd.isna(r["Max_Restor_Days"]) else int(r["Max_Restor_Days"]), + "FunctionalDepth": None if pd.isna(r["FunctionalDepth"]) else float(np.float32(r["FunctionalDepth"]))} + + depths = [*np.arange(-6.0, 30.0, 0.25), -4, 0, 1, 4, 8, 12, 24, np.nan] + for ef_class in [*set(rs["EfClass"]), " EFHS ", "NOPE"]: + for d in depths: + assert f.restoration(ef_class, d) == reference(ef_class, d), (ef_class, d) + + +def max_depth_reference(rasters, xs, ys): + import rasterio + + xs, ys = np.asarray(xs, dtype="float64"), np.asarray(ys, dtype="float64") + depth = np.zeros(len(xs)) + hazard = np.full(len(xs), "R", dtype=object) + for fr in rasters: + with rasterio.open(fr.path) as ds: + rows, cols = rasterio.transform.rowcol(ds.transform, xs, ys, op=np.floor) + rows, cols = np.asarray(rows).astype("int64"), np.asarray(cols).astype("int64") + inside = (rows >= 0) & (cols >= 0) & (rows < ds.height) & (cols < ds.width) + values = np.full(len(xs), np.nan) + for k in np.flatnonzero(inside): + cell = ds.read(1, window=((rows[k], rows[k] + 1), (cols[k], cols[k] + 1)), masked=True) + values[k] = np.ma.filled(cell.astype("float64"), np.nan)[0, 0] + higher = np.nan_to_num(values, nan=-np.inf) > depth + depth[higher] = values[higher] + hazard[higher] = fr.hazard_type.controlling_code + return depth, hazard + + +@pytest.fixture +def facility_grids(tmp_path): + from hazus.flood.hazard import FloodRaster + + res = 0.00001 # 1300 x 700 cells: several 512-cell read windows + bounds = (-77.0, 38.8, -77.0 + 1300 * res, 38.8 + 700 * res) + river = depth_grid(tmp_path / "r.tif", bounds, res, seed=7, nodata=-9999.0) + coast = depth_grid(tmp_path / "c.tif", bounds, res * 1.3, seed=8, offset=(0.004, 0.0), scale=1.4) + return bounds, res, [FloodRaster(river, "riverine", 100), FloodRaster(coast, "coastal", 100)] + + +def test_point_depths_match_cell_reads(facility_grids): + from hazus.flood.hazard import max_depth_at_points + + (left, bottom, right, top), res, rasters = facility_grids + rng = np.random.default_rng(9) + xs = np.r_[rng.uniform(left - 0.002, right + 0.006, 3000), left + res * np.arange(0, 1300, 97)] + ys = np.r_[rng.uniform(bottom - 0.002, top + 0.002, 3000), np.full(14, top - 3 * res)] # cell edges too + depth, hazard = max_depth_at_points(rasters, xs, ys) + ref_depth, ref_hazard = max_depth_reference(rasters, xs, ys) + np.testing.assert_array_equal(depth, ref_depth) + np.testing.assert_array_equal(hazard, ref_hazard) + assert set(hazard) == {"R", "C"} and (depth == 0).any() and (depth > 0).any() + + +def test_facility_kinds_sampled_together_match_each_alone(facility_grids): + (left, bottom, right, top), _, rasters = facility_grids + rng = np.random.default_rng(10) + care, school, water = ESSENTIAL_FACILITIES["care"], ESSENTIAL_FACILITIES["school"], UTILITIES["potable_water"] + region = {} + for kind, n in ((care, 40), (school, 0), (water, 30)): + ids = [f"{kind.id_col[:2]}{i:04d}" for i in range(n)] + cls = {"UtilFcltyClass": "PWTL"} if kind is water else {"EfClass": rng.choice(["EFHS", "EFS1"], n)} + region[kind.base] = pd.DataFrame({kind.id_col: ids, **cls, "Cost": 1000.0, + "Longitude": rng.uniform(left, right, n), + "Latitude": rng.uniform(bottom, top, n)}) + functions = FacilityDamageFunctions() + together = run_facilities(region, rasters, functions=functions) + for kind in (care, school, water): + alone = run_facilities({kind.base: region[kind.base]}, rasters, functions=functions) + pd.testing.assert_frame_equal(together[kind.result], alone[kind.result], check_exact=True) + + +# --- depth-grid preparation --- + +def test_prepare_depth_grid_does_not_depend_on_threads(tmp_path): + rasterio = pytest.importorskip("rasterio") + from rasterio.transform import from_origin + + from hazus.flood.hazard import prepare_depth_grid + + rng = np.random.default_rng(11) + d = (np.linspace(3, -1, 400)[None, :] + rng.normal(0, 0.3, (300, 400))).astype("float32") + d[d <= 0] = -9999.0 + src = tmp_path / "utm.tif" + with rasterio.open(src, "w", driver="GTiff", height=300, width=400, count=1, dtype="float32", crs="EPSG:26918", + transform=from_origin(320000, 4300000, 1.0, 1.0), nodata=-9999.0) as ds: + ds.write(d, 1) + grids = [prepare_depth_grid(src, tmp_path / f"out{t}.tif", "meters", block=128, threads=t) for t in (1, None)] + with rasterio.open(grids[0]) as a, rasterio.open(grids[1]) as b: + assert a.transform == b.transform + np.testing.assert_array_equal(a.read(1), b.read(1))