From ee798c77d16743e57acff44e3cca9db8401a68af Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 31 Aug 2026 20:27:45 +1200 Subject: [PATCH 1/2] Realisation config plumbing, and SW4 station coordinates Groups the behaviour-neutral housekeeping that later layers of this stack depend on, so it can be reviewed in one pass: - `realisations.path_serialiser`, wired into `write_to_realisation` and `migrate.print_diff`. Later config sections carry `Path` values, which `json` cannot serialise natively. - `nshm2022_to_realisation` opens the NSHMDB inside a `with` block, as nshmdb 2026.08.3 requires. Bumps the floor accordingly. - `srf-to-hdf5`, a thin wrapper over `srf.write_sw4_hdf5`. - `dt` moves from `Resolution` onto `SRFConfig`. The SRF timestep is a property of the SRF discretisation, not of the simulation grid, and `realisation-to-srf` no longer reads `Resolution` at all. It is skipped when building the genslip command line, which takes `dt` explicitly. - `stoch_dx`/`stoch_dy` move out of `HFConfig` into a new `StochConfig`. They size the stoch grid, which is an *input* to the high-frequency simulation rather than one of its parameters. - `generate-station-coordinates` gains a `--format` flag. The EMOD3D path is unchanged, only lifted out of the command body into `write_emod3d_station_format` so a second writer can sit beside it, with the shared `.ll` list factored into `write_ascii_station_locations`. The SW4 writer emits the HDF5 station file `rechdf5` expects: one group per station holding an `STLA,STLO,STDP` triple. Stations outside the domain are dropped by the SW4 writer rather than clamped, because unlike EMOD3D's grid-point files there is no nearest gridpoint to snap them to -- a station SW4 cannot place is not a recording. The containment test is done in NZTM against the domain polygon, matching how the domain itself is defined. `h5py` becomes a declared dependency here: this is the first place workflow code imports it directly rather than relying on it arriving through `xarray[io]`. No change to any computed result: the genslip command line and the srf2stoch invocation are byte-identical to before, as is the EMOD3D station output. Co-Authored-By: Claude Opus 5 --- pyproject.toml | 2 + tests/test_hf.py | 2 - tests/test_realisation.py | 2 + tests/test_realisation_to_srf.py | 1 + uv.lock | 2 + .../default_parameters/root/defaults.yaml | 2 + workflow/realisations.py | 45 ++++- workflow/schemas.py | 8 + .../scripts/generate_station_coordinates.py | 174 ++++++++++++++---- workflow/scripts/generate_stoch.py | 8 +- workflow/scripts/migrate.py | 4 +- workflow/scripts/nshm2022_to_realisation.py | 37 ++-- workflow/scripts/realisation_to_srf.py | 15 +- workflow/scripts/srf_to_hdf5.py | 51 +++++ 14 files changed, 281 insertions(+), 72 deletions(-) create mode 100755 workflow/scripts/srf_to_hdf5.py diff --git a/pyproject.toml b/pyproject.toml index c7a65d2d..085ca9dc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -34,6 +34,7 @@ dependencies = [ "schema", # For loading realisations "structlog", # Logging. "psutil", # To get the CPU affinity for jobs + "h5py>=3.15.1", "parse>=1.21.0", "rich>=14.3.2", ] @@ -66,6 +67,7 @@ generate-station-coordinates = "workflow.scripts.generate_station_coordinates:ap generate-model-coordinates = "workflow.scripts.generate_model_coordinates:app" generate-rupture-propagation = "workflow.scripts.generate_rupture_propagation:app" copy-domain-parameters = "workflow.scripts.copy_velocity_model_parameters:app" +srf-to-hdf5 = "workflow.scripts.srf_to_hdf5:app" create-e3d-par = "workflow.scripts.create_e3d_par:app" generate-stoch = "workflow.scripts.generate_stoch:app" merge-ts = "workflow.scripts.merge_ts:app" diff --git a/tests/test_hf.py b/tests/test_hf.py index 0a4a5358..c006ca3c 100644 --- a/tests/test_hf.py +++ b/tests/test_hf.py @@ -48,8 +48,6 @@ def test_build_hf_input_serialisation() -> None: t_sec=0.0, site_specific=False, dpath_pert=0, - stoch_dx=2.0, - stoch_dy=2.0, stress_parameter_adjustment_fault_area=None, stress_parameter_adjustment_target_magnitude=None, stress_parameter_adjustment_tect_type=0, diff --git a/tests/test_realisation.py b/tests/test_realisation.py index 20c0efa8..4a2223d4 100644 --- a/tests/test_realisation.py +++ b/tests/test_realisation.py @@ -115,6 +115,7 @@ def test_srf_config_example(tmp_path: Path) -> None: ) srf_config = realisations.SRFConfig( resolution=0.1, + dt=0.005, point_source_params=schemas.PointSourceParams( stype=schemas.Stype.cos, risetime=0.5, @@ -249,6 +250,7 @@ def test_srf_config_example(tmp_path: Path) -> None: }, "srf": { "resolution": 0.1, + "dt": 0.005, "point_source_params": { "stype": "cos", "risetime": 0.5, diff --git a/tests/test_realisation_to_srf.py b/tests/test_realisation_to_srf.py index d38d38ab..863f56ab 100644 --- a/tests/test_realisation_to_srf.py +++ b/tests/test_realisation_to_srf.py @@ -8,6 +8,7 @@ def test_build_genslip_command_static_args() -> None: srf_config = SRFConfig( resolution=0.1, + dt=0.005, point_source_params=schemas.PointSourceParams( stype=schemas.Stype.cos, risetime=0.5, diff --git a/uv.lock b/uv.lock index 21e90bda..104bbcab 100644 --- a/uv.lock +++ b/uv.lock @@ -3471,6 +3471,7 @@ name = "workflow" source = { editable = "." } dependencies = [ { name = "geopandas" }, + { name = "h5py" }, { name = "im-calculation" }, { name = "nshmdb" }, { name = "numpy" }, @@ -3518,6 +3519,7 @@ requires-dist = [ { name = "coverage", extras = ["toml"], marker = "extra == 'test'" }, { name = "deptry", marker = "extra == 'dev'" }, { name = "geopandas" }, + { name = "h5py", specifier = ">=3.15.1" }, { name = "hypothesis", extras = ["numpy"], marker = "extra == 'test'", specifier = ">=6.0.0" }, { name = "im-calculation", specifier = ">=2025.12.5" }, { name = "nshmdb", specifier = ">=2026.9.1" }, diff --git a/workflow/default_parameters/root/defaults.yaml b/workflow/default_parameters/root/defaults.yaml index 25ff26e7..d34d270b 100644 --- a/workflow/default_parameters/root/defaults.yaml +++ b/workflow/default_parameters/root/defaults.yaml @@ -93,6 +93,7 @@ hf: stress_parameter_adjustment_tect_type: 0 stress_parameter_adjustment_target_magnitude: null stress_parameter_adjustment_fault_area: null +stoch: stoch_dx: 2.0 stoch_dy: 2.0 rupture_velocity: @@ -163,6 +164,7 @@ srf: read_erf: false read_gsf: true resolution: 0.1 + dt: 0.005 risetime_coef: 1.6 risetimedep: 6.5 risetimedep_range: 1.5 diff --git a/workflow/realisations.py b/workflow/realisations.py index ebb4ec75..b4d8b432 100644 --- a/workflow/realisations.py +++ b/workflow/realisations.py @@ -77,6 +77,25 @@ class RealisationParseError(Exception): """Realisation JSON parse error.""" +def path_serialiser(obj: Any) -> Any: + """Serialise `Path` values, which `json` does not handle natively. + + Parameters + ---------- + obj : Any + The object `json` could not serialise. + + Returns + ------- + Any + The string form of `obj` if it is a `Path`, and `obj` unchanged + otherwise. + """ + if isinstance(obj, Path): + return str(obj) + return obj + + @dataclasses.dataclass class RealisationConfiguration(ABC): """Abstract base class for RealisationConfiguration.""" @@ -220,7 +239,12 @@ def write_to_realisation( realisation_configuration = json.load(realisation_file_handle) realisation_configuration.update({self._config_key: self.to_dict()}) with open(realisation_ffp, "w", encoding="utf-8") as realisation_file_handle: - json.dump(realisation_configuration, realisation_file_handle, indent=4) + json.dump( + realisation_configuration, + realisation_file_handle, + indent=4, + default=path_serialiser, + ) @dataclasses.dataclass @@ -373,6 +397,9 @@ class SRFConfig(RealisationConfiguration): resolution: float """The resolution of the SRF discretisation (different, in general, from the simulation resolution).""" + dt: float + """SRF temporal resolution (timestep).""" + point_source_params: schemas.PointSourceParams | None """Parameters for point source approximation, if applicable.""" @@ -1049,11 +1076,23 @@ class HFConfig(RealisationConfiguration): """Target magnitude (or inferred if None)""" stress_parameter_adjustment_fault_area: float | None """Target magnitude (or inferred if None)""" - # these are used in stoch generation, rather than HF invocation + + +@dataclasses.dataclass +class StochConfig(RealisationConfiguration): + """Stoch file generation. + + Not part of :class:`HFConfig`: these size the stoch grid, which is an input to the + high-frequency simulation rather than one of its parameters. + """ + + _config_key: ClassVar[str] = "stoch" + _schema: ClassVar[Schema] = schemas.STOCH_CONFIG_SCHEMA + stoch_dx: float """stoch file resolution in x.""" stoch_dy: float - """stoch file resolution in x.""" + """stoch file resolution in y.""" @dataclasses.dataclass diff --git a/workflow/schemas.py b/workflow/schemas.py index 6f5e54cf..81460cf0 100644 --- a/workflow/schemas.py +++ b/workflow/schemas.py @@ -428,6 +428,9 @@ def _corners_to_array(corners_spec: list[dict[str, float]]) -> np.ndarray: Literal( "resolution", description="The resolution of the SRF discretisation." ): And(NUMBER, _is_positive), + Literal("dt", description="SRF temporal resolution (timestep)."): And( + NUMBER, _is_positive + ), Optional( Literal( "point_source_params", @@ -1069,6 +1072,11 @@ def _corners_to_array(corners_spec: list[dict[str, float]]) -> np.ndarray: Literal( "stress_parameter_adjustment_fault_area", "Fault area (or inferred if null)" ): Or(NUMBER, None), + } +) + +STOCH_CONFIG_SCHEMA = Schema( + { Literal("stoch_dx", description="Stoch file dx"): And(NUMBER, _is_positive), Literal("stoch_dy", description="Stoch file dy"): And(NUMBER, _is_positive), } diff --git a/workflow/scripts/generate_station_coordinates.py b/workflow/scripts/generate_station_coordinates.py index 13ac2535..b4e9d538 100644 --- a/workflow/scripts/generate_station_coordinates.py +++ b/workflow/scripts/generate_station_coordinates.py @@ -28,10 +28,14 @@ See the output of `generate-station-coordinates --help`. """ +from enum import StrEnum, auto from pathlib import Path from typing import Annotated +import h5py +import numpy as np import pandas as pd +import shapely import typer from qcore import cli, coordinates @@ -41,29 +45,66 @@ app = typer.Typer() -@cli.from_docstring(app) -@log_utils.log_call() -def generate_fd_files( - realisation_ffp: Annotated[Path, typer.Argument(readable=True, dir_okay=False)], - stat_file: Annotated[Path, typer.Argument(readable=True, dir_okay=False)], - output_path: Annotated[Path, typer.Argument(file_okay=False, writable=True)], +class Format(StrEnum): + """The solver the station files are written for.""" + + EMOD3D = auto() + """EMOD3D `.statcords` grid-point files.""" + SW4 = auto() + """SW4 HDF5 station file.""" + + +def write_ascii_station_locations( + stations: pd.DataFrame, + ll_out: Path, + lon_col: str = "lon", + lat_col: str = "lat", + name_col: str = "name", ) -> None: - """Generate station coordinate files. + """Write stations to a whitespace-separated `lon lat name` file. + + Both solvers take the same `.ll` station list, so it is written once + here from whichever coordinate columns the caller has. Parameters ---------- - realisation_ffp : Path - Path to realisation json file. - stat_file : Path - The location of the station files. + stations : pd.DataFrame + The stations to write. + ll_out : Path + Path of the file to write. + lon_col, lat_col, name_col : str, optional + Columns of `stations` holding the longitude, latitude and name. + """ + with open(ll_out, "w", encoding="utf-8") as llf: + stations.apply( + lambda station: llf.write( + f"{station[lon_col]:11.5f} {station[lat_col]:11.5f} {station[name_col]}\n" + ), + axis=1, + ) + + +def write_emod3d_station_format( + domain_parameters: DomainParameters, + resolution_parameters: Resolution, + stations: pd.DataFrame, + output_path: Path, +) -> None: + """Write station coordinates in EMOD3D format to two output files. + + Parameters + ---------- + domain_parameters : DomainParameters + Object containing domain definition and methods to compute grid dimensions. + resolution_parameters : Resolution + Object containing the grid resolution. + stations : pd.DataFrame + DataFrame with columns `lat`, `lon`, and `name`. Latitude and longitude + are in degrees. output_path : Path - Output path for station files. + Directory path where the output files will be written. """ - output_path.mkdir(exist_ok=True) - domain_parameters = DomainParameters.read_from_realisation(realisation_ffp) - resolution_parameters = Resolution.read_from_realisation(realisation_ffp) domain = domain_parameters.domain - nx = domain_parameters.nx(resolution_parameters.resolution) ny = domain_parameters.ny(resolution_parameters.resolution) mlat, mlon = domain.origin @@ -74,14 +115,6 @@ def generate_fd_files( gp_out = output_path / "stations.statcords" ll_out = output_path / "stations.ll" - # retrieve in station names, latitudes and longitudes - stations = pd.read_csv( - stat_file, - delimiter=r"\s+", - comment="#", - names=["lon", "lat", "name"], - ) - x, y = proj( lat=stations["lat"].to_numpy(float), lon=stations["lon"].to_numpy(float) ).T @@ -127,13 +160,90 @@ def generate_fd_files( axis=1, ) - # create ll file - with open(ll_out, "w", encoding="utf-8") as llf: - stations.apply( - lambda station: llf.write( - f"{station['grid_lon']:11.5f} {station['grid_lat']:11.5f} {station['name']}\n" - ), - axis=1, - ) + write_ascii_station_locations( + stations, ll_out, lon_col="grid_lon", lat_col="grid_lat" + ) + + +def write_sw4_station_format( + domain_parameters: DomainParameters, stations: pd.DataFrame, output_path: Path +) -> None: + """Write station coordinates in SW4 format. + + Stations outside the domain are dropped rather than clamped to its + edge: unlike EMOD3D's grid-point files there is no nearest gridpoint + to snap to, and a station SW4 cannot place is not a recording. + + Parameters + ---------- + domain_parameters : DomainParameters + Domain definition, used to test which stations fall inside it. + stations : pd.DataFrame + DataFrame with columns `lat`, `lon`, and `name`. Latitude and + longitude are in degrees. + output_path : Path + Directory to write `stations.h5` and `stations.ll` into. + """ + lat_lon = stations[["lat", "lon"]].to_numpy() + nzvm_coordinates = coordinates.wgs_depth_to_nztm(lat_lon) + poly = domain_parameters.domain.polygon + mask = shapely.contains_xy(poly, nzvm_coordinates[:, 0], nzvm_coordinates[:, 1]) + stations = stations.loc[mask] + + with h5py.File(output_path / "stations.h5", "w") as f: + for station_name, position in stations.set_index("name").iterrows(): + station_dset = f.create_group(station_name) + location = station_dset.create_dataset( + "STLA,STLO,STDP", (3,), dtype=np.float64 + ) + location[0] = position["lat"] + location[1] = position["lon"] + location[2] = 0.0 + + write_ascii_station_locations(stations, output_path / "stations.ll") + + +@cli.from_docstring(app) +@log_utils.log_call() +def generate_fd_files( + realisation_ffp: Annotated[Path, typer.Argument(readable=True, dir_okay=False)], + stat_file: Annotated[Path, typer.Argument(readable=True, dir_okay=False)], + output_path: Annotated[Path, typer.Argument(file_okay=False, writable=True)], + format: Format = Format.EMOD3D, +) -> None: + """Generate station coordinate files. + + Parameters + ---------- + realisation_ffp : Path + Path to realisation json file. + stat_file : Path + The location of the station files. + output_path : Path + Output path for station files. + format : Format, optional + The solver to write station files for. EMOD3D writes + `stations.statcords` and `stations.ll`; SW4 writes `stations.h5` + and `stations.ll`. Defaults to EMOD3D. + """ + output_path.mkdir(exist_ok=True) + domain_parameters = DomainParameters.read_from_realisation(realisation_ffp) + resolution_parameters = Resolution.read_from_realisation(realisation_ffp) + + # retrieve in station names, latitudes and longitudes + stations = pd.read_csv( + stat_file, + delimiter=r"\s+", + comment="#", + names=["lon", "lat", "name"], + ) + + match format: + case Format.EMOD3D: + write_emod3d_station_format( + domain_parameters, resolution_parameters, stations, output_path + ) + case Format.SW4: + write_sw4_station_format(domain_parameters, stations, output_path) realisations.append_log_entry(realisation_ffp) diff --git a/workflow/scripts/generate_stoch.py b/workflow/scripts/generate_stoch.py index 667e3a53..a9ebfaaf 100644 --- a/workflow/scripts/generate_stoch.py +++ b/workflow/scripts/generate_stoch.py @@ -33,7 +33,7 @@ from qcore import cli from source_modelling import sources, srf from workflow import log_utils, realisations -from workflow.realisations import HFConfig, RealisationMetadata, SourceConfig +from workflow.realisations import RealisationMetadata, SourceConfig, StochConfig app = typer.Typer() @@ -64,7 +64,7 @@ def generate_stoch( Path to the `srf2stoch` binary used for the conversion. """ metadata = RealisationMetadata.read_from_realisation(realisation_ffp) - hf_config = HFConfig.read_from_realisation_or_defaults( + stoch_config = StochConfig.read_from_realisation_or_defaults( realisation_ffp, metadata.defaults_version ) @@ -87,8 +87,8 @@ def generate_stoch( min_length = min(fault.length for fault in geometries) min_width = min(fault.width for fault in geometries) # If the stoch dx is greater than the length (resp. dy and width), we might get an empty stoch file - dx = min(hf_config.stoch_dx, min_length / 2) - dy = min(hf_config.stoch_dy, min_width / 2) + dx = min(stoch_config.stoch_dx, min_length / 2) + dy = min(stoch_config.stoch_dy, min_width / 2) log_utils.log_check_call( [ diff --git a/workflow/scripts/migrate.py b/workflow/scripts/migrate.py index 1895324e..055f2f72 100644 --- a/workflow/scripts/migrate.py +++ b/workflow/scripts/migrate.py @@ -304,8 +304,8 @@ def print_diff(config_a: dict, config_b: dict) -> None: config_b : dict Dictionary b. """ - config_a_str = json.dumps(config_a, indent=4) - config_b_str = json.dumps(config_b, indent=4) + config_a_str = json.dumps(config_a, indent=4, default=realisations.path_serialiser) + config_b_str = json.dumps(config_b, indent=4, default=realisations.path_serialiser) diff = difflib.context_diff( config_a_str.splitlines(keepends=True), diff --git a/workflow/scripts/nshm2022_to_realisation.py b/workflow/scripts/nshm2022_to_realisation.py index ecf0ece3..a2d40b68 100755 --- a/workflow/scripts/nshm2022_to_realisation.py +++ b/workflow/scripts/nshm2022_to_realisation.py @@ -317,19 +317,19 @@ def generate_realisation( realisation_ffp, defaults_version ) db = nshmdb.NSHMDB(nshmdb_path) - db.connect() - faults = db.get_rupture_faults(fault_system, rupture_id) - faults = { - fault_name: sources.simplify_fault(fault, srf_config.resolution) - for fault_name, fault in faults.items() - } - - if initial_fault and initial_fault not in faults: - print( - f"Initial fault '{initial_fault}' not found in rupture. Options are {', '.join(list(faults))}" - ) - raise typer.Exit(code=1) - faults_info = db.get_rupture_fault_info(fault_system, rupture_id) + with db: + faults = db.get_rupture_faults(fault_system, rupture_id) + faults = { + fault_name: sources.simplify_fault(fault, srf_config.resolution) + for fault_name, fault in faults.items() + } + + if initial_fault and initial_fault not in faults: + print( + f"Initial fault '{initial_fault}' not found in rupture. Options are {', '.join(list(faults))}" + ) + raise typer.Exit(code=1) + faults_info = db.get_rupture_fault_info(fault_system, rupture_id) seeds = Seeds.read_from_realisation_or_random(realisation_ffp) np.random.seed(seed=seeds.nshm_to_realisation_seed) random.seed(seeds.nshm_to_realisation_seed) @@ -361,11 +361,12 @@ def generate_realisation( else: # The ty ignore below can be removed once NSHM2022DB merges the change to # accept dict[str, BoldM] in most_likely_fault (branch support-BoldM-in-workflow). - mfds_rates = db.most_likely_fault( - fault_system, - rupture_id, - magnitudes, # ty: ignore[invalid-argument-type] - ) + with db: + mfds_rates = db.most_likely_fault( + fault_system, + rupture_id, + magnitudes, # ty: ignore[invalid-argument-type] + ) mfds_probabilities = np.array(list(mfds_rates.values())) if np.allclose(mfds_probabilities, 0): mfds_probabilities = np.ones_like(mfds_probabilities) diff --git a/workflow/scripts/realisation_to_srf.py b/workflow/scripts/realisation_to_srf.py index fb9038e9..bc2dfb78 100644 --- a/workflow/scripts/realisation_to_srf.py +++ b/workflow/scripts/realisation_to_srf.py @@ -61,7 +61,6 @@ Magnitudes, Rakes, RealisationMetadata, - Resolution, RupturePropagationConfig, RuptureVelocity, Seeds, @@ -364,8 +363,6 @@ def process_fault(fault_name: str) -> None: class SRFRealisationContext: """Realisation configuration for the entire SRF generation process.""" - resolution: Resolution - """The spatial/temporal resolution""" source_config: SourceConfig """The sources to generate""" rupture_propagation_config: RupturePropagationConfig @@ -486,7 +483,7 @@ def _build_genslip_command( if rupture_velocity.rvfrac_slip_sig is not None: cmd.append(f"rvfrac_slip_sig={rupture_velocity.rvfrac_slip_sig}") - skipped_fields = {"point_source_params"} + skipped_fields = {"point_source_params", "dt"} for field in dataclasses.fields(srf_config): key = field.name value = getattr(srf_config, key) @@ -563,7 +560,7 @@ def generate_fault_srf( shypo=genslip_hypocentre_coords[0], dhypo=genslip_hypocentre_coords[1], magnitude=params.magnitudes.magnitudes[name], - dt=params.resolution.dt, + dt=params.srf_config.dt, srf_config=params.srf_config, rupture_velocity=params.rupture_velocity, ) @@ -713,7 +710,7 @@ def generate_point_source_srf( fault = params.source_config.source_geometries[name] - resolution = params.resolution.resolution + resolution = params.srf_config.resolution # Get magnitude and convert to seismic moment magnitude = params.magnitudes.magnitudes[name] @@ -753,7 +750,7 @@ def generate_point_source_srf( f"outfile={environment.srf_directory / (normalise_name(name) + '.srf')}", "outbin=0", f"stype={params.srf_config.point_source_params.stype}", - f"dt={params.resolution.dt}", + f"dt={params.srf_config.dt}", "plane_header=1", f"risetime={params.srf_config.point_source_params.risetime}", f"risetimefac={params.srf_config.point_source_params.risetimefac}", @@ -850,15 +847,11 @@ def generate_srf( rakes = Rakes.read_from_realisation(realisation_ffp) magnitudes = Magnitudes.read_from_realisation(realisation_ffp) source_config = SourceConfig.read_from_realisation(realisation_ffp) - resolution = Resolution.read_from_realisation_or_defaults( - realisation_ffp, metadata.defaults_version - ) rupture_velocity = RuptureVelocity.read_from_realisation_or_defaults( realisation_ffp, metadata.defaults_version ) params = SRFRealisationContext( - resolution=resolution, source_config=source_config, rupture_propagation_config=rupture_propagation, magnitudes=magnitudes, diff --git a/workflow/scripts/srf_to_hdf5.py b/workflow/scripts/srf_to_hdf5.py new file mode 100755 index 00000000..dd26a7da --- /dev/null +++ b/workflow/scripts/srf_to_hdf5.py @@ -0,0 +1,51 @@ +#!/usr/bin/env python3 +"""SRF to HDF5. + +Description +----------- +Convert an SRF file into SW4's SRF-HDF5 source format. This is a thin wrapper over `source_modelling.srf.SrfFile.write_sw4_hdf5` so the conversion can be run as a workflow stage. + +Inputs +------ +An SRF file, for example one produced by `realisation-to-srf`. + +Outputs +------- +The same rupture written as an [SRF-HDF5](https://github.com/geodynamics/sw4/blob/master/doc/SW4_UsersGuide.pdf) file, which is the source description SW4 reads. + +Environment +----------- +Can be run in the cybershake container. Can also be run from your own computer using the `srf-to-hdf5` command which is installed after running `pip install workflow@git+https://github.com/ucgmsim/workflow`. + +Usage +----- +`srf-to-hdf5 [OPTIONS] SRF_PATH HDF5_PATH` + +For More Help +------------- +See the output of `srf-to-hdf5 --help`. +""" + +from pathlib import Path + +import typer + +from qcore import cli +from source_modelling import srf + +app = typer.Typer() + + +@cli.from_docstring(app) +def srf_to_hdf5(srf_path: Path, hdf5_path: Path) -> None: + """Convert an SRF file to SW4's SRF-HDF5 format. + + Parameters + ---------- + srf_path : Path + Path to the input SRF file. + hdf5_path : Path + Path where the output HDF5 file will be written. + """ + srf_file = srf.read_srf(srf_path) + srf_file.write_sw4_hdf5(hdf5_path) From 0173c25d95102bb85a1dc8845074fb9a852fbe2f Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 3 Sep 2026 09:07:23 +1200 Subject: [PATCH 2/2] Apply batched suggestions from code review Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- workflow/scripts/generate_station_coordinates.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/workflow/scripts/generate_station_coordinates.py b/workflow/scripts/generate_station_coordinates.py index b4e9d538..84389daf 100644 --- a/workflow/scripts/generate_station_coordinates.py +++ b/workflow/scripts/generate_station_coordinates.py @@ -190,6 +190,9 @@ def write_sw4_station_format( mask = shapely.contains_xy(poly, nzvm_coordinates[:, 0], nzvm_coordinates[:, 1]) stations = stations.loc[mask] + if len(stations) == 0: + raise ValueError("No stations in domain.") + with h5py.File(output_path / "stations.h5", "w") as f: for station_name, position in stations.set_index("name").iterrows(): station_dset = f.create_group(station_name)