Read SW4 station recordings, carrying supergrid penetration through - #114
Read SW4 station recordings, carrying supergrid penetration through#114lispandfound wants to merge 1 commit into
Conversation
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
There was a problem hiding this comment.
Code Review
This pull request adds support for converting SW4 HDF5 station recordings to xarray datasets, introducing dask as a dependency to handle chunkwise data loading. The review feedback highlights several robust improvements to the SW4 parsing logic, including adding explicit type checks to avoid processing non-group datasets, using [()] instead of [0] or direct conversion to safely retrieve scalar values from HDF5 datasets, correcting a misleading comment regarding the CMS conversion unit, and adding a guard check to handle cases where no valid station groups are found.
Important
The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.
|
@gemini-code-assist please re-review |
There was a problem hiding this comment.
Code Review
This pull request introduces support for converting SW4 HDF5 station recordings to xarray datasets using Dask and h5py, and adds waveform resampling functionality to the broadband simulation script. Key feedback points out critical bugs: SciPy submodules like signal must be imported directly to avoid AttributeErrors, and h5py.Group type checks are needed to prevent TypeErrors when iterating over HDF5 file items. Additionally, thread-safe access should be ensured by passing lock=True to da.from_array, and floating-point comparisons for dt should use np.isclose to avoid precision issues.
Important
The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.
| for station_name, group in handle.items(): | ||
| if "NPTS" not in group: | ||
| continue |
There was a problem hiding this comment.
The root of an SW4 HDF5 file contains datasets like DELTA, CYCLE, and TIME in addition to station groups. Since h5py.Dataset does not implement contains (and is not a container), checking 'NPTS' not in group when group is a dataset will raise a TypeError. You should first verify that group is an instance of h5py.Group.
| for station_name, group in handle.items(): | |
| if "NPTS" not in group: | |
| continue | |
| for station_name, group in handle.items(): | |
| if not isinstance(group, h5py.Group) or 'NPTS' not in group: | |
| continue |
There was a problem hiding this comment.
This looks like a valid comment from gemini?
| x = da.from_array(group["X"], chunks=-1) | ||
| xs.append(x * CMS) | ||
| y = da.from_array(group["Y"], chunks=-1) | ||
| ys.append(y * CMS) | ||
| z = da.from_array(group["Z"], chunks=-1) | ||
| zs.append(z * CMS) |
There was a problem hiding this comment.
h5py is not thread-safe by default. When Dask executes the graph using its default multi-threaded scheduler, concurrent reads from the same h5py.File can cause race conditions, hangs, or segmentation faults. You should pass lock=True to da.from_array to ensure thread-safe access.
| x = da.from_array(group["X"], chunks=-1) | |
| xs.append(x * CMS) | |
| y = da.from_array(group["Y"], chunks=-1) | |
| ys.append(y * CMS) | |
| z = da.from_array(group["Z"], chunks=-1) | |
| zs.append(z * CMS) | |
| x = da.from_array(group['X'], chunks=-1, lock=True) | |
| xs.append(x * CMS) | |
| y = da.from_array(group['Y'], chunks=-1, lock=True) | |
| ys.append(y * CMS) | |
| z = da.from_array(group['Z'], chunks=-1, lock=True) | |
| zs.append(z * CMS) |
There was a problem hiding this comment.
@lispandfound Please comment why you thought it doesn't matter?
![github-actions[bot]](https://avatars.githubusercontent.com/in/15368?s=60&v=4){kind=small}
📁 workflow/scripts/lf_to_xarray.py (deepseek-v4-flash:cloud)📍 Line 1 ✅ Fix: Use a custom enum class or require Python >= 3.11. |
🤖 AI Code Review (deepseek-v4-flash:cloud)Found 4 issues. See inline comments for details. |
392031c
| ] | ||
| sdist = { url = "https://files.pythonhosted.org/packages/ea/0d/837dbd5d8430fd0f01ed72c4cfb2f548180f4c68c635df84ce87956cff32/cligj-0.7.2.tar.gz", hash = "sha256:a4bc13d623356b373c2c27c53dbd9c68cae5d526270bfa71f6c6fa69669c6b27", size = 9803, upload-time = "2021-05-28T21:23:27.935Z" } | ||
| wheels = [ | ||
| { url = "https://files.pythonhosted.org/packages/73/86/43fa9f15c5b9fb6e82620428827cd3c284aa933431405d1bcf5231ae3d3e/cligj-0.7.2-py3-none-any.whl", hash = "sha256:c1ca117dbce1fe20a5809dc96f01e1c2840f6dcc939b3ddbb1111bf330ba82df", size = 7069, upload-time = "2021-05-28T21:23:26.877Z" }, |
There was a problem hiding this comment.
📍 Line 420
source = { registry = "https://pypi.org/simple" }
dependencies = [
{ name = "click" },
]
sdist = { url = "https://files.pythonhosted.org/packages/ea/0d/837dbd5d8430fd0f01ed72c4cfb2f548180f4c68c635df84ce87956cff32/cligj-0.7.2.tar.gz", hash = "sha256:a4bc13d623356b373c2c27c53dbd9c68cae5d526270bfa71f6c6fa69669c6b27", size = 9803, upload-time = "2021-05-28T21:23:27.935Z" }
wheels = [
👉 { url = "https://files.pythonhosted.org/packages/73/86/43fa9f15c5b9fb6e82620428827cd3c284aa933431405d1bcf5231ae3d3e/cligj-0.7.2-py3-none-any.whl", hash = "sha256:c1ca117dbce1fe20a5809dc96f01e1c2840f6dcc939b3ddbb1111bf330ba82df", size = 7069, upload-time = "2021-05-28T21:23:26.877Z" },
]
[[package]]
name = "cloudpickle"
version = "3.1.2"
source = { registry = "https://pypi.org/simple" }
✅ Fix: Regenerate the lock file with correct package versions after fixing pyproject.toml.
| def convert_lf_to_xarray_dataset(lfseis_directory: Path, output_ffp: Path) -> None: | ||
| def convert_lf_to_xarray_dataset( | ||
| low_frequency_path: Path, output_ffp: Path, format: Format = Format.EMOD3D | ||
| ) -> None: |
There was a problem hiding this comment.
📍 Line 119
@cli.from_docstring(app)
@log_utils.log_call()
def convert_lf_to_xarray_dataset(
low_frequency_path: Path, output_ffp: Path, format: Format = Format.EMOD3D
👉 ) -> None:
"""Merge low-frequency outputs into an xarray dataset.
Parameters
----------
low_frequency_path : Path
Directory containing station seismogram outputs.
"start_sec" attribute, which is expected by downstream code (e.g., resample_signal in bb_sim.py).
✅ Fix: Set the attribute after reading: lf_dataset.attrs['start_sec'] = 0.0 (or appropriate value).
📁 workflow/scripts/lf_to_xarray.py (deepseek-v4-flash:cloud)📍 Line 1
✅ Fix: Use a custom string enum (e.g., |
🤖 AI Code Review (deepseek-v4-flash:cloud)Found 10 issues. See inline comments for details. |
🤖 AI Code Review (deepseek-v4-flash:cloud)❌ No review content |
sungeunbae
left a comment
sungeunbae
left a comment
There was a problem hiding this comment.
Some explanation will be desirable when you dismiss comments
| for station_name, group in handle.items(): | ||
| if "NPTS" not in group: | ||
| continue |
There was a problem hiding this comment.
This looks like a valid comment from gemini?
| x = da.from_array(group["X"], chunks=-1) | ||
| xs.append(x * CMS) | ||
| y = da.from_array(group["Y"], chunks=-1) | ||
| ys.append(y * CMS) | ||
| z = da.from_array(group["Z"], chunks=-1) | ||
| zs.append(z * CMS) |
There was a problem hiding this comment.
@lispandfound Please comment why you thought it doesn't matter?
77c1544 to
8340b60
Compare
lispandfound
changed the title
706a4ed to
87b1c00
Compare
351014c to
b96f1ac
Compare
b96f1ac to
424cbdc
Compare
424cbdc to
e991bd7
Compare
e991bd7 to
6e2f496
Compare
`lf-to-xarray --format sw4` reads SW4's HDF5 station file into the same dataset shape the EMOD3D path produces, so everything downstream is solver-agnostic. Stations are read in batches sized to a target chunk, lazily via `map_blocks`, because a real run's station file does not fit in memory. Two things worth stating explicitly: The datasets carry SW4's displacement-mode names (EW/NS/UP), but for an SRF rupture source the time function SW4 receives is the slip *rate*, so the nominal displacement output is physically velocity. Components are reordered to EMOD3D's convention (x east, y north, z down). SW4 reports, per station, how far into the supergrid sponge that station sits. That number decides whether a trace is a ground motion prediction at all, so it is carried all the way to the intensity measure file as a station-dimension *coordinate*. That is load-bearing rather than incidental: coordinates ride through `bb-sim` and `im-calc` untouched, whereas data variables are dropped when chunks are recombined. A station that reported nothing gets NaN, never 0.0. Zero means "checked and found in the interior" -- a real claim about the run -- and an old station file, or a solver with no absorbing layer, is not entitled to make it. Stored as float32 for the same reason: readers open these files with `mask_and_scale=False`, so an integer sentinel would read back raw and become a plausible penetration depth. `im-calc` attaches the coordinate for every solver, all-NaN where nothing was reported, so the column always exists and its absence is never mistaken for a clean run. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
6e2f496 to
12f9cb5
Compare
Adds a function to read the station HDF5 output from SW4 and convert it into a compliant NetCDF file. This file is compatible with all downstream processing in hf-sim, bb-sim, so this is the only portion of the workflow required to drop-in SW4 instead of EMOD3D (besides the velocity model).