Read SW4 station recordings, carrying supergrid penetration through

tests/test_im_calc.py

@@ -0,0 +1,179 @@
+"""Pure-function tests for the intensity measure calculation.
+
+Nothing here runs OpenQuake or an IM kernel: what is tested is the metadata
+plumbing that carries the SW4 supergrid (absorbing layer) penetration from the
+waveform file into the intensity measure file, which is arithmetic-free and
+where the silent failures live.
+"""
+
+import numpy as np
+import pytest
+import xarray as xr
+
+from workflow.scripts import im_calc
+
+
+def _waveform_dataset(
+    n_stations: int = 3,
+    supergrid_depth: list[float] | None = None,
+    attrs: dict[str, float] | None = None,
+) -> xr.Dataset:
+    """A minimal stand-in for an opened broadband/LF waveform file."""
+    stations = [f"ST{index:02d}" for index in range(n_stations)]
+    coords: dict[str, object] = {"station": stations}
+    if supergrid_depth is not None:
+        coords["supergrid_depth"] = (
+            "station",
+            np.array(supergrid_depth, dtype=np.float32),
+        )
+        coords["supergrid_depth_gp"] = (
+            "station",
+            np.array(supergrid_depth, dtype=np.float32) / 400.0,
+        )
+    return xr.Dataset(
+        {"waveform": (("station",), np.ones(n_stations, dtype=np.float32))},
+        coords=coords,
+        attrs=attrs or {},
+    )
+
+
+def test_a_solver_with_no_absorbing_layer_still_gets_the_coordinate() -> None:
+    """Every IM file carries the coordinate, whatever the solver produced it.
+
+    EMOD3D and the high-frequency simulation have no supergrid at all, so the
+    value has to be NaN -- "not applicable / not reported" -- and never `0.0`,
+    which would assert the station had been checked and found in the interior.
+    """
+    supergrid = im_calc.supergrid_coordinates(_waveform_dataset(n_stations=4))
+
+    assert set(supergrid) == set(im_calc.SUPERGRID_COORDINATES)
+    for values in supergrid.values():
+        assert values.dims == ("station",)
+        assert values.dtype == np.float32
+        assert np.isnan(values.values).all()
+
+
+def test_a_reported_penetration_is_passed_through_unchanged() -> None:
+    """The three states survive verbatim: clean, flagged, and unknown."""
+    dataset = _waveform_dataset(supergrid_depth=[0.0, 5750.0, np.nan])
+
+    supergrid = im_calc.supergrid_coordinates(dataset)
+
+    depth = supergrid["supergrid_depth"]
+    assert depth.values[0] == 0.0
+    assert depth.values[1] == pytest.approx(5750.0)
+    assert np.isnan(depth.values[2])
+    # `> 0` is the documented threshold; NaN must not satisfy it.
+    np.testing.assert_array_equal(depth.values > 0, [False, True, False])
+
+
+def test_the_coordinate_is_read_eagerly() -> None:
+    """`im-calc` opens the waveform file chunked, so the coordinate arrives as
+    a dask array whose "auto" station chunking differs from the waveform's.
+    Loading it here keeps that mismatch out of the attached coordinates, the
+    way `vs30` is loaded eagerly for the same reason.
+    """
+    dataset = _waveform_dataset(supergrid_depth=[0.0, 1.0, 2.0]).chunk({"station": 1})
+
+    supergrid = im_calc.supergrid_coordinates(dataset)
+
+    assert supergrid["supergrid_depth"].chunks is None
+
+
+def test_nothing_is_claimed_about_a_run_that_reported_nothing() -> None:
+    """An all-NaN flag must not put `absorbing_layer` in the root attributes.
+
+    Writing it would claim the run had an absorbing layer that somebody
+    measured, on the strength of a coordinate that says only "unknown".
+    """
+    dataset = _waveform_dataset(attrs={"SGWIDTH": 12000.0})
+    supergrid = im_calc.supergrid_coordinates(dataset)
+
+    assert im_calc.supergrid_attributes(dataset, supergrid) == {}
+
+
+def test_the_sponge_width_comes_from_the_waveform_file() -> None:
+    """Not from the realisation configuration.
+
+    The configuration is editable after a run; the waveform file is what SW4
+    actually wrote. Reading the config here would let the IM file's
+    self-description drift away from the run it describes.
+    """
+    dataset = _waveform_dataset(
+        supergrid_depth=[0.0, 5750.0, 0.0],
+        attrs={"SGWIDTH": 12000.0, "SGWIDTHGP": 30.0},
+    )
+    supergrid = im_calc.supergrid_coordinates(dataset)
+
+    attributes = im_calc.supergrid_attributes(dataset, supergrid)
+
+    assert attributes["absorbing_layer"] == "sw4_supergrid"
+    assert attributes["absorbing_layer_width_m"] == pytest.approx(12000.0)
+    assert attributes["absorbing_layer_width_gp"] == pytest.approx(30.0)
+
+
+def test_a_flag_without_a_width_still_names_the_layer() -> None:
+    """A station file written with penetrations but no file-level width (or a
+    broadband file whose width did not survive) must not lose the layer name.
+    """
+    dataset = _waveform_dataset(supergrid_depth=[0.0, 5750.0, 0.0])
+    supergrid = im_calc.supergrid_coordinates(dataset)
+
+    attributes = im_calc.supergrid_attributes(dataset, supergrid)
+
+    assert attributes == {"absorbing_layer": "sw4_supergrid"}
+
+
+def test_the_flag_is_attached_as_a_coordinate_on_every_leaf() -> None:
+    """A data variable here would die at `bb_sim`'s `combined` dataset on the
+    next run through the pipeline, and cannot be selected alongside an IM.
+    The root deliberately has no station dimension, so it stays untouched.
+    """
+    dataset = _waveform_dataset(supergrid_depth=[0.0, 5750.0, np.nan])
+    supergrid = im_calc.supergrid_coordinates(dataset)
+    dtree = xr.DataTree.from_dict(
+        {
+            "PGA": xr.Dataset(
+                {"rotd50": (("station",), np.ones(3))},
+                coords={"station": dataset.station},
+            ),
+            "pSA": xr.Dataset(
+                {"rotd50": (("station",), np.ones(3))},
+                coords={"station": dataset.station},
+            ),
+        }
+    )
+
+    parameterised = im_calc.add_station_parameters(dtree, supergrid)
+
+    for group in ("PGA", "pSA"):
+        leaf = parameterised[group].dataset
+        for name in im_calc.SUPERGRID_COORDINATES:
+            assert name in leaf.coords
+            assert name not in leaf.data_vars
+    assert "supergrid_depth" not in parameterised.dataset.coords
+
+
+def test_the_recommended_threshold_travels_with_the_data() -> None:
+    """The coordinate's `description` is the only documentation a downstream
+    user ever sees, so it has to carry the threshold and say what the trace
+    is, or every tool invents its own cut.
+    """
+    for name in im_calc.SUPERGRID_COORDINATES:
+        assert name in im_calc.COORDINATE_METADATA
+        description = im_calc.COORDINATE_METADATA[name]["description"]
+        assert "> 0" in description
+        assert "NaN" in description
+
+    dataset = xr.Dataset(
+        {"rotd50": (("station",), np.ones(2))},
+        coords={
+            "station": ["ST00", "ST01"],
+            "supergrid_depth": ("station", np.array([0.0, 1.0], dtype=np.float32)),
+        },
+    )
+    annotated = im_calc.add_units(xr.DataTree.from_dict({"PGA": dataset}))
+
+    attrs = annotated["PGA"].dataset["supergrid_depth"].attrs
+    assert attrs["units"] == "m"
+    assert "> 0" in attrs["description"]

tests/test_lf_to_xarray.py

@@ -0,0 +1,212 @@
+"""Tests for reading SW4 station recordings, and in particular for the
+supergrid (absorbing layer) penetration SW4 reports per station.
+
+The station file fixture below is the first synthetic SW4 recording in the
+suite; it is deliberately written against the layout documented in Section
+12.9 of the SW4 User Guide (a root `DELTA`, one group per station holding
+`NPTS`, `STLA,STLO,STDP` and the three geographic components) so that it is
+reusable for any other SW4-read test.
+"""
+
+from pathlib import Path
+
+import h5py
+import numpy as np
+import pytest
+import xarray as xr
+
+from workflow.scripts import lf_to_xarray
+
+
+def write_sw4_station_file(
+    path: Path,
+    stations: dict[str, dict[str, float] | None],
+    npts: int = 8,
+    dt: float = 0.05,
+    widths: dict[str, float] | None = None,
+) -> Path:
+    """Write a synthetic SW4 HDF5 station recording.
+
+    Parameters
+    ----------
+    path : Path
+        Where to write the file.
+    stations : dict
+        Map from station name to either `None` (an old-style station, with no
+        supergrid datasets at all) or a dict which may hold `SGDEPTH` and
+        `SGDEPTHGP`. A dict holding only one of the two produces a
+        deliberately corrupt station.
+    npts : int
+        Number of samples per component.
+    dt : float
+        Sample spacing, written to the root `DELTA`.
+    widths : dict, optional
+        File-level scalars (`SGWIDTH`, `SGWIDTHGP`) written beside `DELTA`.
+
+    Returns
+    -------
+    Path
+        The path written, for convenience.
+    """
+    with h5py.File(path, "w") as handle:
+        handle.create_dataset("DELTA", data=np.array([dt]))
+        for name, value in (widths or {}).items():
+            handle.create_dataset(name, data=np.array([value]))
+        for index, (station, supergrid) in enumerate(stations.items()):
+            group = handle.create_group(station)
+            group.create_dataset("NPTS", data=np.array([npts]))
+            group.create_dataset(
+                "STLA,STLO,STDP",
+                data=np.array([-43.5 + index, 172.6 + index, 0.0]),
+            )
+            for component in ("EW", "NS", "UP"):
+                group.create_dataset(
+                    component, data=np.arange(npts, dtype=np.float32) + index
+                )
+            for key, value in (supergrid or {}).items():
+                group.create_dataset(key, data=np.array([value]))
+    return path
+
+
+def test_supergrid_penetration_arrives_as_float32_coordinates(tmp_path: Path) -> None:
+    """The flag must be a *coordinate*, and it must be floating point.
+
+    Both halves are load bearing. A station-dimension coordinate rides
+    through `bb-sim` and `im-calc` untouched, whereas a data variable is
+    silently dropped by `bb_sim._process_bb_chunk`, so a data variable here
+    would mean the flag never reaches the intensity measures. And the
+    downstream consumer opens IM files with `mask_and_scale=False`, so an
+    integer with a `_FillValue` would read back raw and become a plausible
+    penetration depth; only a real float NaN survives that.
+    """
+    ffp = write_sw4_station_file(
+        tmp_path / "stations.h5",
+        {
+            "AAAA": {"SGDEPTH": 0.0, "SGDEPTHGP": 0.0},
+            "BBBB": {"SGDEPTH": 5750.0, "SGDEPTHGP": 14.375},
+        },
+    )
+
+    dset = lf_to_xarray.read_station_metadata(ffp)
+
+    for name in ("supergrid_depth", "supergrid_depth_gp"):
+        assert name in dset.coords
+        assert name not in dset.data_vars
+        assert dset.coords[name].dims == ("station",)
+        assert dset.coords[name].dtype == np.float32
+
+    ordered = dset.sortby("station")
+    np.testing.assert_array_equal(ordered["supergrid_depth"].values, [0.0, 5750.0])
+    np.testing.assert_allclose(
+        ordered["supergrid_depth_gp"].values, [0.0, 14.375], rtol=1e-6
+    )
+
+
+def test_an_old_station_file_converts_with_an_all_nan_flag(tmp_path: Path) -> None:
+    """A file written before SW4 reported the supergrid must not raise.
+
+    This is the common case for every recording made so far, so it has to be
+    a no-op rather than an error. The value must be NaN and never `0.0`: `0.0`
+    is the positive claim "this station was checked and is in the interior",
+    which nobody checked here.
+    """
+    ffp = write_sw4_station_file(
+        tmp_path / "old.h5", {"AAAA": None, "BBBB": None, "CCCC": None}
+    )
+
+    dset = lf_to_xarray.read_station_metadata(ffp)
+
+    assert dset.sizes["station"] == 3
+    for name in ("supergrid_depth", "supergrid_depth_gp"):
+        assert name in dset.coords
+        assert dset.coords[name].dtype == np.float32
+        assert np.isnan(dset.coords[name].values).all()
+    assert "SGWIDTH" not in dset.attrs
+    assert "SGWIDTHGP" not in dset.attrs
+
+
+def test_stations_missing_the_flag_are_nan_not_zero(tmp_path: Path) -> None:
+    """Mixed groups: only the stations SW4 reported on get a number."""
+    ffp = write_sw4_station_file(
+        tmp_path / "mixed.h5",
+        {
+            "AAAA": {"SGDEPTH": 0.0, "SGDEPTHGP": 0.0},
+            "BBBB": None,
+            "CCCC": {"SGDEPTH": 1200.0, "SGDEPTHGP": 3.0},
+        },
+    )
+
+    depth = lf_to_xarray.read_station_metadata(ffp).sortby("station")["supergrid_depth"]
+
+    assert depth.values[0] == 0.0
+    assert np.isnan(depth.values[1])
+    assert depth.values[2] == 1200.0
+
+
+def test_one_dataset_without_the_other_is_a_corrupt_file(tmp_path: Path) -> None:
+    """`SGDEPTHGP` missing while `SGDEPTH` is present is corruption, not age.
+
+    The back-compatibility guard is deliberately on `SGDEPTH` alone, so this
+    raises rather than quietly reporting a metre depth with no grid-point
+    depth beside it.
+    """
+    ffp = write_sw4_station_file(tmp_path / "corrupt.h5", {"AAAA": {"SGDEPTH": 900.0}})
+
+    with pytest.raises(KeyError):
+        lf_to_xarray.read_station_metadata(ffp)
+
+
+def test_the_sponge_width_is_lifted_into_the_dataset_attributes(
+    tmp_path: Path,
+) -> None:
+    """`SGWIDTH`/`SGWIDTHGP` make the file self-describing.
+
+    They are what turns the penetration into a severity fraction downstream,
+    and taking them from the file rather than from the realisation
+    configuration is the point: the configuration can be edited after the run.
+    """
+    ffp = write_sw4_station_file(
+        tmp_path / "width.h5",
+        {"AAAA": {"SGDEPTH": 0.0, "SGDEPTHGP": 0.0}},
+        widths={"SGWIDTH": 12000.0, "SGWIDTHGP": 30.0},
+    )
+
+    dset = lf_to_xarray.read_station_metadata(ffp)
+
+    assert dset.attrs["SGWIDTH"] == pytest.approx(12000.0)
+    assert dset.attrs["SGWIDTHGP"] == pytest.approx(30.0)
+    # The pre-existing attributes must survive alongside them.
+    assert dset.attrs["nt"] == 8
+    assert dset.attrs["dt"] == pytest.approx(0.05)
+
+
+def test_the_flag_survives_a_netcdf_round_trip(tmp_path: Path) -> None:
+    """As a coordinate, and as NaN -- checked the way a consumer reads it.
+
+    `mask_and_scale=False` is what `eqvis`'s `open_ims` passes, so this is the
+    exact read path the flag has to survive: no fill-value decoding, NaN read
+    straight off disk.
+    """
+    ffp = write_sw4_station_file(
+        tmp_path / "roundtrip.h5",
+        {
+            "AAAA": {"SGDEPTH": 0.0, "SGDEPTHGP": 0.0},
+            "BBBB": None,
+            "CCCC": {"SGDEPTH": 5750.0, "SGDEPTHGP": 14.0},
+        },
+        widths={"SGWIDTH": 12000.0, "SGWIDTHGP": 30.0},
+    )
+    dset = lf_to_xarray.convert_sw4_station_recording(ffp)
+    output = tmp_path / "lf.nc"
+    # The same engine `lf-to-xarray` itself writes with.
+    dset.to_netcdf(output, engine="h5netcdf")
+
+    with xr.open_dataset(output, mask_and_scale=False) as reopened:
+        assert "supergrid_depth" in reopened.coords
+        assert "supergrid_depth" not in reopened.data_vars
+        assert reopened["supergrid_depth"].dtype == np.float32
+        depth = reopened.sortby("station")["supergrid_depth"].values
+        assert depth[0] == 0.0
+        assert np.isnan(depth[1])
+        assert depth[2] == 5750.0
+        assert reopened.attrs["SGWIDTH"] == pytest.approx(12000.0)