Skip to content

SW4 station post-processing - #114

Open
lispandfound wants to merge 23 commits into
pegasusfrom
sw4_stations
Open

SW4 station post-processing#114
lispandfound wants to merge 23 commits into
pegasusfrom
sw4_stations

Conversation

@lispandfound

Copy link
Copy Markdown
Contributor

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).

Copilot AI review requested due to automatic review settings June 23, 2026 03:21

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

@gemini-code-assist gemini-code-assist Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Comment thread workflow/scripts/lf_to_xarray.py
Comment thread workflow/scripts/lf_to_xarray.py Outdated
Comment thread workflow/scripts/lf_to_xarray.py Outdated
Comment thread workflow/scripts/lf_to_xarray.py Outdated
Comment thread workflow/scripts/lf_to_xarray.py
@lispandfound

Copy link
Copy Markdown
Contributor Author

@gemini-code-assist please re-review

@gemini-code-assist gemini-code-assist Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Comment thread workflow/scripts/bb_sim.py
Comment thread workflow/scripts/bb_sim.py
Comment on lines +74 to +76
for station_name, group in handle.items():
if "NPTS" not in group:
continue

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

high

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.

Suggested change
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

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks like a valid comment from gemini?

Comment on lines +83 to +88
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)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

high

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.

Suggested change
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)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lispandfound Please comment why you thought it doesn't matter?

Comment thread workflow/scripts/bb_sim.py Outdated
joelridden
joelridden previously approved these changes Jun 24, 2026
Comment thread workflow/scripts/bb_sim.py
raise RuntimeError(f"SW4 output is corrupted: {npts=} but {global_npts=}")
global_npts = npts
# Dask arrays here ensure that data is read chunkwise from the HDF5 file
# without putting it all in-memory.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

📍 Line 82

            continue
        npts = int(group["NPTS"][()])
        if global_npts is not None and npts != global_npts:
            raise RuntimeError(f"SW4 output is corrupted: {npts=} but {global_npts=}")
        global_npts = npts
        # Dask arrays here ensure that data is read chunkwise from the HDF5 file
👉         # without putting it all in-memory.
        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)

⚠️ Problem: Using chunks=-1 creates a single chunk, which defeats the purpose of dask for memory efficiency. The comment claims it ensures chunkwise reading, but it does not. This will load the entire array into memory when computed.

✅ Fix: Use a proper chunk size, e.g., chunks='auto' or a specific size like chunks=(1000,) to enable chunkwise processing.

@github-actions

Copy link
Copy Markdown
Contributor

📁 workflow/scripts/lf_to_xarray.py (deepseek-v4-flash:cloud)

📍 Line 1

👉 #!/usr/bin/env python
"""Low-frequency output merger.

Description
-----------
Merges low-frequency outputs into one xarray dataset.

⚠️ Problem: StrEnum is only available in Python 3.11+. If the project supports older Python versions, this will cause an ImportError.

✅ Fix: Use a custom enum class or require Python >= 3.11.

return lf_aligned, hf_aligned, common_time


def resample_signal(dset: xr.Dataset, dt: float) -> xr.Dataset:

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

📍 Line 117

    lf_aligned[:, lf_indices] = lf_waveform
    hf_aligned[:, hf_indices] = hf_waveform

    return lf_aligned, hf_aligned, common_time


👉 def resample_signal(dset: xr.Dataset, dt: float) -> xr.Dataset:
    """Resample waveform dataset to a new time step.

    Parameters
    ----------
    dset : xr.Dataset
        Input dataset with dimensions (component, station, time) and

⚠️ Problem: The function only updates the 'waveform' variable and the 'dt' attribute. If the dataset contains other data variables with a time dimension, their time coordinates will become inconsistent after resampling, leading to incorrect data.

✅ Fix: Either ensure the dataset only contains the waveform variable, or update all variables with time dimension accordingly.

@github-actions

Copy link
Copy Markdown
Contributor

🤖 AI Code Review (deepseek-v4-flash:cloud)

Found 4 issues. See inline comments for details.

Comment thread pyproject.toml
"psutil", # To get the CPU affinity for jobs
"parse>=1.21.0",
"rich>=14.3.2",
"dask>=2026.6.0",

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

📍 Line 39

  "requests",  # For gcmt-to-realisation
  "schema",    # For loading realisations
  "structlog", # Logging.
  "psutil",    # To get the CPU affinity for jobs
  "parse>=1.21.0",
  "rich>=14.3.2",
👉   "dask>=2026.6.0",
  "h5py>=3.15.1",
]

[project.optional-dependencies]
test = [
  "pytest>=6.0.0",            # required for the tool.pytest section

⚠️ Problem: The version specifier 2026.6.0 refers to a future version that does not exist (as of 2025). This will cause dependency resolution failures.

✅ Fix: Use an existing version, e.g., "dask>=2024.12.0".

Comment thread pyproject.toml
"parse>=1.21.0",
"rich>=14.3.2",
"dask>=2026.6.0",
"h5py>=3.15.1",

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

📍 Line 40

  "schema",    # For loading realisations
  "structlog", # Logging.
  "psutil",    # To get the CPU affinity for jobs
  "parse>=1.21.0",
  "rich>=14.3.2",
  "dask>=2026.6.0",
👉   "h5py>=3.15.1",
]

[project.optional-dependencies]
test = [
  "pytest>=6.0.0",            # required for the tool.pytest section
  "hypothesis[numpy]>=6.0.0",

⚠️ Problem: The version specifier 3.15.1 refers to a future version that does not exist. This will cause dependency resolution failures.

✅ Fix: Use an existing version, e.g., "h5py>=3.12.0".

Comment thread uv.lock
@@ -420,6 +420,15 @@ 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" },

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

📍 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" }

⚠️ Problem: The lock file contains entries for packages with future versions (e.g., dask 2026.6.0, cloudpickle 3.1.2) that do not exist, making the lock file invalid and installation impossible.

✅ Fix: Regenerate the lock file with correct package versions after fixing pyproject.toml.

import numpy.typing as npt
import pandas as pd
import scipy as sp
import typer

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

📍 Line 45

from typing import Annotated

import numpy as np
import numpy.typing as npt
import pandas as pd
import scipy as sp
👉 import typer
import xarray as xr

from qcore import cli, siteamp_models, timeseries
from workflow import log_utils, realisations
from workflow.realisations import (
    BroadbandParameters,

⚠️ Problem: scipy is not declared as a dependency in pyproject.toml. The script will fail at runtime if scipy is not installed.

✅ Fix: Add "scipy" to the dependencies in pyproject.toml.

# error not to provide one (no implicit magic behaviour).
new_time = np.arange(nt) * dt + dset.attrs["start_sec"]

resampled_waveform = xr.apply_ufunc(

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

📍 Line 140

    nt = round(duration / dt)

    # NOTE: I am not providing a default start second because we consider it an
    # error not to provide one (no implicit magic behaviour).
    new_time = np.arange(nt) * dt + dset.attrs["start_sec"]

👉     resampled_waveform = xr.apply_ufunc(
        sp.signal.resample,
        dset["waveform"],
        # This tells xarray that resample expects an array with all of the time component intact.
        # So it will be passed arrays of shape (n_component, n_stations, n_time) = (i, j, nt)
        input_core_dims=[["time"]],
        # This tells xarray that the time dimension is going to be returned in

⚠️ Problem: scipy.signal.resample defaults to axis=0, but the time dimension is the last axis in the input array. This will resample along the component dimension instead of time, producing incorrect results.

✅ Fix: Add axis=-1 to kwargs: kwargs=dict(num=nt, axis=-1),

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lispandfound Does this make sense? I am not sure

resampled_waveform = xr.apply_ufunc(
sp.signal.resample,
dset["waveform"],
# This tells xarray that resample expects an array with all of the time component intact.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

📍 Line 143

    # error not to provide one (no implicit magic behaviour).
    new_time = np.arange(nt) * dt + dset.attrs["start_sec"]

    resampled_waveform = xr.apply_ufunc(
        sp.signal.resample,
        dset["waveform"],
👉         # This tells xarray that resample expects an array with all of the time component intact.
        # So it will be passed arrays of shape (n_component, n_stations, n_time) = (i, j, nt)
        input_core_dims=[["time"]],
        # This tells xarray that the time dimension is going to be returned in
        # its entirety by scipy resample.
        output_core_dims=[["time"]],
        # This tells xarray that the time coordinates from the dset dataset are no

⚠️ Problem: The code assumes the dataset has a "start_sec" attribute. If it is missing (e.g., from EMOD3D conversion), a KeyError is raised.

✅ Fix: Use dset.attrs.get("start_sec", 0.0) to provide a default value, or ensure the attribute is always set.

Comment thread workflow/scripts/lf_to_xarray.py
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:

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

📍 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.

⚠️ Problem: The resulting dataset may lack a "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).


class Format(StrEnum):
"""Input low frequency file format."""

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

📍 Line 108

        attrs=attributes,
    )


class Format(StrEnum):
    """Input low frequency file format."""
👉 
    SW4 = auto()
    """SW4 HDF5 station recording."""
    EMOD3D = auto()
    """EMOD3D LFSeis directory."""


⚠️ Problem: The parameter name changed from lfseis_directory to low_frequency_path, breaking backward compatibility for any callers using keyword arguments.

✅ Fix: Keep the old parameter name as an alias or update all callers. Alternatively, add a deprecation warning.

@github-actions

Copy link
Copy Markdown
Contributor

📁 workflow/scripts/lf_to_xarray.py (deepseek-v4-flash:cloud)

📍 Line 1

👉 #!/usr/bin/env python
"""Low-frequency output merger.

Description
-----------
Merges low-frequency outputs into one xarray dataset.

⚠️ Problem: StrEnum is only available in Python 3.11+. If the project supports older Python versions, this will cause a syntax error.

✅ Fix: Use a custom string enum (e.g., class Format(str, Enum)) or require Python 3.11+.

@github-actions

Copy link
Copy Markdown
Contributor

🤖 AI Code Review (deepseek-v4-flash:cloud)

Found 10 issues. See inline comments for details.

@github-actions

Copy link
Copy Markdown
Contributor

🤖 AI Code Review (deepseek-v4-flash:cloud)

❌ No review content

@sungeunbae sungeunbae left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some explanation will be desirable when you dismiss comments

Comment on lines +74 to +76
for station_name, group in handle.items():
if "NPTS" not in group:
continue

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks like a valid comment from gemini?

# error not to provide one (no implicit magic behaviour).
new_time = np.arange(nt) * dt + dset.attrs["start_sec"]

resampled_waveform = xr.apply_ufunc(

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lispandfound Does this make sense? I am not sure

Comment on lines +83 to +88
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)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lispandfound Please comment why you thought it doesn't matter?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants