From 6ec9e495899bc2c924029f0f58fd2baf3e7eb69d Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Fri, 20 Mar 2026 15:30:55 +1300 Subject: [PATCH 01/18] Add write_sw4_hdf5 for writing SRF data in SW4's HDF5 format --- source_modelling/srf.py | 109 ++++++++++++++++++++++++++++++++++++++-- tests/test_srf.py | 57 +++++++++++++++++++++ 2 files changed, 163 insertions(+), 3 deletions(-) diff --git a/source_modelling/srf.py b/source_modelling/srf.py index 4e7377e..bedbdbf 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -47,6 +47,7 @@ from pathlib import Path from typing import Self +import h5py import numpy as np import pandas as pd import scipy as sp @@ -60,6 +61,47 @@ PLANE_COUNT_RE = r"PLANE (\d+)" POINT_COUNT_RE = r"POINTS (\d+)" +SW4_PLANE_DTYPE = np.dtype( + [ + ("ELON", "f4"), + ("ELAT", "f4"), + ("NSTK", "i4"), + ("NDIP", "i4"), + ("LEN", "f4"), + ("WID", "f4"), + ("STK", "f4"), + ("DIP", "f4"), + ("DTOP", "f4"), + ("SHYP", "f4"), + ("DHYP", "f4"), + ] +) + +SW4_POINTS_DTYPE = np.dtype( + [ + ("LON", "f4"), + ("LAT", "f4"), + ("DEP", "f4"), + ("STK", "f4"), + ("DIP", "f4"), + ("AREA", "f4"), + ("TINIT", "f4"), + ("DT", "f4"), + ("VS", "f4"), + ("DEN", "f4"), + ("RAKE", "f4"), + ("SLIP1", "f4"), + ("NT1", "i4"), + ("SLIP2", "f4"), + ("NT2", "i4"), + ("SLIP3", "f4"), + ("NT3", "i4"), + ] +) + +_SW4_POINTS_COL_OVERRIDES = {"SLIP1": "slip"} +_SW4_POINTS_EXTERNAL_FIELDS = {"VS", "DEN", "NT1", "SLIP2", "NT2", "SLIP3", "NT3"} + class Segments(Sequence): """A read-only view for SRF segments. @@ -234,7 +276,7 @@ def from_file(cls, srf_ffp: Path | str) -> Self: point_count = int(points_count_match.group(1)) position = srf_file_handle.tell() - points_metadata, slipt1_array = srf_parser.parse_srf( # type: ignore + points_metadata, slipt1_array = srf_parser.parse_srf( # type: ignore str(srf_ffp), position, point_count ) @@ -290,13 +332,72 @@ def write_srf(self, srf_ffp: Path) -> None: srf_file_handle.write(f"POINTS {len(self.points)}\n") - srf_parser.write_srf_points( # type: ignore + srf_parser.write_srf_points( # type: ignore str(srf_ffp), self.points.values.astype(np.float32), self.slip.indptr, self.slip.data, ) + def write_sw4_hdf5( + self, + output_ffp: Path, + vs: float | np.ndarray = 0.0, + den: float | np.ndarray = 0.0, + include_slip_time_function: bool = True, + ) -> None: + """Write the SRF file in SW4's SRF-HDF5 format. + + Parameters + ---------- + output_ffp : Path + The path to the output HDF5 file. + vs : float or np.ndarray + Shear wave velocity (m/s). A scalar broadcasts to all points, + or provide a per-point array of shape (npoints,). The Version 1.0 + SRF format does not include this field; defaults to 0.0. + den : float or np.ndarray + Density (kg/m^3). A scalar broadcasts to all points, + or provide a per-point array of shape (npoints,). The Version 1.0 + SRF format does not include this field; defaults to 0.0. + include_slip_time_function : bool + If True, include the SR1 slip rate time series dataset. + """ + # Build PLANE structured array (header has one row per fault plane) + plane_data = np.empty(len(self.header), dtype=SW4_PLANE_DTYPE) + for field in SW4_PLANE_DTYPE.names: + plane_data[field] = self.header[field.lower()].values.astype( + SW4_PLANE_DTYPE[field].type + ) + + # Build POINTS structured array + points_data = np.zeros(len(self.points), dtype=SW4_POINTS_DTYPE) + for field in SW4_POINTS_DTYPE.names: + if field in _SW4_POINTS_EXTERNAL_FIELDS: + continue + srf_col = _SW4_POINTS_COL_OVERRIDES.get(field, field.lower()) + points_data[field] = self.points[srf_col].values.astype( + SW4_POINTS_DTYPE[field].type + ) + + points_data["VS"] = np.float32(vs) + points_data["DEN"] = np.float32(den) + + if include_slip_time_function: + points_data["NT1"] = np.diff(self.slipt1_array.indptr).astype(np.int32) + + with h5py.File(output_ffp, "w") as h5file: + h5file.attrs.create("VERSION", np.float32(self.version)) + h5file.attrs.create( + "PLANE", plane_data, plane_data.shape, SW4_PLANE_DTYPE + ) + h5file.create_dataset("POINTS", data=points_data) + + if include_slip_time_function: + h5file.create_dataset( + "SR1", data=self.slipt1_array.data.astype(np.float32) + ) + def write_hdf5( self, hdf5_ffp: Path, include_slip_time_function: bool = True ) -> None: @@ -352,7 +453,9 @@ def from_hdf5(cls, hdf5_ffp: Path) -> Self: points_data = { col: ds[col].values for col in ds.data_vars - if isinstance(col, str) and not col.startswith("plane_") and col not in {"data", "indices", "indptr"} + if isinstance(col, str) + and not col.startswith("plane_") + and col not in {"data", "indices", "indptr"} } points_df = pd.DataFrame(points_data) diff --git a/tests/test_srf.py b/tests/test_srf.py index 847f539..5d65d2b 100644 --- a/tests/test_srf.py +++ b/tests/test_srf.py @@ -503,3 +503,60 @@ def test_hdf5_read_write(): assert (original_srf.slipt1_array != reconstructed_srf.slipt1_array).nnz == 0, ( # type: ignore "slipt1_array content mismatch" ) + + +def test_sw4_hdf5_read_write(): + """Test that write_sw4_hdf5 preserves header, points, and slip data.""" + import h5py + + original_srf = srf.read_srf(SRF_DIR / "3468575.srf") + + with tempfile.NamedTemporaryFile(suffix=".h5") as tmp_file: + output_path = Path(tmp_file.name) + original_srf.write_sw4_hdf5(output_path) + + with h5py.File(output_path, "r") as h5file: + # VERSION + assert h5file.attrs["VERSION"] == np.float32(1.0) + + # PLANE — compare every field against original header + plane = h5file.attrs["PLANE"] + assert plane.shape == (len(original_srf.header),) + for i, (_, row) in enumerate(original_srf.header.iterrows()): + for field in srf.SW4_PLANE_DTYPE.names: + assert plane[i][field] == pytest.approx( + row[field.lower()], abs=1e-3 + ) + + # POINTS — compare each DataFrame-sourced field directly + points = h5file["POINTS"] + assert points.shape == (len(original_srf.points),) + for field in ("LON", "LAT", "DEP", "STK", "DIP", "AREA", "TINIT", "DT", "RAKE"): + np.testing.assert_array_almost_equal( + points[field], + original_srf.points[field.lower()].values, + decimal=3, + ) + np.testing.assert_array_almost_equal( + points["SLIP1"], + original_srf.points["slip"].values, + decimal=3, + ) + + # VS/DEN default to 0.0 for Version 1.0 SRF + np.testing.assert_array_equal(points["VS"], 0.0) + np.testing.assert_array_equal(points["DEN"], 0.0) + + # NT1 from slipt1_array.indptr + expected_nt1 = np.diff(original_srf.slipt1_array.indptr).astype(np.int32) + np.testing.assert_array_equal(points["NT1"], expected_nt1) + + # SR1 slip-time function data + np.testing.assert_array_almost_equal( + h5file["SR1"][...], + original_srf.slipt1_array.data.astype(np.float32), + ) + + # Unused slip components stay zero + for zero_field in ("SLIP2", "NT2", "SLIP3", "NT3"): + assert np.all(points[zero_field] == 0) From ad86fe9dd2e3e20b405e97a6814679d28835e309 Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Mon, 23 Mar 2026 16:11:04 +1300 Subject: [PATCH 02/18] minor code cleanup --- source_modelling/srf.py | 12 ++++-------- tests/test_srf.py | 19 ++++++++++--------- 2 files changed, 14 insertions(+), 17 deletions(-) diff --git a/source_modelling/srf.py b/source_modelling/srf.py index bedbdbf..f77e862 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -99,7 +99,6 @@ ] ) -_SW4_POINTS_COL_OVERRIDES = {"SLIP1": "slip"} _SW4_POINTS_EXTERNAL_FIELDS = {"VS", "DEN", "NT1", "SLIP2", "NT2", "SLIP3", "NT3"} @@ -375,10 +374,9 @@ def write_sw4_hdf5( for field in SW4_POINTS_DTYPE.names: if field in _SW4_POINTS_EXTERNAL_FIELDS: continue - srf_col = _SW4_POINTS_COL_OVERRIDES.get(field, field.lower()) - points_data[field] = self.points[srf_col].values.astype( - SW4_POINTS_DTYPE[field].type - ) + points_data[field] = self.points[ + "slip" if field == "SLIP1" else field.lower() + ].values.astype(SW4_POINTS_DTYPE[field].type) points_data["VS"] = np.float32(vs) points_data["DEN"] = np.float32(den) @@ -388,9 +386,7 @@ def write_sw4_hdf5( with h5py.File(output_ffp, "w") as h5file: h5file.attrs.create("VERSION", np.float32(self.version)) - h5file.attrs.create( - "PLANE", plane_data, plane_data.shape, SW4_PLANE_DTYPE - ) + h5file.attrs.create("PLANE", plane_data) h5file.create_dataset("POINTS", data=points_data) if include_slip_time_function: diff --git a/tests/test_srf.py b/tests/test_srf.py index 5d65d2b..d2a6478 100644 --- a/tests/test_srf.py +++ b/tests/test_srf.py @@ -512,23 +512,22 @@ def test_sw4_hdf5_read_write(): original_srf = srf.read_srf(SRF_DIR / "3468575.srf") with tempfile.NamedTemporaryFile(suffix=".h5") as tmp_file: - output_path = Path(tmp_file.name) - original_srf.write_sw4_hdf5(output_path) + original_srf.write_sw4_hdf5(tmp_file.name) - with h5py.File(output_path, "r") as h5file: + with h5py.File(tmp_file.name, "r") as h5file: # VERSION assert h5file.attrs["VERSION"] == np.float32(1.0) - # PLANE — compare every field against original header + # PLANE -- compare every field against original header plane = h5file.attrs["PLANE"] assert plane.shape == (len(original_srf.header),) - for i, (_, row) in enumerate(original_srf.header.iterrows()): + for idx, row in original_srf.header.iterrows(): for field in srf.SW4_PLANE_DTYPE.names: - assert plane[i][field] == pytest.approx( + assert plane[idx][field] == pytest.approx( row[field.lower()], abs=1e-3 ) - # POINTS — compare each DataFrame-sourced field directly + # POINTS -- compare each DataFrame-sourced field directly points = h5file["POINTS"] assert points.shape == (len(original_srf.points),) for field in ("LON", "LAT", "DEP", "STK", "DIP", "AREA", "TINIT", "DT", "RAKE"): @@ -548,8 +547,10 @@ def test_sw4_hdf5_read_write(): np.testing.assert_array_equal(points["DEN"], 0.0) # NT1 from slipt1_array.indptr - expected_nt1 = np.diff(original_srf.slipt1_array.indptr).astype(np.int32) - np.testing.assert_array_equal(points["NT1"], expected_nt1) + np.testing.assert_array_equal( + points["NT1"], + np.diff(original_srf.slipt1_array.indptr).astype(np.int32), + ) # SR1 slip-time function data np.testing.assert_array_almost_equal( From e213b80b366063b8c94c44e66e927955e2563828 Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Mon, 23 Mar 2026 16:55:10 +1300 Subject: [PATCH 03/18] fix deptry check --- pyproject.toml | 1 + source_modelling/srf.py | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 719e167..6da12f1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,6 +12,7 @@ dynamic = ["version"] dependencies = [ "fiona", "geopandas", + "h5py", "networkx", "numpy", "pandas", diff --git a/source_modelling/srf.py b/source_modelling/srf.py index f77e862..a701196 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -47,7 +47,6 @@ from pathlib import Path from typing import Self -import h5py import numpy as np import pandas as pd import scipy as sp @@ -362,6 +361,8 @@ def write_sw4_hdf5( include_slip_time_function : bool If True, include the SR1 slip rate time series dataset. """ + import h5py + # Build PLANE structured array (header has one row per fault plane) plane_data = np.empty(len(self.header), dtype=SW4_PLANE_DTYPE) for field in SW4_PLANE_DTYPE.names: From 0c6285db7d0ea3b02f138ad806c1dc1684d45a1d Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Mon, 23 Mar 2026 17:19:42 +1300 Subject: [PATCH 04/18] fix(srf): use np.asarray for vs/den to support array inputs np.float32() fails on arrays, but the docstring allows np.ndarray. Co-Authored-By: Claude Opus 4.6 --- source_modelling/srf.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source_modelling/srf.py b/source_modelling/srf.py index a701196..430e7e0 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -379,8 +379,8 @@ def write_sw4_hdf5( "slip" if field == "SLIP1" else field.lower() ].values.astype(SW4_POINTS_DTYPE[field].type) - points_data["VS"] = np.float32(vs) - points_data["DEN"] = np.float32(den) + points_data["VS"] = np.asarray(vs, dtype=np.float32) + points_data["DEN"] = np.asarray(den, dtype=np.float32) if include_slip_time_function: points_data["NT1"] = np.diff(self.slipt1_array.indptr).astype(np.int32) From 1cbc393333fe9ca0ba9b4a3639f37f1f92a932a1 Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Mon, 23 Mar 2026 17:40:21 +1300 Subject: [PATCH 05/18] pass Path to write_sw4_hdf5 in test --- tests/test_srf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_srf.py b/tests/test_srf.py index d2a6478..82ab6c3 100644 --- a/tests/test_srf.py +++ b/tests/test_srf.py @@ -512,7 +512,7 @@ def test_sw4_hdf5_read_write(): original_srf = srf.read_srf(SRF_DIR / "3468575.srf") with tempfile.NamedTemporaryFile(suffix=".h5") as tmp_file: - original_srf.write_sw4_hdf5(tmp_file.name) + original_srf.write_sw4_hdf5(Path(tmp_file.name)) with h5py.File(tmp_file.name, "r") as h5file: # VERSION From 6e4a12142e21b9bc45d42625badfa17b09807872 Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Tue, 24 Mar 2026 14:42:56 +1300 Subject: [PATCH 06/18] fix ty type issues --- pyproject.toml | 3 +++ source_modelling/community_fault_model.py | 2 +- source_modelling/moment.py | 2 +- source_modelling/sources.py | 21 ++++++++++++--------- source_modelling/srf.py | 14 +++++++++++--- tests/test_gsf.py | 2 +- tests/test_magnitude_scaling.py | 6 +++--- tests/test_rupture_propagation.py | 7 ++++--- tests/test_sources.py | 8 ++++---- tests/test_srf.py | 21 +++++++++++++++++---- 10 files changed, 57 insertions(+), 29 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 6da12f1..02dcad3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -47,6 +47,9 @@ source_modelling = ['NZ_CFM/*'] [tool.setuptools.package-dir] source_modelling = "source_modelling" +[tool.ty.src] +exclude = ["source_modelling/ccldpy.py"] + [tool.ruff] extend-exclude = ["source_modelling/ccldpy.py"] diff --git a/source_modelling/community_fault_model.py b/source_modelling/community_fault_model.py index c1698f0..1dd6f3d 100644 --- a/source_modelling/community_fault_model.py +++ b/source_modelling/community_fault_model.py @@ -372,7 +372,7 @@ def community_fault_model_as_geodataframe() -> gpd.GeoDataFrame: # Transform the trace to WGS84 for fault in model: fault.trace = shapely.transform(fault.trace, lambda coord: coord[:, ::-1]) - return gpd.GeoDataFrame( + return gpd.GeoDataFrame( # type: ignore # set_index returns pd.DataFrame which ty flags as not being a GeoDataFrame [vars(fault) for fault in model], geometry="trace", crs="EPSG:2193" ).set_index("name") diff --git a/source_modelling/moment.py b/source_modelling/moment.py index 2f533b8..49e59d3 100644 --- a/source_modelling/moment.py +++ b/source_modelling/moment.py @@ -52,7 +52,7 @@ def find_connected_faults( """ fault_names = list(faults) - fault_components: DisjointSet[str] = DisjointSet(fault_names) # type: ignore + fault_components: DisjointSet[str] = DisjointSet(fault_names) for fault_a_name, fault_b_name in itertools.product(fault_names, repeat=2): if not fault_b_name or fault_components.connected(fault_a_name, fault_b_name): continue diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 3af5e8d..a295ffe 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -21,6 +21,7 @@ import itertools import json import warnings +from collections.abc import Sequence from typing import NamedTuple, Self import networkx as nx @@ -110,8 +111,8 @@ def geometry(self) -> shapely.Point: # numpydoc ignore=RT01 return shapely.Point(self.bounds) @property - def geojson(self) -> dict: # numpydoc ignore=RT01 - """dict: A GeoJSON representation of the fault.""" + def geojson(self) -> str: # numpydoc ignore=RT01 + """str: A GeoJSON representation of the fault.""" return shapely.to_geojson( shapely.transform( self.geometry, @@ -438,8 +439,8 @@ def trace_geometry(self) -> shapely.LineString: # numpydoc ignore=RT01 return shapely.LineString(self.trace) @property - def geojson(self) -> dict: # numpydoc ignore=RT01 - """dict: A GeoJSON representation of the fault.""" + def geojson(self) -> str: # numpydoc ignore=RT01 + """str: A GeoJSON representation of the fault.""" return shapely.to_geojson( shapely.transform( self.geometry, @@ -450,7 +451,7 @@ def geojson(self) -> dict: # numpydoc ignore=RT01 @classmethod def from_nztm_trace( cls, - trace_points_nztm: npt.NDArray[float], + trace_points_nztm: npt.NDArray[np.float64], dtop: float, dbottom: float, dip: float, @@ -523,6 +524,8 @@ def from_nztm_trace( (trace_points_nztm, np.array([dbottom, dbottom])) ) else: + # Non-None guaranteed by validation above; assert narrows type for ty. + assert dip_dir_nztm is not None dip_dir_nztm_rad = np.deg2rad(dip_dir_nztm) proj_width = (dbottom - dtop) / np.tan(np.deg2rad(dip)) @@ -1247,8 +1250,8 @@ def wgs_depth_coordinates_to_fault_coordinates( raise ValueError("Given coordinates are not on fault.") @property - def geojson(self) -> dict: # numpydoc ignore=RT01 - """dict: A GeoJSON representation of the fault.""" + def geojson(self) -> str: # numpydoc ignore=RT01 + """str: A GeoJSON representation of the fault.""" return shapely.to_geojson( shapely.transform( self.geometry, @@ -1334,12 +1337,12 @@ def rjb_distance(self, point: np.ndarray) -> float: IsSource = Plane | Fault | Point -def sources_as_geojson_features(sources: list[IsSource]) -> str: +def sources_as_geojson_features(sources: Sequence[IsSource]) -> str: """Convert a list of sources to a GeoJSON FeatureCollection. Parameters ---------- - sources : list[IsSource] + sources : Sequence[IsSource] The sources to convert. Returns diff --git a/source_modelling/srf.py b/source_modelling/srf.py index 430e7e0..f46b3ad 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -125,7 +125,9 @@ def __init__(self, header: pd.DataFrame, points: pd.DataFrame) -> None: self._header = header self._points = points - def __getitem__(self, index: int) -> pd.DataFrame: + # ty: slice overload missing to satisfy Sequence LSP; fix by adding + # @overload stubs for int and slice once slice support is implemented. + def __getitem__(self, index: int) -> pd.DataFrame: # ty: ignore[invalid-method-override] """Get the nth segment in the SRF. Parameters @@ -302,7 +304,7 @@ def from_file(cls, srf_ffp: Path | str) -> Self: slipt1_array, ) - def write_srf(self, srf_ffp: Path) -> None: + def write_srf(self, srf_ffp: str | Path) -> None: """Write an SRFFile object to a file. Parameters @@ -365,6 +367,9 @@ def write_sw4_hdf5( # Build PLANE structured array (header has one row per fault plane) plane_data = np.empty(len(self.header), dtype=SW4_PLANE_DTYPE) + assert ( + SW4_PLANE_DTYPE.names is not None + ) # always has names as it is a structured dtype, but ty needs a type guard for field in SW4_PLANE_DTYPE.names: plane_data[field] = self.header[field.lower()].values.astype( SW4_PLANE_DTYPE[field].type @@ -372,6 +377,9 @@ def write_sw4_hdf5( # Build POINTS structured array points_data = np.zeros(len(self.points), dtype=SW4_POINTS_DTYPE) + assert ( + SW4_POINTS_DTYPE.names is not None + ) # always has names as it is a structured dtype, but ty needs a type guard for field in SW4_POINTS_DTYPE.names: if field in _SW4_POINTS_EXTERNAL_FIELDS: continue @@ -694,7 +702,7 @@ def read_srf(srf_ffp: Path | str) -> SrfFile: return SrfFile.from_file(srf_ffp) -def write_srf(srf_ffp: Path, srf: SrfFile) -> None: +def write_srf(srf_ffp: str | Path, srf: SrfFile) -> None: """Write an SRF object to a filepath. Parameters diff --git a/tests/test_gsf.py b/tests/test_gsf.py index de9b04e..f15a0e3 100644 --- a/tests/test_gsf.py +++ b/tests/test_gsf.py @@ -56,7 +56,7 @@ def test_plane_gsf(): # Use tmp_path fixture def test_bad_gsf_type(): with pytest.raises(TypeError): - gsf.source_to_gsf_dataframe(1, 0.1) + gsf.source_to_gsf_dataframe(1, 0.1) # ty: ignore[invalid-argument-type] - intentionally passing wrong type to test error handling def test_write_gsf(tmp_path: Path): diff --git a/tests/test_magnitude_scaling.py b/tests/test_magnitude_scaling.py index 58d6041..9c252b5 100644 --- a/tests/test_magnitude_scaling.py +++ b/tests/test_magnitude_scaling.py @@ -86,7 +86,7 @@ def sampler( magnitude = draw(st.floats(min_value=min_magnitude, max_value=max_magnitude)) return scaling_relation, rake, magnitude - return sampler() + return sampler() # ty: ignore[missing-argument] - @st.composite strips the `draw` parameter; sampler() is valid with no args # The contreras slab 2020 relation (which is the same as Strasser 2010) is not invertible. See the coefficients of Table 2 in the paper @@ -116,7 +116,7 @@ def test_inversion( else: mag_to_area = MAGNITUDE_TO_AREA[scaling_relation] area_to_mag = AREA_TO_MAGNITUDE[scaling_relation] - assert area_to_mag(mag_to_area(magnitude)) == pytest.approx(magnitude) + assert area_to_mag(mag_to_area(magnitude)) == pytest.approx(magnitude) # ty: ignore[missing-argument] - rake is handled by functools.partial above for relations that require it @pytest.mark.parametrize( @@ -245,7 +245,7 @@ def test_monotonicity_mag_to_area( mag_to_area = functools.partial(MAGNITUDE_TO_AREA[scaling_relation], rake=rake) else: mag_to_area = MAGNITUDE_TO_AREA[scaling_relation] - assert mag_to_area(magnitude + 0.1) > mag_to_area(magnitude) + assert mag_to_area(magnitude + 0.1) > mag_to_area(magnitude) # ty: ignore[missing-argument] - rake is handled by functools.partial above for relations that require it class RandomFunction(Protocol): diff --git a/tests/test_rupture_propagation.py b/tests/test_rupture_propagation.py index ca1321f..2efc693 100644 --- a/tests/test_rupture_propagation.py +++ b/tests/test_rupture_propagation.py @@ -11,7 +11,7 @@ from source_modelling import rupture_propagation, sources -random_strongly_connected_graph = hnx.graph_builder( +random_strongly_connected_graph = hnx.graph_builder( # ty: ignore[missing-argument] - hypothesis_networkx has no stubs; ty can't resolve ParamSpec through the decorator graph_type=nx.Graph, node_keys=st.text( alphabet=st.characters( @@ -450,7 +450,7 @@ def test_sample_rupture_propagation( n_samples = 100 sampled_trees = [ rupture_propagation.sample_rupture_propagation( - sources_map, + sources_map, # ty: ignore[invalid-argument-type] - dict invariance: Point satisfies IsSource but dict[str, Point] is not assignable to dict[str, IsSource] initial_source=initial_source, initial_source_distribution=initial_source_distribution, jump_impossibility_limit_distance=jump_impossibility_limit_distance, @@ -566,7 +566,8 @@ def test_jump_points_from_rupture_tree( expected_jump_points: dict[str, rupture_propagation.JumpPair], ): result_jump_points = rupture_propagation.jump_points_from_rupture_tree( - source_map, rupture_causality_tree + source_map, # ty: ignore[invalid-argument-type] - dict invariance: Point satisfies IsSource but dict[str, Point] is not assignable to dict[str, IsSource] + rupture_causality_tree, ) # Check if the jump points match the expected values diff --git a/tests/test_sources.py b/tests/test_sources.py index b0c0b05..d5888a9 100644 --- a/tests/test_sources.py +++ b/tests/test_sources.py @@ -357,7 +357,7 @@ def test_general_invalid_input(): def trace( - start_trace_nztm: npt.NDArray[float], length: float, strike: float + start_trace_nztm: npt.NDArray[np.float64], length: float, strike: float ) -> np.ndarray: # Do this in NZTM to prevent any issues with the coordinate system conversions strike_vec = np.array([np.cos(np.radians(strike)), np.sin(np.radians(strike))]) @@ -418,7 +418,7 @@ def valid_trace_definition(draw: st.DrawFn): ) -@given(valid_trace_definition()) +@given(valid_trace_definition()) # ty: ignore[missing-argument] - @st.composite strips the `draw` parameter; valid_trace_definition() is valid with no args def test_plane_from_trace(data: tuple): ( trace_points_nztm, @@ -1227,7 +1227,7 @@ def fault( ) -@given(fault(min_segments=2, max_segments=10, min_length=0.4, max_length=100)) +@given(fault(min_segments=2, max_segments=10, min_length=0.4, max_length=100)) # ty: ignore[missing-argument] - @st.composite strips the `draw` parameter; fault() is valid with no positional args def test_simplify_fault(fault: Fault): tolerance = 0.4 simplified_fault = sources.simplify_fault(fault, tolerance) @@ -1253,7 +1253,7 @@ def test_simplify_fault(fault: Fault): @given( st.lists( - fault(min_segments=2, max_segments=10, min_length=0.4, max_length=100), + fault(min_segments=2, max_segments=10, min_length=0.4, max_length=100), # ty: ignore[missing-argument] - @st.composite strips the `draw` parameter; fault() is valid with no positional args min_size=2, max_size=10, ), diff --git a/tests/test_srf.py b/tests/test_srf.py index 82ab6c3..1d764d7 100644 --- a/tests/test_srf.py +++ b/tests/test_srf.py @@ -407,7 +407,7 @@ def test_planes_nstk_1_ndip_gt_1(): version="1.0", header=header, points=points, - slipt1_array=None, # type: ignore + slipt1_array=None, ) planes = mock_srf.planes @@ -450,7 +450,7 @@ def test_planes_nstk_1_ndip_1(): version="1.0", header=header, points=points, - slipt1_array=None, # type: ignore + slipt1_array=None, ) planes = mock_srf.planes @@ -500,7 +500,7 @@ def test_hdf5_read_write(): original_srf.slipt1_array.indptr, reconstructed_srf.slipt1_array.indptr ), "slipt1_array indptr mismatch" - assert (original_srf.slipt1_array != reconstructed_srf.slipt1_array).nnz == 0, ( # type: ignore + assert (original_srf.slipt1_array != reconstructed_srf.slipt1_array).nnz == 0, ( "slipt1_array content mismatch" ) @@ -521,6 +521,9 @@ def test_sw4_hdf5_read_write(): # PLANE -- compare every field against original header plane = h5file.attrs["PLANE"] assert plane.shape == (len(original_srf.header),) + assert ( + srf.SW4_PLANE_DTYPE.names is not None + ) # always has names as it is a structured dtype, but ty needs a type guard for idx, row in original_srf.header.iterrows(): for field in srf.SW4_PLANE_DTYPE.names: assert plane[idx][field] == pytest.approx( @@ -530,7 +533,17 @@ def test_sw4_hdf5_read_write(): # POINTS -- compare each DataFrame-sourced field directly points = h5file["POINTS"] assert points.shape == (len(original_srf.points),) - for field in ("LON", "LAT", "DEP", "STK", "DIP", "AREA", "TINIT", "DT", "RAKE"): + for field in ( + "LON", + "LAT", + "DEP", + "STK", + "DIP", + "AREA", + "TINIT", + "DT", + "RAKE", + ): np.testing.assert_array_almost_equal( points[field], original_srf.points[field.lower()].values, From 973a027959aafa10321cbc242e5e83bd97b3e7a4 Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Tue, 24 Mar 2026 15:09:43 +1300 Subject: [PATCH 07/18] ty typing fix attempt 2 --- source_modelling/community_fault_model.py | 2 +- source_modelling/rupture_propagation.py | 26 +++++++++++++++-------- source_modelling/sources.py | 5 +++-- tests/test_gsf.py | 4 ++-- tests/test_rupture_propagation.py | 2 +- tests/test_srf.py | 10 ++++----- 6 files changed, 29 insertions(+), 20 deletions(-) diff --git a/source_modelling/community_fault_model.py b/source_modelling/community_fault_model.py index 1dd6f3d..c1698f0 100644 --- a/source_modelling/community_fault_model.py +++ b/source_modelling/community_fault_model.py @@ -372,7 +372,7 @@ def community_fault_model_as_geodataframe() -> gpd.GeoDataFrame: # Transform the trace to WGS84 for fault in model: fault.trace = shapely.transform(fault.trace, lambda coord: coord[:, ::-1]) - return gpd.GeoDataFrame( # type: ignore # set_index returns pd.DataFrame which ty flags as not being a GeoDataFrame + return gpd.GeoDataFrame( [vars(fault) for fault in model], geometry="trace", crs="EPSG:2193" ).set_index("name") diff --git a/source_modelling/rupture_propagation.py b/source_modelling/rupture_propagation.py index 09b9cdc..028bbcc 100644 --- a/source_modelling/rupture_propagation.py +++ b/source_modelling/rupture_propagation.py @@ -16,7 +16,7 @@ import warnings from collections import defaultdict, namedtuple from collections.abc import Generator -from typing import Literal +from typing import Literal, overload import networkx as nx import numpy as np @@ -59,7 +59,7 @@ def spanning_tree_with_probabilities( trees = [] probabilities = [] - for tree in mst.SpanningTreeIterator(graph): + for tree in mst.SpanningTreeIterator(graph): # ty: ignore[invalid-argument-type] - SpanningTreeIterator accepts undirected Graph at runtime; the networkx-stubs incorrectly restrict its __init__ to DiGraph only p_tree = 1.0 for u, v in graph.edges: if tree.has_edge(u, v): @@ -72,6 +72,12 @@ def spanning_tree_with_probabilities( return trees, probabilities +@overload +def sampled_spanning_tree(graph: nx.Graph, n_samples: Literal[1] = ...) -> nx.Graph: ... +@overload +def sampled_spanning_tree( + graph: nx.Graph, n_samples: int +) -> list[nx.Graph] | nx.Graph: ... def sampled_spanning_tree( graph: nx.Graph, n_samples: int = 1 ) -> list[nx.Graph] | nx.Graph: @@ -209,7 +215,7 @@ def select_top_spanning_trees( cumulative_tree_weight = 0.0 spanning_trees = [] - for spanning_tree in mst.SpanningTreeIterator(weighted_graph, minimum=False): + for spanning_tree in mst.SpanningTreeIterator(weighted_graph, minimum=False): # ty: ignore[invalid-argument-type] - SpanningTreeIterator accepts undirected Graph at runtime; the networkx-stubs incorrectly restrict its __init__ to DiGraph only spanning_trees.append(spanning_tree) tree_log_probability = sum( spanning_tree[node_u][node_v]["weight"] @@ -320,7 +326,7 @@ def prune_distance_graph(distances: DistanceGraph, cutoff: float) -> DistanceGra def probability_graph( distances: DistanceGraph, d0: float = 3, delta: float = 1 -) -> nx.DiGraph: +) -> nx.Graph: """ Convert a distance graph into a probability graph. @@ -339,8 +345,8 @@ def probability_graph( Returns ------- - nx.DiGraph - A directed graph where edges are weighted by probabilities of rupture + nx.Graph + An undirected graph where edges are weighted by probabilities of rupture propagation between faults. """ @@ -390,9 +396,11 @@ def distance_between( """ global_point_a = source_a.fault_coordinates_to_wgs_depth_coordinates(source_a_point) global_point_b = source_b.fault_coordinates_to_wgs_depth_coordinates(source_b_point) - return float(coordinates.distance_between_wgs_depth_coordinates( - global_point_a, global_point_b - )) + return float( + coordinates.distance_between_wgs_depth_coordinates( + global_point_a, global_point_b + ) + ) def sample_rupture_propagation( diff --git a/source_modelling/sources.py b/source_modelling/sources.py index a295ffe..2a0344b 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -996,7 +996,8 @@ def __post_init__(self) -> None: # This relation can now be used to identify if the list of planes given is a line. points_into_graph: nx.DiGraph = nx.from_dict_of_lists( - points_into_relation, create_using=nx.DiGraph + points_into_relation, # ty: ignore[invalid-argument-type] - dict[int, list[int]] satisfies dict[Unknown, Iterable[Unknown]]; ty can't resolve types through networkx's _dispatchable wrapper + create_using=nx.DiGraph, ) try: self._validate_fault_plane_connectivity(points_into_graph) @@ -1174,7 +1175,7 @@ def centroid(self) -> np.ndarray: # numpydoc ignore=RT01 return self.fault_coordinates_to_wgs_depth_coordinates(np.array([1 / 2, 1 / 2])) @property - def geometry(self) -> shapely.Geometry: # numpydoc ignore=RT01 + def geometry(self) -> shapely.Polygon | shapely.LineString: # numpydoc ignore=RT01 """shapely.Polygon or LineString: A shapely geometry for the fault (projected onto the surface). Geometry will be LineString if `dip = 90`. diff --git a/tests/test_gsf.py b/tests/test_gsf.py index f15a0e3..2376a4c 100644 --- a/tests/test_gsf.py +++ b/tests/test_gsf.py @@ -49,7 +49,7 @@ def test_plane_gsf(): # Use tmp_path fixture for _, point in gsf_df.iterrows(): assert plane.geometry.contains( shapely.Point( - coordinates.wgs_depth_to_nztm(point[["lat", "lon", "dep"]].values) + coordinates.wgs_depth_to_nztm(point[["lat", "lon", "dep"]].to_numpy()) ) ) @@ -304,7 +304,7 @@ def test_fault_to_gsf(fault: Fault): assert shapely.contains( fault.geometry, shapely.Point( - coordinates.wgs_depth_to_nztm(point[["lat", "lon", "dep"]].values) + coordinates.wgs_depth_to_nztm(point[["lat", "lon", "dep"]].to_numpy()) ), ) diff --git a/tests/test_rupture_propagation.py b/tests/test_rupture_propagation.py index 2efc693..20562c4 100644 --- a/tests/test_rupture_propagation.py +++ b/tests/test_rupture_propagation.py @@ -154,7 +154,7 @@ def test_random_sampling_root(graph: nx.Graph, n: int): random.seed(1) probability_of_spanning_trees = 0.0 tree_probabilities: dict[str, float] = {} - for tree in mst.SpanningTreeIterator(graph): + for tree in mst.SpanningTreeIterator(graph): # ty: ignore[invalid-argument-type] - SpanningTreeIterator accepts undirected Graph at runtime; the networkx-stubs incorrectly restrict its __init__ to DiGraph only p_tree = 1.0 for u, v in graph.edges: if tree.has_edge(u, v): diff --git a/tests/test_srf.py b/tests/test_srf.py index 1d764d7..b6e8d9c 100644 --- a/tests/test_srf.py +++ b/tests/test_srf.py @@ -407,7 +407,7 @@ def test_planes_nstk_1_ndip_gt_1(): version="1.0", header=header, points=points, - slipt1_array=None, + slipt1_array=None, # ty: ignore[invalid-argument-type] - intentionally passing None to construct a minimal mock SrfFile that doesn't need slip time functions ) planes = mock_srf.planes @@ -450,7 +450,7 @@ def test_planes_nstk_1_ndip_1(): version="1.0", header=header, points=points, - slipt1_array=None, + slipt1_array=None, # ty: ignore[invalid-argument-type] - intentionally passing None to construct a minimal mock SrfFile that doesn't need slip time functions ) planes = mock_srf.planes @@ -500,7 +500,7 @@ def test_hdf5_read_write(): original_srf.slipt1_array.indptr, reconstructed_srf.slipt1_array.indptr ), "slipt1_array indptr mismatch" - assert (original_srf.slipt1_array != reconstructed_srf.slipt1_array).nnz == 0, ( + assert (original_srf.slipt1_array != reconstructed_srf.slipt1_array).nnz == 0, ( # ty: ignore[unresolved-attribute] - csr_array.__ne__ returns a sparse array at runtime; stubs incorrectly type it as bool "slipt1_array content mismatch" ) @@ -546,12 +546,12 @@ def test_sw4_hdf5_read_write(): ): np.testing.assert_array_almost_equal( points[field], - original_srf.points[field.lower()].values, + original_srf.points[field.lower()].to_numpy(), decimal=3, ) np.testing.assert_array_almost_equal( points["SLIP1"], - original_srf.points["slip"].values, + original_srf.points["slip"].to_numpy(), decimal=3, ) From a77bbc5f45abb1dd875493758b9a29182b3ae548 Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Tue, 24 Mar 2026 15:13:40 +1300 Subject: [PATCH 08/18] ty type check fix attempt 3 --- source_modelling/sources.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 2a0344b..d0961dc 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -1180,7 +1180,7 @@ def geometry(self) -> shapely.Polygon | shapely.LineString: # numpydoc ignore=R Geometry will be LineString if `dip = 90`. """ - return shapely.normalize( + return shapely.normalize( # ty: ignore[invalid-return-type] - shapely.normalize and union_all are typed to return BaseGeometry, but the result is always Polygon or LineString for plane geometries shapely.union_all([plane.geometry for plane in self.planes]) ) From defb0d3e6cf2c5a0468bb219e52053814bb16a43 Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Tue, 24 Mar 2026 15:17:16 +1300 Subject: [PATCH 09/18] numpydoc-lint fix --- source_modelling/rupture_propagation.py | 10 ++++++++-- source_modelling/srf.py | 4 ++-- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/source_modelling/rupture_propagation.py b/source_modelling/rupture_propagation.py index 028bbcc..7f165c2 100644 --- a/source_modelling/rupture_propagation.py +++ b/source_modelling/rupture_propagation.py @@ -73,11 +73,17 @@ def spanning_tree_with_probabilities( @overload -def sampled_spanning_tree(graph: nx.Graph, n_samples: Literal[1] = ...) -> nx.Graph: ... +def sampled_spanning_tree(graph: nx.Graph, n_samples: Literal[1] = ...) -> nx.Graph: + """Sample a single spanning tree from a graph.""" + ... + + @overload def sampled_spanning_tree( graph: nx.Graph, n_samples: int -) -> list[nx.Graph] | nx.Graph: ... +) -> list[nx.Graph] | nx.Graph: + """Sample one or more spanning trees from a graph.""" + ... def sampled_spanning_tree( graph: nx.Graph, n_samples: int = 1 ) -> list[nx.Graph] | nx.Graph: diff --git a/source_modelling/srf.py b/source_modelling/srf.py index f46b3ad..cfb19f3 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -5,8 +5,8 @@ See https://wiki.canterbury.ac.nz/display/QuakeCore/File+Formats+Used+On+GM for details on the SRF format. -Why Not qcore.srf? ------------------- +**Why not qcore.srf?** + You might use this module instead of the `qcore.srf` module because: 1. The `qcore.srf` module does not support writing SRF files. From b3a26bdb8b76b6498827625a5af261b97cf8e8fd Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Tue, 24 Mar 2026 15:20:53 +1300 Subject: [PATCH 10/18] numpydoc-lint fix attempt 2 --- source_modelling/fsp.py | 13 ++++------- source_modelling/gsf.py | 11 +++------ source_modelling/rupture_propagation.py | 30 +++++++++++++++++++------ source_modelling/sources.py | 12 ++-------- source_modelling/srf.py | 14 +++++------- 5 files changed, 37 insertions(+), 43 deletions(-) diff --git a/source_modelling/fsp.py b/source_modelling/fsp.py index 0567f5a..66d2067 100644 --- a/source_modelling/fsp.py +++ b/source_modelling/fsp.py @@ -5,16 +5,11 @@ See http://equake-rc.info/SRCMOD/fileformats/fsp/ for details on the FSP format. -Classes -------- -- FSPFile: Representation of an FSP file. +Classes: ``FSPFile`` (representation of an FSP file). +Exceptions: ``FSPParseError`` (raised for errors in parsing FSP files). -Exceptions ----------- -- FSPParseError: Exception raised for errors in parsing FSP files. - -Example -------- +Examples +-------- >>> fsp_file = FSPFile.read_from_file(fsp_ffp) >>> (fsp_file.data['trup'] + fsp_file.data['rise']).max() # Get time of final rise for subfaults. """ diff --git a/source_modelling/gsf.py b/source_modelling/gsf.py index 069e786..7dcd260 100644 --- a/source_modelling/gsf.py +++ b/source_modelling/gsf.py @@ -1,13 +1,8 @@ """This module provides functions for working with and generating GSF files. -Functions ---------- -source_to_gsf_dataframe(gsf_filepath, source, resolution) - Generates a pandas DataFrame suitable for writing to a GSF file from a source object. -write_gsf(gsf_df, gsf_filepath) - Writes a pandas DataFrame to a GSF file. -read_gsf(gsf_filepath) - Parses a GSF file into a pandas DataFrame. +Functions: ``source_to_gsf_dataframe`` (generate a DataFrame from a source object), +``write_gsf`` (write a DataFrame to a GSF file), ``read_gsf`` (parse a GSF file into +a DataFrame). References ---------- diff --git a/source_modelling/rupture_propagation.py b/source_modelling/rupture_propagation.py index 7f165c2..2385d80 100644 --- a/source_modelling/rupture_propagation.py +++ b/source_modelling/rupture_propagation.py @@ -4,8 +4,8 @@ This module provides functions for computing likely rupture paths from information about the distances between faults. -Reference ---------- +References +---------- To understand the purpose and implementation of the algorithms in the 'Rupture Propagation' page[0] on the source modelling wiki. @@ -74,16 +74,32 @@ def spanning_tree_with_probabilities( @overload def sampled_spanning_tree(graph: nx.Graph, n_samples: Literal[1] = ...) -> nx.Graph: - """Sample a single spanning tree from a graph.""" + """Sample a single spanning tree from a graph. + + Parameters + ---------- + graph : nx.Graph + The graph to sample from. + n_samples : Literal[1], optional + Number of trees to sample. Default is 1. + """ ... @overload -def sampled_spanning_tree( - graph: nx.Graph, n_samples: int -) -> list[nx.Graph] | nx.Graph: - """Sample one or more spanning trees from a graph.""" +def sampled_spanning_tree(graph: nx.Graph, n_samples: int) -> list[nx.Graph] | nx.Graph: + """Sample one or more spanning trees from a graph. + + Parameters + ---------- + graph : nx.Graph + The graph to sample from. + n_samples : int + Number of trees to sample. + """ ... + + def sampled_spanning_tree( graph: nx.Graph, n_samples: int = 1 ) -> list[nx.Graph] | nx.Graph: diff --git a/source_modelling/sources.py b/source_modelling/sources.py index d0961dc..33fcc0a 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -4,16 +4,8 @@ faults, along with methods for calculating various properties such as dimensions, orientation, and coordinate transformations. -Classes -------- -Point: - A representation of a point source. - -Plane: - A representation of a single plane of a Fault. - -Fault: - A representation of a fault, consisting of one or more Planes. +Classes: ``Point`` (a point source), ``Plane`` (a single fault plane), +``Fault`` (a fault consisting of one or more planes). """ import copy diff --git a/source_modelling/srf.py b/source_modelling/srf.py index cfb19f3..84ab711 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -21,17 +21,13 @@ points of the SRF file (it is memory efficient), or you are working with code that already uses `qcore.srf`. -Classes -------- -- SrfFile: Representation of an SRF file. +Classes: ``SrfFile`` (representation of an SRF file). -Functions ---------- -- read_srf: Read an SRF file into memory. -- write_srf: Write an SRF object to a filepath. +Functions: ``read_srf`` (read an SRF file into memory), ``write_srf`` (write an SRF +object to a filepath). -Example -------- +Examples +-------- >>> srf_file = srf.read_srf('/path/to/srf') >>> srf_file.points['tinit'].max() # get the last time any point in the SRF ruptures >>> srf_file.points['tinit'] += 1 # delay all points by one second From 4dce64030a9224258ffd3dd4d5332dfc504b7bae Mon Sep 17 00:00:00 2001 From: AndrewRidden-Harper <52001209+AndrewRidden-Harper@users.noreply.github.com> Date: Sat, 28 Mar 2026 16:17:55 +1300 Subject: [PATCH 11/18] del comments Co-authored-by: Jake Faulkner --- source_modelling/srf.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/source_modelling/srf.py b/source_modelling/srf.py index 84ab711..fbc9004 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -361,11 +361,10 @@ def write_sw4_hdf5( """ import h5py - # Build PLANE structured array (header has one row per fault plane) plane_data = np.empty(len(self.header), dtype=SW4_PLANE_DTYPE) assert ( SW4_PLANE_DTYPE.names is not None - ) # always has names as it is a structured dtype, but ty needs a type guard + ) for field in SW4_PLANE_DTYPE.names: plane_data[field] = self.header[field.lower()].values.astype( SW4_PLANE_DTYPE[field].type From e08fb131efaaf769e43a27c67e4f8c9450097a5a Mon Sep 17 00:00:00 2001 From: AndrewRidden-Harper <52001209+AndrewRidden-Harper@users.noreply.github.com> Date: Sat, 28 Mar 2026 16:18:50 +1300 Subject: [PATCH 12/18] del comments Co-authored-by: Jake Faulkner --- source_modelling/srf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source_modelling/srf.py b/source_modelling/srf.py index fbc9004..5b8bfa4 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -374,7 +374,7 @@ def write_sw4_hdf5( points_data = np.zeros(len(self.points), dtype=SW4_POINTS_DTYPE) assert ( SW4_POINTS_DTYPE.names is not None - ) # always has names as it is a structured dtype, but ty needs a type guard + ) for field in SW4_POINTS_DTYPE.names: if field in _SW4_POINTS_EXTERNAL_FIELDS: continue From a70b61c2deb06551e29c198947bfbbcc0c8d1b14 Mon Sep 17 00:00:00 2001 From: AndrewRidden-Harper <52001209+AndrewRidden-Harper@users.noreply.github.com> Date: Sat, 28 Mar 2026 16:20:00 +1300 Subject: [PATCH 13/18] del comment Co-authored-by: Jake Faulkner --- tests/test_srf.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_srf.py b/tests/test_srf.py index b6e8d9c..fda9ab9 100644 --- a/tests/test_srf.py +++ b/tests/test_srf.py @@ -518,7 +518,6 @@ def test_sw4_hdf5_read_write(): # VERSION assert h5file.attrs["VERSION"] == np.float32(1.0) - # PLANE -- compare every field against original header plane = h5file.attrs["PLANE"] assert plane.shape == (len(original_srf.header),) assert ( From ca0c46d39dfe17c0d9adba524742cc7c71030274 Mon Sep 17 00:00:00 2001 From: AndrewRidden-Harper <52001209+AndrewRidden-Harper@users.noreply.github.com> Date: Sat, 28 Mar 2026 16:20:48 +1300 Subject: [PATCH 14/18] del comment Co-authored-by: Jake Faulkner --- tests/test_srf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_srf.py b/tests/test_srf.py index fda9ab9..9e4cb38 100644 --- a/tests/test_srf.py +++ b/tests/test_srf.py @@ -522,7 +522,7 @@ def test_sw4_hdf5_read_write(): assert plane.shape == (len(original_srf.header),) assert ( srf.SW4_PLANE_DTYPE.names is not None - ) # always has names as it is a structured dtype, but ty needs a type guard + ) for idx, row in original_srf.header.iterrows(): for field in srf.SW4_PLANE_DTYPE.names: assert plane[idx][field] == pytest.approx( From 4652ea63e44bba0bc9a89e61a1f4e034fd2e0a71 Mon Sep 17 00:00:00 2001 From: AndrewRidden-Harper <52001209+AndrewRidden-Harper@users.noreply.github.com> Date: Sat, 28 Mar 2026 16:21:33 +1300 Subject: [PATCH 15/18] output_ffp: Path | str Co-authored-by: Jake Faulkner --- source_modelling/srf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source_modelling/srf.py b/source_modelling/srf.py index 5b8bfa4..66acefd 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -337,7 +337,7 @@ def write_srf(self, srf_ffp: str | Path) -> None: def write_sw4_hdf5( self, - output_ffp: Path, + output_ffp: Path | str, vs: float | np.ndarray = 0.0, den: float | np.ndarray = 0.0, include_slip_time_function: bool = True, From 9d1dff4e9ffe41f4305c3619d8c2be6302a3a1b0 Mon Sep 17 00:00:00 2001 From: AndrewRidden-Harper <52001209+AndrewRidden-Harper@users.noreply.github.com> Date: Sat, 28 Mar 2026 16:23:22 +1300 Subject: [PATCH 16/18] del comment Co-authored-by: Jake Faulkner --- tests/test_srf.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_srf.py b/tests/test_srf.py index 9e4cb38..ee99264 100644 --- a/tests/test_srf.py +++ b/tests/test_srf.py @@ -529,7 +529,6 @@ def test_sw4_hdf5_read_write(): row[field.lower()], abs=1e-3 ) - # POINTS -- compare each DataFrame-sourced field directly points = h5file["POINTS"] assert points.shape == (len(original_srf.points),) for field in ( From 1c1814c8992dd16cf36644877db5c0cded56c3b1 Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Sat, 28 Mar 2026 16:59:35 +1300 Subject: [PATCH 17/18] adjust comments --- source_modelling/rupture_propagation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source_modelling/rupture_propagation.py b/source_modelling/rupture_propagation.py index 2385d80..e337365 100644 --- a/source_modelling/rupture_propagation.py +++ b/source_modelling/rupture_propagation.py @@ -59,7 +59,7 @@ def spanning_tree_with_probabilities( trees = [] probabilities = [] - for tree in mst.SpanningTreeIterator(graph): # ty: ignore[invalid-argument-type] - SpanningTreeIterator accepts undirected Graph at runtime; the networkx-stubs incorrectly restrict its __init__ to DiGraph only + for tree in mst.SpanningTreeIterator(graph): # ty: ignore[invalid-argument-type] p_tree = 1.0 for u, v in graph.edges: if tree.has_edge(u, v): From 3984534a94f4f671a780eb2c2aaaa9ad6230a445 Mon Sep 17 00:00:00 2001 From: Andrew Ridden-Harper Date: Sat, 28 Mar 2026 17:51:43 +1300 Subject: [PATCH 18/18] PR comments --- source_modelling/fsp.py | 9 +- source_modelling/rupture_propagation.py | 24 +---- source_modelling/sources.py | 4 +- source_modelling/srf.py | 40 ++++----- tests/test_gsf.py | 2 +- tests/test_magnitude_scaling.py | 6 +- tests/test_rupture_propagation.py | 8 +- tests/test_sources.py | 6 +- tests/test_srf.py | 114 +++++++++++------------- 9 files changed, 92 insertions(+), 121 deletions(-) diff --git a/source_modelling/fsp.py b/source_modelling/fsp.py index 66d2067..c0acc0e 100644 --- a/source_modelling/fsp.py +++ b/source_modelling/fsp.py @@ -5,8 +5,13 @@ See http://equake-rc.info/SRCMOD/fileformats/fsp/ for details on the FSP format. -Classes: ``FSPFile`` (representation of an FSP file). -Exceptions: ``FSPParseError`` (raised for errors in parsing FSP files). +Classes +------- +- FSPFile: Representation of an FSP file. + +Exceptions +---------- +- FSPParseError: Exception raised for errors in parsing FSP files. Examples -------- diff --git a/source_modelling/rupture_propagation.py b/source_modelling/rupture_propagation.py index e337365..e39ea12 100644 --- a/source_modelling/rupture_propagation.py +++ b/source_modelling/rupture_propagation.py @@ -73,30 +73,12 @@ def spanning_tree_with_probabilities( @overload -def sampled_spanning_tree(graph: nx.Graph, n_samples: Literal[1] = ...) -> nx.Graph: - """Sample a single spanning tree from a graph. - - Parameters - ---------- - graph : nx.Graph - The graph to sample from. - n_samples : Literal[1], optional - Number of trees to sample. Default is 1. - """ +def sampled_spanning_tree(graph: nx.Graph, n_samples: Literal[1] = ...) -> nx.Graph: # numpydoc ignore=GL08 ... @overload -def sampled_spanning_tree(graph: nx.Graph, n_samples: int) -> list[nx.Graph] | nx.Graph: - """Sample one or more spanning trees from a graph. - - Parameters - ---------- - graph : nx.Graph - The graph to sample from. - n_samples : int - Number of trees to sample. - """ +def sampled_spanning_tree(graph: nx.Graph, n_samples: int) -> list[nx.Graph] | nx.Graph: # numpydoc ignore=GL08 ... @@ -237,7 +219,7 @@ def select_top_spanning_trees( cumulative_tree_weight = 0.0 spanning_trees = [] - for spanning_tree in mst.SpanningTreeIterator(weighted_graph, minimum=False): # ty: ignore[invalid-argument-type] - SpanningTreeIterator accepts undirected Graph at runtime; the networkx-stubs incorrectly restrict its __init__ to DiGraph only + for spanning_tree in mst.SpanningTreeIterator(weighted_graph, minimum=False): # ty: ignore[invalid-argument-type] spanning_trees.append(spanning_tree) tree_log_probability = sum( spanning_tree[node_u][node_v]["weight"] diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 33fcc0a..51b299b 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -988,7 +988,7 @@ def __post_init__(self) -> None: # This relation can now be used to identify if the list of planes given is a line. points_into_graph: nx.DiGraph = nx.from_dict_of_lists( - points_into_relation, # ty: ignore[invalid-argument-type] - dict[int, list[int]] satisfies dict[Unknown, Iterable[Unknown]]; ty can't resolve types through networkx's _dispatchable wrapper + points_into_relation, # ty: ignore[invalid-argument-type] create_using=nx.DiGraph, ) try: @@ -1172,7 +1172,7 @@ def geometry(self) -> shapely.Polygon | shapely.LineString: # numpydoc ignore=R Geometry will be LineString if `dip = 90`. """ - return shapely.normalize( # ty: ignore[invalid-return-type] - shapely.normalize and union_all are typed to return BaseGeometry, but the result is always Polygon or LineString for plane geometries + return shapely.normalize( # ty: ignore[invalid-return-type] shapely.union_all([plane.geometry for plane in self.planes]) ) diff --git a/source_modelling/srf.py b/source_modelling/srf.py index 66acefd..68fe408 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -43,6 +43,7 @@ from pathlib import Path from typing import Self +import h5py import numpy as np import pandas as pd import scipy as sp @@ -338,9 +339,6 @@ def write_srf(self, srf_ffp: str | Path) -> None: def write_sw4_hdf5( self, output_ffp: Path | str, - vs: float | np.ndarray = 0.0, - den: float | np.ndarray = 0.0, - include_slip_time_function: bool = True, ) -> None: """Write the SRF file in SW4's SRF-HDF5 format. @@ -348,19 +346,17 @@ def write_sw4_hdf5( ---------- output_ffp : Path The path to the output HDF5 file. - vs : float or np.ndarray - Shear wave velocity (m/s). A scalar broadcasts to all points, - or provide a per-point array of shape (npoints,). The Version 1.0 - SRF format does not include this field; defaults to 0.0. - den : float or np.ndarray - Density (kg/m^3). A scalar broadcasts to all points, - or provide a per-point array of shape (npoints,). The Version 1.0 - SRF format does not include this field; defaults to 0.0. - include_slip_time_function : bool - If True, include the SR1 slip rate time series dataset. - """ - import h5py + References + ---------- + .. [1] Petersson, N.A. and B. Sjogreen (2017). SW4 v2.0. + Computational Infrastructure of Geodynamics, Davis, CA. + DOI: 10.5281/zenodo.1045297. + .. [2] Petersson, N.A. and B. Sjogreen (2017). User's guide to + SW4, version 2.0. Technical report LLNL-SM-741439, Lawrence + Livermore National Laboratory, Livermore, CA. + https://github.com/geodynamics/sw4/blob/master/doc/SW4_UsersGuide.pdf + """ plane_data = np.empty(len(self.header), dtype=SW4_PLANE_DTYPE) assert ( SW4_PLANE_DTYPE.names is not None @@ -382,21 +378,15 @@ def write_sw4_hdf5( "slip" if field == "SLIP1" else field.lower() ].values.astype(SW4_POINTS_DTYPE[field].type) - points_data["VS"] = np.asarray(vs, dtype=np.float32) - points_data["DEN"] = np.asarray(den, dtype=np.float32) - - if include_slip_time_function: - points_data["NT1"] = np.diff(self.slipt1_array.indptr).astype(np.int32) + points_data["NT1"] = np.diff(self.slipt1_array.indptr).astype(np.int32) with h5py.File(output_ffp, "w") as h5file: h5file.attrs.create("VERSION", np.float32(self.version)) h5file.attrs.create("PLANE", plane_data) h5file.create_dataset("POINTS", data=points_data) - - if include_slip_time_function: - h5file.create_dataset( - "SR1", data=self.slipt1_array.data.astype(np.float32) - ) + h5file.create_dataset( + "SR1", data=self.slipt1_array.data.astype(np.float32) + ) def write_hdf5( self, hdf5_ffp: Path, include_slip_time_function: bool = True diff --git a/tests/test_gsf.py b/tests/test_gsf.py index 2376a4c..a3ec59c 100644 --- a/tests/test_gsf.py +++ b/tests/test_gsf.py @@ -56,7 +56,7 @@ def test_plane_gsf(): # Use tmp_path fixture def test_bad_gsf_type(): with pytest.raises(TypeError): - gsf.source_to_gsf_dataframe(1, 0.1) # ty: ignore[invalid-argument-type] - intentionally passing wrong type to test error handling + gsf.source_to_gsf_dataframe(1, 0.1) # ty: ignore[invalid-argument-type] def test_write_gsf(tmp_path: Path): diff --git a/tests/test_magnitude_scaling.py b/tests/test_magnitude_scaling.py index 9c252b5..feb7100 100644 --- a/tests/test_magnitude_scaling.py +++ b/tests/test_magnitude_scaling.py @@ -86,7 +86,7 @@ def sampler( magnitude = draw(st.floats(min_value=min_magnitude, max_value=max_magnitude)) return scaling_relation, rake, magnitude - return sampler() # ty: ignore[missing-argument] - @st.composite strips the `draw` parameter; sampler() is valid with no args + return sampler() # ty: ignore[missing-argument] # The contreras slab 2020 relation (which is the same as Strasser 2010) is not invertible. See the coefficients of Table 2 in the paper @@ -116,7 +116,7 @@ def test_inversion( else: mag_to_area = MAGNITUDE_TO_AREA[scaling_relation] area_to_mag = AREA_TO_MAGNITUDE[scaling_relation] - assert area_to_mag(mag_to_area(magnitude)) == pytest.approx(magnitude) # ty: ignore[missing-argument] - rake is handled by functools.partial above for relations that require it + assert area_to_mag(mag_to_area(magnitude)) == pytest.approx(magnitude) # ty: ignore[missing-argument] @pytest.mark.parametrize( @@ -245,7 +245,7 @@ def test_monotonicity_mag_to_area( mag_to_area = functools.partial(MAGNITUDE_TO_AREA[scaling_relation], rake=rake) else: mag_to_area = MAGNITUDE_TO_AREA[scaling_relation] - assert mag_to_area(magnitude + 0.1) > mag_to_area(magnitude) # ty: ignore[missing-argument] - rake is handled by functools.partial above for relations that require it + assert mag_to_area(magnitude + 0.1) > mag_to_area(magnitude) # ty: ignore[missing-argument] class RandomFunction(Protocol): diff --git a/tests/test_rupture_propagation.py b/tests/test_rupture_propagation.py index 20562c4..9ba6efc 100644 --- a/tests/test_rupture_propagation.py +++ b/tests/test_rupture_propagation.py @@ -11,7 +11,7 @@ from source_modelling import rupture_propagation, sources -random_strongly_connected_graph = hnx.graph_builder( # ty: ignore[missing-argument] - hypothesis_networkx has no stubs; ty can't resolve ParamSpec through the decorator +random_strongly_connected_graph = hnx.graph_builder( # ty: ignore[missing-argument] graph_type=nx.Graph, node_keys=st.text( alphabet=st.characters( @@ -154,7 +154,7 @@ def test_random_sampling_root(graph: nx.Graph, n: int): random.seed(1) probability_of_spanning_trees = 0.0 tree_probabilities: dict[str, float] = {} - for tree in mst.SpanningTreeIterator(graph): # ty: ignore[invalid-argument-type] - SpanningTreeIterator accepts undirected Graph at runtime; the networkx-stubs incorrectly restrict its __init__ to DiGraph only + for tree in mst.SpanningTreeIterator(graph): # ty: ignore[invalid-argument-type] p_tree = 1.0 for u, v in graph.edges: if tree.has_edge(u, v): @@ -450,7 +450,7 @@ def test_sample_rupture_propagation( n_samples = 100 sampled_trees = [ rupture_propagation.sample_rupture_propagation( - sources_map, # ty: ignore[invalid-argument-type] - dict invariance: Point satisfies IsSource but dict[str, Point] is not assignable to dict[str, IsSource] + sources_map, # ty: ignore[invalid-argument-type] initial_source=initial_source, initial_source_distribution=initial_source_distribution, jump_impossibility_limit_distance=jump_impossibility_limit_distance, @@ -566,7 +566,7 @@ def test_jump_points_from_rupture_tree( expected_jump_points: dict[str, rupture_propagation.JumpPair], ): result_jump_points = rupture_propagation.jump_points_from_rupture_tree( - source_map, # ty: ignore[invalid-argument-type] - dict invariance: Point satisfies IsSource but dict[str, Point] is not assignable to dict[str, IsSource] + source_map, # ty: ignore[invalid-argument-type] rupture_causality_tree, ) diff --git a/tests/test_sources.py b/tests/test_sources.py index d5888a9..6dc78cd 100644 --- a/tests/test_sources.py +++ b/tests/test_sources.py @@ -418,7 +418,7 @@ def valid_trace_definition(draw: st.DrawFn): ) -@given(valid_trace_definition()) # ty: ignore[missing-argument] - @st.composite strips the `draw` parameter; valid_trace_definition() is valid with no args +@given(valid_trace_definition()) # ty: ignore[missing-argument] def test_plane_from_trace(data: tuple): ( trace_points_nztm, @@ -1227,7 +1227,7 @@ def fault( ) -@given(fault(min_segments=2, max_segments=10, min_length=0.4, max_length=100)) # ty: ignore[missing-argument] - @st.composite strips the `draw` parameter; fault() is valid with no positional args +@given(fault(min_segments=2, max_segments=10, min_length=0.4, max_length=100)) # ty: ignore[missing-argument] def test_simplify_fault(fault: Fault): tolerance = 0.4 simplified_fault = sources.simplify_fault(fault, tolerance) @@ -1253,7 +1253,7 @@ def test_simplify_fault(fault: Fault): @given( st.lists( - fault(min_segments=2, max_segments=10, min_length=0.4, max_length=100), # ty: ignore[missing-argument] - @st.composite strips the `draw` parameter; fault() is valid with no positional args + fault(min_segments=2, max_segments=10, min_length=0.4, max_length=100), # ty: ignore[missing-argument] min_size=2, max_size=10, ), diff --git a/tests/test_srf.py b/tests/test_srf.py index ee99264..f89b1c9 100644 --- a/tests/test_srf.py +++ b/tests/test_srf.py @@ -1,6 +1,7 @@ import tempfile from pathlib import Path +import h5py import numpy as np import pandas as pd import pytest @@ -407,7 +408,7 @@ def test_planes_nstk_1_ndip_gt_1(): version="1.0", header=header, points=points, - slipt1_array=None, # ty: ignore[invalid-argument-type] - intentionally passing None to construct a minimal mock SrfFile that doesn't need slip time functions + slipt1_array=None, # ty: ignore[invalid-argument-type] ) planes = mock_srf.planes @@ -450,7 +451,7 @@ def test_planes_nstk_1_ndip_1(): version="1.0", header=header, points=points, - slipt1_array=None, # ty: ignore[invalid-argument-type] - intentionally passing None to construct a minimal mock SrfFile that doesn't need slip time functions + slipt1_array=None, # ty: ignore[invalid-argument-type] ) planes = mock_srf.planes @@ -500,75 +501,68 @@ def test_hdf5_read_write(): original_srf.slipt1_array.indptr, reconstructed_srf.slipt1_array.indptr ), "slipt1_array indptr mismatch" - assert (original_srf.slipt1_array != reconstructed_srf.slipt1_array).nnz == 0, ( # ty: ignore[unresolved-attribute] - csr_array.__ne__ returns a sparse array at runtime; stubs incorrectly type it as bool + assert (original_srf.slipt1_array != reconstructed_srf.slipt1_array).nnz == 0, ( # ty: ignore[unresolved-attribute] "slipt1_array content mismatch" ) -def test_sw4_hdf5_read_write(): +def test_sw4_hdf5_read_write(tmp_path: Path): """Test that write_sw4_hdf5 preserves header, points, and slip data.""" - import h5py original_srf = srf.read_srf(SRF_DIR / "3468575.srf") - with tempfile.NamedTemporaryFile(suffix=".h5") as tmp_file: - original_srf.write_sw4_hdf5(Path(tmp_file.name)) - - with h5py.File(tmp_file.name, "r") as h5file: - # VERSION - assert h5file.attrs["VERSION"] == np.float32(1.0) - - plane = h5file.attrs["PLANE"] - assert plane.shape == (len(original_srf.header),) - assert ( - srf.SW4_PLANE_DTYPE.names is not None - ) - for idx, row in original_srf.header.iterrows(): - for field in srf.SW4_PLANE_DTYPE.names: - assert plane[idx][field] == pytest.approx( - row[field.lower()], abs=1e-3 - ) - - points = h5file["POINTS"] - assert points.shape == (len(original_srf.points),) - for field in ( - "LON", - "LAT", - "DEP", - "STK", - "DIP", - "AREA", - "TINIT", - "DT", - "RAKE", - ): - np.testing.assert_array_almost_equal( - points[field], - original_srf.points[field.lower()].to_numpy(), - decimal=3, + output_path = tmp_path / "test.h5" + original_srf.write_sw4_hdf5(output_path) + + with h5py.File(output_path, "r") as h5file: + # VERSION + assert h5file.attrs["VERSION"] == np.float32(1.0) + + plane = h5file.attrs["PLANE"] + assert plane.shape == (len(original_srf.header),) + assert ( + srf.SW4_PLANE_DTYPE.names is not None + ) + for idx, row in original_srf.header.iterrows(): + for field in srf.SW4_PLANE_DTYPE.names: + assert plane[idx][field] == pytest.approx( + row[field.lower()], abs=1e-3 ) - np.testing.assert_array_almost_equal( - points["SLIP1"], - original_srf.points["slip"].to_numpy(), - decimal=3, + + points = h5file["POINTS"] + assert points.shape == (len(original_srf.points),) + for field in ( + "LON", + "LAT", + "DEP", + "STK", + "DIP", + "AREA", + "TINIT", + "DT", + "RAKE", + ): + assert points[field] == pytest.approx( + original_srf.points[field.lower()].to_numpy(), abs=1e-3 ) + assert points["SLIP1"] == pytest.approx( + original_srf.points["slip"].to_numpy(), abs=1e-3 + ) - # VS/DEN default to 0.0 for Version 1.0 SRF - np.testing.assert_array_equal(points["VS"], 0.0) - np.testing.assert_array_equal(points["DEN"], 0.0) + # VS/DEN default to 0.0 for Version 1.0 SRF + assert points["VS"] == pytest.approx(0.0) + assert points["DEN"] == pytest.approx(0.0) - # NT1 from slipt1_array.indptr - np.testing.assert_array_equal( - points["NT1"], - np.diff(original_srf.slipt1_array.indptr).astype(np.int32), - ) + # NT1 from slipt1_array.indptr + assert points["NT1"] == pytest.approx( + np.diff(original_srf.slipt1_array.indptr).astype(np.int32) + ) - # SR1 slip-time function data - np.testing.assert_array_almost_equal( - h5file["SR1"][...], - original_srf.slipt1_array.data.astype(np.float32), - ) + # SR1 slip-time function data + assert h5file["SR1"][...] == pytest.approx( + original_srf.slipt1_array.data.astype(np.float32) + ) - # Unused slip components stay zero - for zero_field in ("SLIP2", "NT2", "SLIP3", "NT3"): - assert np.all(points[zero_field] == 0) + # Unused slip components stay zero + for zero_field in ("SLIP2", "NT2", "SLIP3", "NT3"): + assert points[zero_field] == pytest.approx(0)