From 67f0cbad5b8efc2a6bdad2ea7a938ad86a4a29c6 Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Mon, 21 Sep 2026 16:20:55 +1200 Subject: [PATCH 1/8] Ingest the cs_nshm_2022 run3 campaign into an IMDB Adapted from nz_sim_validation_ingest.py for run3's intensity_measures.h5 layout (one netCDF4 group per IM, one dataset per component), trimmed to geom and rotd50, 25 pSA periods and no FAS or GMM records. Magnitudes stay BoldM, the group convention, rather than converting to Mw. The events table takes the NSHM catalogue magnitude from a CSV exported from nshmdb. The build stops on inconsistent input: NaNs, site values that differ between ruptures, distances that differ between realisations, or a magnitude or hypocentre that disagrees with the h5. It writes a .partial file and renames it only after validate(), exact row counts and a 1,000-record spot check against the h5 all pass. Co-Authored-By: Claude Opus 5 --- scripts/cs_nshm_2022_ingest.py | 636 ++++++++++++++++++++++ scripts/tests/test_cs_nshm_2022_ingest.py | 469 ++++++++++++++++ 2 files changed, 1105 insertions(+) create mode 100644 scripts/cs_nshm_2022_ingest.py create mode 100644 scripts/tests/test_cs_nshm_2022_ingest.py diff --git a/scripts/cs_nshm_2022_ingest.py b/scripts/cs_nshm_2022_ingest.py new file mode 100644 index 0000000..653b2b4 --- /dev/null +++ b/scripts/cs_nshm_2022_ingest.py @@ -0,0 +1,636 @@ +#!/usr/bin/env -S uv run --script +# /// script +# requires-python = ">=3.12,<3.14" # fiona (via source-modelling) has no 3.14 wheel yet +# dependencies = [ +# "numpy>=2", +# "pandas>=3", +# "h5py>=3.11", +# "shapely", +# "source-modelling>=2026.8.3", +# "ucgmsim-imdb>=2026.9.2", +# ] +# +# [tool.uv] +# python-preference = "only-managed" # avoid a broken/non-standard system Python on PATH +# /// + +"""Ingest the cs_nshm_2022 run3 simulation campaign into an IMDB. + +One event per NSHM 2022 crustal rupture and one realisation per simulated +`R`. The event_id is the campaign's rupture id, which is the rupture's +crustal `nshm_id` in nshmdb, NOT its `rupture_id`; the rel_id is +`_R`. + +- IMs, distances, vs30 and z1.0/z2.5 come from each realisation's + `intensity_measures.h5`, as written by the workflow's im-calc: netCDF4, one + group per IM, one dataset per component. +- Source and rupture metadata come from its `realisation.json`. +- The rupture's NSHM magnitude and annual rate come from a CSV exported from + nshmdb, which is not available where this runs. +- Canonical site coordinates come from the campaign's `stations_input.ll`. + +Trimmed for size: only the geom and rotd50 components, pSA on 25 of im-calc's +111 periods, and no FAS or empirical-GMM records. + +Every magnitude is BoldM (Hanks & Kanamori 1979 eq. 7), the group convention. +nz_sim_validation_ingest.py converts to Mw instead. +""" + +import argparse +import json +import re +import time +from collections.abc import Iterable, Sequence +from dataclasses import dataclass +from pathlib import Path + +import h5py +import numpy as np +import pandas as pd +from imdb import IMDB, schema +from shapely.geometry import LineString, MultiLineString, MultiPolygon, Polygon +from source_modelling import moment +from source_modelling.sources import Fault + +# 0.1 s steps to 1 s, then 1 s steps. im-calc computed no 16-19 s; 20 s is its last. +PSA_PERIODS = ( + 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, + 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, + 11.0, 12.0, 13.0, 14.0, 15.0, 20.0, +) # fmt: skip +COMPONENTS = ("geom", "rotd50") +N_IM_CALC_PERIODS = 111 +SITE_FIELDS = ("vs30", "z1pt0", "z2pt5") +DISTANCE_FIELDS = ("rrup", "rjb", "rx", "ry") +# How closely realisation.json must reproduce the h5's own attributes. Measured +# over all 221 realisations done on 2026-09-21: magnitude identical, hypocentre +# within 3e-14 degrees. +TOLERANCES = {"magnitude": 1e-9, "hypo_lat": 1e-9, "hypo_lon": 1e-9, "hypo_depth": 1e-6} + +DB_META = { + "dataset_id": "cs_nshm_2022_run3", + "magnitude_convention": ( + "BoldM (Hanks & Kanamori 1979 eq. 7). events.magnitude is the NSHM 2022 " + "catalogue magnitude; realisations.magnitude is the realisation's own " + "moment-summed total." + ), + "psa_period_selection": ( + "25 of im-calc's 111 periods: 0.1-1.0 s by 0.1 s, 2-15 s by 1 s, and 20 s. " + "im-calc computed no 16-19 s." + ), + "station_coordinates": ( + "sites.lat/lon are the canonical stations_input.ll coordinates. The " + "simulation snapped each station to its rupture's grid (~0.1 km away, so " + "up to ~0.25 km apart between ruptures); those coordinates are not stored." + ), +} + + +# ---- inputs ------------------------------------------------------------------ + + +def select_period_indices( + available: np.ndarray, wanted: Sequence[float] = PSA_PERIODS +) -> np.ndarray: + """Column of each wanted period in im-calc's period grid. + + Raises ValueError unless every wanted period matches exactly one available one. + """ + indices = [] + for period in wanted: + (matches,) = np.nonzero(np.isclose(available, period, rtol=1e-9, atol=0.0)) + if len(matches) != 1: + raise ValueError( + f"pSA period {period} s matches {len(matches)} of im-calc's " + "periods, expected exactly 1" + ) + indices.append(int(matches[0])) + return np.array(indices) + + +def parse_realisation_ids(lines: Iterable[str]) -> list[tuple[str, int]]: + """Parse `/R` lines into (rupture, n) pairs, in ingest order. + + Blank lines are skipped and duplicates dropped; the order is by integer + rupture id, then realisation. + """ + ids = set() + for line in lines: + line = line.strip() + if not line: + continue + match = re.fullmatch(r"(\d+)/R(\d+)", line) + if match is None: + raise ValueError(f"not a /R id: {line!r}") + ids.add((str(int(match[1])), int(match[2]))) + return sorted(ids, key=lambda rel: (int(rel[0]), rel[1])) + + +def group_by_rupture(realisations: Sequence[tuple[str, int]]) -> dict[str, list[int]]: + """Realisation numbers per rupture, keeping the input order.""" + grouped: dict[str, list[int]] = {} + for rupture, n in realisations: + grouped.setdefault(rupture, []).append(n) + return grouped + + +def rel_id(rupture: str, n: int) -> str: + """The realisation's stable id.""" + return f"{rupture}_R{n}" + + +def load_stations_input(path: Path) -> pd.DataFrame: + """Canonical station coordinates (`lon lat name` per line), indexed by name.""" + df = pd.read_csv( + path, sep=r"\s+", header=None, names=["lon", "lat", "site_id"], dtype={"site_id": str} + ) + duplicated = df["site_id"][df["site_id"].duplicated()] + if len(duplicated): + raise ValueError(f"{path}: duplicate station names, e.g. {list(duplicated[:5])}") + return df.set_index("site_id") + + +def load_nshm_attributes(path: Path) -> pd.DataFrame: + """Rows of export_nshm_rupture_attributes.py's CSV, indexed by rupture id.""" + return pd.read_csv(path, dtype={"rupture_id": str}).set_index("rupture_id") + + +def is_real_station(site_id: str) -> bool: + """GeoNet station codes are 3-4 characters; the virtual grid's codes are 7.""" + return len(site_id) != 7 + + +# ---- intensity_measures.h5 --------------------------------------------------- + + +@dataclass +class RealisationIMs: + """One realisation's intensity_measures.h5, reduced to what this database keeps.""" + + stations: np.ndarray + """Station names, str.""" + site: dict[str, np.ndarray] + """SITE_FIELDS per station.""" + distances: dict[str, np.ndarray] + """DISTANCE_FIELDS per station, km.""" + psa: dict[str, np.ndarray] + """Component -> (n_stations, len(PSA_PERIODS)) float32.""" + scalars: dict[str, dict[str, np.ndarray]] + """Scalar IM -> component -> (n_stations,) float32; rotd-undefined IMs omitted.""" + attrs: dict[str, float] + """The file's magnitude and hypocentre attributes (TOLERANCES keys).""" + + +def _decode(values: np.ndarray) -> np.ndarray: + return np.array([v.decode() if isinstance(v, bytes) else v for v in values], dtype=object) + + +def read_ims(h5_path: Path) -> RealisationIMs: + """Read the parts of one intensity_measures.h5 this database keeps.""" + with h5py.File(h5_path, "r") as h5: + psa_group = h5["pSA"] + periods = psa_group["period"][:] + if len(periods) != N_IM_CALC_PERIODS: + raise ValueError( + f"{h5_path}: {len(periods)} pSA periods, expected {N_IM_CALC_PERIODS}" + ) + columns = select_period_indices(periods) + raw_stations = psa_group["station"][:] + scalars: dict[str, dict[str, np.ndarray]] = {} + for im in schema.SCALAR_IMS: + if not np.array_equal(h5[im]["station"][:], raw_stations): + raise ValueError(f"{h5_path}: {im} lists stations in a different order from pSA") + scalars[im] = {} + for component in COMPONENTS: + if component in h5[im]: + scalars[im][component] = h5[im][component][:].astype(np.float32) + elif not (component.startswith("rotd") and im in schema.ROTD_UNDEFINED): + raise ValueError(f"{h5_path}: {im} has no {component} component") + return RealisationIMs( + stations=_decode(raw_stations), + site={key: psa_group[key][:].astype(np.float64) for key in SITE_FIELDS}, + distances={key: psa_group[key][:].astype(np.float64) for key in DISTANCE_FIELDS}, + psa={c: psa_group[c][:][:, columns].astype(np.float32) for c in COMPONENTS}, + scalars=scalars, + attrs={key: float(np.atleast_1d(h5.attrs[key])[0]) for key in TOLERANCES}, + ) + + +def check_finite(rel: str, ims: RealisationIMs) -> None: + """Raise if any kept IM, site or distance value is NaN or infinite.""" + arrays = { + **{f"pSA/{c}": values for c, values in ims.psa.items()}, + **{f"{im}/{c}": values for im, by_c in ims.scalars.items() for c, values in by_c.items()}, + **ims.site, + **ims.distances, + } + for name, values in arrays.items(): + bad = ~np.isfinite(values) + if bad.any(): + station = ims.stations[np.argwhere(bad)[0][0]] + raise ValueError( + f"{rel}: {name} has {int(bad.sum())} non-finite values, the first at {station}" + ) + + +def check_matches_first_realisation( + rel: str, first_rel: str, ims: RealisationIMs, first: RealisationIMs +) -> None: + """Raise unless stations, site fields and distances equal the rupture's first realisation's.""" + if not np.array_equal(ims.stations, first.stations): + raise ValueError(f"{rel}: station list differs from {first_rel}") + for mine, theirs in ((ims.site, first.site), (ims.distances, first.distances)): + for key, values in mine.items(): + if not np.array_equal(values, theirs[key]): + raise ValueError( + f"{rel}: {key} differs from {first_rel}; site and distance fields " + "must be identical across a rupture's realisations" + ) + + +def build_records(rel: str, ims: RealisationIMs) -> tuple[pd.DataFrame, np.ndarray]: + """IMDB records (one per station and component) and their pSA rows.""" + n = len(ims.stations) + frames = [] + for component in COMPONENTS: + frame = {"rel_id": rel, "site_id": ims.stations, "component": component, "kind": "simulated"} + for im, by_component in ims.scalars.items(): + frame[im] = by_component.get(component, np.full(n, np.nan, dtype=np.float32)) + frames.append(pd.DataFrame(frame)) + records = pd.concat(frames, ignore_index=True) + return records, np.concatenate([ims.psa[c] for c in COMPONENTS]) + + +def build_site_event(rupture: str, ims: RealisationIMs) -> pd.DataFrame: + """site_event rows for one rupture.""" + return pd.DataFrame({"site_id": ims.stations, "event_id": rupture, **ims.distances}) + + +class SiteRegistry: + """The sites added so far, with the vs30/z1pt0/z2pt5 first seen for each.""" + + def __init__(self, stations_input: pd.DataFrame) -> None: + self.stations_input = stations_input + self.seen: pd.DataFrame | None = None + + def new_sites(self, rupture: str, ims: RealisationIMs) -> pd.DataFrame: + """Rows for stations not added yet; raise if a known one disagrees.""" + here = pd.DataFrame(ims.site, index=pd.Index(ims.stations, name="site_id")) + if self.seen is None: + new = here + else: + known = here.index.isin(self.seen.index) + if known.any(): + earlier = self.seen.loc[here.index[known]] + mismatch = (here[known] != earlier).any(axis=1) + if mismatch.any(): + site = mismatch[mismatch].index[0] + raise ValueError( + f"rupture {rupture}: {site} has {here.loc[site].to_dict()}, but an " + f"earlier rupture gave {earlier.loc[site].to_dict()}" + ) + new = here[~known] + missing = new.index.difference(self.stations_input.index) + if len(missing): + raise ValueError( + f"rupture {rupture}: {len(missing)} stations missing from " + f"stations_input.ll, e.g. {list(missing[:5])}" + ) + self.seen = new if self.seen is None else pd.concat([self.seen, new]) + coords = self.stations_input.loc[new.index] + return pd.DataFrame( + { + "site_id": new.index.to_numpy(), + "lat": coords["lat"].to_numpy(), + "lon": coords["lon"].to_numpy(), + "vs30": new["vs30"].to_numpy(), + "z1p0": new["z1pt0"].to_numpy(), + "z2p5": new["z2pt5"].to_numpy(), + "is_real": [is_real_station(site) for site in new.index], + } + ) + + +# ---- realisation.json (geometry helpers as in nz_sim_validation_ingest.py) ---- + + +def fault_geometries(realisation: dict) -> dict[str, Fault]: + return { + name: Fault.from_corners( + np.array([[c["latitude"], c["longitude"], c["depth"]] for c in entry["corners"]]).reshape(-1, 4, 3) + ) + for name, entry in realisation["sources"]["source_geometries"].items() + } + + +def initial_fault_name(causality_tree: dict[str, str | None]) -> str: + return next(name for name, parent in causality_tree.items() if parent is None) + + +def source_wkt(geometries: dict[str, Fault]) -> str: + planes = [plane for f in geometries.values() for plane in f.planes] + return MultiPolygon([Polygon(p.corners[:, [1, 0]]) for p in planes]).wkt + + +def trace_wkt(geometries: dict[str, Fault]) -> str: + planes = [plane for f in geometries.values() for plane in f.planes] + return MultiLineString([LineString(p.corners[:2, [1, 0]]) for p in planes]).wkt + + +def domain_wkt(realisation: dict) -> str: + corners = realisation["domain"]["domain"] + return Polygon([(c["longitude"], c["latitude"]) for c in corners]).wkt + + +def total_magnitude(realisation: dict) -> float: + """The realisation's total magnitude, BoldM: segment moments summed, converted back.""" + magnitudes = realisation["magnitudes"]["magnitudes"] + total_moment = sum(moment.magnitude_to_moment(m, bold_m=True) for m in magnitudes.values()) + return float(moment.moment_to_magnitude(total_moment, bold_m=True)) + + +def build_event(event_id: str, realisation: dict, nshm: pd.Series) -> tuple[dict, dict[str, Fault]]: + """The events row for a rupture, from one of its realisations and its NSHM attributes. + + Geometry is identical across a rupture's realisations, but magnitude and + causality tree are not, so those live on the realisation instead. + """ + geometries = fault_geometries(realisation) + rakes = realisation["rakes"]["rakes"] + tect_type = realisation.get("empirical", {}).get("tect_type", "").upper() + if tect_type not in schema.TECT_TYPES: + raise ValueError(f"rupture {event_id}: tect_type {tect_type!r} is not one of {schema.TECT_TYPES}") + + if len(geometries) == 1: + (f,) = geometries.values() + dip, dip_dir = f.dip, f.dip_dir + dtop, dbottom, length = f.top_m / 1000, f.bottom_m / 1000, f.length + else: + # Cascading multi-fault rupture: no single dip/length fits, so the scalar + # columns stay null and the per-segment breakdown goes in metadata. + dip = dip_dir = dtop = dbottom = length = None + + segments = { + name: { + "dip_deg": f.dip, + "dip_dir_deg": f.dip_dir, + "dtop_km": f.top_m / 1000, + "dbottom_km": f.bottom_m / 1000, + "length_km": f.length, + "rake": rakes[name], + } + for name, f in geometries.items() + } + event = { + "event_id": event_id, + "magnitude": float(nshm["magnitude_boldm"]), + "tect_type": tect_type, + "dip": dip, + "dip_dir": dip_dir, + "dtop": dtop, + "dbottom": dbottom, + "length": length, + "source_wkt": source_wkt(geometries), + "trace_wkt": trace_wkt(geometries), + "domain_wkt": domain_wkt(realisation), + "metadata": json.dumps( + { + "nshm_rupture_name": realisation["metadata"]["name"], + "nshmdb_rupture_id": int(nshm["nshmdb_rupture_id"]), + "annual_rate": float(nshm["annual_rate"]), + "nshm_area_km2": float(nshm["area_km2"]), + "nshm_length_km": float(nshm["length_km"]), + "segments": segments, + } + ), + } + return event, geometries + + +def build_realisation(rupture: str, n: int, geometries: dict[str, Fault], realisation: dict) -> dict: + """The realisations row for `/R`.""" + causality_tree = realisation["rupture_propagation"]["rupture_causality_tree"] + initial_fault = initial_fault_name(causality_tree) + hypocentre_sd = realisation["rupture_propagation"]["hypocentre"] + hypo_lat, hypo_lon, hypo_depth_m = geometries[initial_fault].fault_coordinates_to_wgs_depth_coordinates( + np.array([hypocentre_sd["s"], hypocentre_sd["d"]]) + ) + return { + "rel_id": rel_id(rupture, n), + "event_id": rupture, + "magnitude": total_magnitude(realisation), + "rake": realisation["rakes"]["rakes"][initial_fault], + "hypo_lat": float(hypo_lat), + "hypo_lon": float(hypo_lon), + "hypo_depth": float(hypo_depth_m) / 1000, + "metadata": json.dumps( + { + "seeds": realisation["seeds"], + "segment_magnitudes_boldm": realisation["magnitudes"]["magnitudes"], + "rupture_causality_tree": causality_tree, + "duration_s": realisation["domain"]["duration"], + } + ), + } + + +def check_against_h5(rel: str, realisation_row: dict, attrs: dict[str, float]) -> None: + """Raise unless realisation.json reproduces the h5's magnitude and hypocentre.""" + for key, tolerance in TOLERANCES.items(): + ours, theirs = realisation_row[key], attrs[key] + if not abs(ours - theirs) <= tolerance: + raise ValueError( + f"{rel}: {key} from realisation.json is {ours!r}, but " + f"intensity_measures.h5 says {theirs!r}" + ) + + +# ---- build and verify ---------------------------------------------------------- + + +def check_row_counts(db: IMDB, expected: dict[str, int]) -> None: + """Raise unless every table has exactly the expected number of rows.""" + for table, n in expected.items(): + found = int(db.con.table(table).count().to_pandas()) + if found != n: + raise RuntimeError(f"{table}: {found} rows, expected {n}") + + +def spot_check( + db_path: Path, share_dir: Path, n_realisations: int = 25, per_realisation: int = 40, seed: int = 0 +) -> int: + """Re-read random records from their h5 and compare exactly; return how many.""" + rng = np.random.default_rng(seed) + checked = 0 + with IMDB(db_path) as db: + rels = sorted(db.get_realisations().index) + for rel in rng.choice(rels, size=min(n_realisations, len(rels)), replace=False): + rupture, n = rel.rsplit("_R", 1) + records = db.get_records(rel_ids=[rel]) + sample = records.sample(n=min(per_realisation, len(records)), random_state=seed) + ids = sample.index.tolist() + psa = db.get_psa(record_int_ids=ids) + scalars = db.get_scalars(record_int_ids=ids) + ims = read_ims(share_dir / rupture / f"R{n}" / "intensity_measures.h5") + row_of = {site: i for i, site in enumerate(ims.stations)} + for record_int_id, record in sample.iterrows(): + i, component = row_of[record["site_id"]], record["component"] + stored = psa.loc[record_int_id].to_numpy(dtype=np.float32) + if not np.array_equal(stored, ims.psa[component][i]): + raise AssertionError(f"{rel} {record['site_id']} {component}: pSA differs from the h5") + for im in schema.SCALAR_IMS: + value = scalars.loc[record_int_id, im] + expected = ims.scalars[im].get(component) + if expected is None: + if not pd.isna(value): + raise AssertionError(f"{rel} {record['site_id']} {component}: {im} should be NULL") + elif np.float32(value) != expected[i]: + raise AssertionError(f"{rel} {record['site_id']} {component}: {im} differs from the h5") + checked += 1 + return checked + + +def build_database( + share_dir: Path, + events_dir: Path, + stations_input: Path, + nshm_csv: Path, + realisations: Sequence[tuple[str, int]], + out_path: Path, + source: str, + script_commit: str, + memory_limit: str | None = None, +) -> dict[str, int]: + """Build the database at out_path and return its table row counts. + + Writes `.partial` and renames it only once validation, the row + counts and the spot check have all passed. A failure leaves the .partial + behind for inspection, and the next run refuses to start until it is removed. + """ + if out_path.exists(): + raise FileExistsError(f"{out_path} already exists") + partial = out_path.with_name(out_path.name + ".partial") + if partial.exists(): + raise FileExistsError(f"{partial} exists: a previous run died; remove it first") + + stations = load_stations_input(stations_input) + nshm = load_nshm_attributes(nshm_csv) + by_rupture = group_by_rupture(realisations) + missing = sorted(set(by_rupture) - set(nshm.index), key=int) + if missing: + raise ValueError(f"ruptures missing from {nshm_csv}: {missing}") + (nshmdb_file,) = nshm["nshmdb_file"].unique() + (nshmdb_sha256,) = nshm["nshmdb_sha256"].unique() + + db = IMDB.create( + partial, + periods=list(PSA_PERIODS), + components=COMPONENTS, + db_meta={ + **DB_META, + "source": source, + "n_realisations": str(len(realisations)), + "nshmdb": f"{nshmdb_file} sha256 {nshmdb_sha256}", + "ingest_script_commit": script_commit, + }, + ) + counts = dict.fromkeys(("events", "realisations", "sites", "site_event", "records"), 0) + started = time.monotonic() + try: + if memory_limit: + db.con.raw_sql(f"SET memory_limit = '{memory_limit}'") + sites = SiteRegistry(stations) + for index, (rupture, ns) in enumerate(by_rupture.items(), start=1): + loaded = [] + for n in ns: + rel = rel_id(rupture, n) + realisation = json.loads((events_dir / rupture / f"R{n}" / "realisation.json").read_text()) + h5_path = share_dir / rupture / f"R{n}" / "intensity_measures.h5" + if not h5_path.exists(): + raise FileNotFoundError(f"{rel}: no {h5_path}") + ims = read_ims(h5_path) + check_finite(rel, ims) + loaded.append((n, realisation, ims)) + + first_n, first_realisation, first_ims = loaded[0] + event, geometries = build_event(rupture, first_realisation, nshm.loc[rupture]) + rows = [] + for n, realisation, ims in loaded: + row = build_realisation(rupture, n, geometries, realisation) + check_against_h5(row["rel_id"], row, ims.attrs) + check_matches_first_realisation(row["rel_id"], rel_id(rupture, first_n), ims, first_ims) + rows.append(row) + new_sites = sites.new_sites(rupture, first_ims) + + db.add_events(pd.DataFrame([event])) + db.add_realisations(pd.DataFrame(rows)) + if len(new_sites): + db.add_sites(new_sites) + db.add_site_event(build_site_event(rupture, first_ims)) + for n, _, ims in loaded: + records, psa = build_records(rel_id(rupture, n), ims) + db.add_records(records, pSA=psa) + counts["records"] += len(records) + + counts["events"] += 1 + counts["realisations"] += len(loaded) + counts["sites"] += len(new_sites) + counts["site_event"] += len(first_ims.stations) + if index % 25 == 0 or index == len(by_rupture): + print( + f"{index}/{len(by_rupture)} ruptures, {counts['records']} records, " + f"{time.monotonic() - started:.0f} s", + flush=True, + ) + + problems = db.validate() + if problems: + raise RuntimeError(f"validate(): {problems}") + expected = {**counts, "psa_ims": counts["records"], "scalars_ims": counts["records"]} + check_row_counts(db, expected) + finally: + db.close() + + checked = spot_check(partial, share_dir) + print(f"validate() clean; row counts match; spot check: {checked} records identical to their h5") + partial.rename(out_path) + return expected + + +def main() -> None: + parser = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter + ) + parser.add_argument("share_dir", type=Path, help="run3 share/: /R/intensity_measures.h5") + parser.add_argument("events_dir", type=Path, help="run3 events/: /R/realisation.json") + parser.add_argument("stations_input", type=Path, help="the campaign's stations_input.ll (lon lat name)") + parser.add_argument("nshm_csv", type=Path, help="campaign_nshm_ruptures.csv from export_nshm_rupture_attributes.py") + parser.add_argument("realisations_file", type=Path, help="one /R per line") + parser.add_argument("out", type=Path, help="output .duckdb; must not exist") + parser.add_argument("--source", required=True, help="provenance for db_meta.source") + parser.add_argument("--script-commit", required=True, help="git commit of this script, for db_meta") + parser.add_argument( + "--memory-limit", + help="DuckDB memory_limit, e.g. 12GB. Without it DuckDB sizes itself from " + "the node's RAM, not the job's memory allocation.", + ) + args = parser.parse_args() + + realisations = parse_realisation_ids(args.realisations_file.read_text().splitlines()) + print(f"ingesting {len(realisations)} realisations into {args.out}", flush=True) + counts = build_database( + args.share_dir, + args.events_dir, + args.stations_input, + args.nshm_csv, + realisations, + args.out, + source=args.source, + script_commit=args.script_commit, + memory_limit=args.memory_limit, + ) + print("row counts:", counts) + + +if __name__ == "__main__": + main() diff --git a/scripts/tests/test_cs_nshm_2022_ingest.py b/scripts/tests/test_cs_nshm_2022_ingest.py new file mode 100644 index 0000000..42a8a28 --- /dev/null +++ b/scripts/tests/test_cs_nshm_2022_ingest.py @@ -0,0 +1,469 @@ +"""Tests for scripts/cs_nshm_2022_ingest.py, on a synthetic campaign in run3's layout. + +imdb's own CI installs neither h5py nor source-modelling, so this module skips +there. Run it inside the imdb_tools image: + apptainer exec imdb_tools_.sif python3 -m pytest -q scripts/tests +""" + +import importlib.util +import json +import sys +from dataclasses import dataclass +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + +h5py = pytest.importorskip("h5py") +pytest.importorskip("source_modelling") + +from imdb import IMDB # noqa: E402 +from source_modelling import moment, sources # noqa: E402 + +SCRIPT = Path(__file__).resolve().parents[1] / "cs_nshm_2022_ingest.py" +_spec = importlib.util.spec_from_file_location("cs_nshm_2022_ingest", SCRIPT) +ingest = importlib.util.module_from_spec(_spec) +sys.modules[_spec.name] = ingest # dataclasses resolve their module through sys.modules +_spec.loader.exec_module(ingest) + +# im-calc's period grid, as in every run3 realisation.json `im.valid_periods` +VALID_PERIODS = ( + 0.01, 0.02, 0.022, 0.025, 0.029, 0.03, 0.032, 0.035, 0.036, 0.04, 0.042, 0.044, 0.045, + 0.046, 0.048, 0.05, 0.055, 0.06, 0.065, 0.067, 0.07, 0.075, 0.08, 0.085, 0.09, 0.095, + 0.1, 0.11, 0.12, 0.13, 0.133, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.22, 0.24, 0.25, + 0.26, 0.28, 0.29, 0.3, 0.32, 0.34, 0.35, 0.36, 0.38, 0.4, 0.42, 0.44, 0.45, 0.46, 0.48, + 0.5, 0.55, 0.6, 0.65, 0.667, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, + 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.5, 2.6, 2.8, 3.0, 3.2, 3.4, 3.5, 3.6, 3.8, 4.0, + 4.2, 4.4, 4.6, 4.8, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 11.0, 12.0, + 13.0, 14.0, 15.0, 20.0, +) # fmt: skip +ALL_COMPONENTS = ("000", "090", "ver", "geom", "rotd0", "rotd50", "rotd100") +AS_RECORDED = ("000", "090", "ver", "geom") +H5_GROUPS = { + "pSA": ALL_COMPONENTS, + "PGA": ALL_COMPONENTS, + "PGV": ALL_COMPONENTS, + "PGD": ALL_COMPONENTS, + "CAV": AS_RECORDED, + "AI": AS_RECORDED, + "Ds575": AS_RECORDED, + "Ds595": AS_RECORDED, + "FAS": (*AS_RECORDED, "eas"), +} + +# Canonical coordinates, deliberately not what the h5 files carry. +STATIONS_INPUT = { + "AMBC": (172.1, -43.1), + "2O8aUSL": (172.2, -43.2), + "2O8bSSL": (172.3, -43.3), + "2TfEMSL": (172.4, -43.4), + "2TfFKSL": (172.5, -43.5), +} +# vs30, z1pt0, z2pt5 per station: global, like the real campaign's. +SITE_VALUES = { + "AMBC": (300.0, 0.3, 1.3), + "2O8aUSL": (400.0, 0.2, 1.2), + "2O8bSSL": (500.0, 0.15, 1.1), + "2TfEMSL": (600.0, 0.1, 1.0), + "2TfFKSL": (700.0, 0.05, 0.9), +} + + +def _plane(lat: float, lon: float, strike: float) -> sources.Plane: + return sources.Plane.from_centroid_strike_dip( + np.array([lat, lon, 8.0]), 60.0, 10.0, 6.0, strike=strike + ) + + +def _realisation( + name: str, + faults: dict[str, sources.Plane], + magnitudes: dict[str, float], + tree: dict[str, str | None], + hypocentre: tuple[float, float], + duration: float, +) -> dict: + return { + "metadata": {"name": f"Rupture {name}"}, + "sources": { + "source_geometries": { + fault: { + "type": "fault", + "corners": [ + {"latitude": c[0], "longitude": c[1], "depth": c[2]} for c in plane.corners + ], + } + for fault, plane in faults.items() + } + }, + "magnitudes": {"magnitudes": magnitudes}, + "rakes": {"rakes": dict.fromkeys(faults, 110.0)}, + "rupture_propagation": { + "rupture_causality_tree": tree, + "hypocentre": {"s": hypocentre[0], "d": hypocentre[1]}, + }, + "seeds": {"hf_seed": 7}, + "domain": { + "domain": [ + {"latitude": -43.0, "longitude": 172.0}, + {"latitude": -43.0, "longitude": 173.0}, + {"latitude": -44.0, "longitude": 173.0}, + {"latitude": -44.0, "longitude": 172.0}, + ], + "duration": duration, + "depth": 40.0, + }, + "empirical": {"tect_type": "active_shallow", "models": ["NSHM2022"]}, + } + + +def _boldm_total(realisation: dict) -> float: + magnitudes = realisation["magnitudes"]["magnitudes"].values() + return moment.moment_to_magnitude( + sum(moment.magnitude_to_moment(m, bold_m=True) for m in magnitudes), bold_m=True + ) + + +def _hypocentre(realisation: dict) -> tuple[float, float, float]: + tree = realisation["rupture_propagation"]["rupture_causality_tree"] + initial = next(name for name, parent in tree.items() if parent is None) + entry = realisation["sources"]["source_geometries"][initial] + corners = np.array([[c["latitude"], c["longitude"], c["depth"]] for c in entry["corners"]]) + fault = sources.Fault.from_corners(corners.reshape(-1, 4, 3)) + h = realisation["rupture_propagation"]["hypocentre"] + lat, lon, depth_m = fault.fault_coordinates_to_wgs_depth_coordinates(np.array([h["s"], h["d"]])) + return float(lat), float(lon), float(depth_m) / 1000 + + +def _write_h5(path: Path, stations: list[str], realisation: dict, distances: np.ndarray, seed: int) -> None: + """An intensity_measures.h5 in im-calc's layout, attributes matching `realisation`.""" + rng = np.random.default_rng(seed) + n = len(stations) + hypo_lat, hypo_lon, hypo_depth = _hypocentre(realisation) + path.parent.mkdir(parents=True, exist_ok=True) + with h5py.File(path, "w") as h5: + for key, value in { + "magnitude": _boldm_total(realisation), + "hypo_lat": hypo_lat, + "hypo_lon": hypo_lon, + "hypo_depth": hypo_depth, + "rake": 45.0, + }.items(): + h5.attrs[key] = np.array([value]) + for group_name, components in H5_GROUPS.items(): + group = h5.create_group(group_name) + group.create_dataset("station", data=stations, dtype=h5py.string_dtype()) + group["latitude"] = np.float32([STATIONS_INPUT[s][1] + 0.001 for s in stations]) + group["longitude"] = np.float32([STATIONS_INPUT[s][0] + 0.001 for s in stations]) + for i, key in enumerate(("vs30", "z1pt0", "z2pt5")): + group[key] = np.array([SITE_VALUES[s][i] for s in stations]) + for i, key in enumerate(("rrup", "rjb", "rx", "ry")): + group[key] = distances[:, i] + group["epi"] = rng.uniform(1, 100, n) + group["hyp"] = rng.uniform(1, 100, n) + if group_name == "pSA": + group["period"] = np.array(VALID_PERIODS) + for c in components: + group[c] = rng.uniform(1e-3, 2.0, (n, len(VALID_PERIODS))) + elif group_name == "FAS": + group["frequency"] = np.array([0.1, 1.0, 10.0]) + for c in components: + group[c] = rng.uniform(size=(n, 3)) + else: + for c in components: + group[c] = rng.uniform(0.01, 5.0, n) + + +@dataclass +class Campaign: + share: Path + events: Path + stations_input: Path + nshm_csv: Path + realisations: list[tuple[str, int]] + out: Path + + def h5(self, rel: str) -> Path: + rupture, n = rel.split("/R") + return self.share / rupture / f"R{n}" / "intensity_measures.h5" + + def realisation_json(self, rel: str) -> Path: + rupture, n = rel.split("/R") + return self.events / rupture / f"R{n}" / "realisation.json" + + def build(self) -> dict[str, int]: + return ingest.build_database( + self.share, + self.events, + self.stations_input, + self.nshm_csv, + self.realisations, + self.out, + source="pytest", + script_commit="abc123", + ) + + +RUPTURE_STATIONS = { + "161984": ["AMBC", "2O8aUSL", "2O8bSSL", "2TfEMSL"], + "288271": ["2O8bSSL", "2TfEMSL", "2TfFKSL"], +} +NSHM_MAGNITUDE = {"161984": 7.480455658305042, "288271": 7.13806256610197} + + +@pytest.fixture +def campaign(tmp_path: Path) -> Campaign: + """Two ruptures: 161984 (two faults, R2 and R3) and 288271 (one fault, R2).""" + share, events = tmp_path / "share", tmp_path / "events" + two_faults = {"Clarence": _plane(-43.5, 172.5, 45.0), "Kekerengu": _plane(-43.45, 172.58, 50.0)} + one_fault = {"Monowai": _plane(-43.6, 172.4, 30.0)} + realisations = { + "161984/R2": _realisation( + "161984", two_faults, {"Clarence": 7.11, "Kekerengu": 6.59}, + {"Kekerengu": "Clarence", "Clarence": None}, (0.17, 0.19), 161.7, + ), + "161984/R3": _realisation( + "161984", two_faults, {"Clarence": 6.54, "Kekerengu": 6.02}, + {"Clarence": "Kekerengu", "Kekerengu": None}, (0.59, 0.92), 156.5, + ), + "288271/R2": _realisation( + "288271", one_fault, {"Monowai": 7.32}, {"Monowai": None}, (0.5, 0.6), 120.0, + ), + } # fmt: skip + rng = np.random.default_rng(1) + distances = {rupture: rng.uniform(0, 80, (len(s), 4)) for rupture, s in RUPTURE_STATIONS.items()} + for seed, (rel, realisation) in enumerate(realisations.items()): + rupture = rel.split("/")[0] + path = events / rel / "realisation.json" + path.parent.mkdir(parents=True) + path.write_text(json.dumps(realisation)) + _write_h5(share / rel / "intensity_measures.h5", RUPTURE_STATIONS[rupture], realisation, distances[rupture], seed) + + stations_input = tmp_path / "stations_input.ll" + stations_input.write_text("".join(f"{lon} {lat} {s}\n" for s, (lon, lat) in STATIONS_INPUT.items())) + nshm_csv = tmp_path / "campaign_nshm_ruptures.csv" + pd.DataFrame( + { + "rupture_id": ["607", "161984", "288271"], + "nshmdb_rupture_id": [40083, 201460, 327747], + "magnitude_boldm": [7.47, NSHM_MAGNITUDE["161984"], NSHM_MAGNITUDE["288271"]], + "area_km2": [1872.5, 1907.5, 867.1], + "length_km": [59.9, 80.2, 34.6], + "annual_rate": [6.1e-06, 1.86e-04, 7.97e-05], + "nshmdb_file": "nshmdb_v2026.08.3.db", + "nshmdb_sha256": "3fe692cb", + } + ).to_csv(nshm_csv, index=False) + return Campaign( + share, events, stations_input, nshm_csv, + [("161984", 2), ("161984", 3), ("288271", 2)], tmp_path / "out.duckdb", + ) # fmt: skip + + +# ---- helpers ------------------------------------------------------------------- + + +def test_select_period_indices_picks_the_exact_columns(): + columns = ingest.select_period_indices(np.array(VALID_PERIODS)) + assert tuple(np.array(VALID_PERIODS)[columns]) == ingest.PSA_PERIODS + assert len(ingest.PSA_PERIODS) == 25 + + +def test_select_period_indices_rejects_a_missing_period(): + without = np.array([p for p in VALID_PERIODS if p != 0.3]) + with pytest.raises(ValueError, match="0.3 s matches 0"): + ingest.select_period_indices(without) + + +def test_parse_realisation_ids(): + lines = ["161984/R3", "", "607/R2", "161984/R2", "607/R2"] + assert ingest.parse_realisation_ids(lines) == [("607", 2), ("161984", 2), ("161984", 3)] + with pytest.raises(ValueError, match="161984_R2"): + ingest.parse_realisation_ids(["161984_R2"]) + + +def test_is_real_station(): + assert ingest.is_real_station("AMBC") + assert ingest.is_real_station("BFZ") + assert not ingest.is_real_station("2O8aUSL") + + +# ---- the whole build ------------------------------------------------------------- + + +def test_row_counts_and_meta(campaign: Campaign): + counts = campaign.build() + assert counts == { + "events": 2, + "realisations": 3, + "sites": 5, + "site_event": 7, + "records": 22, + "psa_ims": 22, + "scalars_ims": 22, + } + assert campaign.out.exists() + assert not campaign.out.with_name("out.duckdb.partial").exists() + with IMDB(campaign.out) as db: + meta = db.db_meta + assert db.validate() == [] + assert meta["dataset_id"] == "cs_nshm_2022_run3" + assert meta["components"] == "geom,rotd50" + assert meta["n_periods"] == "25" + assert meta["n_frequencies"] == "0" + assert meta["magnitude_convention"].startswith("BoldM") + assert meta["nshmdb"] == "nshmdb_v2026.08.3.db sha256 3fe692cb" + assert meta["ingest_script_commit"] == "abc123" + assert meta["source"] == "pytest" + assert meta["n_realisations"] == "3" + + +def test_events(campaign: Campaign): + campaign.build() + with IMDB(campaign.out) as db: + events = db.get_events() + assert events.loc["161984", "magnitude"] == pytest.approx(NSHM_MAGNITUDE["161984"], abs=1e-6) + assert events.loc["288271", "magnitude"] == pytest.approx(NSHM_MAGNITUDE["288271"], abs=1e-6) + assert (events["tect_type"] == "ACTIVE_SHALLOW").all() + assert pd.isna(events.loc["161984", "dip"]) # multi-fault: no single dip + assert events.loc["288271", "dip"] == pytest.approx(60.0, abs=1e-3) + assert events.loc["288271", "length"] == pytest.approx(10.0, abs=1e-3) + metadata = json.loads(events.loc["161984", "metadata"]) + assert metadata["nshmdb_rupture_id"] == 201460 + assert metadata["annual_rate"] == 1.86e-04 + assert set(metadata["segments"]) == {"Clarence", "Kekerengu"} + assert events.loc["161984", "source_wkt"].startswith("MULTIPOLYGON") + + +def test_realisations_are_boldm_and_match_their_h5(campaign: Campaign): + campaign.build() + with IMDB(campaign.out) as db: + realisations = db.get_realisations() + for rel in ("161984/R2", "161984/R3", "288271/R2"): + realisation = json.loads(campaign.realisation_json(rel).read_text()) + row = realisations.loc[rel.replace("/", "_")] + assert row["magnitude"] == pytest.approx(_boldm_total(realisation), abs=1e-5) + assert (row["hypo_lat"], row["hypo_lon"], row["hypo_depth"]) == pytest.approx( + _hypocentre(realisation), abs=1e-4 + ) + metadata = json.loads(row["metadata"]) + assert metadata["duration_s"] == realisation["domain"]["duration"] + # R3 starts on the other fault, so its causality tree is its own + assert json.loads(realisations.loc["161984_R3", "metadata"])["rupture_causality_tree"]["Kekerengu"] is None + + +def test_sites_use_canonical_coordinates(campaign: Campaign): + campaign.build() + with IMDB(campaign.out) as db: + sites = db.get_sites() + for site, (lon, lat) in STATIONS_INPUT.items(): + assert (sites.loc[site, "lat"], sites.loc[site, "lon"]) == pytest.approx((lat, lon), abs=1e-5) + assert (sites.loc[site, "vs30"], sites.loc[site, "z1p0"], sites.loc[site, "z2p5"]) == pytest.approx( + SITE_VALUES[site], abs=1e-5 + ) + assert sites["is_real"].to_dict() == {s: s == "AMBC" for s in STATIONS_INPUT} + + +def test_site_event_comes_from_the_h5(campaign: Campaign): + campaign.build() + with h5py.File(campaign.h5("288271/R2")) as h5: + rrup = h5["pSA/rrup"][:] + with IMDB(campaign.out) as db: + site_event = db.get_site_event(event_ids=["288271"]).set_index("site_id") + assert site_event.loc[RUPTURE_STATIONS["288271"], "rrup"].to_numpy() == pytest.approx(rrup, abs=1e-4) + + +def test_ims_are_the_h5_values(campaign: Campaign): + campaign.build() + with h5py.File(campaign.h5("161984/R3")) as h5: + columns = ingest.select_period_indices(h5["pSA/period"][:]) + station = RUPTURE_STATIONS["161984"].index("2TfEMSL") + psa_rotd50 = h5["pSA/rotd50"][station, columns].astype(np.float32) + cav_geom = np.float32(h5["CAV/geom"][station]) + pga_rotd50 = np.float32(h5["PGA/rotd50"][station]) + with IMDB(campaign.out) as db: + rotd50 = db.get_records(rel_ids=["161984_R3"], site_ids=["2TfEMSL"], component="rotd50").index + geom = db.get_records(rel_ids=["161984_R3"], site_ids=["2TfEMSL"], component="geom").index + psa = db.get_psa(record_int_ids=list(rotd50)) + scalars = db.get_scalars(record_int_ids=[*rotd50, *geom]) + np.testing.assert_array_equal(psa.iloc[0].to_numpy(dtype=np.float32), psa_rotd50) + assert list(psa.columns) == [str(p) for p in ingest.PSA_PERIODS] + assert np.float32(scalars.loc[rotd50[0], "PGA"]) == pga_rotd50 + assert pd.isna(scalars.loc[rotd50[0], "CAV"]) # undefined on rotd + assert np.float32(scalars.loc[geom[0], "CAV"]) == cav_geom + + +# ---- refusals and consistency checks ------------------------------------------- + + +def test_refuses_to_overwrite(campaign: Campaign): + campaign.out.write_text("") + with pytest.raises(FileExistsError): + campaign.build() + + +def test_refuses_a_leftover_partial(campaign: Campaign): + campaign.out.with_name("out.duckdb.partial").write_text("") + with pytest.raises(FileExistsError, match="previous run died"): + campaign.build() + + +KEPT_COLUMN = VALID_PERIODS.index(0.1) # 26: the first period the database keeps +DROPPED_COLUMN = VALID_PERIODS.index(0.025) # a period the database drops + + +def _edit_h5(path: Path, dataset: str, index, value, groups: tuple[str, ...] | None = None) -> None: + """Overwrite one element of `dataset` in `groups` (default: every group that has it).""" + with h5py.File(path, "r+") as h5: + for name, group in h5.items(): + if dataset in group and (groups is None or name in groups): + data = group[dataset][...] + data[index] = value + group[dataset][...] = data + + +def _set_attr(path: Path, key: str, delta: float) -> None: + with h5py.File(path, "r+") as h5: + h5.attrs[key] = h5.attrs[key] + delta + + +@pytest.mark.parametrize( + "corrupt, error, match", + [ + (lambda c: _edit_h5(c.h5("288271/R2"), "geom", (0, KEPT_COLUMN), np.nan, ("pSA",)), ValueError, "pSA/geom has 1 non-finite"), + (lambda c: _edit_h5(c.h5("288271/R2"), "vs30", 1, 999.0), ValueError, "earlier rupture"), + (lambda c: _edit_h5(c.h5("161984/R3"), "rrup", 2, 1.5), ValueError, "rrup differs from 161984_R2"), + (lambda c: _set_attr(c.h5("161984/R3"), "magnitude", 0.01), ValueError, "magnitude from realisation.json"), + (lambda c: _set_attr(c.h5("288271/R2"), "hypo_lat", 0.001), ValueError, "hypo_lat from realisation.json"), + (lambda c: c.stations_input.write_text(c.stations_input.read_text().replace("2TfFKSL", "XXXXXXX")), + ValueError, "missing from stations_input.ll"), + (lambda c: c.h5("288271/R2").unlink(), FileNotFoundError, "288271_R2"), + ], + ids=["nan", "site-mismatch", "distance-mismatch", "magnitude", "hypocentre", "unknown-station", "no-h5"], +) # fmt: skip +def test_bad_input_stops_the_build(campaign: Campaign, corrupt, error, match): + corrupt(campaign) + with pytest.raises(error, match=match): + campaign.build() + assert not campaign.out.exists() + + +def test_unknown_tect_type_stops_the_build(campaign: Campaign): + path = campaign.realisation_json("288271/R2") + realisation = json.loads(path.read_text()) + realisation["empirical"]["tect_type"] = "subduction_mystery" + path.write_text(json.dumps(realisation)) + with pytest.raises(ValueError, match="tect_type"): + campaign.build() + + +def test_rupture_missing_from_the_nshm_csv(campaign: Campaign): + pd.read_csv(campaign.nshm_csv).iloc[:2].to_csv(campaign.nshm_csv, index=False) + with pytest.raises(ValueError, match=r"missing from .*\['288271'\]"): + campaign.build() + + +def test_nan_in_a_dropped_period_is_ignored(campaign: Campaign): + _edit_h5(campaign.h5("288271/R2"), "geom", (0, DROPPED_COLUMN), np.nan, ("pSA",)) + assert campaign.build()["records"] == 22 From 33a9746d44ebccc950cf5672fa8a2425bc853893 Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Wed, 23 Sep 2026 20:22:32 +1200 Subject: [PATCH 2/8] Ingest the pilot and the main campaign into one cs_nshm_2022 IMDB A manifest now lists every realisation with its rel_id, pilot flag and files, so the script knows neither campaign's layout. An event takes its geometry from its first main realisation and its domain from the union of all of them. site_event keeps the main campaign's distances wherever it has them and the pilot's elsewhere, flagged {"fault_geometry": "pilot"}. Each hypocentre is placed on its own realisation's geometry, which differs for the pilot. db_meta now speaks to researchers, not to our run layout. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UgVr1htJMAqx8tABJLTnD5 --- scripts/cs_nshm_2022_ingest.py | 507 +++++++++++++++++----- scripts/tests/test_cs_nshm_2022_ingest.py | 398 ++++++++++++----- 2 files changed, 670 insertions(+), 235 deletions(-) diff --git a/scripts/cs_nshm_2022_ingest.py b/scripts/cs_nshm_2022_ingest.py index 653b2b4..7029458 100644 --- a/scripts/cs_nshm_2022_ingest.py +++ b/scripts/cs_nshm_2022_ingest.py @@ -14,12 +14,22 @@ # python-preference = "only-managed" # avoid a broken/non-standard system Python on PATH # /// -"""Ingest the cs_nshm_2022 run3 simulation campaign into an IMDB. +"""Ingest the cs_nshm_2022 simulations into one IMDB. -One event per NSHM 2022 crustal rupture and one realisation per simulated -`R`. The event_id is the campaign's rupture id, which is the rupture's -crustal `nshm_id` in nshmdb, NOT its `rupture_id`; the rel_id is -`_R`. +Two simulation campaigns become one dataset. Each event is one NSHM 2022 +crustal rupture and each realisation is one simulation of it: + +- the main campaign: up to two realisations per rupture, each with its own + randomly drawn magnitude, hypocentre, rupture propagation and slip; +- the pilot: one realisation for each of 293 ruptures. Its magnitude is the + scaling relation's central value, and its fault geometry comes from an + earlier NSHM fault database release. + +A manifest CSV (cs_nshm_2022's scripts/build_imdb_manifest.py) lists every +realisation with its rel_id, pilot flag and file paths, so this script knows +nothing about either campaign's layout. Each event's realisations are R1..Rn +with no gaps and the pilot's last. The event_id is the rupture's crustal +`nshm_id` in nshmdb, NOT its `rupture_id`. - IMs, distances, vs30 and z1.0/z2.5 come from each realisation's `intensity_measures.h5`, as written by the workflow's im-calc: netCDF4, one @@ -29,6 +39,10 @@ nshmdb, which is not available where this runs. - Canonical site coordinates come from the campaign's `stations_input.ll`. +An event's geometry, and the distances stored for it, come from its first main +realisation, or from the pilot where the main campaign has none. At a site only +the pilot covers, site_event holds the pilot's distances, flagged in its metadata. + Trimmed for size: only the geom and rotd50 components, pSA on 25 of im-calc's 111 periods, and no FAS or empirical-GMM records. @@ -37,10 +51,11 @@ """ import argparse +import csv +import hashlib import json -import re import time -from collections.abc import Iterable, Sequence +from collections.abc import Sequence from dataclasses import dataclass from pathlib import Path @@ -48,6 +63,7 @@ import numpy as np import pandas as pd from imdb import IMDB, schema +from shapely import unary_union from shapely.geometry import LineString, MultiLineString, MultiPolygon, Polygon from source_modelling import moment from source_modelling.sources import Fault @@ -66,9 +82,55 @@ # over all 221 realisations done on 2026-09-21: magnitude identical, hypocentre # within 3e-14 degrees. TOLERANCES = {"magnitude": 1e-9, "hypo_lat": 1e-9, "hypo_lon": 1e-9, "hypo_depth": 1e-6} +MANIFEST_FIELDS = ("rel_id", "event_id", "realisation", "pilot", "original_id", "im_path", "realisation_path") +# site_event.metadata on a row whose distances come from the pilot's geometry +# although the event's geometry is the main campaign's. +PILOT_DISTANCES = json.dumps({"fault_geometry": "pilot"}) +# The NSHM fault database releases behind each campaign's geometry. +MAIN_RELEASE = "v2026.08.3" +MAIN_NSHMDB_SHA256 = "3fe692cb9b22c769b6a9baa6a526b4eed989bd61ec054a344f0e74d82855bc4e" +PILOT_RELEASE = "pre-v2026.08" +PILOT_NSHMDB_SHA256 = "00e256480618cd15e11fbf744037d037bf3fc2d523fb977ee30e0b84a640bc57" +N_TOP_PEAKS = 20 DB_META = { - "dataset_id": "cs_nshm_2022_run3", + "dataset_id": "cs_nshm_2022", + "description": ( + "Physics-based ground-motion simulations of New Zealand NSHM 2022 crustal " + "ruptures. Each event is one rupture and each realisation is one simulation of " + "it, with its own randomly drawn magnitude, hypocentre, rupture propagation and " + "slip. Most realisations come from the main campaign, up to two per rupture; " + "those from an earlier pilot campaign are marked realisations.metadata.pilot = " + "true." + ), + "realisation_numbering": ( + "Each event's realisations are R1..Rn with no gaps. Where an event has a pilot " + "realisation, it is always the last. The seeds in realisations.metadata identify " + "each simulation uniquely." + ), + "pilot_realisations": ( + "realisations.metadata.pilot = true marks the pilot's realisations, at most one " + "per event. Their magnitude is the magnitude-area scaling relation's central value " + "rather than a random draw, so leave them out of between-realisation variability " + "estimates. Their fault geometry comes from an earlier release of the NSHM fault " + "database: the same trace and depths, with the bottom edge placed differently (by " + "a median of 226 m, at most 1.8 km)." + ), + "distances": ( + "site_event rrup, rjb, rx and ry are computed from the event's own fault geometry " + "(events.source_wkt): the main campaign's where the event has a main-campaign " + "realisation, else the pilot's. At a site that only the pilot's realisation " + f"covers, they come from the pilot's geometry and site_event.metadata is " + f"{PILOT_DISTANCES}. Epicentral and hypocentral distances are not stored; derive " + "them from realisations.hypo_* and the site coordinates." + ), + "nshm_fault_database": ( + "events.metadata.fault_geometry_release names the NSHM2022DB release the event's " + f"geometry comes from: {MAIN_RELEASE} (sha256 {MAIN_NSHMDB_SHA256}), or " + f"{PILOT_RELEASE}, the earlier release the pilot used (sha256 " + f"{PILOT_NSHMDB_SHA256}). events.magnitude and the NSHM attributes in " + f"events.metadata come from {MAIN_RELEASE} for every event." + ), "magnitude_convention": ( "BoldM (Hanks & Kanamori 1979 eq. 7). events.magnitude is the NSHM 2022 " "catalogue magnitude; realisations.magnitude is the realisation's own " @@ -108,35 +170,90 @@ def select_period_indices( return np.array(indices) -def parse_realisation_ids(lines: Iterable[str]) -> list[tuple[str, int]]: - """Parse `/R` lines into (rupture, n) pairs, in ingest order. - - Blank lines are skipped and duplicates dropped; the order is by integer - rupture id, then realisation. - """ - ids = set() - for line in lines: - line = line.strip() - if not line: - continue - match = re.fullmatch(r"(\d+)/R(\d+)", line) - if match is None: - raise ValueError(f"not a /R id: {line!r}") - ids.add((str(int(match[1])), int(match[2]))) - return sorted(ids, key=lambda rel: (int(rel[0]), rel[1])) - - -def group_by_rupture(realisations: Sequence[tuple[str, int]]) -> dict[str, list[int]]: - """Realisation numbers per rupture, keeping the input order.""" - grouped: dict[str, list[int]] = {} - for rupture, n in realisations: - grouped.setdefault(rupture, []).append(n) +def rel_id(event_id: str, n: int) -> str: + """The realisation's stable id.""" + return f"{event_id}_R{n}" + + +@dataclass(frozen=True) +class ManifestRow: + """One realisation to ingest, as listed in the build manifest.""" + + rel_id: str + event_id: str + realisation: int + pilot: bool + original_id: str + """Where the result came from; for messages only, never stored.""" + im_path: Path + realisation_path: Path + + @property + def label(self) -> str: + return f"{self.rel_id} ({self.original_id})" + + +def read_manifest(path: Path) -> list[ManifestRow]: + """The manifest's rows in ingest order: events by integer id, then realisations.""" + with path.open(newline="") as f: + reader = csv.DictReader(f) + missing = set(MANIFEST_FIELDS) - set(reader.fieldnames or ()) + if missing: + raise ValueError(f"{path}: missing columns {sorted(missing)}") + rows = [] + for raw in reader: + if raw["pilot"] not in ("true", "false"): + raise ValueError(f"{raw['rel_id']}: pilot must be true or false, not {raw['pilot']!r}") + rows.append( + ManifestRow( + rel_id=raw["rel_id"], + event_id=raw["event_id"], + realisation=int(raw["realisation"]), + pilot=raw["pilot"] == "true", + original_id=raw["original_id"], + im_path=Path(raw["im_path"]), + realisation_path=Path(raw["realisation_path"]), + ) + ) + check_manifest(rows) + return sorted(rows, key=lambda row: (int(row.event_id), row.realisation)) + + +def check_manifest(rows: Sequence[ManifestRow]) -> None: + """Raise unless ids are consistent and unique, and every event is R1..Rn with at most one pilot, last.""" + seen: set[str] = set() + by_event: dict[str, list[ManifestRow]] = {} + for row in rows: + if not row.event_id.isdigit() or str(int(row.event_id)) != row.event_id: + raise ValueError(f"{row.rel_id}: event_id {row.event_id!r} is not an unpadded integer") + if row.rel_id != rel_id(row.event_id, row.realisation): + raise ValueError(f"{row.rel_id}: expected {rel_id(row.event_id, row.realisation)} from event_id and realisation") + if row.rel_id in seen: + raise ValueError(f"{row.rel_id}: listed twice") + seen.add(row.rel_id) + by_event.setdefault(row.event_id, []).append(row) + for event_id, group in by_event.items(): + numbers = sorted(row.realisation for row in group) + if numbers != list(range(1, len(group) + 1)): + raise ValueError(f"event {event_id}: realisations {numbers}, expected 1..{len(group)}") + pilots = [row for row in group if row.pilot] + if len(pilots) > 1: + raise ValueError(f"event {event_id}: {len(pilots)} pilot realisations, expected at most 1") + if pilots and pilots[0].realisation != len(group): + raise ValueError(f"event {event_id}: the pilot is R{pilots[0].realisation}, but it must be the last, R{len(group)}") + + +def group_by_event(rows: Sequence[ManifestRow]) -> dict[str, list[ManifestRow]]: + """Rows per event, keeping the input order.""" + grouped: dict[str, list[ManifestRow]] = {} + for row in rows: + grouped.setdefault(row.event_id, []).append(row) return grouped -def rel_id(rupture: str, n: int) -> str: - """The realisation's stable id.""" - return f"{rupture}_R{n}" +def md5_of(path: Path) -> str: + """Hex md5 of a file.""" + return hashlib.md5(path.read_bytes()).hexdigest() def load_stations_input(path: Path) -> pd.DataFrame: @@ -151,8 +268,17 @@ def load_stations_input(path: Path) -> pd.DataFrame: def load_nshm_attributes(path: Path) -> pd.DataFrame: - """Rows of export_nshm_rupture_attributes.py's CSV, indexed by rupture id.""" - return pd.read_csv(path, dtype={"rupture_id": str}).set_index("rupture_id") + """Rows of export_nshm_rupture_attributes.py's CSV, indexed by rupture id. + + Raises unless they all come from the release DB_META documents. + """ + nshm = pd.read_csv(path, dtype={"rupture_id": str}).set_index("rupture_id") + if set(nshm["nshmdb_sha256"]) != {MAIN_NSHMDB_SHA256}: + raise ValueError( + f"{path}: NSHM attributes must all come from {MAIN_RELEASE} " + f"(sha256 {MAIN_NSHMDB_SHA256}), found {sorted(set(nshm['nshmdb_sha256']))}" + ) + return nshm def is_real_station(site_id: str) -> bool: @@ -179,12 +305,20 @@ class RealisationIMs: """Scalar IM -> component -> (n_stations,) float32; rotd-undefined IMs omitted.""" attrs: dict[str, float] """The file's magnitude and hypocentre attributes (TOLERANCES keys).""" + peaks: dict[str, tuple[float, str]] + """PGA and PGV (geom) -> (largest value, its station), whatever the database keeps.""" def _decode(values: np.ndarray) -> np.ndarray: return np.array([v.decode() if isinstance(v, bytes) else v for v in values], dtype=object) +def _peak(values: np.ndarray, stations: np.ndarray) -> tuple[float, str]: + i = int(np.nanargmax(values)) + station = stations[i] + return float(values[i]), station.decode() if isinstance(station, bytes) else str(station) + + def read_ims(h5_path: Path) -> RealisationIMs: """Read the parts of one intensity_measures.h5 this database keeps.""" with h5py.File(h5_path, "r") as h5: @@ -213,6 +347,7 @@ def read_ims(h5_path: Path) -> RealisationIMs: psa={c: psa_group[c][:][:, columns].astype(np.float32) for c in COMPONENTS}, scalars=scalars, attrs={key: float(np.atleast_1d(h5.attrs[key])[0]) for key in TOLERANCES}, + peaks={im: _peak(h5[im]["geom"][:], raw_stations) for im in ("PGA", "PGV")}, ) @@ -236,7 +371,7 @@ def check_finite(rel: str, ims: RealisationIMs) -> None: def check_matches_first_realisation( rel: str, first_rel: str, ims: RealisationIMs, first: RealisationIMs ) -> None: - """Raise unless stations, site fields and distances equal the rupture's first realisation's.""" + """Raise unless stations, site fields and distances equal the event's first main realisation's.""" if not np.array_equal(ims.stations, first.stations): raise ValueError(f"{rel}: station list differs from {first_rel}") for mine, theirs in ((ims.site, first.site), (ims.distances, first.distances)): @@ -244,7 +379,7 @@ def check_matches_first_realisation( if not np.array_equal(values, theirs[key]): raise ValueError( f"{rel}: {key} differs from {first_rel}; site and distance fields " - "must be identical across a rupture's realisations" + "must be identical across an event's main realisations" ) @@ -261,9 +396,36 @@ def build_records(rel: str, ims: RealisationIMs) -> tuple[pd.DataFrame, np.ndarr return records, np.concatenate([ims.psa[c] for c in COMPONENTS]) -def build_site_event(rupture: str, ims: RealisationIMs) -> pd.DataFrame: - """site_event rows for one rupture.""" - return pd.DataFrame({"site_id": ims.stations, "event_id": rupture, **ims.distances}) +def _site_event_rows(event_id: str, stations: np.ndarray, distances: dict[str, np.ndarray], metadata: str | None) -> pd.DataFrame: + return pd.DataFrame({"site_id": stations, "event_id": event_id, **distances, "metadata": metadata}) + + +def build_site_event( + event_id: str, main: RealisationIMs | None, pilot: RealisationIMs | None +) -> tuple[pd.DataFrame, int]: + """site_event rows for one event, and how many carry the pilot's distances. + + The main campaign's distances wherever its realisations cover the site: they + all share one station list and one set of distances, so the first covers + them all. The pilot's only where nothing else exists, flagged when the + event's geometry is the main campaign's and unflagged when the event is the + pilot's alone. + """ + if main is None: + if pilot is None: + raise ValueError(f"event {event_id}: no realisations") + return _site_event_rows(event_id, pilot.stations, pilot.distances, None), 0 + rows = _site_event_rows(event_id, main.stations, main.distances, None) + if pilot is None: + return rows, 0 + gap = ~pd.Index(pilot.stations).isin(main.stations) + flagged = _site_event_rows( + event_id, + pilot.stations[gap], + {key: values[gap] for key, values in pilot.distances.items()}, + PILOT_DISTANCES, + ) + return pd.concat([rows, flagged], ignore_index=True), int(gap.sum()) class SiteRegistry: @@ -273,7 +435,7 @@ def __init__(self, stations_input: pd.DataFrame) -> None: self.stations_input = stations_input self.seen: pd.DataFrame | None = None - def new_sites(self, rupture: str, ims: RealisationIMs) -> pd.DataFrame: + def new_sites(self, label: str, ims: RealisationIMs) -> pd.DataFrame: """Rows for stations not added yet; raise if a known one disagrees.""" here = pd.DataFrame(ims.site, index=pd.Index(ims.stations, name="site_id")) if self.seen is None: @@ -286,14 +448,14 @@ def new_sites(self, rupture: str, ims: RealisationIMs) -> pd.DataFrame: if mismatch.any(): site = mismatch[mismatch].index[0] raise ValueError( - f"rupture {rupture}: {site} has {here.loc[site].to_dict()}, but an " - f"earlier rupture gave {earlier.loc[site].to_dict()}" + f"{label}: {site} has {here.loc[site].to_dict()}, but an " + f"earlier realisation gave {earlier.loc[site].to_dict()}" ) new = here[~known] missing = new.index.difference(self.stations_input.index) if len(missing): raise ValueError( - f"rupture {rupture}: {len(missing)} stations missing from " + f"{label}: {len(missing)} stations missing from " f"stations_input.ll, e.g. {list(missing[:5])}" ) self.seen = new if self.seen is None else pd.concat([self.seen, new]) @@ -337,9 +499,14 @@ def trace_wkt(geometries: dict[str, Fault]) -> str: return MultiLineString([LineString(p.corners[:2, [1, 0]]) for p in planes]).wkt -def domain_wkt(realisation: dict) -> str: +def domain_polygon(realisation: dict) -> Polygon: corners = realisation["domain"]["domain"] - return Polygon([(c["longitude"], c["latitude"]) for c in corners]).wkt + return Polygon([(c["longitude"], c["latitude"]) for c in corners]) + + +def domain_wkt(realisations: Sequence[dict]) -> str: + """The union of the realisations' simulation domains, so it contains every site with data.""" + return unary_union([domain_polygon(r) for r in realisations]).wkt def total_magnitude(realisation: dict) -> float: @@ -349,15 +516,19 @@ def total_magnitude(realisation: dict) -> float: return float(moment.moment_to_magnitude(total_moment, bold_m=True)) -def build_event(event_id: str, realisation: dict, nshm: pd.Series) -> tuple[dict, dict[str, Fault]]: - """The events row for a rupture, from one of its realisations and its NSHM attributes. +def build_event( + event_id: str, defining: dict, realisations: Sequence[dict], nshm: pd.Series, release: str +) -> dict: + """The events row for a rupture. - Geometry is identical across a rupture's realisations, but magnitude and - causality tree are not, so those live on the realisation instead. + Geometry comes from `defining`, the realisation whose geometry the event's + unflagged distances share; the domain is the union over `realisations`. + Magnitude and causality tree differ between realisations, so those live on + the realisation instead. """ - geometries = fault_geometries(realisation) - rakes = realisation["rakes"]["rakes"] - tect_type = realisation.get("empirical", {}).get("tect_type", "").upper() + geometries = fault_geometries(defining) + rakes = defining["rakes"]["rakes"] + tect_type = defining.get("empirical", {}).get("tect_type", "").upper() if tect_type not in schema.TECT_TYPES: raise ValueError(f"rupture {event_id}: tect_type {tect_type!r} is not one of {schema.TECT_TYPES}") @@ -381,7 +552,7 @@ def build_event(event_id: str, realisation: dict, nshm: pd.Series) -> tuple[dict } for name, f in geometries.items() } - event = { + return { "event_id": event_id, "magnitude": float(nshm["magnitude_boldm"]), "tect_type": tect_type, @@ -392,23 +563,24 @@ def build_event(event_id: str, realisation: dict, nshm: pd.Series) -> tuple[dict "length": length, "source_wkt": source_wkt(geometries), "trace_wkt": trace_wkt(geometries), - "domain_wkt": domain_wkt(realisation), + "domain_wkt": domain_wkt(realisations), "metadata": json.dumps( { - "nshm_rupture_name": realisation["metadata"]["name"], + "nshm_rupture_name": defining["metadata"]["name"], "nshmdb_rupture_id": int(nshm["nshmdb_rupture_id"]), "annual_rate": float(nshm["annual_rate"]), "nshm_area_km2": float(nshm["area_km2"]), "nshm_length_km": float(nshm["length_km"]), + "fault_geometry_release": release, "segments": segments, } ), } - return event, geometries -def build_realisation(rupture: str, n: int, geometries: dict[str, Fault], realisation: dict) -> dict: - """The realisations row for `/R`.""" +def build_realisation(row: ManifestRow, realisation: dict) -> dict: + """The realisations row, from the realisation's own realisation.json and geometry.""" + geometries = fault_geometries(realisation) causality_tree = realisation["rupture_propagation"]["rupture_causality_tree"] initial_fault = initial_fault_name(causality_tree) hypocentre_sd = realisation["rupture_propagation"]["hypocentre"] @@ -416,8 +588,8 @@ def build_realisation(rupture: str, n: int, geometries: dict[str, Fault], realis np.array([hypocentre_sd["s"], hypocentre_sd["d"]]) ) return { - "rel_id": rel_id(rupture, n), - "event_id": rupture, + "rel_id": row.rel_id, + "event_id": row.event_id, "magnitude": total_magnitude(realisation), "rake": realisation["rakes"]["rakes"][initial_fault], "hypo_lat": float(hypo_lat), @@ -425,6 +597,7 @@ def build_realisation(rupture: str, n: int, geometries: dict[str, Fault], realis "hypo_depth": float(hypo_depth_m) / 1000, "metadata": json.dumps( { + "pilot": row.pilot, "seeds": realisation["seeds"], "segment_magnitudes_boldm": realisation["magnitudes"]["magnitudes"], "rupture_causality_tree": causality_tree, @@ -448,6 +621,15 @@ def check_against_h5(rel: str, realisation_row: dict, attrs: dict[str, float]) - # ---- build and verify ---------------------------------------------------------- +@dataclass +class Loaded: + """One realisation read from disk.""" + + row: ManifestRow + realisation: dict + ims: RealisationIMs + + def check_row_counts(db: IMDB, expected: dict[str, int]) -> None: """Raise unless every table has exactly the expected number of rows.""" for table, n in expected.items(): @@ -456,22 +638,72 @@ def check_row_counts(db: IMDB, expected: dict[str, int]) -> None: raise RuntimeError(f"{table}: {found} rows, expected {n}") +# Counts the flagged site_event rows three ways, from the tables alone: the rows +# flagged; the (site, event) pairs where only a pilot realisation has records in +# an event that also has main realisations; and the overlap of the two. +FLAG_QUERY = """ +WITH rel AS ( + SELECT rel_int_id, event_int_id, + coalesce(json_extract_string(metadata, '$.pilot') = 'true', false) AS pilot + FROM realisations +), site_rel AS ( + SELECT DISTINCT records.site_int_id, records.event_int_id, rel.pilot + FROM records JOIN rel USING (rel_int_id) +), pilot_only AS ( + SELECT site_int_id, event_int_id FROM site_rel + GROUP BY site_int_id, event_int_id HAVING bool_and(pilot) +), shared AS ( + SELECT event_int_id FROM rel GROUP BY event_int_id HAVING bool_or(pilot) AND bool_or(NOT pilot) +) +SELECT + (SELECT count(*) FROM site_event WHERE metadata IS NOT NULL), + (SELECT count(*) FROM pilot_only WHERE event_int_id IN (SELECT event_int_id FROM shared)), + (SELECT count(*) FROM site_event JOIN pilot_only USING (site_int_id, event_int_id) + WHERE site_event.metadata IS NOT NULL + AND site_event.event_int_id IN (SELECT event_int_id FROM shared)) +""" + + +def check_flags(db: IMDB, n_flagged: int) -> None: + """Raise unless the flagged site_event rows are exactly the pilot-only sites of shared events.""" + flagged, pilot_only, both = db.con.raw_sql(FLAG_QUERY).fetchone() + if not flagged == pilot_only == both == n_flagged: + raise RuntimeError( + f"site_event flags: {flagged} flagged rows, {pilot_only} pilot-only sites in " + f"shared events, {both} in both, {n_flagged} written" + ) + + def spot_check( - db_path: Path, share_dir: Path, n_realisations: int = 25, per_realisation: int = 40, seed: int = 0 + db_path: Path, + manifest: dict[str, ManifestRow], + n_realisations: int = 25, + per_realisation: int = 40, + min_pilot: int = 5, + seed: int = 0, ) -> int: - """Re-read random records from their h5 and compare exactly; return how many.""" + """Re-read random records from their h5 and compare exactly; return how many. + + At least `min_pilot` of the sampled realisations are pilot ones, where there are any. + """ rng = np.random.default_rng(seed) + pilot_rels = sorted(rel for rel, row in manifest.items() if row.pilot) + main_rels = sorted(rel for rel, row in manifest.items() if not row.pilot) + n_pilot = min(len(pilot_rels), max(min_pilot, n_realisations - len(main_rels))) + n_main = min(len(main_rels), n_realisations - n_pilot) + chosen = [ + *rng.choice(pilot_rels, size=n_pilot, replace=False), + *rng.choice(main_rels, size=n_main, replace=False), + ] checked = 0 with IMDB(db_path) as db: - rels = sorted(db.get_realisations().index) - for rel in rng.choice(rels, size=min(n_realisations, len(rels)), replace=False): - rupture, n = rel.rsplit("_R", 1) + for rel in chosen: records = db.get_records(rel_ids=[rel]) sample = records.sample(n=min(per_realisation, len(records)), random_state=seed) ids = sample.index.tolist() psa = db.get_psa(record_int_ids=ids) scalars = db.get_scalars(record_int_ids=ids) - ims = read_ims(share_dir / rupture / f"R{n}" / "intensity_measures.h5") + ims = read_ims(manifest[rel].im_path) row_of = {site: i for i, site in enumerate(ims.stations)} for record_int_id, record in sample.iterrows(): i, component = row_of[record["site_id"]], record["component"] @@ -490,22 +722,35 @@ def spot_check( return checked +def peak_line(row: ManifestRow, peak: dict[str, tuple[float, str]]) -> str: + """One realisation's largest PGA and PGV (geom), with where they are.""" + (pga, pga_site), (pgv, pgv_site) = peak["PGA"], peak["PGV"] + return f"{row.label}: PGA {pga:.3f} g at {pga_site}; PGV {pgv:.1f} cm/s at {pgv_site}" + + +def print_peaks(peaks: Sequence[tuple[ManifestRow, dict[str, tuple[float, str]]]]) -> None: + """The realisations with the largest PGA, for a human to judge before the database is shared.""" + ranked = sorted(peaks, key=lambda item: item[1]["PGA"][0], reverse=True) + print(f"top {min(N_TOP_PEAKS, len(ranked))} realisations by max PGA (geom):", flush=True) + for row, peak in ranked[:N_TOP_PEAKS]: + print(f" {peak_line(row, peak)}", flush=True) + + def build_database( - share_dir: Path, - events_dir: Path, + manifest_path: Path, stations_input: Path, nshm_csv: Path, - realisations: Sequence[tuple[str, int]], out_path: Path, - source: str, + harvested_at: str, script_commit: str, memory_limit: str | None = None, ) -> dict[str, int]: """Build the database at out_path and return its table row counts. Writes `.partial` and renames it only once validation, the row - counts and the spot check have all passed. A failure leaves the .partial - behind for inspection, and the next run refuses to start until it is removed. + counts, the site_event flags and the spot check have all passed. A failure + leaves the .partial behind for inspection, and the next run refuses to + start until it is removed. """ if out_path.exists(): raise FileExistsError(f"{out_path} already exists") @@ -513,14 +758,13 @@ def build_database( if partial.exists(): raise FileExistsError(f"{partial} exists: a previous run died; remove it first") + rows = read_manifest(manifest_path) + by_event = group_by_event(rows) stations = load_stations_input(stations_input) nshm = load_nshm_attributes(nshm_csv) - by_rupture = group_by_rupture(realisations) - missing = sorted(set(by_rupture) - set(nshm.index), key=int) + missing = sorted(set(by_event) - set(nshm.index), key=int) if missing: raise ValueError(f"ruptures missing from {nshm_csv}: {missing}") - (nshmdb_file,) = nshm["nshmdb_file"].unique() - (nshmdb_sha256,) = nshm["nshmdb_sha256"].unique() db = IMDB.create( partial, @@ -528,57 +772,79 @@ def build_database( components=COMPONENTS, db_meta={ **DB_META, - "source": source, - "n_realisations": str(len(realisations)), - "nshmdb": f"{nshmdb_file} sha256 {nshmdb_sha256}", + "n_realisations": str(len(rows)), + "n_pilot_realisations": str(sum(row.pilot for row in rows)), + "im_calc_harvested_at": harvested_at, + "build_manifest_md5": md5_of(manifest_path), "ingest_script_commit": script_commit, }, ) counts = dict.fromkeys(("events", "realisations", "sites", "site_event", "records"), 0) + n_flagged = 0 + seeds_seen: dict[str, str] = {} + peaks: list[tuple[ManifestRow, dict[str, tuple[float, str]]]] = [] started = time.monotonic() try: if memory_limit: db.con.raw_sql(f"SET memory_limit = '{memory_limit}'") sites = SiteRegistry(stations) - for index, (rupture, ns) in enumerate(by_rupture.items(), start=1): + for index, (event_id, group) in enumerate(by_event.items(), start=1): loaded = [] - for n in ns: - rel = rel_id(rupture, n) - realisation = json.loads((events_dir / rupture / f"R{n}" / "realisation.json").read_text()) - h5_path = share_dir / rupture / f"R{n}" / "intensity_measures.h5" - if not h5_path.exists(): - raise FileNotFoundError(f"{rel}: no {h5_path}") - ims = read_ims(h5_path) - check_finite(rel, ims) - loaded.append((n, realisation, ims)) - - first_n, first_realisation, first_ims = loaded[0] - event, geometries = build_event(rupture, first_realisation, nshm.loc[rupture]) - rows = [] - for n, realisation, ims in loaded: - row = build_realisation(rupture, n, geometries, realisation) - check_against_h5(row["rel_id"], row, ims.attrs) - check_matches_first_realisation(row["rel_id"], rel_id(rupture, first_n), ims, first_ims) - rows.append(row) - new_sites = sites.new_sites(rupture, first_ims) + for row in group: + for path in (row.im_path, row.realisation_path): + if not path.exists(): + raise FileNotFoundError(f"{row.label}: no {path}") + realisation = json.loads(row.realisation_path.read_text()) + seeds = json.dumps(realisation["seeds"], sort_keys=True) + if seeds in seeds_seen: + raise ValueError(f"{row.label}: the same seeds as {seeds_seen[seeds]}") + seeds_seen[seeds] = row.rel_id + ims = read_ims(row.im_path) + check_finite(row.label, ims) + loaded.append(Loaded(row, realisation, ims)) + peaks.append((row, ims.peaks)) + print(f"peak {peak_line(row, ims.peaks)}", flush=True) + + mains = [x for x in loaded if not x.row.pilot] + pilots = [x for x in loaded if x.row.pilot] + defining = mains[0] if mains else pilots[0] + event = build_event( + event_id, + defining.realisation, + [x.realisation for x in loaded], + nshm.loc[event_id], + MAIN_RELEASE if mains else PILOT_RELEASE, + ) + realisation_rows = [] + for x in loaded: + realisation_row = build_realisation(x.row, x.realisation) + check_against_h5(x.row.label, realisation_row, x.ims.attrs) + realisation_rows.append(realisation_row) + for x in mains[1:]: + check_matches_first_realisation(x.row.label, mains[0].row.label, x.ims, mains[0].ims) + new_sites = pd.concat([sites.new_sites(x.row.label, x.ims) for x in loaded], ignore_index=True) + site_event, flagged = build_site_event( + event_id, mains[0].ims if mains else None, pilots[0].ims if pilots else None + ) db.add_events(pd.DataFrame([event])) - db.add_realisations(pd.DataFrame(rows)) + db.add_realisations(pd.DataFrame(realisation_rows)) if len(new_sites): db.add_sites(new_sites) - db.add_site_event(build_site_event(rupture, first_ims)) - for n, _, ims in loaded: - records, psa = build_records(rel_id(rupture, n), ims) + db.add_site_event(site_event) + for x in loaded: + records, psa = build_records(x.row.rel_id, x.ims) db.add_records(records, pSA=psa) counts["records"] += len(records) counts["events"] += 1 counts["realisations"] += len(loaded) counts["sites"] += len(new_sites) - counts["site_event"] += len(first_ims.stations) - if index % 25 == 0 or index == len(by_rupture): + counts["site_event"] += len(site_event) + n_flagged += flagged + if index % 25 == 0 or index == len(by_event): print( - f"{index}/{len(by_rupture)} ruptures, {counts['records']} records, " + f"{index}/{len(by_event)} events, {counts['records']} records, " f"{time.monotonic() - started:.0f} s", flush=True, ) @@ -588,11 +854,17 @@ def build_database( raise RuntimeError(f"validate(): {problems}") expected = {**counts, "psa_ims": counts["records"], "scalars_ims": counts["records"]} check_row_counts(db, expected) + check_flags(db, n_flagged) finally: db.close() - checked = spot_check(partial, share_dir) - print(f"validate() clean; row counts match; spot check: {checked} records identical to their h5") + checked = spot_check(partial, {row.rel_id: row for row in rows}) + print( + f"validate() clean; row counts match; {n_flagged} site_event rows carry the " + f"pilot's distances; spot check: {checked} records identical to their h5", + flush=True, + ) + print_peaks(peaks) partial.rename(out_path) return expected @@ -601,13 +873,11 @@ def main() -> None: parser = argparse.ArgumentParser( description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter ) - parser.add_argument("share_dir", type=Path, help="run3 share/: /R/intensity_measures.h5") - parser.add_argument("events_dir", type=Path, help="run3 events/: /R/realisation.json") + parser.add_argument("manifest", type=Path, help="build manifest CSV from build_imdb_manifest.py") parser.add_argument("stations_input", type=Path, help="the campaign's stations_input.ll (lon lat name)") parser.add_argument("nshm_csv", type=Path, help="campaign_nshm_ruptures.csv from export_nshm_rupture_attributes.py") - parser.add_argument("realisations_file", type=Path, help="one /R per line") parser.add_argument("out", type=Path, help="output .duckdb; must not exist") - parser.add_argument("--source", required=True, help="provenance for db_meta.source") + parser.add_argument("--harvested-at", required=True, help="when the manifest was built from cylc, UTC ISO 8601") parser.add_argument("--script-commit", required=True, help="git commit of this script, for db_meta") parser.add_argument( "--memory-limit", @@ -616,16 +886,13 @@ def main() -> None: ) args = parser.parse_args() - realisations = parse_realisation_ids(args.realisations_file.read_text().splitlines()) - print(f"ingesting {len(realisations)} realisations into {args.out}", flush=True) + print(f"ingesting {args.manifest} into {args.out}", flush=True) counts = build_database( - args.share_dir, - args.events_dir, + args.manifest, args.stations_input, args.nshm_csv, - realisations, args.out, - source=args.source, + harvested_at=args.harvested_at, script_commit=args.script_commit, memory_limit=args.memory_limit, ) diff --git a/scripts/tests/test_cs_nshm_2022_ingest.py b/scripts/tests/test_cs_nshm_2022_ingest.py index 42a8a28..1db8714 100644 --- a/scripts/tests/test_cs_nshm_2022_ingest.py +++ b/scripts/tests/test_cs_nshm_2022_ingest.py @@ -1,10 +1,12 @@ -"""Tests for scripts/cs_nshm_2022_ingest.py, on a synthetic campaign in run3's layout. +"""Tests for scripts/cs_nshm_2022_ingest.py, on a synthetic two-campaign dataset. imdb's own CI installs neither h5py nor source-modelling, so this module skips there. Run it inside the imdb_tools image: apptainer exec imdb_tools_.sif python3 -m pytest -q scripts/tests """ +import csv +import hashlib import importlib.util import json import sys @@ -19,6 +21,7 @@ pytest.importorskip("source_modelling") from imdb import IMDB # noqa: E402 +from shapely import unary_union, wkt # noqa: E402 from source_modelling import moment, sources # noqa: E402 SCRIPT = Path(__file__).resolve().parents[1] / "cs_nshm_2022_ingest.py" @@ -38,6 +41,7 @@ 4.2, 4.4, 4.6, 4.8, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 20.0, ) # fmt: skip +FAS_FREQUENCIES = (0.1, 1.0, 10.0) ALL_COMPONENTS = ("000", "090", "ver", "geom", "rotd0", "rotd50", "rotd100") AS_RECORDED = ("000", "090", "ver", "geom") H5_GROUPS = { @@ -60,7 +64,7 @@ "2TfEMSL": (172.4, -43.4), "2TfFKSL": (172.5, -43.5), } -# vs30, z1pt0, z2pt5 per station: global, like the real campaign's. +# vs30, z1pt0, z2pt5 per station: global, like the real campaigns'. SITE_VALUES = { "AMBC": (300.0, 0.3, 1.3), "2O8aUSL": (400.0, 0.2, 1.2), @@ -68,6 +72,7 @@ "2TfEMSL": (600.0, 0.1, 1.0), "2TfFKSL": (700.0, 0.05, 0.9), } +MAIN_SHA256 = "3fe692cb9b22c769b6a9baa6a526b4eed989bd61ec054a344f0e74d82855bc4e" def _plane(lat: float, lon: float, strike: float) -> sources.Plane: @@ -83,6 +88,8 @@ def _realisation( tree: dict[str, str | None], hypocentre: tuple[float, float], duration: float, + seed: int, + domain_shift: float = 0.0, ) -> dict: return { "metadata": {"name": f"Rupture {name}"}, @@ -103,13 +110,13 @@ def _realisation( "rupture_causality_tree": tree, "hypocentre": {"s": hypocentre[0], "d": hypocentre[1]}, }, - "seeds": {"hf_seed": 7}, + "seeds": {"hf_seed": seed, "genslip_seed": 1000 + seed}, "domain": { "domain": [ - {"latitude": -43.0, "longitude": 172.0}, - {"latitude": -43.0, "longitude": 173.0}, - {"latitude": -44.0, "longitude": 173.0}, - {"latitude": -44.0, "longitude": 172.0}, + {"latitude": -43.0 + domain_shift, "longitude": 172.0 + domain_shift}, + {"latitude": -43.0 + domain_shift, "longitude": 173.0 + domain_shift}, + {"latitude": -44.0 + domain_shift, "longitude": 173.0 + domain_shift}, + {"latitude": -44.0 + domain_shift, "longitude": 172.0 + domain_shift}, ], "duration": duration, "depth": 40.0, @@ -167,98 +174,143 @@ def _write_h5(path: Path, stations: list[str], realisation: dict, distances: np. for c in components: group[c] = rng.uniform(1e-3, 2.0, (n, len(VALID_PERIODS))) elif group_name == "FAS": - group["frequency"] = np.array([0.1, 1.0, 10.0]) + group["frequency"] = np.array(FAS_FREQUENCIES) for c in components: - group[c] = rng.uniform(size=(n, 3)) + group[c] = rng.uniform(size=(n, len(FAS_FREQUENCIES))) else: for c in components: group[c] = rng.uniform(0.01, 5.0, n) +def _write_manifest(path: Path, rows: list[dict[str, str]]) -> None: + with path.open("w", newline="") as f: + writer = csv.DictWriter(f, fieldnames=list(rows[0])) + writer.writeheader() + writer.writerows(rows) + + @dataclass class Campaign: - share: Path - events: Path + manifest: Path + files: dict[str, tuple[Path, Path]] + """original_id -> (intensity_measures.h5, realisation.json)""" stations_input: Path nshm_csv: Path - realisations: list[tuple[str, int]] out: Path - def h5(self, rel: str) -> Path: - rupture, n = rel.split("/R") - return self.share / rupture / f"R{n}" / "intensity_measures.h5" + def h5(self, original_id: str) -> Path: + return self.files[original_id][0] - def realisation_json(self, rel: str) -> Path: - rupture, n = rel.split("/R") - return self.events / rupture / f"R{n}" / "realisation.json" + def realisation_json(self, original_id: str) -> Path: + return self.files[original_id][1] + + def realisation(self, original_id: str) -> dict: + return json.loads(self.realisation_json(original_id).read_text()) def build(self) -> dict[str, int]: return ingest.build_database( - self.share, - self.events, + self.manifest, self.stations_input, self.nshm_csv, - self.realisations, self.out, - source="pytest", + harvested_at="2026-09-23T07:35:00Z", script_commit="abc123", ) -RUPTURE_STATIONS = { - "161984": ["AMBC", "2O8aUSL", "2O8bSSL", "2TfEMSL"], - "288271": ["2O8bSSL", "2TfEMSL", "2TfFKSL"], -} -NSHM_MAGNITUDE = {"161984": 7.480455658305042, "288271": 7.13806256610197} +# (original_id, rel_id, pilot, stations) in manifest order +REALISATIONS = [ + ("3", "3_R1", True, ["AMBC", "2O8aUSL"]), + ("161984/R2", "161984_R1", False, ["AMBC", "2O8aUSL", "2O8bSSL", "2TfEMSL"]), + ("161984/R3", "161984_R2", False, ["AMBC", "2O8aUSL", "2O8bSSL", "2TfEMSL"]), + ("288271/R2", "288271_R1", False, ["2O8bSSL", "2TfEMSL", "2TfFKSL"]), + ("288271", "288271_R2", True, ["2TfEMSL", "2TfFKSL", "AMBC", "2O8aUSL"]), +] +STATIONS = {original_id: stations for original_id, _, _, stations in REALISATIONS} +NSHM_MAGNITUDE = {"3": 6.9, "161984": 7.480455658305042, "288271": 7.13806256610197} @pytest.fixture def campaign(tmp_path: Path) -> Campaign: - """Two ruptures: 161984 (two faults, R2 and R3) and 288271 (one fault, R2).""" - share, events = tmp_path / "share", tmp_path / "events" + """Three events. + + - 161984: two main realisations on two faults. + - 288271: one main realisation and a pilot whose fault's bottom edge sits + elsewhere, whose domain is shifted, and whose stations only partly overlap. + - 3: a pilot realisation alone. + """ two_faults = {"Clarence": _plane(-43.5, 172.5, 45.0), "Kekerengu": _plane(-43.45, 172.58, 50.0)} - one_fault = {"Monowai": _plane(-43.6, 172.4, 30.0)} realisations = { + "3": _realisation("3", {"Acton": _plane(-43.2, 172.2, 10.0)}, {"Acton": 6.8}, {"Acton": None}, (0.4, 0.5), 90.0, 1), "161984/R2": _realisation( "161984", two_faults, {"Clarence": 7.11, "Kekerengu": 6.59}, - {"Kekerengu": "Clarence", "Clarence": None}, (0.17, 0.19), 161.7, + {"Kekerengu": "Clarence", "Clarence": None}, (0.17, 0.19), 161.7, 2, ), "161984/R3": _realisation( "161984", two_faults, {"Clarence": 6.54, "Kekerengu": 6.02}, - {"Clarence": "Kekerengu", "Kekerengu": None}, (0.59, 0.92), 156.5, + {"Clarence": "Kekerengu", "Kekerengu": None}, (0.59, 0.92), 156.5, 3, ), "288271/R2": _realisation( - "288271", one_fault, {"Monowai": 7.32}, {"Monowai": None}, (0.5, 0.6), 120.0, + "288271", {"Monowai": _plane(-43.6, 172.4, 30.0)}, {"Monowai": 7.32}, {"Monowai": None}, (0.5, 0.6), 120.0, 4, + ), + "288271": _realisation( + "288271", {"Monowai": _plane(-43.6, 172.4, 33.0)}, {"Monowai": 7.14}, {"Monowai": None}, (0.1, 0.9), 110.0, 5, + domain_shift=0.2, ), } # fmt: skip rng = np.random.default_rng(1) - distances = {rupture: rng.uniform(0, 80, (len(s), 4)) for rupture, s in RUPTURE_STATIONS.items()} - for seed, (rel, realisation) in enumerate(realisations.items()): - rupture = rel.split("/")[0] - path = events / rel / "realisation.json" - path.parent.mkdir(parents=True) - path.write_text(json.dumps(realisation)) - _write_h5(share / rel / "intensity_measures.h5", RUPTURE_STATIONS[rupture], realisation, distances[rupture], seed) + distances = { + "3": rng.uniform(0, 80, (2, 4)), + "161984": rng.uniform(0, 80, (4, 4)), + "288271/R2": rng.uniform(0, 80, (3, 4)), + "288271": rng.uniform(0, 80, (4, 4)), + } + files, rows = {}, [] + for seed, (original_id, rel, pilot, stations) in enumerate(REALISATIONS): + base = tmp_path / ("pilot" if pilot else "main") / original_id + h5_path, json_path = base / "intensity_measures.h5", base / "realisation.json" + json_path.parent.mkdir(parents=True, exist_ok=True) + json_path.write_text(json.dumps(realisations[original_id])) + rupture_distances = distances["161984"] if original_id.startswith("161984") else distances[original_id] + _write_h5(h5_path, stations, realisations[original_id], rupture_distances, seed) + files[original_id] = (h5_path, json_path) + event_id, n = rel.split("_R") + rows.append( + { + "rel_id": rel, + "event_id": event_id, + "realisation": n, + "pilot": "true" if pilot else "false", + "campaign": "run2" if pilot else "run3", + "original_id": original_id, + "im_path": str(h5_path), + "realisation_path": str(json_path), + } + ) + manifest = tmp_path / "manifest.csv" + _write_manifest(manifest, rows) stations_input = tmp_path / "stations_input.ll" stations_input.write_text("".join(f"{lon} {lat} {s}\n" for s, (lon, lat) in STATIONS_INPUT.items())) nshm_csv = tmp_path / "campaign_nshm_ruptures.csv" pd.DataFrame( { - "rupture_id": ["607", "161984", "288271"], - "nshmdb_rupture_id": [40083, 201460, 327747], - "magnitude_boldm": [7.47, NSHM_MAGNITUDE["161984"], NSHM_MAGNITUDE["288271"]], - "area_km2": [1872.5, 1907.5, 867.1], - "length_km": [59.9, 80.2, 34.6], - "annual_rate": [6.1e-06, 1.86e-04, 7.97e-05], + "rupture_id": ["3", "607", "161984", "288271"], + "nshmdb_rupture_id": [3001, 40083, 201460, 327747], + "magnitude_boldm": [NSHM_MAGNITUDE["3"], 7.47, NSHM_MAGNITUDE["161984"], NSHM_MAGNITUDE["288271"]], + "area_km2": [300.0, 1872.5, 1907.5, 867.1], + "length_km": [20.0, 59.9, 80.2, 34.6], + "annual_rate": [1e-05, 6.1e-06, 1.86e-04, 7.97e-05], "nshmdb_file": "nshmdb_v2026.08.3.db", - "nshmdb_sha256": "3fe692cb", + "nshmdb_sha256": MAIN_SHA256, } ).to_csv(nshm_csv, index=False) - return Campaign( - share, events, stations_input, nshm_csv, - [("161984", 2), ("161984", 3), ("288271", 2)], tmp_path / "out.duckdb", - ) # fmt: skip + return Campaign(manifest, files, stations_input, nshm_csv, tmp_path / "out.duckdb") + + +def _manifest_rows(campaign: Campaign) -> list[dict[str, str]]: + with campaign.manifest.open(newline="") as f: + return list(csv.DictReader(f)) # ---- helpers ------------------------------------------------------------------- @@ -276,81 +328,145 @@ def test_select_period_indices_rejects_a_missing_period(): ingest.select_period_indices(without) -def test_parse_realisation_ids(): - lines = ["161984/R3", "", "607/R2", "161984/R2", "607/R2"] - assert ingest.parse_realisation_ids(lines) == [("607", 2), ("161984", 2), ("161984", 3)] - with pytest.raises(ValueError, match="161984_R2"): - ingest.parse_realisation_ids(["161984_R2"]) - - def test_is_real_station(): assert ingest.is_real_station("AMBC") assert ingest.is_real_station("BFZ") assert not ingest.is_real_station("2O8aUSL") +def test_read_manifest_returns_ingest_order(campaign: Campaign): + rows = _manifest_rows(campaign) + _write_manifest(campaign.manifest, rows[::-1]) + parsed = ingest.read_manifest(campaign.manifest) + assert [row.rel_id for row in parsed] == ["3_R1", "161984_R1", "161984_R2", "288271_R1", "288271_R2"] + assert [row.pilot for row in parsed] == [True, False, False, False, True] + assert parsed[1].original_id == "161984/R2" + assert parsed[1].im_path == campaign.h5("161984/R2") + + +def _set(rows: list[dict[str, str]], rel: str, **changes: str) -> list[dict[str, str]]: + return [{**row, **changes} if row["rel_id"] == rel else row for row in rows] + + +@pytest.mark.parametrize( + "edit, match", + [ + (lambda rows: _set(rows, "161984_R2", rel_id="161984_R3", realisation="3"), r"realisations \[1, 3\]"), + (lambda rows: _set(_set(rows, "288271_R1", pilot="true"), "288271_R2", pilot="false"), "the pilot is R1"), + (lambda rows: _set(rows, "288271_R1", pilot="true"), "2 pilot realisations"), + (lambda rows: [*rows, rows[0]], "3_R1: listed twice"), + (lambda rows: _set(rows, "3_R1", rel_id="3_R9"), "3_R9: expected 3_R1"), + (lambda rows: _set(rows, "3_R1", event_id="03", rel_id="03_R1"), "not an unpadded integer"), + (lambda rows: _set(rows, "3_R1", pilot="yes"), "pilot must be true or false"), + (lambda rows: [{k: v for k, v in row.items() if k != "original_id"} for row in rows], "missing columns"), + ], + ids=["gap", "pilot-not-last", "two-pilots", "duplicate", "rel-id-mismatch", "padded-event", "bad-pilot", "no-column"], +) # fmt: skip +def test_bad_manifests_are_refused(campaign: Campaign, edit, match): + _write_manifest(campaign.manifest, edit(_manifest_rows(campaign))) + with pytest.raises(ValueError, match=match): + ingest.read_manifest(campaign.manifest) + + # ---- the whole build ------------------------------------------------------------- def test_row_counts_and_meta(campaign: Campaign): counts = campaign.build() assert counts == { - "events": 2, - "realisations": 3, + "events": 3, + "realisations": 5, "sites": 5, - "site_event": 7, - "records": 22, - "psa_ims": 22, - "scalars_ims": 22, + "site_event": 11, # 2 + 4 + (3 main + 2 only the pilot has) + "records": 34, + "psa_ims": 34, + "scalars_ims": 34, } assert campaign.out.exists() assert not campaign.out.with_name("out.duckdb.partial").exists() with IMDB(campaign.out) as db: meta = db.db_meta assert db.validate() == [] - assert meta["dataset_id"] == "cs_nshm_2022_run3" + assert meta["dataset_id"] == "cs_nshm_2022" assert meta["components"] == "geom,rotd50" assert meta["n_periods"] == "25" assert meta["n_frequencies"] == "0" - assert meta["magnitude_convention"].startswith("BoldM") - assert meta["nshmdb"] == "nshmdb_v2026.08.3.db sha256 3fe692cb" + assert meta["n_realisations"] == "5" + assert meta["n_pilot_realisations"] == "2" + assert meta["im_calc_harvested_at"] == "2026-09-23T07:35:00Z" + assert meta["build_manifest_md5"] == hashlib.md5(campaign.manifest.read_bytes()).hexdigest() assert meta["ingest_script_commit"] == "abc123" - assert meta["source"] == "pytest" - assert meta["n_realisations"] == "3" + assert meta["magnitude_convention"].startswith("BoldM") + for key in ("description", "realisation_numbering", "pilot_realisations", "distances", "nshm_fault_database"): + assert meta[key] + assert "source" not in meta and "nshmdb" not in meta + for key, value in meta.items(): # nothing researcher-facing names our internals + for word in ("run2", "run3", "/gpfs", "BSC", "envelope"): + assert word not in value, f"db_meta {key} mentions {word}" def test_events(campaign: Campaign): campaign.build() with IMDB(campaign.out) as db: events = db.get_events() - assert events.loc["161984", "magnitude"] == pytest.approx(NSHM_MAGNITUDE["161984"], abs=1e-6) - assert events.loc["288271", "magnitude"] == pytest.approx(NSHM_MAGNITUDE["288271"], abs=1e-6) + for event_id, magnitude in NSHM_MAGNITUDE.items(): + assert events.loc[event_id, "magnitude"] == pytest.approx(magnitude, abs=1e-6) assert (events["tect_type"] == "ACTIVE_SHALLOW").all() assert pd.isna(events.loc["161984", "dip"]) # multi-fault: no single dip assert events.loc["288271", "dip"] == pytest.approx(60.0, abs=1e-3) - assert events.loc["288271", "length"] == pytest.approx(10.0, abs=1e-3) - metadata = json.loads(events.loc["161984", "metadata"]) - assert metadata["nshmdb_rupture_id"] == 201460 - assert metadata["annual_rate"] == 1.86e-04 - assert set(metadata["segments"]) == {"Clarence", "Kekerengu"} - assert events.loc["161984", "source_wkt"].startswith("MULTIPOLYGON") + metadata = {e: json.loads(events.loc[e, "metadata"]) for e in NSHM_MAGNITUDE} + assert metadata["161984"]["nshmdb_rupture_id"] == 201460 + assert metadata["161984"]["annual_rate"] == 1.86e-04 + assert set(metadata["161984"]["segments"]) == {"Clarence", "Kekerengu"} + assert {e: m["fault_geometry_release"] for e, m in metadata.items()} == { + "3": "pre-v2026.08", + "161984": "v2026.08.3", + "288271": "v2026.08.3", + } + + +def test_event_geometry_comes_from_the_main_realisation(campaign: Campaign): + campaign.build() + with IMDB(campaign.out) as db: + events = db.get_events() + + def geometry_of(original_id: str) -> str: + return ingest.source_wkt(ingest.fault_geometries(campaign.realisation(original_id))) + + assert events.loc["288271", "source_wkt"] == geometry_of("288271/R2") + assert events.loc["288271", "source_wkt"] != geometry_of("288271") + assert events.loc["3", "source_wkt"] == geometry_of("3") + + +def test_domain_is_the_union_of_the_realisations_domains(campaign: Campaign): + campaign.build() + with IMDB(campaign.out) as db: + events = db.get_events() + main, pilot = (ingest.domain_polygon(campaign.realisation(o)) for o in ("288271/R2", "288271")) + stored = wkt.loads(events.loc["288271", "domain_wkt"]) + assert stored.equals(unary_union([main, pilot])) + assert stored.area > main.area + assert wkt.loads(events.loc["161984", "domain_wkt"]).equals(ingest.domain_polygon(campaign.realisation("161984/R2"))) -def test_realisations_are_boldm_and_match_their_h5(campaign: Campaign): +def test_realisations_match_their_own_realisation_json(campaign: Campaign): campaign.build() with IMDB(campaign.out) as db: realisations = db.get_realisations() - for rel in ("161984/R2", "161984/R3", "288271/R2"): - realisation = json.loads(campaign.realisation_json(rel).read_text()) - row = realisations.loc[rel.replace("/", "_")] + for original_id, rel, pilot, _ in REALISATIONS: + realisation = campaign.realisation(original_id) + row = realisations.loc[rel] assert row["magnitude"] == pytest.approx(_boldm_total(realisation), abs=1e-5) + # the pilot's hypocentre sits on the pilot's own fault, not the event's assert (row["hypo_lat"], row["hypo_lon"], row["hypo_depth"]) == pytest.approx( _hypocentre(realisation), abs=1e-4 ) metadata = json.loads(row["metadata"]) + assert metadata["pilot"] is pilot + assert metadata["seeds"] == realisation["seeds"] assert metadata["duration_s"] == realisation["domain"]["duration"] - # R3 starts on the other fault, so its causality tree is its own - assert json.loads(realisations.loc["161984_R3", "metadata"])["rupture_causality_tree"]["Kekerengu"] is None + # 161984_R2 starts on the other fault, so its causality tree is its own + assert json.loads(realisations.loc["161984_R2", "metadata"])["rupture_causality_tree"]["Kekerengu"] is None def test_sites_use_canonical_coordinates(campaign: Campaign): @@ -365,33 +481,63 @@ def test_sites_use_canonical_coordinates(campaign: Campaign): assert sites["is_real"].to_dict() == {s: s == "AMBC" for s in STATIONS_INPUT} -def test_site_event_comes_from_the_h5(campaign: Campaign): +def _h5_distances(path: Path) -> pd.DataFrame: + with h5py.File(path) as h5: + stations = [s.decode() for s in h5["pSA/station"][:]] + return pd.DataFrame({key: h5[f"pSA/{key}"][:] for key in ingest.DISTANCE_FIELDS}, index=stations) + + +def test_site_event_prefers_the_main_campaign_and_flags_the_pilots(campaign: Campaign): campaign.build() - with h5py.File(campaign.h5("288271/R2")) as h5: - rrup = h5["pSA/rrup"][:] with IMDB(campaign.out) as db: - site_event = db.get_site_event(event_ids=["288271"]).set_index("site_id") - assert site_event.loc[RUPTURE_STATIONS["288271"], "rrup"].to_numpy() == pytest.approx(rrup, abs=1e-4) + site_event = db.get_site_event() + by_event = {e: df.set_index("site_id") for e, df in site_event.groupby("event_id")} + main, pilot = _h5_distances(campaign.h5("288271/R2")), _h5_distances(campaign.h5("288271")) + shared = by_event["288271"] + assert sorted(shared.index) == sorted(STATIONS_INPUT) + for site in STATIONS["288271/R2"]: # the main realisation's stations, 2TfEMSL included + assert shared.loc[site, list(ingest.DISTANCE_FIELDS)].to_numpy() == pytest.approx(main.loc[site].to_numpy()) + assert pd.isna(shared.loc[site, "metadata"]) + for site in ("AMBC", "2O8aUSL"): # only the pilot has these + assert shared.loc[site, list(ingest.DISTANCE_FIELDS)].to_numpy() == pytest.approx(pilot.loc[site].to_numpy()) + assert json.loads(shared.loc[site, "metadata"]) == {"fault_geometry": "pilot"} + # the pilot-only event keeps the pilot's distances, unflagged + alone = _h5_distances(campaign.h5("3")) + for site in STATIONS["3"]: + assert by_event["3"].loc[site, list(ingest.DISTANCE_FIELDS)].to_numpy() == pytest.approx(alone.loc[site].to_numpy()) + assert by_event["3"]["metadata"].isna().all() + assert by_event["161984"]["metadata"].isna().all() def test_ims_are_the_h5_values(campaign: Campaign): campaign.build() - with h5py.File(campaign.h5("161984/R3")) as h5: - columns = ingest.select_period_indices(h5["pSA/period"][:]) - station = RUPTURE_STATIONS["161984"].index("2TfEMSL") - psa_rotd50 = h5["pSA/rotd50"][station, columns].astype(np.float32) - cav_geom = np.float32(h5["CAV/geom"][station]) - pga_rotd50 = np.float32(h5["PGA/rotd50"][station]) - with IMDB(campaign.out) as db: - rotd50 = db.get_records(rel_ids=["161984_R3"], site_ids=["2TfEMSL"], component="rotd50").index - geom = db.get_records(rel_ids=["161984_R3"], site_ids=["2TfEMSL"], component="geom").index - psa = db.get_psa(record_int_ids=list(rotd50)) - scalars = db.get_scalars(record_int_ids=[*rotd50, *geom]) - np.testing.assert_array_equal(psa.iloc[0].to_numpy(dtype=np.float32), psa_rotd50) - assert list(psa.columns) == [str(p) for p in ingest.PSA_PERIODS] - assert np.float32(scalars.loc[rotd50[0], "PGA"]) == pga_rotd50 - assert pd.isna(scalars.loc[rotd50[0], "CAV"]) # undefined on rotd - assert np.float32(scalars.loc[geom[0], "CAV"]) == cav_geom + for original_id, rel in (("161984/R3", "161984_R2"), ("288271", "288271_R2")): + with h5py.File(campaign.h5(original_id)) as h5: + columns = ingest.select_period_indices(h5["pSA/period"][:]) + station = STATIONS[original_id].index("2TfEMSL") + psa_rotd50 = h5["pSA/rotd50"][station, columns].astype(np.float32) + cav_geom = np.float32(h5["CAV/geom"][station]) + pga_rotd50 = np.float32(h5["PGA/rotd50"][station]) + with IMDB(campaign.out) as db: + rotd50 = db.get_records(rel_ids=[rel], site_ids=["2TfEMSL"], component="rotd50").index + geom = db.get_records(rel_ids=[rel], site_ids=["2TfEMSL"], component="geom").index + psa = db.get_psa(record_int_ids=list(rotd50)) + scalars = db.get_scalars(record_int_ids=[*rotd50, *geom]) + np.testing.assert_array_equal(psa.iloc[0].to_numpy(dtype=np.float32), psa_rotd50) + assert list(psa.columns) == [str(p) for p in ingest.PSA_PERIODS] + assert np.float32(scalars.loc[rotd50[0], "PGA"]) == pga_rotd50 + assert pd.isna(scalars.loc[rotd50[0], "CAV"]) # undefined on rotd + assert np.float32(scalars.loc[geom[0], "CAV"]) == cav_geom + + +def test_peaks_are_logged(campaign: Campaign, capsys: pytest.CaptureFixture): + campaign.build() + out = capsys.readouterr().out + with h5py.File(campaign.h5("288271")) as h5: + pga = h5["PGA/geom"][:] + site = h5["PGA/station"][int(np.argmax(pga))].decode() + assert f"peak 288271_R2 (288271): PGA {pga.max():.3f} g at {site}" in out + assert "top 5 realisations by max PGA (geom):" in out # ---- refusals and consistency checks ------------------------------------------- @@ -428,19 +574,29 @@ def _set_attr(path: Path, key: str, delta: float) -> None: h5.attrs[key] = h5.attrs[key] + delta +def _edit_json(path: Path, edit) -> None: + realisation = json.loads(path.read_text()) + edit(realisation) + path.write_text(json.dumps(realisation)) + + @pytest.mark.parametrize( "corrupt, error, match", [ (lambda c: _edit_h5(c.h5("288271/R2"), "geom", (0, KEPT_COLUMN), np.nan, ("pSA",)), ValueError, "pSA/geom has 1 non-finite"), - (lambda c: _edit_h5(c.h5("288271/R2"), "vs30", 1, 999.0), ValueError, "earlier rupture"), - (lambda c: _edit_h5(c.h5("161984/R3"), "rrup", 2, 1.5), ValueError, "rrup differs from 161984_R2"), + (lambda c: _edit_h5(c.h5("288271/R2"), "vs30", 1, 999.0), ValueError, "earlier realisation"), + (lambda c: _edit_h5(c.h5("288271"), "vs30", 0, 999.0), ValueError, r"288271_R2 \(288271\): 2TfEMSL"), + (lambda c: _edit_h5(c.h5("161984/R3"), "rrup", 2, 1.5), ValueError, r"rrup differs from 161984_R1 \(161984/R2\)"), (lambda c: _set_attr(c.h5("161984/R3"), "magnitude", 0.01), ValueError, "magnitude from realisation.json"), - (lambda c: _set_attr(c.h5("288271/R2"), "hypo_lat", 0.001), ValueError, "hypo_lat from realisation.json"), + (lambda c: _set_attr(c.h5("288271"), "hypo_lat", 0.001), ValueError, "hypo_lat from realisation.json"), (lambda c: c.stations_input.write_text(c.stations_input.read_text().replace("2TfFKSL", "XXXXXXX")), ValueError, "missing from stations_input.ll"), - (lambda c: c.h5("288271/R2").unlink(), FileNotFoundError, "288271_R2"), + (lambda c: c.h5("288271/R2").unlink(), FileNotFoundError, r"288271_R1 \(288271/R2\)"), + (lambda c: _edit_json(c.realisation_json("3"), lambda r: r.update(seeds=c.realisation("161984/R3")["seeds"])), + ValueError, "the same seeds as 3_R1"), ], - ids=["nan", "site-mismatch", "distance-mismatch", "magnitude", "hypocentre", "unknown-station", "no-h5"], + ids=["nan", "site-mismatch", "pilot-site-mismatch", "distance-mismatch", "magnitude", "hypocentre", + "unknown-station", "no-h5", "same-seeds"], ) # fmt: skip def test_bad_input_stops_the_build(campaign: Campaign, corrupt, error, match): corrupt(campaign) @@ -449,21 +605,33 @@ def test_bad_input_stops_the_build(campaign: Campaign, corrupt, error, match): assert not campaign.out.exists() +def test_the_pilots_distances_are_not_compared_with_the_main_campaigns(campaign: Campaign): + main, pilot = _h5_distances(campaign.h5("288271/R2")), _h5_distances(campaign.h5("288271")) + assert not np.allclose(main.loc["2TfEMSL"], pilot.loc["2TfEMSL"]) + campaign.build() + + def test_unknown_tect_type_stops_the_build(campaign: Campaign): - path = campaign.realisation_json("288271/R2") - realisation = json.loads(path.read_text()) - realisation["empirical"]["tect_type"] = "subduction_mystery" - path.write_text(json.dumps(realisation)) + _edit_json(campaign.realisation_json("288271/R2"), lambda r: r["empirical"].update(tect_type="subduction_mystery")) with pytest.raises(ValueError, match="tect_type"): campaign.build() def test_rupture_missing_from_the_nshm_csv(campaign: Campaign): - pd.read_csv(campaign.nshm_csv).iloc[:2].to_csv(campaign.nshm_csv, index=False) + nshm = pd.read_csv(campaign.nshm_csv) + nshm[nshm["rupture_id"] != 288271].to_csv(campaign.nshm_csv, index=False) with pytest.raises(ValueError, match=r"missing from .*\['288271'\]"): campaign.build() +def test_nshm_attributes_from_another_release_are_refused(campaign: Campaign): + nshm = pd.read_csv(campaign.nshm_csv) + nshm.loc[0, "nshmdb_sha256"] = "00e25648" + nshm.to_csv(campaign.nshm_csv, index=False) + with pytest.raises(ValueError, match="must all come from v2026.08.3"): + campaign.build() + + def test_nan_in_a_dropped_period_is_ignored(campaign: Campaign): _edit_h5(campaign.h5("288271/R2"), "geom", (0, DROPPED_COLUMN), np.nan, ("pSA",)) - assert campaign.build()["records"] == 22 + assert campaign.build()["records"] == 34 From 7e3443bb1ddc7dcfe5684b8cf7a88724145554c2 Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Wed, 23 Sep 2026 20:23:27 +1200 Subject: [PATCH 3/8] Add --components, --psa-periods all and --fas to the cs_nshm_2022 ingest The default build stays lean (geom and rotd50, 25 periods, no FAS), but a full IMDB can now be built from the same script. eas carries FAS only, so its records go in their own batch with no pSA or scalar rows, and rotd records get no FAS row. Every realisation's frequency grid must match the first one's, as must the period grid under --psa-periods all. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UgVr1htJMAqx8tABJLTnD5 --- scripts/cs_nshm_2022_ingest.py | 230 ++++++++++++++++++---- scripts/tests/test_cs_nshm_2022_ingest.py | 89 ++++++++- 2 files changed, 277 insertions(+), 42 deletions(-) diff --git a/scripts/cs_nshm_2022_ingest.py b/scripts/cs_nshm_2022_ingest.py index 7029458..2eb68c5 100644 --- a/scripts/cs_nshm_2022_ingest.py +++ b/scripts/cs_nshm_2022_ingest.py @@ -43,8 +43,9 @@ realisation, or from the pilot where the main campaign has none. At a site only the pilot covers, site_event holds the pilot's distances, flagged in its metadata. -Trimmed for size: only the geom and rotd50 components, pSA on 25 of im-calc's -111 periods, and no FAS or empirical-GMM records. +Trimmed for size by default: only the geom and rotd50 components, pSA on 25 of +im-calc's 111 periods, and no FAS. --components, --psa-periods all and --fas +put the rest in. Empirical-GMM records are never included. Every magnitude is BoldM (Hanks & Kanamori 1979 eq. 7), the group convention. nz_sim_validation_ingest.py converts to Mw instead. @@ -74,7 +75,10 @@ 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 20.0, ) # fmt: skip -COMPONENTS = ("geom", "rotd50") +DEFAULT_COMPONENTS = ("geom", "rotd50") +# Components im-calc writes for pSA, PGA, PGV and PGD, and for FAS. +PSA_COMPONENTS = ("000", "090", "ver", "geom", "rotd0", "rotd50", "rotd100") +FAS_COMPONENTS = ("000", "090", "ver", "geom", "eas") N_IM_CALC_PERIODS = 111 SITE_FIELDS = ("vs30", "z1pt0", "z2pt5") DISTANCE_FIELDS = ("rrup", "rjb", "rx", "ry") @@ -136,10 +140,6 @@ "catalogue magnitude; realisations.magnitude is the realisation's own " "moment-summed total." ), - "psa_period_selection": ( - "25 of im-calc's 111 periods: 0.1-1.0 s by 0.1 s, 2-15 s by 1 s, and 20 s. " - "im-calc computed no 16-19 s." - ), "station_coordinates": ( "sites.lat/lon are the canonical stations_input.ll coordinates. The " "simulation snapped each station to its rupture's grid (~0.1 km away, so " @@ -170,6 +170,78 @@ def select_period_indices( return np.array(indices) +@dataclass(frozen=True) +class Content: + """Which components, pSA periods and spectra go into the database.""" + + components: tuple[str, ...] = DEFAULT_COMPONENTS + all_periods: bool = False + """All 111 of im-calc's pSA periods, instead of the 25 in PSA_PERIODS.""" + fas: bool = False + """FAS, on im-calc's full frequency grid, for the components that have it.""" + + def __post_init__(self) -> None: + if not self.components: + raise ValueError("no components chosen") + if len(set(self.components)) != len(self.components): + raise ValueError(f"components listed twice: {self.components}") + unknown = set(self.components) - set(PSA_COMPONENTS) - set(FAS_COMPONENTS) + if unknown: + raise ValueError(f"unknown components {sorted(unknown)}; choose from {schema.COMPONENTS}") + if "eas" in self.components and not self.fas: + raise ValueError("eas has only FAS, so it needs --fas") + + @property + def psa_components(self) -> tuple[str, ...]: + """The chosen components that have pSA and the scalar IMs: all but eas.""" + return tuple(c for c in self.components if c in PSA_COMPONENTS) + + @property + def fas_components(self) -> tuple[str, ...]: + """The chosen components that get FAS: none without --fas.""" + return tuple(c for c in self.components if self.fas and c in FAS_COMPONENTS) + + +@dataclass(frozen=True) +class Grids: + """The pSA periods and FAS frequencies the database is created with.""" + + periods: tuple[float, ...] + frequencies: tuple[float, ...] + + +DEFAULT_CONTENT = Content() +DEFAULT_GRIDS = Grids(periods=PSA_PERIODS, frequencies=()) + + +def read_grids(h5_path: Path, content: Content) -> Grids: + """The database's grids, taken from one realisation's h5 where the content needs its full grid.""" + with h5py.File(h5_path, "r") as h5: + periods = tuple(float(p) for p in h5["pSA"]["period"][:]) if content.all_periods else PSA_PERIODS + frequencies = tuple(float(f) for f in h5["FAS"]["frequency"][:]) if content.fas else () + return Grids(periods=periods, frequencies=frequencies) + + +def content_meta(content: Content, grids: Grids) -> dict[str, str]: + """db_meta entries describing what the options put in.""" + if content.all_periods: + periods = f"all {len(grids.periods)} of im-calc's periods, {min(grids.periods):g}-{max(grids.periods):g} s." + else: + periods = ( + "25 of im-calc's 111 periods: 0.1-1.0 s by 0.1 s, 2-15 s by 1 s, and 20 s. " + "im-calc computed no 16-19 s." + ) + if content.fas_components: + fas = ( + f"all {len(grids.frequencies)} of im-calc's FAS frequencies, " + f"{min(grids.frequencies):g}-{max(grids.frequencies):g} Hz, for " + f"{', '.join(content.fas_components)}. rotd components have no FAS." + ) + else: + fas = "not included" + return {"psa_period_selection": periods, "fas": fas} + + def rel_id(event_id: str, n: int) -> str: """The realisation's stable id.""" return f"{event_id}_R{n}" @@ -300,7 +372,11 @@ class RealisationIMs: distances: dict[str, np.ndarray] """DISTANCE_FIELDS per station, km.""" psa: dict[str, np.ndarray] - """Component -> (n_stations, len(PSA_PERIODS)) float32.""" + """Component -> (n_stations, n_periods) float32, for Content.psa_components.""" + fas: dict[str, np.ndarray] + """Component -> (n_stations, n_frequencies) float32, for Content.fas_components.""" + n_frequencies: int + """Length of the database's FAS grid; 0 without FAS.""" scalars: dict[str, dict[str, np.ndarray]] """Scalar IM -> component -> (n_stations,) float32; rotd-undefined IMs omitted.""" attrs: dict[str, float] @@ -319,7 +395,9 @@ def _peak(values: np.ndarray, stations: np.ndarray) -> tuple[float, str]: return float(values[i]), station.decode() if isinstance(station, bytes) else str(station) -def read_ims(h5_path: Path) -> RealisationIMs: +def read_ims( + h5_path: Path, content: Content = DEFAULT_CONTENT, grids: Grids = DEFAULT_GRIDS +) -> RealisationIMs: """Read the parts of one intensity_measures.h5 this database keeps.""" with h5py.File(h5_path, "r") as h5: psa_group = h5["pSA"] @@ -328,23 +406,33 @@ def read_ims(h5_path: Path) -> RealisationIMs: raise ValueError( f"{h5_path}: {len(periods)} pSA periods, expected {N_IM_CALC_PERIODS}" ) - columns = select_period_indices(periods) + columns = select_period_indices(periods, grids.periods) raw_stations = psa_group["station"][:] scalars: dict[str, dict[str, np.ndarray]] = {} for im in schema.SCALAR_IMS: if not np.array_equal(h5[im]["station"][:], raw_stations): raise ValueError(f"{h5_path}: {im} lists stations in a different order from pSA") scalars[im] = {} - for component in COMPONENTS: + for component in content.psa_components: if component in h5[im]: scalars[im][component] = h5[im][component][:].astype(np.float32) elif not (component.startswith("rotd") and im in schema.ROTD_UNDEFINED): raise ValueError(f"{h5_path}: {im} has no {component} component") + fas: dict[str, np.ndarray] = {} + if content.fas: + fas_group = h5["FAS"] + if not np.array_equal(fas_group["frequency"][:], np.array(grids.frequencies)): + raise ValueError(f"{h5_path}: FAS frequency grid differs from the database's") + if not np.array_equal(fas_group["station"][:], raw_stations): + raise ValueError(f"{h5_path}: FAS lists stations in a different order from pSA") + fas = {c: fas_group[c][:].astype(np.float32) for c in content.fas_components} return RealisationIMs( stations=_decode(raw_stations), site={key: psa_group[key][:].astype(np.float64) for key in SITE_FIELDS}, distances={key: psa_group[key][:].astype(np.float64) for key in DISTANCE_FIELDS}, - psa={c: psa_group[c][:][:, columns].astype(np.float32) for c in COMPONENTS}, + psa={c: psa_group[c][:][:, columns].astype(np.float32) for c in content.psa_components}, + fas=fas, + n_frequencies=len(grids.frequencies), scalars=scalars, attrs={key: float(np.atleast_1d(h5.attrs[key])[0]) for key in TOLERANCES}, peaks={im: _peak(h5[im]["geom"][:], raw_stations) for im in ("PGA", "PGV")}, @@ -355,6 +443,7 @@ def check_finite(rel: str, ims: RealisationIMs) -> None: """Raise if any kept IM, site or distance value is NaN or infinite.""" arrays = { **{f"pSA/{c}": values for c, values in ims.psa.items()}, + **{f"FAS/{c}": values for c, values in ims.fas.items()}, **{f"{im}/{c}": values for im, by_c in ims.scalars.items() for c, values in by_c.items()}, **ims.site, **ims.distances, @@ -383,17 +472,35 @@ def check_matches_first_realisation( ) -def build_records(rel: str, ims: RealisationIMs) -> tuple[pd.DataFrame, np.ndarray]: - """IMDB records (one per station and component) and their pSA rows.""" +Batch = tuple[pd.DataFrame, np.ndarray | None, np.ndarray | None] + + +def build_records(rel: str, ims: RealisationIMs, content: Content = DEFAULT_CONTENT) -> list[Batch]: + """IMDB records (one per station and component), as add_records batches of (records, pSA, FAS). + + `eas` has only FAS, so its records form a batch of their own with no pSA + and no scalar columns, and add_records writes no psa_ims or scalars_ims + rows for them. A rotd record's FAS row is all NaN, so it gets no fas_ims row. + """ n = len(ims.stations) frames = [] - for component in COMPONENTS: + for component in content.psa_components: frame = {"rel_id": rel, "site_id": ims.stations, "component": component, "kind": "simulated"} for im, by_component in ims.scalars.items(): frame[im] = by_component.get(component, np.full(n, np.nan, dtype=np.float32)) frames.append(pd.DataFrame(frame)) - records = pd.concat(frames, ignore_index=True) - return records, np.concatenate([ims.psa[c] for c in COMPONENTS]) + batches: list[Batch] = [] + if frames: + psa = np.concatenate([ims.psa[c] for c in content.psa_components]) + fas = None + if content.fas: + no_fas = np.full((n, ims.n_frequencies), np.nan, dtype=np.float32) + fas = np.concatenate([ims.fas.get(c, no_fas) for c in content.psa_components]) + batches.append((pd.concat(frames, ignore_index=True), psa, fas)) + if "eas" in content.fas_components: + eas = pd.DataFrame({"rel_id": rel, "site_id": ims.stations, "component": "eas", "kind": "simulated"}) + batches.append((eas, None, ims.fas["eas"])) + return batches def _site_event_rows(event_id: str, stations: np.ndarray, distances: dict[str, np.ndarray], metadata: str | None) -> pd.DataFrame: @@ -677,6 +784,8 @@ def check_flags(db: IMDB, n_flagged: int) -> None: def spot_check( db_path: Path, manifest: dict[str, ManifestRow], + content: Content = DEFAULT_CONTENT, + grids: Grids = DEFAULT_GRIDS, n_realisations: int = 25, per_realisation: int = 40, min_pilot: int = 5, @@ -703,21 +812,33 @@ def spot_check( ids = sample.index.tolist() psa = db.get_psa(record_int_ids=ids) scalars = db.get_scalars(record_int_ids=ids) - ims = read_ims(manifest[rel].im_path) + fas = db.get_fas(record_int_ids=ids) if content.fas else None + ims = read_ims(manifest[rel].im_path, content, grids) row_of = {site: i for i, site in enumerate(ims.stations)} for record_int_id, record in sample.iterrows(): i, component = row_of[record["site_id"]], record["component"] - stored = psa.loc[record_int_id].to_numpy(dtype=np.float32) - if not np.array_equal(stored, ims.psa[component][i]): - raise AssertionError(f"{rel} {record['site_id']} {component}: pSA differs from the h5") - for im in schema.SCALAR_IMS: - value = scalars.loc[record_int_id, im] - expected = ims.scalars[im].get(component) - if expected is None: - if not pd.isna(value): - raise AssertionError(f"{rel} {record['site_id']} {component}: {im} should be NULL") - elif np.float32(value) != expected[i]: - raise AssertionError(f"{rel} {record['site_id']} {component}: {im} differs from the h5") + where = f"{rel} {record['site_id']} {component}" + if component in content.psa_components: + stored = psa.loc[record_int_id].to_numpy(dtype=np.float32) + if not np.array_equal(stored, ims.psa[component][i]): + raise AssertionError(f"{where}: pSA differs from the h5") + for im in schema.SCALAR_IMS: + value = scalars.loc[record_int_id, im] + expected = ims.scalars[im].get(component) + if expected is None: + if not pd.isna(value): + raise AssertionError(f"{where}: {im} should be NULL") + elif np.float32(value) != expected[i]: + raise AssertionError(f"{where}: {im} differs from the h5") + elif record_int_id in psa.index or record_int_id in scalars.index: + raise AssertionError(f"{where}: should have no pSA or scalar IMs") + if fas is not None: + if component in content.fas_components: + stored = fas.loc[record_int_id].to_numpy(dtype=np.float32) + if not np.array_equal(stored, ims.fas[component][i]): + raise AssertionError(f"{where}: FAS differs from the h5") + elif record_int_id in fas.index: + raise AssertionError(f"{where}: should have no FAS") checked += 1 return checked @@ -743,6 +864,7 @@ def build_database( out_path: Path, harvested_at: str, script_commit: str, + content: Content = DEFAULT_CONTENT, memory_limit: str | None = None, ) -> dict[str, int]: """Build the database at out_path and return its table row counts. @@ -766,12 +888,15 @@ def build_database( if missing: raise ValueError(f"ruptures missing from {nshm_csv}: {missing}") + grids = read_grids(rows[0].im_path, content) db = IMDB.create( partial, - periods=list(PSA_PERIODS), - components=COMPONENTS, + periods=list(grids.periods), + frequencies=list(grids.frequencies), + components=content.components, db_meta={ **DB_META, + **content_meta(content, grids), "n_realisations": str(len(rows)), "n_pilot_realisations": str(sum(row.pilot for row in rows)), "im_calc_harvested_at": harvested_at, @@ -779,7 +904,9 @@ def build_database( "ingest_script_commit": script_commit, }, ) - counts = dict.fromkeys(("events", "realisations", "sites", "site_event", "records"), 0) + counts = dict.fromkeys( + ("events", "realisations", "sites", "site_event", "records", "psa_ims", "scalars_ims", "fas_ims"), 0 + ) n_flagged = 0 seeds_seen: dict[str, str] = {} peaks: list[tuple[ManifestRow, dict[str, tuple[float, str]]]] = [] @@ -799,7 +926,7 @@ def build_database( if seeds in seeds_seen: raise ValueError(f"{row.label}: the same seeds as {seeds_seen[seeds]}") seeds_seen[seeds] = row.rel_id - ims = read_ims(row.im_path) + ims = read_ims(row.im_path, content, grids) check_finite(row.label, ims) loaded.append(Loaded(row, realisation, ims)) peaks.append((row, ims.peaks)) @@ -833,9 +960,13 @@ def build_database( db.add_sites(new_sites) db.add_site_event(site_event) for x in loaded: - records, psa = build_records(x.row.rel_id, x.ims) - db.add_records(records, pSA=psa) - counts["records"] += len(records) + for records, psa, fas in build_records(x.row.rel_id, x.ims, content): + db.add_records(records, pSA=psa, FAS=fas) + n = len(x.ims.stations) + counts["records"] += n * len(content.components) + counts["psa_ims"] += n * len(content.psa_components) + counts["scalars_ims"] += n * len(content.psa_components) + counts["fas_ims"] += n * len(content.fas_components) counts["events"] += 1 counts["realisations"] += len(loaded) @@ -852,13 +983,12 @@ def build_database( problems = db.validate() if problems: raise RuntimeError(f"validate(): {problems}") - expected = {**counts, "psa_ims": counts["records"], "scalars_ims": counts["records"]} - check_row_counts(db, expected) + check_row_counts(db, counts) check_flags(db, n_flagged) finally: db.close() - checked = spot_check(partial, {row.rel_id: row for row in rows}) + checked = spot_check(partial, {row.rel_id: row for row in rows}, content, grids) print( f"validate() clean; row counts match; {n_flagged} site_event rows carry the " f"pilot's distances; spot check: {checked} records identical to their h5", @@ -866,7 +996,7 @@ def build_database( ) print_peaks(peaks) partial.rename(out_path) - return expected + return counts def main() -> None: @@ -879,14 +1009,31 @@ def main() -> None: parser.add_argument("out", type=Path, help="output .duckdb; must not exist") parser.add_argument("--harvested-at", required=True, help="when the manifest was built from cylc, UTC ISO 8601") parser.add_argument("--script-commit", required=True, help="git commit of this script, for db_meta") + parser.add_argument( + "--components", + default=",".join(DEFAULT_COMPONENTS), + help=f"comma-separated, from {','.join(schema.COMPONENTS)}; eas needs --fas (default: %(default)s)", + ) + parser.add_argument( + "--psa-periods", + choices=("default", "all"), + default="default", + help="default: the 25 in PSA_PERIODS; all: every one of im-calc's 111", + ) + parser.add_argument("--fas", action="store_true", help="add FAS on im-calc's full frequency grid") parser.add_argument( "--memory-limit", help="DuckDB memory_limit, e.g. 12GB. Without it DuckDB sizes itself from " "the node's RAM, not the job's memory allocation.", ) args = parser.parse_args() + content = Content( + components=tuple(args.components.split(",")), + all_periods=args.psa_periods == "all", + fas=args.fas, + ) - print(f"ingesting {args.manifest} into {args.out}", flush=True) + print(f"ingesting {args.manifest} into {args.out}: {content}", flush=True) counts = build_database( args.manifest, args.stations_input, @@ -894,6 +1041,7 @@ def main() -> None: args.out, harvested_at=args.harvested_at, script_commit=args.script_commit, + content=content, memory_limit=args.memory_limit, ) print("row counts:", counts) diff --git a/scripts/tests/test_cs_nshm_2022_ingest.py b/scripts/tests/test_cs_nshm_2022_ingest.py index 1db8714..e1949ff 100644 --- a/scripts/tests/test_cs_nshm_2022_ingest.py +++ b/scripts/tests/test_cs_nshm_2022_ingest.py @@ -207,7 +207,7 @@ def realisation_json(self, original_id: str) -> Path: def realisation(self, original_id: str) -> dict: return json.loads(self.realisation_json(original_id).read_text()) - def build(self) -> dict[str, int]: + def build(self, content=None) -> dict[str, int]: return ingest.build_database( self.manifest, self.stations_input, @@ -215,6 +215,7 @@ def build(self) -> dict[str, int]: self.out, harvested_at="2026-09-23T07:35:00Z", script_commit="abc123", + content=content or ingest.Content(), ) @@ -381,6 +382,7 @@ def test_row_counts_and_meta(campaign: Campaign): "records": 34, "psa_ims": 34, "scalars_ims": 34, + "fas_ims": 0, } assert campaign.out.exists() assert not campaign.out.with_name("out.duckdb.partial").exists() @@ -391,6 +393,8 @@ def test_row_counts_and_meta(campaign: Campaign): assert meta["components"] == "geom,rotd50" assert meta["n_periods"] == "25" assert meta["n_frequencies"] == "0" + assert meta["psa_period_selection"].startswith("25 of im-calc's 111 periods") + assert meta["fas"] == "not included" assert meta["n_realisations"] == "5" assert meta["n_pilot_realisations"] == "2" assert meta["im_calc_harvested_at"] == "2026-09-23T07:35:00Z" @@ -635,3 +639,86 @@ def test_nshm_attributes_from_another_release_are_refused(campaign: Campaign): def test_nan_in_a_dropped_period_is_ignored(campaign: Campaign): _edit_h5(campaign.h5("288271/R2"), "geom", (0, DROPPED_COLUMN), np.nan, ("pSA",)) assert campaign.build()["records"] == 34 + + +# ---- content options ------------------------------------------------------------- + +EVERYTHING = ingest.Content( + components=("000", "090", "ver", "geom", "rotd0", "rotd50", "rotd100", "eas"), + all_periods=True, + fas=True, +) + + +def test_default_content(): + content = ingest.Content() + assert content.components == ("geom", "rotd50") + assert content.psa_components == ("geom", "rotd50") + assert content.fas_components == () + assert EVERYTHING.psa_components == ("000", "090", "ver", "geom", "rotd0", "rotd50", "rotd100") + assert EVERYTHING.fas_components == ("000", "090", "ver", "geom", "eas") + + +@pytest.mark.parametrize( + "kwargs, match", + [ + ({"components": ("geom", "eas")}, "eas has only FAS"), + ({"components": ("geom", "rotd90")}, "unknown components"), + ({"components": ("geom", "geom")}, "listed twice"), + ({"components": ()}, "no components"), + ], +) +def test_bad_content_is_refused(kwargs, match): + with pytest.raises(ValueError, match=match): + ingest.Content(**kwargs) + + +def test_everything_on(campaign: Campaign): + n_stations = sum(len(stations) for stations in STATIONS.values()) # 17 + counts = campaign.build(EVERYTHING) + assert counts == { + "events": 3, + "realisations": 5, + "sites": 5, + "site_event": 11, + "records": 8 * n_stations, + "psa_ims": 7 * n_stations, + "scalars_ims": 7 * n_stations, + "fas_ims": 5 * n_stations, + } + with IMDB(campaign.out) as db: + meta = db.db_meta + assert db.validate() == [] + rel, site = "288271_R2", "2TfEMSL" + ids = {c: db.get_records(rel_ids=[rel], site_ids=[site], component=c).index[0] for c in EVERYTHING.components} + psa = db.get_psa(record_int_ids=list(ids.values())) + fas = db.get_fas(record_int_ids=list(ids.values())) + scalars = db.get_scalars(record_int_ids=list(ids.values())) + assert meta["components"] == "000,090,ver,geom,rotd0,rotd50,rotd100,eas" + assert meta["n_periods"] == "111" + assert meta["n_frequencies"] == "3" + assert meta["psa_period_selection"] == "all 111 of im-calc's periods, 0.01-20 s." + assert meta["fas"].startswith("all 3 of im-calc's FAS frequencies, 0.1-10 Hz, for 000, 090, ver, geom, eas") + with h5py.File(campaign.h5("288271")) as h5: + i = STATIONS["288271"].index(site) + np.testing.assert_array_equal(psa.loc[ids["090"]].to_numpy(dtype=np.float32), h5["pSA/090"][i].astype(np.float32)) + np.testing.assert_array_equal(fas.loc[ids["000"]].to_numpy(dtype=np.float32), h5["FAS/000"][i].astype(np.float32)) + np.testing.assert_array_equal(fas.loc[ids["eas"]].to_numpy(dtype=np.float32), h5["FAS/eas"][i].astype(np.float32)) + assert np.float32(scalars.loc[ids["ver"], "CAV"]) == np.float32(h5["CAV/ver"][i]) + assert ids["rotd50"] not in fas.index # no rotd FAS + assert ids["eas"] not in psa.index and ids["eas"] not in scalars.index # eas is FAS only + assert list(psa.columns) == [str(p) for p in VALID_PERIODS] + + +def test_a_different_fas_grid_stops_the_build(campaign: Campaign): + with h5py.File(campaign.h5("288271"), "r+") as h5: + h5["FAS/frequency"][...] = np.array([0.1, 1.0, 20.0]) + with pytest.raises(ValueError, match="FAS frequency grid differs"): + campaign.build(EVERYTHING) + + +def test_a_different_period_grid_stops_an_all_periods_build(campaign: Campaign): + with h5py.File(campaign.h5("288271"), "r+") as h5: + h5["pSA/period"][0] = 0.015 + with pytest.raises(ValueError, match="0.01 s matches 0"): + campaign.build(ingest.Content(all_periods=True)) From f07bd87e0eabbd59d97e651ad4f9df60c54f4c11 Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Wed, 23 Sep 2026 21:03:04 +1200 Subject: [PATCH 4/8] Address the review of the merged cs_nshm_2022 ingest - Release ibis's in-memory tables after each event. ibis 12 registers every inserted DataFrame with DuckDB and frees them only at exit, so the whole build stayed in memory outside DuckDB's memory_limit. That was ~17 GiB for the lean build, and enough to OOM any full-content build. - Refuse a realisation.json that describes another event, and a manifest whose pilot flags do not follow its campaigns: either would otherwise build a wrong database silently. - Check site_event against the (site, event) pairs that have records, and optionally the flagged count (--expect-flagged). - Find missing files and an empty manifest before building, and name the realisation in every per-realisation error. - Refuse --fas without a FAS component, and store no periods for eas alone. - Correct db_meta: pilot magnitudes are not drawn, pilot-only events are unflagged, the pilot's bottom edges moved up to 2.6 km (not 1.8), and station coordinates need no internal file names. Add a clean source and drop nshm_rupture_name, which only ever held "Rupture ". Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UgVr1htJMAqx8tABJLTnD5 --- scripts/cs_nshm_2022_ingest.py | 136 ++++++++++++++++++---- scripts/tests/test_cs_nshm_2022_ingest.py | 88 +++++++++++++- 2 files changed, 195 insertions(+), 29 deletions(-) diff --git a/scripts/cs_nshm_2022_ingest.py b/scripts/cs_nshm_2022_ingest.py index 2eb68c5..5832f3f 100644 --- a/scripts/cs_nshm_2022_ingest.py +++ b/scripts/cs_nshm_2022_ingest.py @@ -86,7 +86,9 @@ # over all 221 realisations done on 2026-09-21: magnitude identical, hypocentre # within 3e-14 degrees. TOLERANCES = {"magnitude": 1e-9, "hypo_lat": 1e-9, "hypo_lon": 1e-9, "hypo_depth": 1e-6} -MANIFEST_FIELDS = ("rel_id", "event_id", "realisation", "pilot", "original_id", "im_path", "realisation_path") +MANIFEST_FIELDS = ( + "rel_id", "event_id", "realisation", "pilot", "campaign", "original_id", "im_path", "realisation_path", +) # fmt: skip # site_event.metadata on a row whose distances come from the pilot's geometry # although the event's geometry is the main campaign's. PILOT_DISTANCES = json.dumps({"fault_geometry": "pilot"}) @@ -102,11 +104,13 @@ "description": ( "Physics-based ground-motion simulations of New Zealand NSHM 2022 crustal " "ruptures. Each event is one rupture and each realisation is one simulation of " - "it, with its own randomly drawn magnitude, hypocentre, rupture propagation and " - "slip. Most realisations come from the main campaign, up to two per rupture; " - "those from an earlier pilot campaign are marked realisations.metadata.pilot = " - "true." + "it, with its own hypocentre, rupture propagation and slip. Most realisations " + "come from the main campaign, up to two per rupture, which also draws each " + "realisation's magnitude at random. The rest come from an earlier pilot " + "campaign, are marked realisations.metadata.pilot = true, and use the scaling " + "relation's central magnitude instead." ), + "source": "cs_nshm_2022 physics-based simulations; see description.", "realisation_numbering": ( "Each event's realisations are R1..Rn with no gaps. Where an event has a pilot " "realisation, it is always the last. The seeds in realisations.metadata identify " @@ -118,15 +122,17 @@ "rather than a random draw, so leave them out of between-realisation variability " "estimates. Their fault geometry comes from an earlier release of the NSHM fault " "database: the same trace and depths, with the bottom edge placed differently (by " - "a median of 226 m, at most 1.8 km)." + "a median of about 210 m, at most 2.6 km)." ), "distances": ( - "site_event rrup, rjb, rx and ry are computed from the event's own fault geometry " - "(events.source_wkt): the main campaign's where the event has a main-campaign " - "realisation, else the pilot's. At a site that only the pilot's realisation " - f"covers, they come from the pilot's geometry and site_event.metadata is " - f"{PILOT_DISTANCES}. Epicentral and hypocentral distances are not stored; derive " - "them from realisations.hypo_* and the site coordinates." + "site_event rrup, rjb, rx and ry are computed from the event's own fault geometry, " + "events.source_wkt (the main campaign's where the event has a main-campaign " + "realisation, else the pilot's), and have NULL metadata. One exception: in an " + "event that also has main-campaign realisations, at a site only its pilot " + "realisation covers, the distances come from the pilot's own geometry, which is " + f"not events.source_wkt, and site_event.metadata is {PILOT_DISTANCES}. Epicentral " + "and hypocentral distances are not stored; derive them from realisations.hypo_* " + "and the site coordinates." ), "nshm_fault_database": ( "events.metadata.fault_geometry_release names the NSHM2022DB release the event's " @@ -141,9 +147,10 @@ "moment-summed total." ), "station_coordinates": ( - "sites.lat/lon are the canonical stations_input.ll coordinates. The " - "simulation snapped each station to its rupture's grid (~0.1 km away, so " - "up to ~0.25 km apart between ruptures); those coordinates are not stored." + "sites.lat/lon are each station's canonical location. Each simulation moved " + "its stations onto its own computational grid (~0.1 km away, so up to ~0.25 km " + "apart between simulations) and computed the distances there; those " + "coordinates are not stored." ), } @@ -190,6 +197,8 @@ def __post_init__(self) -> None: raise ValueError(f"unknown components {sorted(unknown)}; choose from {schema.COMPONENTS}") if "eas" in self.components and not self.fas: raise ValueError("eas has only FAS, so it needs --fas") + if self.fas and not self.fas_components: + raise ValueError(f"--fas needs a component that has FAS: one of {', '.join(FAS_COMPONENTS)}") @property def psa_components(self) -> tuple[str, ...]: @@ -217,14 +226,21 @@ class Grids: def read_grids(h5_path: Path, content: Content) -> Grids: """The database's grids, taken from one realisation's h5 where the content needs its full grid.""" with h5py.File(h5_path, "r") as h5: - periods = tuple(float(p) for p in h5["pSA"]["period"][:]) if content.all_periods else PSA_PERIODS + if not content.psa_components: + periods: tuple[float, ...] = () + elif content.all_periods: + periods = tuple(float(p) for p in h5["pSA"]["period"][:]) + else: + periods = PSA_PERIODS frequencies = tuple(float(f) for f in h5["FAS"]["frequency"][:]) if content.fas else () return Grids(periods=periods, frequencies=frequencies) def content_meta(content: Content, grids: Grids) -> dict[str, str]: """db_meta entries describing what the options put in.""" - if content.all_periods: + if not content.psa_components: + periods = "not included" + elif content.all_periods: periods = f"all {len(grids.periods)} of im-calc's periods, {min(grids.periods):g}-{max(grids.periods):g} s." else: periods = ( @@ -255,6 +271,8 @@ class ManifestRow: event_id: str realisation: int pilot: bool + campaign: str + """Which simulation campaign; for checks only, never stored.""" original_id: str """Where the result came from; for messages only, never stored.""" im_path: Path @@ -282,17 +300,35 @@ def read_manifest(path: Path) -> list[ManifestRow]: event_id=raw["event_id"], realisation=int(raw["realisation"]), pilot=raw["pilot"] == "true", + campaign=raw["campaign"], original_id=raw["original_id"], im_path=Path(raw["im_path"]), realisation_path=Path(raw["realisation_path"]), ) ) check_manifest(rows) + missing = [ + f"{row.label}: {path}" for row in rows for path in (row.im_path, row.realisation_path) if not path.is_file() + ] + if missing: + raise FileNotFoundError(f"{len(missing)} listed files are missing, e.g. " + "; ".join(missing[:5])) return sorted(rows, key=lambda row: (int(row.event_id), row.realisation)) def check_manifest(rows: Sequence[ManifestRow]) -> None: - """Raise unless ids are consistent and unique, and every event is R1..Rn with at most one pilot, last.""" + """Raise unless ids are consistent and unique, and every event is R1..Rn with at most one pilot, last. + + The pilot flag must also follow the campaign: no campaign may supply both + pilot and non-pilot rows. + """ + if not rows: + raise ValueError("the manifest lists no realisations") + pilot_campaigns = {row.campaign for row in rows if row.pilot} + main_campaigns = {row.campaign for row in rows if not row.pilot} + if pilot_campaigns & main_campaigns: + raise ValueError( + f"campaign {sorted(pilot_campaigns & main_campaigns)} has both pilot and non-pilot rows" + ) seen: set[str] = set() by_event: dict[str, list[ManifestRow]] = {} for row in rows: @@ -673,7 +709,6 @@ def build_event( "domain_wkt": domain_wkt(realisations), "metadata": json.dumps( { - "nshm_rupture_name": defining["metadata"]["name"], "nshmdb_rupture_id": int(nshm["nshmdb_rupture_id"]), "annual_rate": float(nshm["annual_rate"]), "nshm_area_km2": float(nshm["area_km2"]), @@ -728,6 +763,21 @@ def check_against_h5(rel: str, realisation_row: dict, attrs: dict[str, float]) - # ---- build and verify ---------------------------------------------------------- +def load_realisation(row: ManifestRow, content: Content, grids: Grids) -> tuple[dict, RealisationIMs]: + """Read one realisation's realisation.json and h5, naming the realisation in any error. + + Raises unless the realisation.json describes the manifest row's event. + """ + try: + realisation = json.loads(row.realisation_path.read_text()) + name = realisation["metadata"]["name"] + if name != f"Rupture {row.event_id}": + raise ValueError(f"realisation.json describes {name!r}, not event {row.event_id}") + return realisation, read_ims(row.im_path, content, grids) + except (KeyError, OSError, ValueError) as error: + raise type(error)(f"{row.label}: {error}") from error + + @dataclass class Loaded: """One realisation read from disk.""" @@ -737,6 +787,20 @@ class Loaded: ims: RealisationIMs +def release_memtables(db: IMDB) -> None: + """Drop the in-memory tables ibis registers for every insert. + + ibis (12.0) registers each inserted DataFrame with DuckDB as an Arrow table + and releases them only at interpreter exit, so without this every row + written stays in memory, outside DuckDB's memory_limit, for the whole build. + """ + raw = db.con.con + for (name,) in raw.execute( + "SELECT view_name FROM duckdb_views() WHERE temporary AND view_name LIKE 'ibis_pandas_memtable_%'" + ).fetchall(): + raw.unregister(name) + + def check_row_counts(db: IMDB, expected: dict[str, int]) -> None: """Raise unless every table has exactly the expected number of rows.""" for table, n in expected.items(): @@ -771,6 +835,23 @@ def check_row_counts(db: IMDB, expected: dict[str, int]) -> None: """ +# site_event must hold exactly the (site, event) pairs that have records: the +# pairs in site_event but not in records, and those in records but not in site_event. +SITE_EVENT_QUERY = """ +WITH pairs AS (SELECT DISTINCT site_int_id, event_int_id FROM records) +SELECT + (SELECT count(*) FROM site_event ANTI JOIN pairs USING (site_int_id, event_int_id)), + (SELECT count(*) FROM pairs ANTI JOIN site_event USING (site_int_id, event_int_id)) +""" + + +def check_site_event_matches_records(db: IMDB) -> None: + """Raise unless site_event has a row for exactly the (site, event) pairs that have records.""" + extra, missing = db.con.raw_sql(SITE_EVENT_QUERY).fetchone() + if extra or missing: + raise RuntimeError(f"site_event: {extra} rows with no records, {missing} record pairs with no row") + + def check_flags(db: IMDB, n_flagged: int) -> None: """Raise unless the flagged site_event rows are exactly the pilot-only sites of shared events.""" flagged, pilot_only, both = db.con.raw_sql(FLAG_QUERY).fetchone() @@ -866,6 +947,7 @@ def build_database( script_commit: str, content: Content = DEFAULT_CONTENT, memory_limit: str | None = None, + expect_flagged: int | None = None, ) -> dict[str, int]: """Build the database at out_path and return its table row counts. @@ -918,15 +1000,11 @@ def build_database( for index, (event_id, group) in enumerate(by_event.items(), start=1): loaded = [] for row in group: - for path in (row.im_path, row.realisation_path): - if not path.exists(): - raise FileNotFoundError(f"{row.label}: no {path}") - realisation = json.loads(row.realisation_path.read_text()) + realisation, ims = load_realisation(row, content, grids) seeds = json.dumps(realisation["seeds"], sort_keys=True) if seeds in seeds_seen: raise ValueError(f"{row.label}: the same seeds as {seeds_seen[seeds]}") seeds_seen[seeds] = row.rel_id - ims = read_ims(row.im_path, content, grids) check_finite(row.label, ims) loaded.append(Loaded(row, realisation, ims)) peaks.append((row, ims.peaks)) @@ -968,6 +1046,7 @@ def build_database( counts["scalars_ims"] += n * len(content.psa_components) counts["fas_ims"] += n * len(content.fas_components) + release_memtables(db) counts["events"] += 1 counts["realisations"] += len(loaded) counts["sites"] += len(new_sites) @@ -984,7 +1063,10 @@ def build_database( if problems: raise RuntimeError(f"validate(): {problems}") check_row_counts(db, counts) + check_site_event_matches_records(db) check_flags(db, n_flagged) + if expect_flagged is not None and n_flagged != expect_flagged: + raise RuntimeError(f"{n_flagged} site_event rows carry the pilot's distances, expected {expect_flagged}") finally: db.close() @@ -1021,6 +1103,11 @@ def main() -> None: help="default: the 25 in PSA_PERIODS; all: every one of im-calc's 111", ) parser.add_argument("--fas", action="store_true", help="add FAS on im-calc's full frequency grid") + parser.add_argument( + "--expect-flagged", + type=int, + help="fail unless exactly this many site_event rows carry the pilot's distances", + ) parser.add_argument( "--memory-limit", help="DuckDB memory_limit, e.g. 12GB. Without it DuckDB sizes itself from " @@ -1043,6 +1130,7 @@ def main() -> None: script_commit=args.script_commit, content=content, memory_limit=args.memory_limit, + expect_flagged=args.expect_flagged, ) print("row counts:", counts) diff --git a/scripts/tests/test_cs_nshm_2022_ingest.py b/scripts/tests/test_cs_nshm_2022_ingest.py index e1949ff..9f53603 100644 --- a/scripts/tests/test_cs_nshm_2022_ingest.py +++ b/scripts/tests/test_cs_nshm_2022_ingest.py @@ -207,7 +207,7 @@ def realisation_json(self, original_id: str) -> Path: def realisation(self, original_id: str) -> dict: return json.loads(self.realisation_json(original_id).read_text()) - def build(self, content=None) -> dict[str, int]: + def build(self, content=None, **kwargs) -> dict[str, int]: return ingest.build_database( self.manifest, self.stations_input, @@ -216,6 +216,7 @@ def build(self, content=None) -> dict[str, int]: harvested_at="2026-09-23T07:35:00Z", script_commit="abc123", content=content or ingest.Content(), + **kwargs, ) @@ -353,15 +354,18 @@ def _set(rows: list[dict[str, str]], rel: str, **changes: str) -> list[dict[str, "edit, match", [ (lambda rows: _set(rows, "161984_R2", rel_id="161984_R3", realisation="3"), r"realisations \[1, 3\]"), - (lambda rows: _set(_set(rows, "288271_R1", pilot="true"), "288271_R2", pilot="false"), "the pilot is R1"), - (lambda rows: _set(rows, "288271_R1", pilot="true"), "2 pilot realisations"), + (lambda rows: _set(_set(rows, "288271_R1", pilot="true", campaign="run2"), "288271_R2", pilot="false", campaign="run3"), + "the pilot is R1"), + (lambda rows: _set(rows, "288271_R1", pilot="true", campaign="run2"), "2 pilot realisations"), (lambda rows: [*rows, rows[0]], "3_R1: listed twice"), (lambda rows: _set(rows, "3_R1", rel_id="3_R9"), "3_R9: expected 3_R1"), (lambda rows: _set(rows, "3_R1", event_id="03", rel_id="03_R1"), "not an unpadded integer"), (lambda rows: _set(rows, "3_R1", pilot="yes"), "pilot must be true or false"), (lambda rows: [{k: v for k, v in row.items() if k != "original_id"} for row in rows], "missing columns"), + (lambda rows: _set(rows, "3_R1", pilot="false"), r"campaign \['run2'\] has both pilot and non-pilot rows"), ], - ids=["gap", "pilot-not-last", "two-pilots", "duplicate", "rel-id-mismatch", "padded-event", "bad-pilot", "no-column"], + ids=["gap", "pilot-not-last", "two-pilots", "duplicate", "rel-id-mismatch", "padded-event", "bad-pilot", "no-column", + "pilot-flag-vs-campaign"], ) # fmt: skip def test_bad_manifests_are_refused(campaign: Campaign, edit, match): _write_manifest(campaign.manifest, edit(_manifest_rows(campaign))) @@ -369,6 +373,18 @@ def test_bad_manifests_are_refused(campaign: Campaign, edit, match): ingest.read_manifest(campaign.manifest) +def test_an_empty_manifest_is_refused(campaign: Campaign): + campaign.manifest.write_text(campaign.manifest.read_text().splitlines()[0] + "\n") + with pytest.raises(ValueError, match="lists no realisations"): + ingest.read_manifest(campaign.manifest) + + +def test_missing_files_are_found_before_the_build_starts(campaign: Campaign): + campaign.realisation_json("161984/R3").unlink() + with pytest.raises(FileNotFoundError, match=r"161984_R2 \(161984/R3\)"): + ingest.read_manifest(campaign.manifest) + + # ---- the whole build ------------------------------------------------------------- @@ -403,7 +419,8 @@ def test_row_counts_and_meta(campaign: Campaign): assert meta["magnitude_convention"].startswith("BoldM") for key in ("description", "realisation_numbering", "pilot_realisations", "distances", "nshm_fault_database"): assert meta[key] - assert "source" not in meta and "nshmdb" not in meta + assert "nshmdb" not in meta + assert meta["source"] == "cs_nshm_2022 physics-based simulations; see description." for key, value in meta.items(): # nothing researcher-facing names our internals for word in ("run2", "run3", "/gpfs", "BSC", "envelope"): assert word not in value, f"db_meta {key} mentions {word}" @@ -422,6 +439,7 @@ def test_events(campaign: Campaign): assert metadata["161984"]["nshmdb_rupture_id"] == 201460 assert metadata["161984"]["annual_rate"] == 1.86e-04 assert set(metadata["161984"]["segments"]) == {"Clarence", "Kekerengu"} + assert "nshm_rupture_name" not in metadata["161984"] # it only ever held "Rupture " assert {e: m["fault_geometry_release"] for e, m in metadata.items()} == { "3": "pre-v2026.08", "161984": "v2026.08.3", @@ -609,6 +627,53 @@ def test_bad_input_stops_the_build(campaign: Campaign, corrupt, error, match): assert not campaign.out.exists() +def test_files_of_another_event_stop_the_build(campaign: Campaign): + rows = _manifest_rows(campaign) + three = next(row for row in rows if row["rel_id"] == "3_R1") + other = next(row for row in rows if row["rel_id"] == "288271_R2") + for key in ("im_path", "realisation_path"): + three[key], other[key] = other[key], three[key] + _write_manifest(campaign.manifest, rows) + with pytest.raises(ValueError, match=r"3_R1 \(3\): realisation.json describes 'Rupture 288271', not event 3"): + campaign.build() + + +def test_main_realisations_must_share_their_stations(campaign: Campaign): + _write_h5(campaign.h5("161984/R3"), STATIONS["161984/R3"][:3], campaign.realisation("161984/R3"), np.zeros((3, 4)), 9) + with pytest.raises(ValueError, match=r"161984_R2 \(161984/R3\): station list differs from 161984_R1"): + campaign.build() + + +def test_expected_flag_count(campaign: Campaign): + with pytest.raises(RuntimeError, match="2 site_event rows carry the pilot's distances, expected 3"): + campaign.build(expect_flagged=3) + campaign.out.with_name("out.duckdb.partial").unlink() + assert campaign.build(expect_flagged=2)["site_event"] == 11 + + +def test_the_table_checks_catch_a_tampered_database(campaign: Campaign): + campaign.build() + with IMDB(campaign.out, read_only=False) as db: + ingest.check_site_event_matches_records(db) + ingest.check_flags(db, 2) + with pytest.raises(RuntimeError, match="2 flagged rows"): + ingest.check_flags(db, 5) + db.con.raw_sql("DELETE FROM site_event WHERE metadata IS NOT NULL") + with pytest.raises(RuntimeError, match="0 rows with no records, 2 record pairs with no row"): + ingest.check_site_event_matches_records(db) + + +def test_release_memtables(tmp_path: Path): + db = IMDB.create(tmp_path / "m.duckdb", periods=[0.1], components=("geom",)) + db.add_events(pd.DataFrame({"event_id": ["1"]})) + memtables = "SELECT count(*) FROM duckdb_views() WHERE temporary AND view_name LIKE 'ibis_pandas_memtable_%'" + assert db.con.raw_sql(memtables).fetchone()[0] > 0 + ingest.release_memtables(db) + assert db.con.raw_sql(memtables).fetchone()[0] == 0 + db.add_realisations(pd.DataFrame({"rel_id": ["1_R1"], "event_id": ["1"]})) # still writable + db.close() + + def test_the_pilots_distances_are_not_compared_with_the_main_campaigns(campaign: Campaign): main, pilot = _h5_distances(campaign.h5("288271/R2")), _h5_distances(campaign.h5("288271")) assert not np.allclose(main.loc["2TfEMSL"], pilot.loc["2TfEMSL"]) @@ -666,6 +731,7 @@ def test_default_content(): ({"components": ("geom", "rotd90")}, "unknown components"), ({"components": ("geom", "geom")}, "listed twice"), ({"components": ()}, "no components"), + ({"components": ("rotd50",), "fas": True}, "--fas needs a component that has FAS"), ], ) def test_bad_content_is_refused(kwargs, match): @@ -722,3 +788,15 @@ def test_a_different_period_grid_stops_an_all_periods_build(campaign: Campaign): h5["pSA/period"][0] = 0.015 with pytest.raises(ValueError, match="0.01 s matches 0"): campaign.build(ingest.Content(all_periods=True)) + + +def test_eas_alone(campaign: Campaign): + n_stations = sum(len(stations) for stations in STATIONS.values()) + counts = campaign.build(ingest.Content(components=("eas",), fas=True)) + assert (counts["records"], counts["psa_ims"], counts["scalars_ims"], counts["fas_ims"]) == (n_stations, 0, 0, n_stations) + with IMDB(campaign.out) as db: + meta = db.db_meta + assert db.validate() == [] + assert meta["n_periods"] == "0" + assert meta["psa_period_selection"] == "not included" + assert meta["fas"].startswith("all 3 of im-calc's FAS frequencies, 0.1-10 Hz, for eas") From d39293cfd38de3698a58f600effb644fb6fa305e Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Wed, 23 Sep 2026 21:59:30 +1200 Subject: [PATCH 5/8] Store each pilot realisation's own fault geometry in its metadata The pilot's geometry comes from an earlier NSHM database release, so in the 75 events the main campaign also simulated, events.source_wkt is not the geometry the pilot's hypocentre and flagged distances sit on. Every pilot realisation now carries realisations.metadata.fault_geometry_wkt: a MULTIPOLYGON Z with one polygon per plane, (lon, lat, depth km), straight from its realisation.json. It has depth, so distances can be recomputed from it alone. About 2 MB over all 293 pilot realisations. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UgVr1htJMAqx8tABJLTnD5 --- scripts/cs_nshm_2022_ingest.py | 44 ++++++++++++++++------- scripts/tests/test_cs_nshm_2022_ingest.py | 25 ++++++++++++- 2 files changed, 56 insertions(+), 13 deletions(-) diff --git a/scripts/cs_nshm_2022_ingest.py b/scripts/cs_nshm_2022_ingest.py index 5832f3f..5f2a9d4 100644 --- a/scripts/cs_nshm_2022_ingest.py +++ b/scripts/cs_nshm_2022_ingest.py @@ -122,15 +122,18 @@ "rather than a random draw, so leave them out of between-realisation variability " "estimates. Their fault geometry comes from an earlier release of the NSHM fault " "database: the same trace and depths, with the bottom edge placed differently (by " - "a median of about 210 m, at most 2.6 km)." + "a median of about 210 m, at most 2.6 km). Each pilot realisation carries its own " + "fault planes in realisations.metadata.fault_geometry_wkt: a MULTIPOLYGON Z with " + "one polygon per plane, corners as (longitude, latitude, depth in km)." ), "distances": ( "site_event rrup, rjb, rx and ry are computed from the event's own fault geometry, " "events.source_wkt (the main campaign's where the event has a main-campaign " "realisation, else the pilot's), and have NULL metadata. One exception: in an " "event that also has main-campaign realisations, at a site only its pilot " - "realisation covers, the distances come from the pilot's own geometry, which is " - f"not events.source_wkt, and site_event.metadata is {PILOT_DISTANCES}. Epicentral " + "realisation covers, the distances come from the pilot's own geometry " + "(realisations.metadata.fault_geometry_wkt), which is not events.source_wkt, and " + f"site_event.metadata is {PILOT_DISTANCES}. Epicentral " "and hypocentral distances are not stored; derive them from realisations.hypo_* " "and the site coordinates." ), @@ -642,6 +645,20 @@ def trace_wkt(geometries: dict[str, Fault]) -> str: return MultiLineString([LineString(p.corners[:2, [1, 0]]) for p in planes]).wkt +def fault_geometry_wkt(realisation: dict) -> str: + """The realisation's own fault planes, straight from realisation.json. + + A MULTIPOLYGON Z with one polygon per plane, corners as (lon, lat, depth + in km). Unlike events.source_wkt it carries depth, so distances can be + recomputed from it alone. + """ + planes = [] + for entry in realisation["sources"]["source_geometries"].values(): + corners = [[c["longitude"], c["latitude"], c["depth"] / 1000] for c in entry["corners"]] + planes.extend(Polygon(plane) for plane in np.array(corners).reshape(-1, 4, 3)) + return MultiPolygon(planes).wkt + + def domain_polygon(realisation: dict) -> Polygon: corners = realisation["domain"]["domain"] return Polygon([(c["longitude"], c["latitude"]) for c in corners]) @@ -729,6 +746,17 @@ def build_realisation(row: ManifestRow, realisation: dict) -> dict: hypo_lat, hypo_lon, hypo_depth_m = geometries[initial_fault].fault_coordinates_to_wgs_depth_coordinates( np.array([hypocentre_sd["s"], hypocentre_sd["d"]]) ) + metadata = { + "pilot": row.pilot, + "seeds": realisation["seeds"], + "segment_magnitudes_boldm": realisation["magnitudes"]["magnitudes"], + "rupture_causality_tree": causality_tree, + "duration_s": realisation["domain"]["duration"], + } + if row.pilot: + # The pilot's geometry comes from an earlier NSHM database release, so + # it is not the event's wherever the main campaign also simulated it. + metadata["fault_geometry_wkt"] = fault_geometry_wkt(realisation) return { "rel_id": row.rel_id, "event_id": row.event_id, @@ -737,15 +765,7 @@ def build_realisation(row: ManifestRow, realisation: dict) -> dict: "hypo_lat": float(hypo_lat), "hypo_lon": float(hypo_lon), "hypo_depth": float(hypo_depth_m) / 1000, - "metadata": json.dumps( - { - "pilot": row.pilot, - "seeds": realisation["seeds"], - "segment_magnitudes_boldm": realisation["magnitudes"]["magnitudes"], - "rupture_causality_tree": causality_tree, - "duration_s": realisation["domain"]["duration"], - } - ), + "metadata": json.dumps(metadata), } diff --git a/scripts/tests/test_cs_nshm_2022_ingest.py b/scripts/tests/test_cs_nshm_2022_ingest.py index 9f53603..7c0c1cc 100644 --- a/scripts/tests/test_cs_nshm_2022_ingest.py +++ b/scripts/tests/test_cs_nshm_2022_ingest.py @@ -21,7 +21,7 @@ pytest.importorskip("source_modelling") from imdb import IMDB # noqa: E402 -from shapely import unary_union, wkt # noqa: E402 +from shapely import force_2d, unary_union, wkt # noqa: E402 from source_modelling import moment, sources # noqa: E402 SCRIPT = Path(__file__).resolve().parents[1] / "cs_nshm_2022_ingest.py" @@ -491,6 +491,29 @@ def test_realisations_match_their_own_realisation_json(campaign: Campaign): assert json.loads(realisations.loc["161984_R2", "metadata"])["rupture_causality_tree"]["Kekerengu"] is None +def test_pilot_realisations_carry_their_own_fault_geometry(campaign: Campaign): + campaign.build() + with IMDB(campaign.out) as db: + realisations, events, meta = db.get_realisations(), db.get_events(), db.db_meta + for original_id, rel, pilot, _ in REALISATIONS: + metadata = json.loads(realisations.loc[rel, "metadata"]) + if not pilot: + assert "fault_geometry_wkt" not in metadata + continue + stored = wkt.loads(metadata["fault_geometry_wkt"]) + assert stored.geom_type == "MultiPolygon" and stored.has_z + # one polygon per plane, corners as (lon, lat, depth km), straight from realisation.json + geometries = campaign.realisation(original_id)["sources"]["source_geometries"].values() + corners = [c for g in geometries for c in g["corners"]] + expected = np.array([[c["longitude"], c["latitude"], c["depth"] / 1000] for c in corners]).reshape(-1, 4, 3) + got = np.array([np.array(polygon.exterior.coords)[:4] for polygon in stored.geoms]) + np.testing.assert_allclose(got, expected, rtol=0, atol=1e-12) + # In the shared event the pilot's planes are not the event's (the main campaign's) + pilot = wkt.loads(json.loads(realisations.loc["288271_R2", "metadata"])["fault_geometry_wkt"]) + assert not wkt.loads(events.loc["288271", "source_wkt"]).equals(force_2d(pilot)) + assert "realisations.metadata.fault_geometry_wkt" in meta["pilot_realisations"] + + def test_sites_use_canonical_coordinates(campaign: Campaign): campaign.build() with IMDB(campaign.out) as db: From 9adbb0bc3ed7bb3213c6d1eed834db97986281a5 Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Thu, 24 Sep 2026 10:54:16 +1200 Subject: [PATCH 6/8] Reorganise the cs_nshm_2022 ingest to read as a guide The script is meant to sit beside the other ingests in scripts/ as a guide, so it now reads top down: - The module docstring lists the inputs, what is specific to cs_nshm_2022, and a map of the file's six sections. - Every part that exists only because of cs_nshm_2022 (the two campaigns, the pilot's distances, the manifest, grid-snapped stations, the NSHM CSV, the Slurm job's provenance) is marked "# cs_nshm_2022:", so a grep lists what another dataset can drop. - build_database is split into read_event, build_event_rows, write_event and verify_database, and all the checks sit together in section 5. Behaviour is unchanged. Apart from build_database and check_manifest, which now reuses group_by_event, every definition parses to the same code as before, and the 56 tests pass. Building the same three-realisation manifest with this commit and with d39293c, both with the default content and with every component, all periods and FAS, gives identical tables; only db_meta.created_at differs. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01Fo9vmUtWnieoTVWuEaYVii --- scripts/cs_nshm_2022_ingest.py | 1164 ++++++++++++++++++-------------- 1 file changed, 662 insertions(+), 502 deletions(-) diff --git a/scripts/cs_nshm_2022_ingest.py b/scripts/cs_nshm_2022_ingest.py index 5f2a9d4..a6ac5dd 100644 --- a/scripts/cs_nshm_2022_ingest.py +++ b/scripts/cs_nshm_2022_ingest.py @@ -16,39 +16,70 @@ """Ingest the cs_nshm_2022 simulations into one IMDB. -Two simulation campaigns become one dataset. Each event is one NSHM 2022 -crustal rupture and each realisation is one simulation of it: - -- the main campaign: up to two realisations per rupture, each with its own - randomly drawn magnitude, hypocentre, rupture propagation and slip; -- the pilot: one realisation for each of 293 ruptures. Its magnitude is the - scaling relation's central value, and its fault geometry comes from an - earlier NSHM fault database release. - -A manifest CSV (cs_nshm_2022's scripts/build_imdb_manifest.py) lists every -realisation with its rel_id, pilot flag and file paths, so this script knows -nothing about either campaign's layout. Each event's realisations are R1..Rn -with no gaps and the pilot's last. The event_id is the rupture's crustal -`nshm_id` in nshmdb, NOT its `rupture_id`. - -- IMs, distances, vs30 and z1.0/z2.5 come from each realisation's - `intensity_measures.h5`, as written by the workflow's im-calc: netCDF4, one - group per IM, one dataset per component. -- Source and rupture metadata come from its `realisation.json`. -- The rupture's NSHM magnitude and annual rate come from a CSV exported from - nshmdb, which is not available where this runs. -- Canonical site coordinates come from the campaign's `stations_input.ll`. - -An event's geometry, and the distances stored for it, come from its first main -realisation, or from the pilot where the main campaign has none. At a site only -the pilot covers, site_event holds the pilot's distances, flagged in its metadata. - -Trimmed for size by default: only the geom and rotd50 components, pSA on 25 of -im-calc's 111 periods, and no FAS. --components, --psa-periods all and --fas -put the rest in. Empirical-GMM records are never included. +cs_nshm_2022 simulated New Zealand NSHM 2022 crustal ruptures with the ucgmsim +workflow on BSC's MareNostrum 5. In the database each event is one rupture, +and each realisation is one simulation of it. + +This script is longer than most ingests for two reasons. Two simulation +campaigns become one dataset, and it ran unattended as a Slurm job, so it +checks its own output before publishing it. Sections 2 to 4 are the ingest +itself; most of section 5 is that checking. + +Inputs +------ +- A build manifest CSV listing every realisation to ingest, with its rel_id, + pilot flag and file paths. cs_nshm_2022's scripts/build_imdb_manifest.py + writes it (https://github.com/ucgmsim/cs_nshm_2022). +- For each realisation, the `intensity_measures.h5` written by the workflow's + im-calc (IMs, distances, vs30 and z1.0/z2.5 at each station: netCDF4, one + group per IM, one dataset per component) and its `realisation.json` + (source, rupture propagation and simulation domain). +- The campaign's `stations_input.ll`, for each station's canonical location. +- Each rupture's NSHM magnitude and annual rate, in a CSV exported from + nshmdb by cs_nshm_2022's scripts/export_nshm_rupture_attributes.py. + +What is specific to cs_nshm_2022 +-------------------------------- +Code that exists only because of the points below is marked `# cs_nshm_2022:`, +so `grep -n "cs_nshm_2022:" cs_nshm_2022_ingest.py` lists what another dataset +can drop. + +- Two campaigns. The main campaign simulated up to two realisations per + rupture, each with its own randomly drawn magnitude, hypocentre, rupture + propagation and slip. A pilot simulated one realisation for each of 293 + ruptures, at the scaling relation's central magnitude and on fault geometry + from an earlier release of the NSHM fault database. The manifest numbers each + event's realisations R1..Rn with no gaps, the pilot's last. +- One geometry per event. An event's geometry, and the distances stored for + it, come from its first main realisation, or from the pilot where the main + campaign has none. At a site only the pilot reached, site_event holds the + pilot's own distances, flagged in its metadata. Each pilot realisation also + carries its own fault planes. +- Event ids are the rupture's crustal `nshm_id` in nshmdb, NOT its `rupture_id`. +- Stations. Each simulation moved its stations onto its own computational + grid, so the coordinates stored are the canonical ones from + stations_input.ll. Real stations have 3-4 character codes; the virtual + grid's have 7. +- im-calc's pSA grid has 111 periods, with none between 15 s and 20 s. +- nshmdb is not available on BSC, hence the NSHM CSV. + +Layout +------ +1. What goes into the database: the content options. +2. Reading the inputs. +3. Building the rows of each table. +4. The build: each event is read, turned into rows and written, then the + database is checked and published. +5. Checks. +6. Command line. + +By default the database is trimmed for size: only the geom and rotd50 +components, pSA on 25 of im-calc's 111 periods, and no FAS. --components, +--psa-periods all and --fas put the rest in. Empirical-GMM records are never +included. Every magnitude is BoldM (Hanks & Kanamori 1979 eq. 7), the group convention. -nz_sim_validation_ingest.py converts to Mw instead. +Check which convention a script you copy from uses; some convert to Mw. """ import argparse @@ -69,7 +100,10 @@ from source_modelling import moment from source_modelling.sources import Fault -# 0.1 s steps to 1 s, then 1 s steps. im-calc computed no 16-19 s; 20 s is its last. +# ---- Constants ----------------------------------------------------------------------- + +# pSA periods kept by default: 0.1 s steps to 1 s, then 1 s steps. +# cs_nshm_2022: im-calc computed no periods from 16 to 19 s, and 20 s is its last. PSA_PERIODS = ( 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, @@ -79,105 +113,37 @@ # Components im-calc writes for pSA, PGA, PGV and PGD, and for FAS. PSA_COMPONENTS = ("000", "090", "ver", "geom", "rotd0", "rotd50", "rotd100") FAS_COMPONENTS = ("000", "090", "ver", "geom", "eas") +# cs_nshm_2022: every realisation's im-calc used this pSA grid, so any other +# length means a file from a different configuration. N_IM_CALC_PERIODS = 111 SITE_FIELDS = ("vs30", "z1pt0", "z2pt5") DISTANCE_FIELDS = ("rrup", "rjb", "rx", "ry") +TABLES = ("events", "realisations", "sites", "site_event", "records", "psa_ims", "scalars_ims", "fas_ims") # How closely realisation.json must reproduce the h5's own attributes. Measured # over all 221 realisations done on 2026-09-21: magnitude identical, hypocentre # within 3e-14 degrees. TOLERANCES = {"magnitude": 1e-9, "hypo_lat": 1e-9, "hypo_lon": 1e-9, "hypo_depth": 1e-6} + +# cs_nshm_2022: the build manifest's columns. MANIFEST_FIELDS = ( "rel_id", "event_id", "realisation", "pilot", "campaign", "original_id", "im_path", "realisation_path", ) # fmt: skip -# site_event.metadata on a row whose distances come from the pilot's geometry -# although the event's geometry is the main campaign's. +# cs_nshm_2022: site_event.metadata on a row whose distances come from the +# pilot's geometry although the event's geometry is the main campaign's. PILOT_DISTANCES = json.dumps({"fault_geometry": "pilot"}) -# The NSHM fault database releases behind each campaign's geometry. +# cs_nshm_2022: the NSHM fault database releases behind each campaign's geometry. MAIN_RELEASE = "v2026.08.3" MAIN_NSHMDB_SHA256 = "3fe692cb9b22c769b6a9baa6a526b4eed989bd61ec054a344f0e74d82855bc4e" PILOT_RELEASE = "pre-v2026.08" PILOT_NSHMDB_SHA256 = "00e256480618cd15e11fbf744037d037bf3fc2d523fb977ee30e0b84a640bc57" +# cs_nshm_2022: how many realisations the build log lists by largest PGA. N_TOP_PEAKS = 20 -DB_META = { - "dataset_id": "cs_nshm_2022", - "description": ( - "Physics-based ground-motion simulations of New Zealand NSHM 2022 crustal " - "ruptures. Each event is one rupture and each realisation is one simulation of " - "it, with its own hypocentre, rupture propagation and slip. Most realisations " - "come from the main campaign, up to two per rupture, which also draws each " - "realisation's magnitude at random. The rest come from an earlier pilot " - "campaign, are marked realisations.metadata.pilot = true, and use the scaling " - "relation's central magnitude instead." - ), - "source": "cs_nshm_2022 physics-based simulations; see description.", - "realisation_numbering": ( - "Each event's realisations are R1..Rn with no gaps. Where an event has a pilot " - "realisation, it is always the last. The seeds in realisations.metadata identify " - "each simulation uniquely." - ), - "pilot_realisations": ( - "realisations.metadata.pilot = true marks the pilot's realisations, at most one " - "per event. Their magnitude is the magnitude-area scaling relation's central value " - "rather than a random draw, so leave them out of between-realisation variability " - "estimates. Their fault geometry comes from an earlier release of the NSHM fault " - "database: the same trace and depths, with the bottom edge placed differently (by " - "a median of about 210 m, at most 2.6 km). Each pilot realisation carries its own " - "fault planes in realisations.metadata.fault_geometry_wkt: a MULTIPOLYGON Z with " - "one polygon per plane, corners as (longitude, latitude, depth in km)." - ), - "distances": ( - "site_event rrup, rjb, rx and ry are computed from the event's own fault geometry, " - "events.source_wkt (the main campaign's where the event has a main-campaign " - "realisation, else the pilot's), and have NULL metadata. One exception: in an " - "event that also has main-campaign realisations, at a site only its pilot " - "realisation covers, the distances come from the pilot's own geometry " - "(realisations.metadata.fault_geometry_wkt), which is not events.source_wkt, and " - f"site_event.metadata is {PILOT_DISTANCES}. Epicentral " - "and hypocentral distances are not stored; derive them from realisations.hypo_* " - "and the site coordinates." - ), - "nshm_fault_database": ( - "events.metadata.fault_geometry_release names the NSHM2022DB release the event's " - f"geometry comes from: {MAIN_RELEASE} (sha256 {MAIN_NSHMDB_SHA256}), or " - f"{PILOT_RELEASE}, the earlier release the pilot used (sha256 " - f"{PILOT_NSHMDB_SHA256}). events.magnitude and the NSHM attributes in " - f"events.metadata come from {MAIN_RELEASE} for every event." - ), - "magnitude_convention": ( - "BoldM (Hanks & Kanamori 1979 eq. 7). events.magnitude is the NSHM 2022 " - "catalogue magnitude; realisations.magnitude is the realisation's own " - "moment-summed total." - ), - "station_coordinates": ( - "sites.lat/lon are each station's canonical location. Each simulation moved " - "its stations onto its own computational grid (~0.1 km away, so up to ~0.25 km " - "apart between simulations) and computed the distances there; those " - "coordinates are not stored." - ), -} - - -# ---- inputs ------------------------------------------------------------------ - - -def select_period_indices( - available: np.ndarray, wanted: Sequence[float] = PSA_PERIODS -) -> np.ndarray: - """Column of each wanted period in im-calc's period grid. - Raises ValueError unless every wanted period matches exactly one available one. - """ - indices = [] - for period in wanted: - (matches,) = np.nonzero(np.isclose(available, period, rtol=1e-9, atol=0.0)) - if len(matches) != 1: - raise ValueError( - f"pSA period {period} s matches {len(matches)} of im-calc's " - "periods, expected exactly 1" - ) - indices.append(int(matches[0])) - return np.array(indices) +# ---- 1. What goes into the database -------------------------------------------------- +# +# Content is what the command-line options choose; Grids are the pSA periods and +# FAS frequencies the database is created with, which every record then shares. @dataclass(frozen=True) @@ -261,11 +227,22 @@ def content_meta(content: Content, grids: Grids) -> dict[str, str]: return {"psa_period_selection": periods, "fas": fas} +# ---- 2. Reading the inputs ----------------------------------------------------------- +# +# The manifest says what to ingest. Each realisation's intensity_measures.h5 and +# realisation.json hold its data; the station list and the NSHM CSV are shared +# by the whole campaign. The checks run on what is read are in section 5. + + def rel_id(event_id: str, n: int) -> str: """The realisation's stable id.""" return f"{event_id}_R{n}" +# cs_nshm_2022: the two campaigns' results sit in different run directories with +# different layouts, and their realisations are renumbered to be one dataset. A +# manifest names each realisation's files, so this script knows nothing about +# either layout. @dataclass(frozen=True) class ManifestRow: """One realisation to ingest, as listed in the build manifest.""" @@ -318,42 +295,6 @@ def read_manifest(path: Path) -> list[ManifestRow]: return sorted(rows, key=lambda row: (int(row.event_id), row.realisation)) -def check_manifest(rows: Sequence[ManifestRow]) -> None: - """Raise unless ids are consistent and unique, and every event is R1..Rn with at most one pilot, last. - - The pilot flag must also follow the campaign: no campaign may supply both - pilot and non-pilot rows. - """ - if not rows: - raise ValueError("the manifest lists no realisations") - pilot_campaigns = {row.campaign for row in rows if row.pilot} - main_campaigns = {row.campaign for row in rows if not row.pilot} - if pilot_campaigns & main_campaigns: - raise ValueError( - f"campaign {sorted(pilot_campaigns & main_campaigns)} has both pilot and non-pilot rows" - ) - seen: set[str] = set() - by_event: dict[str, list[ManifestRow]] = {} - for row in rows: - if not row.event_id.isdigit() or str(int(row.event_id)) != row.event_id: - raise ValueError(f"{row.rel_id}: event_id {row.event_id!r} is not an unpadded integer") - if row.rel_id != rel_id(row.event_id, row.realisation): - raise ValueError(f"{row.rel_id}: expected {rel_id(row.event_id, row.realisation)} from event_id and realisation") - if row.rel_id in seen: - raise ValueError(f"{row.rel_id}: listed twice") - seen.add(row.rel_id) - by_event.setdefault(row.event_id, []).append(row) - for event_id, group in by_event.items(): - numbers = sorted(row.realisation for row in group) - if numbers != list(range(1, len(group) + 1)): - raise ValueError(f"event {event_id}: realisations {numbers}, expected 1..{len(group)}") - pilots = [row for row in group if row.pilot] - if len(pilots) > 1: - raise ValueError(f"event {event_id}: {len(pilots)} pilot realisations, expected at most 1") - if pilots and pilots[0].realisation != len(group): - raise ValueError(f"event {event_id}: the pilot is R{pilots[0].realisation}, but it must be the last, R{len(group)}") - - def group_by_event(rows: Sequence[ManifestRow]) -> dict[str, list[ManifestRow]]: """Rows per event, keeping the input order.""" grouped: dict[str, list[ManifestRow]] = {} @@ -362,11 +303,9 @@ def group_by_event(rows: Sequence[ManifestRow]) -> dict[str, list[ManifestRow]]: return grouped -def md5_of(path: Path) -> str: - """Hex md5 of a file.""" - return hashlib.md5(path.read_bytes()).hexdigest() - - +# cs_nshm_2022: each simulation moved its stations onto its own computational +# grid (about 0.1 km) and wrote those coordinates into its h5, so they differ +# between simulations. sites stores each station's canonical location instead. def load_stations_input(path: Path) -> pd.DataFrame: """Canonical station coordinates (`lon lat name` per line), indexed by name.""" df = pd.read_csv( @@ -378,6 +317,17 @@ def load_stations_input(path: Path) -> pd.DataFrame: return df.set_index("site_id") +def is_real_station(site_id: str) -> bool: + """GeoNet station codes are 3-4 characters; the virtual grid's codes are 7.""" + # cs_nshm_2022: nothing in the inputs marks a station as real, so the name + # length decides. Check how your own station list names its stations. + return len(site_id) != 7 + + +# cs_nshm_2022: nshmdb is not available on BSC, where this ran, so each +# rupture's NSHM attributes come from a CSV exported beforehand. The CSV is +# indexed by the crustal nshm_id, which is this dataset's event_id; nshmdb's own +# rupture_id is a different number and goes into events.metadata. def load_nshm_attributes(path: Path) -> pd.DataFrame: """Rows of export_nshm_rupture_attributes.py's CSV, indexed by rupture id. @@ -392,14 +342,6 @@ def load_nshm_attributes(path: Path) -> pd.DataFrame: return nshm -def is_real_station(site_id: str) -> bool: - """GeoNet station codes are 3-4 characters; the virtual grid's codes are 7.""" - return len(site_id) != 7 - - -# ---- intensity_measures.h5 --------------------------------------------------- - - @dataclass class RealisationIMs: """One realisation's intensity_measures.h5, reduced to what this database keeps.""" @@ -434,6 +376,25 @@ def _peak(values: np.ndarray, stations: np.ndarray) -> tuple[float, str]: return float(values[i]), station.decode() if isinstance(station, bytes) else str(station) +def select_period_indices( + available: np.ndarray, wanted: Sequence[float] = PSA_PERIODS +) -> np.ndarray: + """Column of each wanted period in im-calc's period grid. + + Raises ValueError unless every wanted period matches exactly one available one. + """ + indices = [] + for period in wanted: + (matches,) = np.nonzero(np.isclose(available, period, rtol=1e-9, atol=0.0)) + if len(matches) != 1: + raise ValueError( + f"pSA period {period} s matches {len(matches)} of im-calc's " + "periods, expected exactly 1" + ) + indices.append(int(matches[0])) + return np.array(indices) + + def read_ims( h5_path: Path, content: Content = DEFAULT_CONTENT, grids: Grids = DEFAULT_GRIDS ) -> RealisationIMs: @@ -478,161 +439,49 @@ def read_ims( ) -def check_finite(rel: str, ims: RealisationIMs) -> None: - """Raise if any kept IM, site or distance value is NaN or infinite.""" - arrays = { - **{f"pSA/{c}": values for c, values in ims.psa.items()}, - **{f"FAS/{c}": values for c, values in ims.fas.items()}, - **{f"{im}/{c}": values for im, by_c in ims.scalars.items() for c, values in by_c.items()}, - **ims.site, - **ims.distances, - } - for name, values in arrays.items(): - bad = ~np.isfinite(values) - if bad.any(): - station = ims.stations[np.argwhere(bad)[0][0]] - raise ValueError( - f"{rel}: {name} has {int(bad.sum())} non-finite values, the first at {station}" - ) +@dataclass +class Loaded: + """One realisation read from disk.""" + row: ManifestRow + realisation: dict + ims: RealisationIMs -def check_matches_first_realisation( - rel: str, first_rel: str, ims: RealisationIMs, first: RealisationIMs -) -> None: - """Raise unless stations, site fields and distances equal the event's first main realisation's.""" - if not np.array_equal(ims.stations, first.stations): - raise ValueError(f"{rel}: station list differs from {first_rel}") - for mine, theirs in ((ims.site, first.site), (ims.distances, first.distances)): - for key, values in mine.items(): - if not np.array_equal(values, theirs[key]): - raise ValueError( - f"{rel}: {key} differs from {first_rel}; site and distance fields " - "must be identical across an event's main realisations" - ) +def load_realisation(row: ManifestRow, content: Content, grids: Grids) -> tuple[dict, RealisationIMs]: + """Read one realisation's realisation.json and h5, naming the realisation in any error. -Batch = tuple[pd.DataFrame, np.ndarray | None, np.ndarray | None] + Raises unless the realisation.json describes the manifest row's event. + """ + try: + realisation = json.loads(row.realisation_path.read_text()) + # cs_nshm_2022: the workflow names each realisation "Rupture ", + # which catches a manifest row that points at another event's files. + name = realisation["metadata"]["name"] + if name != f"Rupture {row.event_id}": + raise ValueError(f"realisation.json describes {name!r}, not event {row.event_id}") + return realisation, read_ims(row.im_path, content, grids) + except (KeyError, OSError, ValueError) as error: + raise type(error)(f"{row.label}: {error}") from error -def build_records(rel: str, ims: RealisationIMs, content: Content = DEFAULT_CONTENT) -> list[Batch]: - """IMDB records (one per station and component), as add_records batches of (records, pSA, FAS). +# ---- 3. Building the rows ------------------------------------------------------------ +# +# One builder per table. The geometry helpers are the same as in the other +# ingest scripts that read realisation.json. - `eas` has only FAS, so its records form a batch of their own with no pSA - and no scalar columns, and add_records writes no psa_ims or scalars_ims - rows for them. A rotd record's FAS row is all NaN, so it gets no fas_ims row. - """ - n = len(ims.stations) - frames = [] - for component in content.psa_components: - frame = {"rel_id": rel, "site_id": ims.stations, "component": component, "kind": "simulated"} - for im, by_component in ims.scalars.items(): - frame[im] = by_component.get(component, np.full(n, np.nan, dtype=np.float32)) - frames.append(pd.DataFrame(frame)) - batches: list[Batch] = [] - if frames: - psa = np.concatenate([ims.psa[c] for c in content.psa_components]) - fas = None - if content.fas: - no_fas = np.full((n, ims.n_frequencies), np.nan, dtype=np.float32) - fas = np.concatenate([ims.fas.get(c, no_fas) for c in content.psa_components]) - batches.append((pd.concat(frames, ignore_index=True), psa, fas)) - if "eas" in content.fas_components: - eas = pd.DataFrame({"rel_id": rel, "site_id": ims.stations, "component": "eas", "kind": "simulated"}) - batches.append((eas, None, ims.fas["eas"])) - return batches +def fault_geometries(realisation: dict) -> dict[str, Fault]: + return { + name: Fault.from_corners( + np.array([[c["latitude"], c["longitude"], c["depth"]] for c in entry["corners"]]).reshape(-1, 4, 3) + ) + for name, entry in realisation["sources"]["source_geometries"].items() + } -def _site_event_rows(event_id: str, stations: np.ndarray, distances: dict[str, np.ndarray], metadata: str | None) -> pd.DataFrame: - return pd.DataFrame({"site_id": stations, "event_id": event_id, **distances, "metadata": metadata}) - -def build_site_event( - event_id: str, main: RealisationIMs | None, pilot: RealisationIMs | None -) -> tuple[pd.DataFrame, int]: - """site_event rows for one event, and how many carry the pilot's distances. - - The main campaign's distances wherever its realisations cover the site: they - all share one station list and one set of distances, so the first covers - them all. The pilot's only where nothing else exists, flagged when the - event's geometry is the main campaign's and unflagged when the event is the - pilot's alone. - """ - if main is None: - if pilot is None: - raise ValueError(f"event {event_id}: no realisations") - return _site_event_rows(event_id, pilot.stations, pilot.distances, None), 0 - rows = _site_event_rows(event_id, main.stations, main.distances, None) - if pilot is None: - return rows, 0 - gap = ~pd.Index(pilot.stations).isin(main.stations) - flagged = _site_event_rows( - event_id, - pilot.stations[gap], - {key: values[gap] for key, values in pilot.distances.items()}, - PILOT_DISTANCES, - ) - return pd.concat([rows, flagged], ignore_index=True), int(gap.sum()) - - -class SiteRegistry: - """The sites added so far, with the vs30/z1pt0/z2pt5 first seen for each.""" - - def __init__(self, stations_input: pd.DataFrame) -> None: - self.stations_input = stations_input - self.seen: pd.DataFrame | None = None - - def new_sites(self, label: str, ims: RealisationIMs) -> pd.DataFrame: - """Rows for stations not added yet; raise if a known one disagrees.""" - here = pd.DataFrame(ims.site, index=pd.Index(ims.stations, name="site_id")) - if self.seen is None: - new = here - else: - known = here.index.isin(self.seen.index) - if known.any(): - earlier = self.seen.loc[here.index[known]] - mismatch = (here[known] != earlier).any(axis=1) - if mismatch.any(): - site = mismatch[mismatch].index[0] - raise ValueError( - f"{label}: {site} has {here.loc[site].to_dict()}, but an " - f"earlier realisation gave {earlier.loc[site].to_dict()}" - ) - new = here[~known] - missing = new.index.difference(self.stations_input.index) - if len(missing): - raise ValueError( - f"{label}: {len(missing)} stations missing from " - f"stations_input.ll, e.g. {list(missing[:5])}" - ) - self.seen = new if self.seen is None else pd.concat([self.seen, new]) - coords = self.stations_input.loc[new.index] - return pd.DataFrame( - { - "site_id": new.index.to_numpy(), - "lat": coords["lat"].to_numpy(), - "lon": coords["lon"].to_numpy(), - "vs30": new["vs30"].to_numpy(), - "z1p0": new["z1pt0"].to_numpy(), - "z2p5": new["z2pt5"].to_numpy(), - "is_real": [is_real_station(site) for site in new.index], - } - ) - - -# ---- realisation.json (geometry helpers as in nz_sim_validation_ingest.py) ---- - - -def fault_geometries(realisation: dict) -> dict[str, Fault]: - return { - name: Fault.from_corners( - np.array([[c["latitude"], c["longitude"], c["depth"]] for c in entry["corners"]]).reshape(-1, 4, 3) - ) - for name, entry in realisation["sources"]["source_geometries"].items() - } - - -def initial_fault_name(causality_tree: dict[str, str | None]) -> str: - return next(name for name, parent in causality_tree.items() if parent is None) +def initial_fault_name(causality_tree: dict[str, str | None]) -> str: + return next(name for name, parent in causality_tree.items() if parent is None) def source_wkt(geometries: dict[str, Fault]) -> str: @@ -685,6 +534,11 @@ def build_event( unflagged distances share; the domain is the union over `realisations`. Magnitude and causality tree differ between realisations, so those live on the realisation instead. + + cs_nshm_2022: with one campaign, `defining` would be any realisation and + the domain its own. Here the pilot's geometry and domain differ from the + main campaign's, so the caller picks `defining`, and the union keeps every + site with data inside the event's domain. """ geometries = fault_geometries(defining) rakes = defining["rakes"]["rakes"] @@ -754,8 +608,9 @@ def build_realisation(row: ManifestRow, realisation: dict) -> dict: "duration_s": realisation["domain"]["duration"], } if row.pilot: - # The pilot's geometry comes from an earlier NSHM database release, so - # it is not the event's wherever the main campaign also simulated it. + # cs_nshm_2022: the pilot's geometry comes from an earlier NSHM database + # release, so it is not the event's wherever the main campaign also + # simulated the rupture. metadata["fault_geometry_wkt"] = fault_geometry_wkt(realisation) return { "rel_id": row.rel_id, @@ -769,42 +624,188 @@ def build_realisation(row: ManifestRow, realisation: dict) -> dict: } -def check_against_h5(rel: str, realisation_row: dict, attrs: dict[str, float]) -> None: - """Raise unless realisation.json reproduces the h5's magnitude and hypocentre.""" - for key, tolerance in TOLERANCES.items(): - ours, theirs = realisation_row[key], attrs[key] - if not abs(ours - theirs) <= tolerance: +class SiteRegistry: + """The sites added so far, with the vs30/z1pt0/z2pt5 first seen for each.""" + + def __init__(self, stations_input: pd.DataFrame) -> None: + self.stations_input = stations_input + self.seen: pd.DataFrame | None = None + + def new_sites(self, label: str, ims: RealisationIMs) -> pd.DataFrame: + """Rows for stations not added yet; raise if a known one disagrees.""" + here = pd.DataFrame(ims.site, index=pd.Index(ims.stations, name="site_id")) + if self.seen is None: + new = here + else: + known = here.index.isin(self.seen.index) + if known.any(): + earlier = self.seen.loc[here.index[known]] + mismatch = (here[known] != earlier).any(axis=1) + if mismatch.any(): + site = mismatch[mismatch].index[0] + raise ValueError( + f"{label}: {site} has {here.loc[site].to_dict()}, but an " + f"earlier realisation gave {earlier.loc[site].to_dict()}" + ) + new = here[~known] + missing = new.index.difference(self.stations_input.index) + if len(missing): raise ValueError( - f"{rel}: {key} from realisation.json is {ours!r}, but " - f"intensity_measures.h5 says {theirs!r}" + f"{label}: {len(missing)} stations missing from " + f"stations_input.ll, e.g. {list(missing[:5])}" ) + self.seen = new if self.seen is None else pd.concat([self.seen, new]) + coords = self.stations_input.loc[new.index] + return pd.DataFrame( + { + "site_id": new.index.to_numpy(), + "lat": coords["lat"].to_numpy(), + "lon": coords["lon"].to_numpy(), + "vs30": new["vs30"].to_numpy(), + "z1p0": new["z1pt0"].to_numpy(), + "z2p5": new["z2pt5"].to_numpy(), + "is_real": [is_real_station(site) for site in new.index], + } + ) -# ---- build and verify ---------------------------------------------------------- +def _site_event_rows( + event_id: str, stations: np.ndarray, distances: dict[str, np.ndarray], metadata: str | None +) -> pd.DataFrame: + return pd.DataFrame({"site_id": stations, "event_id": event_id, **distances, "metadata": metadata}) -def load_realisation(row: ManifestRow, content: Content, grids: Grids) -> tuple[dict, RealisationIMs]: - """Read one realisation's realisation.json and h5, naming the realisation in any error. +def build_site_event( + event_id: str, main: RealisationIMs | None, pilot: RealisationIMs | None +) -> tuple[pd.DataFrame, int]: + """site_event rows for one event, and how many carry the pilot's distances. - Raises unless the realisation.json describes the manifest row's event. + The main campaign's distances wherever its realisations cover the site: they + all share one station list and one set of distances, so the first covers + them all. The pilot's only where nothing else exists, flagged when the + event's geometry is the main campaign's and unflagged when the event is the + pilot's alone. + + cs_nshm_2022: with one campaign, this would be one realisation's distances. """ - try: - realisation = json.loads(row.realisation_path.read_text()) - name = realisation["metadata"]["name"] - if name != f"Rupture {row.event_id}": - raise ValueError(f"realisation.json describes {name!r}, not event {row.event_id}") - return realisation, read_ims(row.im_path, content, grids) - except (KeyError, OSError, ValueError) as error: - raise type(error)(f"{row.label}: {error}") from error + if main is None: + if pilot is None: + raise ValueError(f"event {event_id}: no realisations") + return _site_event_rows(event_id, pilot.stations, pilot.distances, None), 0 + rows = _site_event_rows(event_id, main.stations, main.distances, None) + if pilot is None: + return rows, 0 + gap = ~pd.Index(pilot.stations).isin(main.stations) + flagged = _site_event_rows( + event_id, + pilot.stations[gap], + {key: values[gap] for key, values in pilot.distances.items()}, + PILOT_DISTANCES, + ) + return pd.concat([rows, flagged], ignore_index=True), int(gap.sum()) -@dataclass -class Loaded: - """One realisation read from disk.""" +Batch = tuple[pd.DataFrame, np.ndarray | None, np.ndarray | None] - row: ManifestRow - realisation: dict - ims: RealisationIMs + +def build_records(rel: str, ims: RealisationIMs, content: Content = DEFAULT_CONTENT) -> list[Batch]: + """IMDB records (one per station and component), as add_records batches of (records, pSA, FAS). + + `eas` has only FAS, so its records form a batch of their own with no pSA + and no scalar columns, and add_records writes no psa_ims or scalars_ims + rows for them. A rotd record's FAS row is all NaN, so it gets no fas_ims row. + """ + n = len(ims.stations) + frames = [] + for component in content.psa_components: + frame = {"rel_id": rel, "site_id": ims.stations, "component": component, "kind": "simulated"} + for im, by_component in ims.scalars.items(): + frame[im] = by_component.get(component, np.full(n, np.nan, dtype=np.float32)) + frames.append(pd.DataFrame(frame)) + batches: list[Batch] = [] + if frames: + psa = np.concatenate([ims.psa[c] for c in content.psa_components]) + fas = None + if content.fas: + no_fas = np.full((n, ims.n_frequencies), np.nan, dtype=np.float32) + fas = np.concatenate([ims.fas.get(c, no_fas) for c in content.psa_components]) + batches.append((pd.concat(frames, ignore_index=True), psa, fas)) + if "eas" in content.fas_components: + eas = pd.DataFrame({"rel_id": rel, "site_id": ims.stations, "component": "eas", "kind": "simulated"}) + batches.append((eas, None, ims.fas["eas"])) + return batches + + +# ---- 4. The build -------------------------------------------------------------------- +# +# One event at a time: read its realisations, build its rows, write them. The +# database is written as .partial and renamed only after section 5's +# checks pass, so a failed or killed job never leaves a finished-looking file. + +# What the database says about itself. Every dataset needs its own; most of this +# one explains cs_nshm_2022's two campaigns to the people who will query it. +DB_META = { + "dataset_id": "cs_nshm_2022", + "description": ( + "Physics-based ground-motion simulations of New Zealand NSHM 2022 crustal " + "ruptures. Each event is one rupture and each realisation is one simulation of " + "it, with its own hypocentre, rupture propagation and slip. Most realisations " + "come from the main campaign, up to two per rupture, which also draws each " + "realisation's magnitude at random. The rest come from an earlier pilot " + "campaign, are marked realisations.metadata.pilot = true, and use the scaling " + "relation's central magnitude instead." + ), + "source": "cs_nshm_2022 physics-based simulations; see description.", + "realisation_numbering": ( + "Each event's realisations are R1..Rn with no gaps. Where an event has a pilot " + "realisation, it is always the last. The seeds in realisations.metadata identify " + "each simulation uniquely." + ), + "pilot_realisations": ( + "realisations.metadata.pilot = true marks the pilot's realisations, at most one " + "per event. Their magnitude is the magnitude-area scaling relation's central value " + "rather than a random draw, so leave them out of between-realisation variability " + "estimates. Their fault geometry comes from an earlier release of the NSHM fault " + "database: the same trace and depths, with the bottom edge placed differently (by " + "a median of about 210 m, at most 2.6 km). Each pilot realisation carries its own " + "fault planes in realisations.metadata.fault_geometry_wkt: a MULTIPOLYGON Z with " + "one polygon per plane, corners as (longitude, latitude, depth in km)." + ), + "distances": ( + "site_event rrup, rjb, rx and ry are computed from the event's own fault geometry, " + "events.source_wkt (the main campaign's where the event has a main-campaign " + "realisation, else the pilot's), and have NULL metadata. One exception: in an " + "event that also has main-campaign realisations, at a site only its pilot " + "realisation covers, the distances come from the pilot's own geometry " + "(realisations.metadata.fault_geometry_wkt), which is not events.source_wkt, and " + f"site_event.metadata is {PILOT_DISTANCES}. Epicentral " + "and hypocentral distances are not stored; derive them from realisations.hypo_* " + "and the site coordinates." + ), + "nshm_fault_database": ( + "events.metadata.fault_geometry_release names the NSHM2022DB release the event's " + f"geometry comes from: {MAIN_RELEASE} (sha256 {MAIN_NSHMDB_SHA256}), or " + f"{PILOT_RELEASE}, the earlier release the pilot used (sha256 " + f"{PILOT_NSHMDB_SHA256}). events.magnitude and the NSHM attributes in " + f"events.metadata come from {MAIN_RELEASE} for every event." + ), + "magnitude_convention": ( + "BoldM (Hanks & Kanamori 1979 eq. 7). events.magnitude is the NSHM 2022 " + "catalogue magnitude; realisations.magnitude is the realisation's own " + "moment-summed total." + ), + "station_coordinates": ( + "sites.lat/lon are each station's canonical location. Each simulation moved " + "its stations onto its own computational grid (~0.1 km away, so up to ~0.25 km " + "apart between simulations) and computed the distances there; those " + "coordinates are not stored." + ), +} + + +def md5_of(path: Path) -> str: + """Hex md5 of a file.""" + return hashlib.md5(path.read_bytes()).hexdigest() def release_memtables(db: IMDB) -> None: @@ -813,6 +814,7 @@ def release_memtables(db: IMDB) -> None: ibis (12.0) registers each inserted DataFrame with DuckDB as an Arrow table and releases them only at interpreter exit, so without this every row written stays in memory, outside DuckDB's memory_limit, for the whole build. + Any large ingest needs it until imdb releases them itself. """ raw = db.con.con for (name,) in raw.execute( @@ -821,17 +823,310 @@ def release_memtables(db: IMDB) -> None: raw.unregister(name) -def check_row_counts(db: IMDB, expected: dict[str, int]) -> None: - """Raise unless every table has exactly the expected number of rows.""" - for table, n in expected.items(): - found = int(db.con.table(table).count().to_pandas()) - if found != n: - raise RuntimeError(f"{table}: {found} rows, expected {n}") +def read_event( + group: Sequence[ManifestRow], content: Content, grids: Grids, seeds_seen: dict[str, str] +) -> list[Loaded]: + """Read one event's realisations, checking each on its own. + `seeds_seen` maps each set of seeds read so far to its rel_id, which + catches the same simulation listed twice under different ids. + """ + loaded = [] + for row in group: + realisation, ims = load_realisation(row, content, grids) + seeds = json.dumps(realisation["seeds"], sort_keys=True) + if seeds in seeds_seen: + raise ValueError(f"{row.label}: the same seeds as {seeds_seen[seeds]}") + seeds_seen[seeds] = row.rel_id + check_finite(row.label, ims) + loaded.append(Loaded(row, realisation, ims)) + # cs_nshm_2022: logged so extreme near-fault values could be judged + # before the database was shared. + print(f"peak {peak_line(row, ims.peaks)}", flush=True) + return loaded + + +@dataclass +class EventRows: + """Everything one event adds to the database apart from its records.""" + + event: dict + realisations: list[dict] + sites: pd.DataFrame + """Sites no earlier event reached.""" + site_event: pd.DataFrame + n_flagged: int + """How many site_event rows carry the pilot's distances.""" + + +def build_event_rows(event_id: str, loaded: Sequence[Loaded], nshm: pd.Series, sites: SiteRegistry) -> EventRows: + """The events, realisations, sites and site_event rows for one event.""" + # cs_nshm_2022: the event's geometry and unflagged distances come from its + # first main realisation, or from the pilot where the main campaign never + # simulated the rupture. + mains = [sim for sim in loaded if not sim.row.pilot] + pilots = [sim for sim in loaded if sim.row.pilot] + defining = mains[0] if mains else pilots[0] + event = build_event( + event_id, + defining.realisation, + [sim.realisation for sim in loaded], + nshm, + MAIN_RELEASE if mains else PILOT_RELEASE, + ) + realisation_rows = [] + for sim in loaded: + realisation_row = build_realisation(sim.row, sim.realisation) + check_against_h5(sim.row.label, realisation_row, sim.ims.attrs) + realisation_rows.append(realisation_row) + for sim in mains[1:]: + check_matches_first_realisation(sim.row.label, mains[0].row.label, sim.ims, mains[0].ims) + new_sites = pd.concat([sites.new_sites(sim.row.label, sim.ims) for sim in loaded], ignore_index=True) + site_event, n_flagged = build_site_event( + event_id, mains[0].ims if mains else None, pilots[0].ims if pilots else None + ) + return EventRows(event, realisation_rows, new_sites, site_event, n_flagged) + + +def write_event(db: IMDB, rows: EventRows, loaded: Sequence[Loaded], content: Content) -> dict[str, int]: + """Write one event, and return how many rows each table should have gained.""" + db.add_events(pd.DataFrame([rows.event])) + db.add_realisations(pd.DataFrame(rows.realisations)) + if len(rows.sites): + db.add_sites(rows.sites) + db.add_site_event(rows.site_event) + for sim in loaded: + for records, psa, fas in build_records(sim.row.rel_id, sim.ims, content): + db.add_records(records, pSA=psa, FAS=fas) + release_memtables(db) + + n_stations = sum(len(sim.ims.stations) for sim in loaded) + return { + "events": 1, + "realisations": len(rows.realisations), + "sites": len(rows.sites), + "site_event": len(rows.site_event), + "records": n_stations * len(content.components), + "psa_ims": n_stations * len(content.psa_components), + "scalars_ims": n_stations * len(content.psa_components), + "fas_ims": n_stations * len(content.fas_components), + } + + +def build_database( + manifest_path: Path, + stations_input: Path, + nshm_csv: Path, + out_path: Path, + harvested_at: str, + script_commit: str, + content: Content = DEFAULT_CONTENT, + memory_limit: str | None = None, + expect_flagged: int | None = None, +) -> dict[str, int]: + """Build the database at out_path and return its table row counts. + + Writes `.partial` and renames it only once validation, the row + counts, the site_event flags and the spot check have all passed. A failure + leaves the .partial behind for inspection, and the next run refuses to + start until it is removed. + """ + if out_path.exists(): + raise FileExistsError(f"{out_path} already exists") + partial = out_path.with_name(out_path.name + ".partial") + if partial.exists(): + raise FileExistsError(f"{partial} exists: a previous run died; remove it first") + + rows = read_manifest(manifest_path) + by_event = group_by_event(rows) + stations = load_stations_input(stations_input) + nshm = load_nshm_attributes(nshm_csv) + missing = sorted(set(by_event) - set(nshm.index), key=int) + if missing: + raise ValueError(f"ruptures missing from {nshm_csv}: {missing}") + + grids = read_grids(rows[0].im_path, content) + db = IMDB.create( + partial, + periods=list(grids.periods), + frequencies=list(grids.frequencies), + components=content.components, + db_meta={ + **DB_META, + **content_meta(content, grids), + "n_realisations": str(len(rows)), + "n_pilot_realisations": str(sum(row.pilot for row in rows)), + # cs_nshm_2022: provenance the Slurm job supplies: when the manifest + # was taken from the workflow, and this script's git commit. + "im_calc_harvested_at": harvested_at, + "build_manifest_md5": md5_of(manifest_path), + "ingest_script_commit": script_commit, + }, + ) + counts = dict.fromkeys(TABLES, 0) + n_flagged = 0 + seeds_seen: dict[str, str] = {} + peaks: list[tuple[ManifestRow, dict[str, tuple[float, str]]]] = [] + started = time.monotonic() + try: + if memory_limit: + db.con.raw_sql(f"SET memory_limit = '{memory_limit}'") + sites = SiteRegistry(stations) + for index, (event_id, group) in enumerate(by_event.items(), start=1): + loaded = read_event(group, content, grids, seeds_seen) + event_rows = build_event_rows(event_id, loaded, nshm.loc[event_id], sites) + for table, n in write_event(db, event_rows, loaded, content).items(): + counts[table] += n + n_flagged += event_rows.n_flagged + peaks.extend((sim.row, sim.ims.peaks) for sim in loaded) + if index % 25 == 0 or index == len(by_event): + print( + f"{index}/{len(by_event)} events, {counts['records']} records, " + f"{time.monotonic() - started:.0f} s", + flush=True, + ) + verify_database(db, counts, n_flagged, expect_flagged) + finally: + db.close() + + checked = spot_check(partial, {row.rel_id: row for row in rows}, content, grids) + print( + f"validate() clean; row counts match; {n_flagged} site_event rows carry the " + f"pilot's distances; spot check: {checked} records identical to their h5", + flush=True, + ) + print_peaks(peaks) + partial.rename(out_path) + return counts + + +# ---- 5. Checks ----------------------------------------------------------------------- +# +# The ingest ran unattended, so it checks its inputs as it reads them and its +# output before publishing it. Most of these are worth keeping in any ingest; +# those marked cs_nshm_2022 concern the manifest and the pilot's distances. + + +def check_manifest(rows: Sequence[ManifestRow]) -> None: + """Raise unless ids are consistent and unique, and every event is R1..Rn with at most one pilot, last. + + The pilot flag must also follow the campaign: no campaign may supply both + pilot and non-pilot rows. + """ + if not rows: + raise ValueError("the manifest lists no realisations") + # cs_nshm_2022: the pilot flag is decided per campaign, never per realisation. + pilot_campaigns = {row.campaign for row in rows if row.pilot} + main_campaigns = {row.campaign for row in rows if not row.pilot} + if pilot_campaigns & main_campaigns: + raise ValueError( + f"campaign {sorted(pilot_campaigns & main_campaigns)} has both pilot and non-pilot rows" + ) + seen: set[str] = set() + for row in rows: + if not row.event_id.isdigit() or str(int(row.event_id)) != row.event_id: + raise ValueError(f"{row.rel_id}: event_id {row.event_id!r} is not an unpadded integer") + if row.rel_id != rel_id(row.event_id, row.realisation): + raise ValueError( + f"{row.rel_id}: expected {rel_id(row.event_id, row.realisation)} " + "from event_id and realisation" + ) + if row.rel_id in seen: + raise ValueError(f"{row.rel_id}: listed twice") + seen.add(row.rel_id) + # cs_nshm_2022: renumbering the two campaigns leaves every event R1..Rn, + # with the pilot's realisation last. + for event_id, group in group_by_event(rows).items(): + numbers = sorted(row.realisation for row in group) + if numbers != list(range(1, len(group) + 1)): + raise ValueError(f"event {event_id}: realisations {numbers}, expected 1..{len(group)}") + pilots = [row for row in group if row.pilot] + if len(pilots) > 1: + raise ValueError(f"event {event_id}: {len(pilots)} pilot realisations, expected at most 1") + if pilots and pilots[0].realisation != len(group): + raise ValueError( + f"event {event_id}: the pilot is R{pilots[0].realisation}, " + f"but it must be the last, R{len(group)}" + ) + + +def check_finite(rel: str, ims: RealisationIMs) -> None: + """Raise if any kept IM, site or distance value is NaN or infinite.""" + arrays = { + **{f"pSA/{c}": values for c, values in ims.psa.items()}, + **{f"FAS/{c}": values for c, values in ims.fas.items()}, + **{f"{im}/{c}": values for im, by_c in ims.scalars.items() for c, values in by_c.items()}, + **ims.site, + **ims.distances, + } + for name, values in arrays.items(): + bad = ~np.isfinite(values) + if bad.any(): + station = ims.stations[np.argwhere(bad)[0][0]] + raise ValueError( + f"{rel}: {name} has {int(bad.sum())} non-finite values, the first at {station}" + ) + + +def check_against_h5(rel: str, realisation_row: dict, attrs: dict[str, float]) -> None: + """Raise unless realisation.json reproduces the h5's magnitude and hypocentre.""" + for key, tolerance in TOLERANCES.items(): + ours, theirs = realisation_row[key], attrs[key] + if not abs(ours - theirs) <= tolerance: + raise ValueError( + f"{rel}: {key} from realisation.json is {ours!r}, but " + f"intensity_measures.h5 says {theirs!r}" + ) + + +def check_matches_first_realisation( + rel: str, first_rel: str, ims: RealisationIMs, first: RealisationIMs +) -> None: + """Raise unless stations, site fields and distances equal the event's first main realisation's. + + site_event keeps one set of distances per event, so an event's realisations + must agree on them. Nothing in the workflow enforces that; they agree only + because they share one fault database and one code version. + """ + if not np.array_equal(ims.stations, first.stations): + raise ValueError(f"{rel}: station list differs from {first_rel}") + for mine, theirs in ((ims.site, first.site), (ims.distances, first.distances)): + for key, values in mine.items(): + if not np.array_equal(values, theirs[key]): + raise ValueError( + f"{rel}: {key} differs from {first_rel}; site and distance fields " + "must be identical across an event's main realisations" + ) + + +def check_row_counts(db: IMDB, expected: dict[str, int]) -> None: + """Raise unless every table has exactly the expected number of rows.""" + for table, n in expected.items(): + found = int(db.con.table(table).count().to_pandas()) + if found != n: + raise RuntimeError(f"{table}: {found} rows, expected {n}") + + +# site_event must hold exactly the (site, event) pairs that have records: the +# pairs in site_event but not in records, and those in records but not in site_event. +SITE_EVENT_QUERY = """ +WITH pairs AS (SELECT DISTINCT site_int_id, event_int_id FROM records) +SELECT + (SELECT count(*) FROM site_event ANTI JOIN pairs USING (site_int_id, event_int_id)), + (SELECT count(*) FROM pairs ANTI JOIN site_event USING (site_int_id, event_int_id)) +""" + + +def check_site_event_matches_records(db: IMDB) -> None: + """Raise unless site_event has a row for exactly the (site, event) pairs that have records.""" + extra, missing = db.con.raw_sql(SITE_EVENT_QUERY).fetchone() + if extra or missing: + raise RuntimeError(f"site_event: {extra} rows with no records, {missing} record pairs with no row") -# Counts the flagged site_event rows three ways, from the tables alone: the rows -# flagged; the (site, event) pairs where only a pilot realisation has records in -# an event that also has main realisations; and the overlap of the two. + +# cs_nshm_2022: counts the flagged site_event rows three ways, from the tables +# alone: the rows flagged; the (site, event) pairs where only a pilot realisation +# has records in an event that also has main realisations; and the overlap of the two. FLAG_QUERY = """ WITH rel AS ( SELECT rel_int_id, event_int_id, @@ -855,23 +1150,6 @@ def check_row_counts(db: IMDB, expected: dict[str, int]) -> None: """ -# site_event must hold exactly the (site, event) pairs that have records: the -# pairs in site_event but not in records, and those in records but not in site_event. -SITE_EVENT_QUERY = """ -WITH pairs AS (SELECT DISTINCT site_int_id, event_int_id FROM records) -SELECT - (SELECT count(*) FROM site_event ANTI JOIN pairs USING (site_int_id, event_int_id)), - (SELECT count(*) FROM pairs ANTI JOIN site_event USING (site_int_id, event_int_id)) -""" - - -def check_site_event_matches_records(db: IMDB) -> None: - """Raise unless site_event has a row for exactly the (site, event) pairs that have records.""" - extra, missing = db.con.raw_sql(SITE_EVENT_QUERY).fetchone() - if extra or missing: - raise RuntimeError(f"site_event: {extra} rows with no records, {missing} record pairs with no row") - - def check_flags(db: IMDB, n_flagged: int) -> None: """Raise unless the flagged site_event rows are exactly the pilot-only sites of shared events.""" flagged, pilot_only, both = db.con.raw_sql(FLAG_QUERY).fetchone() @@ -882,6 +1160,20 @@ def check_flags(db: IMDB, n_flagged: int) -> None: ) +def verify_database(db: IMDB, counts: dict[str, int], n_flagged: int, expect_flagged: int | None) -> None: + """Raise unless the finished database passes every check that needs only its tables.""" + problems = db.validate() + if problems: + raise RuntimeError(f"validate(): {problems}") + check_row_counts(db, counts) + check_site_event_matches_records(db) + check_flags(db, n_flagged) + # cs_nshm_2022: the Slurm job knows how many gap sites the campaigns leave, + # counted from their station lists before the build. + if expect_flagged is not None and n_flagged != expect_flagged: + raise RuntimeError(f"{n_flagged} site_event rows carry the pilot's distances, expected {expect_flagged}") + + def spot_check( db_path: Path, manifest: dict[str, ManifestRow], @@ -894,7 +1186,8 @@ def spot_check( ) -> int: """Re-read random records from their h5 and compare exactly; return how many. - At least `min_pilot` of the sampled realisations are pilot ones, where there are any. + At least `min_pilot` of the sampled realisations are pilot ones, where there + are any (cs_nshm_2022: so both campaigns are always sampled). """ rng = np.random.default_rng(seed) pilot_rels = sorted(rel for rel, row in manifest.items() if row.pilot) @@ -944,6 +1237,9 @@ def spot_check( return checked +# cs_nshm_2022: some realisations reached physically implausible PGA near the +# fault, so the build log lists the largest for a human to judge before the +# database is shared. Another dataset can drop this, but it costs little to keep. def peak_line(row: ManifestRow, peak: dict[str, tuple[float, str]]) -> str: """One realisation's largest PGA and PGV (geom), with where they are.""" (pga, pga_site), (pgv, pgv_site) = peak["PGA"], peak["PGV"] @@ -958,147 +1254,7 @@ def print_peaks(peaks: Sequence[tuple[ManifestRow, dict[str, tuple[float, str]]] print(f" {peak_line(row, peak)}", flush=True) -def build_database( - manifest_path: Path, - stations_input: Path, - nshm_csv: Path, - out_path: Path, - harvested_at: str, - script_commit: str, - content: Content = DEFAULT_CONTENT, - memory_limit: str | None = None, - expect_flagged: int | None = None, -) -> dict[str, int]: - """Build the database at out_path and return its table row counts. - - Writes `.partial` and renames it only once validation, the row - counts, the site_event flags and the spot check have all passed. A failure - leaves the .partial behind for inspection, and the next run refuses to - start until it is removed. - """ - if out_path.exists(): - raise FileExistsError(f"{out_path} already exists") - partial = out_path.with_name(out_path.name + ".partial") - if partial.exists(): - raise FileExistsError(f"{partial} exists: a previous run died; remove it first") - - rows = read_manifest(manifest_path) - by_event = group_by_event(rows) - stations = load_stations_input(stations_input) - nshm = load_nshm_attributes(nshm_csv) - missing = sorted(set(by_event) - set(nshm.index), key=int) - if missing: - raise ValueError(f"ruptures missing from {nshm_csv}: {missing}") - - grids = read_grids(rows[0].im_path, content) - db = IMDB.create( - partial, - periods=list(grids.periods), - frequencies=list(grids.frequencies), - components=content.components, - db_meta={ - **DB_META, - **content_meta(content, grids), - "n_realisations": str(len(rows)), - "n_pilot_realisations": str(sum(row.pilot for row in rows)), - "im_calc_harvested_at": harvested_at, - "build_manifest_md5": md5_of(manifest_path), - "ingest_script_commit": script_commit, - }, - ) - counts = dict.fromkeys( - ("events", "realisations", "sites", "site_event", "records", "psa_ims", "scalars_ims", "fas_ims"), 0 - ) - n_flagged = 0 - seeds_seen: dict[str, str] = {} - peaks: list[tuple[ManifestRow, dict[str, tuple[float, str]]]] = [] - started = time.monotonic() - try: - if memory_limit: - db.con.raw_sql(f"SET memory_limit = '{memory_limit}'") - sites = SiteRegistry(stations) - for index, (event_id, group) in enumerate(by_event.items(), start=1): - loaded = [] - for row in group: - realisation, ims = load_realisation(row, content, grids) - seeds = json.dumps(realisation["seeds"], sort_keys=True) - if seeds in seeds_seen: - raise ValueError(f"{row.label}: the same seeds as {seeds_seen[seeds]}") - seeds_seen[seeds] = row.rel_id - check_finite(row.label, ims) - loaded.append(Loaded(row, realisation, ims)) - peaks.append((row, ims.peaks)) - print(f"peak {peak_line(row, ims.peaks)}", flush=True) - - mains = [x for x in loaded if not x.row.pilot] - pilots = [x for x in loaded if x.row.pilot] - defining = mains[0] if mains else pilots[0] - event = build_event( - event_id, - defining.realisation, - [x.realisation for x in loaded], - nshm.loc[event_id], - MAIN_RELEASE if mains else PILOT_RELEASE, - ) - realisation_rows = [] - for x in loaded: - realisation_row = build_realisation(x.row, x.realisation) - check_against_h5(x.row.label, realisation_row, x.ims.attrs) - realisation_rows.append(realisation_row) - for x in mains[1:]: - check_matches_first_realisation(x.row.label, mains[0].row.label, x.ims, mains[0].ims) - new_sites = pd.concat([sites.new_sites(x.row.label, x.ims) for x in loaded], ignore_index=True) - site_event, flagged = build_site_event( - event_id, mains[0].ims if mains else None, pilots[0].ims if pilots else None - ) - - db.add_events(pd.DataFrame([event])) - db.add_realisations(pd.DataFrame(realisation_rows)) - if len(new_sites): - db.add_sites(new_sites) - db.add_site_event(site_event) - for x in loaded: - for records, psa, fas in build_records(x.row.rel_id, x.ims, content): - db.add_records(records, pSA=psa, FAS=fas) - n = len(x.ims.stations) - counts["records"] += n * len(content.components) - counts["psa_ims"] += n * len(content.psa_components) - counts["scalars_ims"] += n * len(content.psa_components) - counts["fas_ims"] += n * len(content.fas_components) - - release_memtables(db) - counts["events"] += 1 - counts["realisations"] += len(loaded) - counts["sites"] += len(new_sites) - counts["site_event"] += len(site_event) - n_flagged += flagged - if index % 25 == 0 or index == len(by_event): - print( - f"{index}/{len(by_event)} events, {counts['records']} records, " - f"{time.monotonic() - started:.0f} s", - flush=True, - ) - - problems = db.validate() - if problems: - raise RuntimeError(f"validate(): {problems}") - check_row_counts(db, counts) - check_site_event_matches_records(db) - check_flags(db, n_flagged) - if expect_flagged is not None and n_flagged != expect_flagged: - raise RuntimeError(f"{n_flagged} site_event rows carry the pilot's distances, expected {expect_flagged}") - finally: - db.close() - - checked = spot_check(partial, {row.rel_id: row for row in rows}, content, grids) - print( - f"validate() clean; row counts match; {n_flagged} site_event rows carry the " - f"pilot's distances; spot check: {checked} records identical to their h5", - flush=True, - ) - print_peaks(peaks) - partial.rename(out_path) - return counts +# ---- 6. Command line ----------------------------------------------------------------- def main() -> None: @@ -1109,6 +1265,7 @@ def main() -> None: parser.add_argument("stations_input", type=Path, help="the campaign's stations_input.ll (lon lat name)") parser.add_argument("nshm_csv", type=Path, help="campaign_nshm_ruptures.csv from export_nshm_rupture_attributes.py") parser.add_argument("out", type=Path, help="output .duckdb; must not exist") + # cs_nshm_2022: provenance for db_meta, supplied by the Slurm job. parser.add_argument("--harvested-at", required=True, help="when the manifest was built from cylc, UTC ISO 8601") parser.add_argument("--script-commit", required=True, help="git commit of this script, for db_meta") parser.add_argument( @@ -1123,11 +1280,14 @@ def main() -> None: help="default: the 25 in PSA_PERIODS; all: every one of im-calc's 111", ) parser.add_argument("--fas", action="store_true", help="add FAS on im-calc's full frequency grid") + # cs_nshm_2022: the number of pilot-only sites in shared events, from the Slurm job. parser.add_argument( "--expect-flagged", type=int, help="fail unless exactly this many site_event rows carry the pilot's distances", ) + # DuckDB sizes its memory from the node's RAM, not from the Slurm job's + # allocation, and a job that outgrows its allocation is killed. parser.add_argument( "--memory-limit", help="DuckDB memory_limit, e.g. 12GB. Without it DuckDB sizes itself from " From 76f99428c3edc01eb6db0d386836ebccfb68ed66 Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Thu, 24 Sep 2026 10:54:16 +1200 Subject: [PATCH 7/8] Say how to run the cs_nshm_2022 ingest tests with uv The module docstring pointed at the cs_nshm_2022 container, which imdb's contributors do not have. The uv command it now gives passes all 56 tests. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01Fo9vmUtWnieoTVWuEaYVii --- scripts/tests/test_cs_nshm_2022_ingest.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/scripts/tests/test_cs_nshm_2022_ingest.py b/scripts/tests/test_cs_nshm_2022_ingest.py index 7c0c1cc..6e7ad93 100644 --- a/scripts/tests/test_cs_nshm_2022_ingest.py +++ b/scripts/tests/test_cs_nshm_2022_ingest.py @@ -1,8 +1,9 @@ """Tests for scripts/cs_nshm_2022_ingest.py, on a synthetic two-campaign dataset. imdb's own CI installs neither h5py nor source-modelling, so this module skips -there. Run it inside the imdb_tools image: - apptainer exec imdb_tools_.sif python3 -m pytest -q scripts/tests +there. To run it, add them from the repository root (fiona, which +source-modelling needs, has no Python 3.14 wheel yet): + uv run --python 3.13 --group test --with h5py --with source-modelling pytest scripts/tests """ import csv From 9141ad2e77744192937c5d52129c3d01f976716b Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Thu, 24 Sep 2026 11:17:58 +1200 Subject: [PATCH 8/8] Write the cs_nshm_2022 ingest's docstrings in numpydoc form Every function, method and class now has a numpydoc docstring with the sections imdb's numpydoc_validation checks ask for: Parameters, Returns and, where a function raises its own errors, Raises. Dataclass fields are documented in Attributes sections instead of attribute docstrings, and the module docstring's guide sits under Notes, since numpydoc rejects custom section headings. `numpydoc lint` now reports nothing (it reported 77 findings). Only docstrings changed: every definition parses to the same code with docstrings stripped, every comment is unchanged, and the 56 tests pass. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01Fo9vmUtWnieoTVWuEaYVii --- scripts/cs_nshm_2022_ingest.py | 907 ++++++++++++++++++++++++++++++--- 1 file changed, 849 insertions(+), 58 deletions(-) diff --git a/scripts/cs_nshm_2022_ingest.py b/scripts/cs_nshm_2022_ingest.py index a6ac5dd..3f74948 100644 --- a/scripts/cs_nshm_2022_ingest.py +++ b/scripts/cs_nshm_2022_ingest.py @@ -25,8 +25,10 @@ checks its own output before publishing it. Sections 2 to 4 are the ingest itself; most of section 5 is that checking. -Inputs ------- +Notes +----- +Inputs: + - A build manifest CSV listing every realisation to ingest, with its rel_id, pilot flag and file paths. cs_nshm_2022's scripts/build_imdb_manifest.py writes it (https://github.com/ucgmsim/cs_nshm_2022). @@ -38,8 +40,8 @@ - Each rupture's NSHM magnitude and annual rate, in a CSV exported from nshmdb by cs_nshm_2022's scripts/export_nshm_rupture_attributes.py. -What is specific to cs_nshm_2022 --------------------------------- +What is specific to cs_nshm_2022: + Code that exists only because of the points below is marked `# cs_nshm_2022:`, so `grep -n "cs_nshm_2022:" cs_nshm_2022_ingest.py` lists what another dataset can drop. @@ -63,8 +65,8 @@ - im-calc's pSA grid has 111 periods, with none between 15 s and 20 s. - nshmdb is not available on BSC, hence the NSHM CSV. -Layout ------- +Layout: + 1. What goes into the database: the content options. 2. Reading the inputs. 3. Building the rows of each table. @@ -148,15 +150,33 @@ @dataclass(frozen=True) class Content: - """Which components, pSA periods and spectra go into the database.""" + """Which components, pSA periods and spectra go into the database. + + Attributes + ---------- + components : tuple of str + Which components to include, e.g. `"geom"`, `"rotd50"`, `"eas"`. + all_periods : bool + All 111 of im-calc's pSA periods, instead of the 25 in PSA_PERIODS. + fas : bool + FAS, on im-calc's full frequency grid, for the components that have it. + """ components: tuple[str, ...] = DEFAULT_COMPONENTS all_periods: bool = False - """All 111 of im-calc's pSA periods, instead of the 25 in PSA_PERIODS.""" fas: bool = False - """FAS, on im-calc's full frequency grid, for the components that have it.""" def __post_init__(self) -> None: + """Raise unless the chosen components are non-empty, unique, known and consistent with fas. + + Raises + ------ + ValueError + If `components` is empty, lists a component twice, includes a + component that is not one of PSA_COMPONENTS or FAS_COMPONENTS, + includes `eas` without `fas`, or `fas` is set without any + component that has FAS. + """ if not self.components: raise ValueError("no components chosen") if len(set(self.components)) != len(self.components): @@ -171,18 +191,38 @@ def __post_init__(self) -> None: @property def psa_components(self) -> tuple[str, ...]: - """The chosen components that have pSA and the scalar IMs: all but eas.""" + """The chosen components that have pSA and the scalar IMs: all but eas. + + Returns + ------- + tuple of str + The chosen components other than `eas`. + """ return tuple(c for c in self.components if c in PSA_COMPONENTS) @property def fas_components(self) -> tuple[str, ...]: - """The chosen components that get FAS: none without --fas.""" + """The chosen components that get FAS: none without --fas. + + Returns + ------- + tuple of str + The chosen components with FAS, empty unless `fas` is set. + """ return tuple(c for c in self.components if self.fas and c in FAS_COMPONENTS) @dataclass(frozen=True) class Grids: - """The pSA periods and FAS frequencies the database is created with.""" + """The pSA periods and FAS frequencies the database is created with. + + Attributes + ---------- + periods : tuple of float + Response spectral periods, in seconds, that pSA records are indexed by. + frequencies : tuple of float + Frequencies, in Hz, that FAS records are indexed by; empty without FAS. + """ periods: tuple[float, ...] frequencies: tuple[float, ...] @@ -193,7 +233,20 @@ class Grids: def read_grids(h5_path: Path, content: Content) -> Grids: - """The database's grids, taken from one realisation's h5 where the content needs its full grid.""" + """The database's grids, taken from one realisation's h5 where the content needs its full grid. + + Parameters + ---------- + h5_path : Path + Path to one realisation's intensity_measures.h5. + content : Content + Which components, periods and spectra to include. + + Returns + ------- + Grids + The periods and frequencies the database should be created with. + """ with h5py.File(h5_path, "r") as h5: if not content.psa_components: periods: tuple[float, ...] = () @@ -206,7 +259,20 @@ def read_grids(h5_path: Path, content: Content) -> Grids: def content_meta(content: Content, grids: Grids) -> dict[str, str]: - """db_meta entries describing what the options put in.""" + """db_meta entries describing what the options put in. + + Parameters + ---------- + content : Content + Which components, periods and spectra were chosen. + grids : Grids + The database's pSA periods and FAS frequencies. + + Returns + ------- + dict of str to str + The `psa_period_selection` and `fas` entries for db_meta. + """ if not content.psa_components: periods = "not included" elif content.all_periods: @@ -235,7 +301,20 @@ def content_meta(content: Content, grids: Grids) -> dict[str, str]: def rel_id(event_id: str, n: int) -> str: - """The realisation's stable id.""" + """The realisation's stable id. + + Parameters + ---------- + event_id : str + The rupture's nshm_id. + n : int + The realisation number within the event, starting at 1. + + Returns + ------- + str + The event id and realisation number joined as `_R`. + """ return f"{event_id}_R{n}" @@ -245,26 +324,70 @@ def rel_id(event_id: str, n: int) -> str: # either layout. @dataclass(frozen=True) class ManifestRow: - """One realisation to ingest, as listed in the build manifest.""" + """One realisation to ingest, as listed in the build manifest. + + Attributes + ---------- + rel_id : str + The realisation's stable id, `rel_id(event_id, realisation)`. + event_id : str + The rupture's nshm_id, as an unpadded integer string. + realisation : int + The realisation number within the event, starting at 1. + pilot : bool + Whether this realisation is from the pilot campaign. + campaign : str + Which simulation campaign; for checks only, never stored. + original_id : str + Where the result came from; for messages only, never stored. + im_path : Path + Path to the realisation's intensity_measures.h5. + realisation_path : Path + Path to the realisation's realisation.json. + """ rel_id: str event_id: str realisation: int pilot: bool campaign: str - """Which simulation campaign; for checks only, never stored.""" original_id: str - """Where the result came from; for messages only, never stored.""" im_path: Path realisation_path: Path @property def label(self) -> str: + """The realisation's id and original id, for messages. + + Returns + ------- + str + The realisation formatted as ` ()`. + """ return f"{self.rel_id} ({self.original_id})" def read_manifest(path: Path) -> list[ManifestRow]: - """The manifest's rows in ingest order: events by integer id, then realisations.""" + """The manifest's rows in ingest order: events by integer id, then realisations. + + Parameters + ---------- + path : Path + Path to the build manifest CSV. + + Returns + ------- + list of ManifestRow + The manifest's rows, sorted by event id and realisation number. + + Raises + ------ + ValueError + If the CSV is missing a required column, or a row's `pilot` field is + not "true" or "false". + FileNotFoundError + If any listed `im_path` or `realisation_path` does not exist. + """ with path.open(newline="") as f: reader = csv.DictReader(f) missing = set(MANIFEST_FIELDS) - set(reader.fieldnames or ()) @@ -296,7 +419,18 @@ def read_manifest(path: Path) -> list[ManifestRow]: def group_by_event(rows: Sequence[ManifestRow]) -> dict[str, list[ManifestRow]]: - """Rows per event, keeping the input order.""" + """Rows per event, keeping the input order. + + Parameters + ---------- + rows : Sequence of ManifestRow + Manifest rows to group. + + Returns + ------- + dict of str to list of ManifestRow + Each event id mapped to its rows, in the order they appear in `rows`. + """ grouped: dict[str, list[ManifestRow]] = {} for row in rows: grouped.setdefault(row.event_id, []).append(row) @@ -307,7 +441,23 @@ def group_by_event(rows: Sequence[ManifestRow]) -> dict[str, list[ManifestRow]]: # grid (about 0.1 km) and wrote those coordinates into its h5, so they differ # between simulations. sites stores each station's canonical location instead. def load_stations_input(path: Path) -> pd.DataFrame: - """Canonical station coordinates (`lon lat name` per line), indexed by name.""" + """Canonical station coordinates (`lon lat name` per line), indexed by name. + + Parameters + ---------- + path : Path + Path to the campaign's stations_input.ll. + + Returns + ------- + pd.DataFrame + Its `lon` and `lat` columns, indexed by `site_id`. + + Raises + ------ + ValueError + If the file lists the same station name twice. + """ df = pd.read_csv( path, sep=r"\s+", header=None, names=["lon", "lat", "site_id"], dtype={"site_id": str} ) @@ -318,7 +468,18 @@ def load_stations_input(path: Path) -> pd.DataFrame: def is_real_station(site_id: str) -> bool: - """GeoNet station codes are 3-4 characters; the virtual grid's codes are 7.""" + """GeoNet station codes are 3-4 characters; the virtual grid's codes are 7. + + Parameters + ---------- + site_id : str + A station's canonical name. + + Returns + ------- + bool + Whether the name's length marks it as a real, not virtual-grid, station. + """ # cs_nshm_2022: nothing in the inputs marks a station as real, so the name # length decides. Check how your own station list names its stations. return len(site_id) != 7 @@ -331,7 +492,21 @@ def is_real_station(site_id: str) -> bool: def load_nshm_attributes(path: Path) -> pd.DataFrame: """Rows of export_nshm_rupture_attributes.py's CSV, indexed by rupture id. - Raises unless they all come from the release DB_META documents. + Parameters + ---------- + path : Path + Path to the exported NSHM attributes CSV. + + Returns + ------- + pd.DataFrame + One row per rupture, indexed by `rupture_id`. + + Raises + ------ + ValueError + If the CSV's `nshmdb_sha256` values are not all MAIN_NSHMDB_SHA256, + i.e. they do not all come from the release DB_META documents. """ nshm = pd.read_csv(path, dtype={"rupture_id": str}).set_index("rupture_id") if set(nshm["nshmdb_sha256"]) != {MAIN_NSHMDB_SHA256}: @@ -344,33 +519,74 @@ def load_nshm_attributes(path: Path) -> pd.DataFrame: @dataclass class RealisationIMs: - """One realisation's intensity_measures.h5, reduced to what this database keeps.""" + """One realisation's intensity_measures.h5, reduced to what this database keeps. + + Attributes + ---------- + stations : np.ndarray + Station names, str. + site : dict of str to np.ndarray + SITE_FIELDS per station. + distances : dict of str to np.ndarray + DISTANCE_FIELDS per station, km. + psa : dict of str to np.ndarray + Component -> (n_stations, n_periods) float32, for Content.psa_components. + fas : dict of str to np.ndarray + Component -> (n_stations, n_frequencies) float32, for Content.fas_components. + n_frequencies : int + Length of the database's FAS grid; 0 without FAS. + scalars : dict of str to dict of str to np.ndarray + Scalar IM -> component -> (n_stations,) float32; rotd-undefined IMs omitted. + attrs : dict of str to float + The file's magnitude and hypocentre attributes (TOLERANCES keys). + peaks : dict of str to tuple of (float, str) + PGA and PGV (geom) -> (largest value, its station), whatever the database keeps. + """ stations: np.ndarray - """Station names, str.""" site: dict[str, np.ndarray] - """SITE_FIELDS per station.""" distances: dict[str, np.ndarray] - """DISTANCE_FIELDS per station, km.""" psa: dict[str, np.ndarray] - """Component -> (n_stations, n_periods) float32, for Content.psa_components.""" fas: dict[str, np.ndarray] - """Component -> (n_stations, n_frequencies) float32, for Content.fas_components.""" n_frequencies: int - """Length of the database's FAS grid; 0 without FAS.""" scalars: dict[str, dict[str, np.ndarray]] - """Scalar IM -> component -> (n_stations,) float32; rotd-undefined IMs omitted.""" attrs: dict[str, float] - """The file's magnitude and hypocentre attributes (TOLERANCES keys).""" peaks: dict[str, tuple[float, str]] - """PGA and PGV (geom) -> (largest value, its station), whatever the database keeps.""" def _decode(values: np.ndarray) -> np.ndarray: + """Decode any bytes elements of an array to str, leaving others unchanged. + + Parameters + ---------- + values : np.ndarray + Array whose elements may be bytes, as h5py returns for HDF5 strings. + + Returns + ------- + np.ndarray + Object array with every bytes element decoded to str. + """ return np.array([v.decode() if isinstance(v, bytes) else v for v in values], dtype=object) def _peak(values: np.ndarray, stations: np.ndarray) -> tuple[float, str]: + """The largest value in `values` and the station it belongs to, ignoring NaN. + + Parameters + ---------- + values : np.ndarray + One value per station. + stations : np.ndarray + Station names in the same order as `values`, str or bytes. + + Returns + ------- + value : float + The largest value. + station : str + The station it belongs to, decoded to str. + """ i = int(np.nanargmax(values)) station = stations[i] return float(values[i]), station.decode() if isinstance(station, bytes) else str(station) @@ -381,7 +597,22 @@ def select_period_indices( ) -> np.ndarray: """Column of each wanted period in im-calc's period grid. - Raises ValueError unless every wanted period matches exactly one available one. + Parameters + ---------- + available : np.ndarray + The full period grid read from one realisation's h5. + wanted : Sequence of float, optional + Periods to find columns for (default PSA_PERIODS). + + Returns + ------- + np.ndarray + Column index in `available` of each period in `wanted`, in order. + + Raises + ------ + ValueError + If any wanted period does not match exactly one available one. """ indices = [] for period in wanted: @@ -398,7 +629,30 @@ def select_period_indices( def read_ims( h5_path: Path, content: Content = DEFAULT_CONTENT, grids: Grids = DEFAULT_GRIDS ) -> RealisationIMs: - """Read the parts of one intensity_measures.h5 this database keeps.""" + """Read the parts of one intensity_measures.h5 this database keeps. + + Parameters + ---------- + h5_path : Path + Path to the realisation's intensity_measures.h5. + content : Content, optional + Which components, periods and spectra to keep (default DEFAULT_CONTENT). + grids : Grids, optional + The database's pSA periods and FAS frequencies (default DEFAULT_GRIDS). + + Returns + ------- + RealisationIMs + The stations, site fields, distances, IMs and peaks this database keeps. + + Raises + ------ + ValueError + If the file's pSA period count is not `N_IM_CALC_PERIODS`, an IM or + FAS group lists stations in a different order from pSA, a wanted + component is missing where it is not rotd-undefined, or the file's + FAS frequency grid differs from the database's. + """ with h5py.File(h5_path, "r") as h5: psa_group = h5["pSA"] periods = psa_group["period"][:] @@ -441,7 +695,17 @@ def read_ims( @dataclass class Loaded: - """One realisation read from disk.""" + """One realisation read from disk. + + Attributes + ---------- + row : ManifestRow + The manifest row this realisation was loaded from. + realisation : dict + The parsed realisation.json. + ims : RealisationIMs + The parts of intensity_measures.h5 this database keeps. + """ row: ManifestRow realisation: dict @@ -451,7 +715,28 @@ class Loaded: def load_realisation(row: ManifestRow, content: Content, grids: Grids) -> tuple[dict, RealisationIMs]: """Read one realisation's realisation.json and h5, naming the realisation in any error. - Raises unless the realisation.json describes the manifest row's event. + Parameters + ---------- + row : ManifestRow + The manifest row to load. + content : Content + Which components, periods and spectra to keep. + grids : Grids + The database's pSA periods and FAS frequencies. + + Returns + ------- + realisation : dict + The parsed realisation.json. + ims : RealisationIMs + The parts of intensity_measures.h5 this database keeps. + + Raises + ------ + KeyError, OSError, ValueError + Any of these from reading either file, re-raised with the + realisation's label at the front of the message; ValueError also if + realisation.json's `metadata.name` is not `"Rupture "`. """ try: realisation = json.loads(row.realisation_path.read_text()) @@ -472,6 +757,19 @@ def load_realisation(row: ManifestRow, content: Content, grids: Grids) -> tuple[ def fault_geometries(realisation: dict) -> dict[str, Fault]: + """Each named source geometry in realisation.json, as a Fault. + + Parameters + ---------- + realisation : dict + The parsed realisation.json. + + Returns + ------- + dict of str to Fault + Each `sources.source_geometries` entry's name mapped to its Fault, + built from its corners' latitude, longitude and depth. + """ return { name: Fault.from_corners( np.array([[c["latitude"], c["longitude"], c["depth"]] for c in entry["corners"]]).reshape(-1, 4, 3) @@ -481,15 +779,52 @@ def fault_geometries(realisation: dict) -> dict[str, Fault]: def initial_fault_name(causality_tree: dict[str, str | None]) -> str: + """The name of the fault the rupture starts on: the one entry with no parent. + + Parameters + ---------- + causality_tree : dict of str to str or None + Each fault name mapped to its parent's name, or None for the fault + the rupture starts on. + + Returns + ------- + str + The initiating fault's name. + """ return next(name for name, parent in causality_tree.items() if parent is None) def source_wkt(geometries: dict[str, Fault]) -> str: + """The geometries' fault planes as a 2D MULTIPOLYGON WKT, for events.source_wkt. + + Parameters + ---------- + geometries : dict of str to Fault + Each source's fault geometry, as from `fault_geometries`. + + Returns + ------- + str + A MULTIPOLYGON WKT with one polygon per plane, corners as (lon, lat). + """ planes = [plane for f in geometries.values() for plane in f.planes] return MultiPolygon([Polygon(p.corners[:, [1, 0]]) for p in planes]).wkt def trace_wkt(geometries: dict[str, Fault]) -> str: + """The geometries' fault plane top edges as a MULTILINESTRING WKT, for events.trace_wkt. + + Parameters + ---------- + geometries : dict of str to Fault + Each source's fault geometry, as from `fault_geometries`. + + Returns + ------- + str + A MULTILINESTRING WKT with one line per plane's top edge, as (lon, lat). + """ planes = [plane for f in geometries.values() for plane in f.planes] return MultiLineString([LineString(p.corners[:2, [1, 0]]) for p in planes]).wkt @@ -500,6 +835,16 @@ def fault_geometry_wkt(realisation: dict) -> str: A MULTIPOLYGON Z with one polygon per plane, corners as (lon, lat, depth in km). Unlike events.source_wkt it carries depth, so distances can be recomputed from it alone. + + Parameters + ---------- + realisation : dict + The parsed realisation.json. + + Returns + ------- + str + A MULTIPOLYGON Z WKT built directly from `sources.source_geometries`. """ planes = [] for entry in realisation["sources"]["source_geometries"].values(): @@ -509,17 +854,51 @@ def fault_geometry_wkt(realisation: dict) -> str: def domain_polygon(realisation: dict) -> Polygon: + """The realisation's simulation domain, from realisation.json. + + Parameters + ---------- + realisation : dict + The parsed realisation.json. + + Returns + ------- + Polygon + The domain, with its corners as (lon, lat). + """ corners = realisation["domain"]["domain"] return Polygon([(c["longitude"], c["latitude"]) for c in corners]) def domain_wkt(realisations: Sequence[dict]) -> str: - """The union of the realisations' simulation domains, so it contains every site with data.""" + """The union of the realisations' simulation domains, so it contains every site with data. + + Parameters + ---------- + realisations : Sequence of dict + Parsed realisation.json documents to union. + + Returns + ------- + str + The union of `domain_polygon` over `realisations`, as WKT. + """ return unary_union([domain_polygon(r) for r in realisations]).wkt def total_magnitude(realisation: dict) -> float: - """The realisation's total magnitude, BoldM: segment moments summed, converted back.""" + """The realisation's total magnitude, BoldM: segment moments summed, converted back. + + Parameters + ---------- + realisation : dict + The parsed realisation.json. + + Returns + ------- + float + The BoldM magnitude equivalent to the sum of the segments' moments. + """ magnitudes = realisation["magnitudes"]["magnitudes"] total_moment = sum(moment.magnitude_to_moment(m, bold_m=True) for m in magnitudes.values()) return float(moment.moment_to_magnitude(total_moment, bold_m=True)) @@ -539,6 +918,29 @@ def build_event( the domain its own. Here the pilot's geometry and domain differ from the main campaign's, so the caller picks `defining`, and the union keeps every site with data inside the event's domain. + + Parameters + ---------- + event_id : str + The rupture's nshm_id. + defining : dict + The parsed realisation.json whose geometry and tect_type the event uses. + realisations : Sequence of dict + Every realisation of this event, parsed, for the domain union. + nshm : pd.Series + The rupture's row from `load_nshm_attributes`. + release : str + The NSHM fault database release `defining`'s geometry came from. + + Returns + ------- + dict + One row for the `events` table. + + Raises + ------ + ValueError + If `defining`'s `empirical.tect_type` is not one of `schema.TECT_TYPES`. """ geometries = fault_geometries(defining) rakes = defining["rakes"]["rakes"] @@ -592,7 +994,20 @@ def build_event( def build_realisation(row: ManifestRow, realisation: dict) -> dict: - """The realisations row, from the realisation's own realisation.json and geometry.""" + """The realisations row, from the realisation's own realisation.json and geometry. + + Parameters + ---------- + row : ManifestRow + The manifest row this realisation was loaded from. + realisation : dict + The parsed realisation.json. + + Returns + ------- + dict + One row for the `realisations` table. + """ geometries = fault_geometries(realisation) causality_tree = realisation["rupture_propagation"]["rupture_causality_tree"] initial_fault = initial_fault_name(causality_tree) @@ -625,14 +1040,47 @@ def build_realisation(row: ManifestRow, realisation: dict) -> dict: class SiteRegistry: - """The sites added so far, with the vs30/z1pt0/z2pt5 first seen for each.""" + """The sites added so far, with the vs30/z1pt0/z2pt5 first seen for each. + + Parameters + ---------- + stations_input : pd.DataFrame + Canonical station coordinates, as from `load_stations_input`. + """ def __init__(self, stations_input: pd.DataFrame) -> None: + """Start with no sites seen yet. + + Parameters + ---------- + stations_input : pd.DataFrame + Canonical station coordinates, as from `load_stations_input`. + """ self.stations_input = stations_input self.seen: pd.DataFrame | None = None def new_sites(self, label: str, ims: RealisationIMs) -> pd.DataFrame: - """Rows for stations not added yet; raise if a known one disagrees.""" + """Rows for stations not added yet; raise if a known one disagrees. + + Parameters + ---------- + label : str + The realisation's label, for the error message. + ims : RealisationIMs + The realisation's IMs, whose stations and site fields are checked + against and added to the registry. + + Returns + ------- + pd.DataFrame + One `sites` row per station in `ims` not already registered. + + Raises + ------ + ValueError + If a station already seen has different site fields this time, or + a new station is missing from `stations_input`. + """ here = pd.DataFrame(ims.site, index=pd.Index(ims.stations, name="site_id")) if self.seen is None: new = here @@ -672,6 +1120,24 @@ def new_sites(self, label: str, ims: RealisationIMs) -> pd.DataFrame: def _site_event_rows( event_id: str, stations: np.ndarray, distances: dict[str, np.ndarray], metadata: str | None ) -> pd.DataFrame: + """Assemble site_event rows from parallel station and distance arrays. + + Parameters + ---------- + event_id : str + The event these rows belong to. + stations : np.ndarray + Station names, one per row. + distances : dict of str to np.ndarray + DISTANCE_FIELDS arrays, one value per row, in the same order as `stations`. + metadata : str or None + `site_event.metadata` value for every row. + + Returns + ------- + pd.DataFrame + One `site_event` row per station. + """ return pd.DataFrame({"site_id": stations, "event_id": event_id, **distances, "metadata": metadata}) @@ -687,6 +1153,27 @@ def build_site_event( pilot's alone. cs_nshm_2022: with one campaign, this would be one realisation's distances. + + Parameters + ---------- + event_id : str + The event these rows belong to. + main : RealisationIMs or None + The event's first main-campaign realisation's IMs, or None if it has none. + pilot : RealisationIMs or None + The event's pilot realisation's IMs, or None if it has none. + + Returns + ------- + rows : pd.DataFrame + One `site_event` row per site the event reaches. + n_flagged : int + How many rows carry the pilot's distances. + + Raises + ------ + ValueError + If both `main` and `pilot` are None. """ if main is None: if pilot is None: @@ -714,6 +1201,21 @@ def build_records(rel: str, ims: RealisationIMs, content: Content = DEFAULT_CONT `eas` has only FAS, so its records form a batch of their own with no pSA and no scalar columns, and add_records writes no psa_ims or scalars_ims rows for them. A rotd record's FAS row is all NaN, so it gets no fas_ims row. + + Parameters + ---------- + rel : str + The realisation's rel_id. + ims : RealisationIMs + The realisation's IMs to turn into records. + content : Content, optional + Which components, periods and spectra to include (default DEFAULT_CONTENT). + + Returns + ------- + list of Batch + One (records, pSA, FAS) batch per `IMDB.add_records` call needed: the + pSA components' batch, if any, then eas's, if included. """ n = len(ims.stations) frames = [] @@ -804,7 +1306,18 @@ def build_records(rel: str, ims: RealisationIMs, content: Content = DEFAULT_CONT def md5_of(path: Path) -> str: - """Hex md5 of a file.""" + """Hex md5 of a file. + + Parameters + ---------- + path : Path + File to hash. + + Returns + ------- + str + The file's contents' hex-encoded md5 digest. + """ return hashlib.md5(path.read_bytes()).hexdigest() @@ -815,6 +1328,11 @@ def release_memtables(db: IMDB) -> None: and releases them only at interpreter exit, so without this every row written stays in memory, outside DuckDB's memory_limit, for the whole build. Any large ingest needs it until imdb releases them itself. + + Parameters + ---------- + db : IMDB + The database to release memtables from. """ raw = db.con.con for (name,) in raw.execute( @@ -828,8 +1346,28 @@ def read_event( ) -> list[Loaded]: """Read one event's realisations, checking each on its own. - `seeds_seen` maps each set of seeds read so far to its rel_id, which - catches the same simulation listed twice under different ids. + Parameters + ---------- + group : Sequence of ManifestRow + The event's manifest rows. + content : Content + Which components, periods and spectra to keep. + grids : Grids + The database's pSA periods and FAS frequencies. + seeds_seen : dict of str to str + Every realisation's seeds read so far, mapped to its rel_id, which + catches the same simulation listed twice under different ids. + Updated in place with this event's realisations. + + Returns + ------- + list of Loaded + The event's realisations, read and checked. + + Raises + ------ + ValueError + If a realisation's seeds match one already in `seeds_seen`. """ loaded = [] for row in group: @@ -848,19 +1386,49 @@ def read_event( @dataclass class EventRows: - """Everything one event adds to the database apart from its records.""" + """Everything one event adds to the database apart from its records. + + Attributes + ---------- + event : dict + The event's `events` row. + realisations : list of dict + The event's `realisations` rows. + sites : pd.DataFrame + Sites no earlier event reached. + site_event : pd.DataFrame + The event's `site_event` rows. + n_flagged : int + How many site_event rows carry the pilot's distances. + """ event: dict realisations: list[dict] sites: pd.DataFrame - """Sites no earlier event reached.""" site_event: pd.DataFrame n_flagged: int - """How many site_event rows carry the pilot's distances.""" def build_event_rows(event_id: str, loaded: Sequence[Loaded], nshm: pd.Series, sites: SiteRegistry) -> EventRows: - """The events, realisations, sites and site_event rows for one event.""" + """The events, realisations, sites and site_event rows for one event. + + Parameters + ---------- + event_id : str + The rupture's nshm_id. + loaded : Sequence of Loaded + The event's realisations, read and checked. + nshm : pd.Series + The rupture's row from `load_nshm_attributes`. + sites : SiteRegistry + The sites added by earlier events; updated in place with this event's + new ones. + + Returns + ------- + EventRows + The event's rows for every table but records. + """ # cs_nshm_2022: the event's geometry and unflagged distances come from its # first main realisation, or from the pilot where the main campaign never # simulated the rupture. @@ -889,7 +1457,24 @@ def build_event_rows(event_id: str, loaded: Sequence[Loaded], nshm: pd.Series, s def write_event(db: IMDB, rows: EventRows, loaded: Sequence[Loaded], content: Content) -> dict[str, int]: - """Write one event, and return how many rows each table should have gained.""" + """Write one event, and return how many rows each table should have gained. + + Parameters + ---------- + db : IMDB + The database to write to, open for writing. + rows : EventRows + The event's rows for every table but records. + loaded : Sequence of Loaded + The event's realisations, for building and writing their records. + content : Content + Which components, periods and spectra to write. + + Returns + ------- + dict of str to int + How many rows each table should have gained. + """ db.add_events(pd.DataFrame([rows.event])) db.add_realisations(pd.DataFrame(rows.realisations)) if len(rows.sites): @@ -930,6 +1515,42 @@ def build_database( counts, the site_event flags and the spot check have all passed. A failure leaves the .partial behind for inspection, and the next run refuses to start until it is removed. + + Parameters + ---------- + manifest_path : Path + Path to the build manifest CSV. + stations_input : Path + Path to the campaign's stations_input.ll. + nshm_csv : Path + Path to export_nshm_rupture_attributes.py's CSV. + out_path : Path + Path to write the finished database to; must not exist. + harvested_at : str + When the manifest was built from cylc, UTC ISO 8601, stored in db_meta. + script_commit : str + This script's git commit, stored in db_meta. + content : Content, optional + Which components, periods and spectra to include (default DEFAULT_CONTENT). + memory_limit : str, optional + DuckDB `memory_limit`, e.g. `"12GB"`. Without it DuckDB sizes itself + from the node's RAM, not the job's memory allocation. + expect_flagged : int, optional + If given, fail unless exactly this many site_event rows carry the + pilot's distances. + + Returns + ------- + dict of str to int + Each table's row count in the finished database. + + Raises + ------ + FileExistsError + If `out_path` already exists, or a `.partial` file from a previous, + unfinished run is still there. + ValueError + If the manifest lists a rupture not found in `nshm_csv`. """ if out_path.exists(): raise FileExistsError(f"{out_path} already exists") @@ -1012,6 +1633,21 @@ def check_manifest(rows: Sequence[ManifestRow]) -> None: The pilot flag must also follow the campaign: no campaign may supply both pilot and non-pilot rows. + + Parameters + ---------- + rows : Sequence of ManifestRow + The manifest rows to check. + + Raises + ------ + ValueError + If `rows` is empty, a campaign supplies both pilot and non-pilot + rows, an `event_id` is not an unpadded integer, a `rel_id` does not + match `rel_id(event_id, realisation)`, a `rel_id` is listed twice, an + event's realisation numbers are not `1..n` with no gaps, an event has + more than one pilot realisation, or an event's pilot realisation is + not last. """ if not rows: raise ValueError("the manifest lists no realisations") @@ -1051,7 +1687,20 @@ def check_manifest(rows: Sequence[ManifestRow]) -> None: def check_finite(rel: str, ims: RealisationIMs) -> None: - """Raise if any kept IM, site or distance value is NaN or infinite.""" + """Raise if any kept IM, site or distance value is NaN or infinite. + + Parameters + ---------- + rel : str + The realisation's label, for the error message. + ims : RealisationIMs + The realisation's IMs to check. + + Raises + ------ + ValueError + If any pSA, FAS, scalar IM, site or distance array has a non-finite value. + """ arrays = { **{f"pSA/{c}": values for c, values in ims.psa.items()}, **{f"FAS/{c}": values for c, values in ims.fas.items()}, @@ -1069,7 +1718,23 @@ def check_finite(rel: str, ims: RealisationIMs) -> None: def check_against_h5(rel: str, realisation_row: dict, attrs: dict[str, float]) -> None: - """Raise unless realisation.json reproduces the h5's magnitude and hypocentre.""" + """Raise unless realisation.json reproduces the h5's magnitude and hypocentre. + + Parameters + ---------- + rel : str + The realisation's label, for the error message. + realisation_row : dict + The realisation's `realisations` row, built from realisation.json. + attrs : dict of str to float + The h5's own magnitude and hypocentre attributes (TOLERANCES keys). + + Raises + ------ + ValueError + If any TOLERANCES field differs between `realisation_row` and `attrs` + by more than its tolerance. + """ for key, tolerance in TOLERANCES.items(): ours, theirs = realisation_row[key], attrs[key] if not abs(ours - theirs) <= tolerance: @@ -1087,6 +1752,23 @@ def check_matches_first_realisation( site_event keeps one set of distances per event, so an event's realisations must agree on them. Nothing in the workflow enforces that; they agree only because they share one fault database and one code version. + + Parameters + ---------- + rel : str + This realisation's label, for the error message. + first_rel : str + The event's first main realisation's label, for the error message. + ims : RealisationIMs + This realisation's IMs. + first : RealisationIMs + The event's first main realisation's IMs, to compare against. + + Raises + ------ + ValueError + If the station list differs from `first`'s, or any site or distance + field differs from `first`'s. """ if not np.array_equal(ims.stations, first.stations): raise ValueError(f"{rel}: station list differs from {first_rel}") @@ -1100,7 +1782,20 @@ def check_matches_first_realisation( def check_row_counts(db: IMDB, expected: dict[str, int]) -> None: - """Raise unless every table has exactly the expected number of rows.""" + """Raise unless every table has exactly the expected number of rows. + + Parameters + ---------- + db : IMDB + The database to check, open. + expected : dict of str to int + Each table's expected row count. + + Raises + ------ + RuntimeError + If a table's row count differs from `expected`. + """ for table, n in expected.items(): found = int(db.con.table(table).count().to_pandas()) if found != n: @@ -1118,7 +1813,19 @@ def check_row_counts(db: IMDB, expected: dict[str, int]) -> None: def check_site_event_matches_records(db: IMDB) -> None: - """Raise unless site_event has a row for exactly the (site, event) pairs that have records.""" + """Raise unless site_event has a row for exactly the (site, event) pairs that have records. + + Parameters + ---------- + db : IMDB + The database to check, open. + + Raises + ------ + RuntimeError + If any site_event row has no matching records, or any (site, event) + pair with records has no site_event row. + """ extra, missing = db.con.raw_sql(SITE_EVENT_QUERY).fetchone() if extra or missing: raise RuntimeError(f"site_event: {extra} rows with no records, {missing} record pairs with no row") @@ -1151,7 +1858,21 @@ def check_site_event_matches_records(db: IMDB) -> None: def check_flags(db: IMDB, n_flagged: int) -> None: - """Raise unless the flagged site_event rows are exactly the pilot-only sites of shared events.""" + """Raise unless the flagged site_event rows are exactly the pilot-only sites of shared events. + + Parameters + ---------- + db : IMDB + The database to check, open. + n_flagged : int + How many site_event rows the build wrote with the pilot's distances. + + Raises + ------ + RuntimeError + If the flagged row count, the pilot-only-site count and `n_flagged` + are not all equal. + """ flagged, pilot_only, both = db.con.raw_sql(FLAG_QUERY).fetchone() if not flagged == pilot_only == both == n_flagged: raise RuntimeError( @@ -1161,7 +1882,25 @@ def check_flags(db: IMDB, n_flagged: int) -> None: def verify_database(db: IMDB, counts: dict[str, int], n_flagged: int, expect_flagged: int | None) -> None: - """Raise unless the finished database passes every check that needs only its tables.""" + """Raise unless the finished database passes every check that needs only its tables. + + Parameters + ---------- + db : IMDB + The database to check, open. + counts : dict of str to int + Each table's expected row count. + n_flagged : int + How many site_event rows the build wrote with the pilot's distances. + expect_flagged : int or None + If given, the exact `n_flagged` the caller expects. + + Raises + ------ + RuntimeError + If `db.validate()` reports any problems, or `n_flagged` does not + equal `expect_flagged` when one is given. + """ problems = db.validate() if problems: raise RuntimeError(f"validate(): {problems}") @@ -1188,6 +1927,38 @@ def spot_check( At least `min_pilot` of the sampled realisations are pilot ones, where there are any (cs_nshm_2022: so both campaigns are always sampled). + + Parameters + ---------- + db_path : Path + Path to the finished (or `.partial`) database. + manifest : dict of str to ManifestRow + Every ingested realisation's manifest row, by rel_id. + content : Content, optional + Which components, periods and spectra the database holds (default DEFAULT_CONTENT). + grids : Grids, optional + The database's pSA periods and FAS frequencies (default DEFAULT_GRIDS). + n_realisations : int, optional + How many realisations to sample (default 25). + per_realisation : int, optional + How many records to sample per realisation (default 40). + min_pilot : int, optional + The minimum number of sampled realisations that must be pilot ones, + where any exist (default 5). + seed : int, optional + Random seed for reproducible sampling (default 0). + + Returns + ------- + int + How many records were checked. + + Raises + ------ + AssertionError + If a sampled record's pSA, scalar IM or FAS values differ from the + h5, are present for a component that should have none, or are not + NULL where the h5 has none. """ rng = np.random.default_rng(seed) pilot_rels = sorted(rel for rel, row in manifest.items() if row.pilot) @@ -1241,13 +2012,32 @@ def spot_check( # fault, so the build log lists the largest for a human to judge before the # database is shared. Another dataset can drop this, but it costs little to keep. def peak_line(row: ManifestRow, peak: dict[str, tuple[float, str]]) -> str: - """One realisation's largest PGA and PGV (geom), with where they are.""" + """One realisation's largest PGA and PGV (geom), with where they are. + + Parameters + ---------- + row : ManifestRow + The realisation, for its label. + peak : dict of str to tuple of (float, str) + `PGA` and `PGV` mapped to (largest value, its station). + + Returns + ------- + str + A one-line summary naming the realisation, its peak PGA and its peak PGV. + """ (pga, pga_site), (pgv, pgv_site) = peak["PGA"], peak["PGV"] return f"{row.label}: PGA {pga:.3f} g at {pga_site}; PGV {pgv:.1f} cm/s at {pgv_site}" def print_peaks(peaks: Sequence[tuple[ManifestRow, dict[str, tuple[float, str]]]]) -> None: - """The realisations with the largest PGA, for a human to judge before the database is shared.""" + """The realisations with the largest PGA, for a human to judge before the database is shared. + + Parameters + ---------- + peaks : Sequence of (ManifestRow, dict of str to tuple of (float, str)) + Each realisation with its `PGA`/`PGV` peaks, as from `read_event`. + """ ranked = sorted(peaks, key=lambda item: item[1]["PGA"][0], reverse=True) print(f"top {min(N_TOP_PEAKS, len(ranked))} realisations by max PGA (geom):", flush=True) for row, peak in ranked[:N_TOP_PEAKS]: @@ -1258,6 +2048,7 @@ def print_peaks(peaks: Sequence[tuple[ManifestRow, dict[str, tuple[float, str]]] def main() -> None: + """Parse the command line, build the database, and print its row counts.""" parser = argparse.ArgumentParser( description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter )