diff --git a/pyproject.toml b/pyproject.toml index 719e167..02dcad3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,6 +12,7 @@ dynamic = ["version"] dependencies = [ "fiona", "geopandas", + "h5py", "networkx", "numpy", "pandas", @@ -46,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/fsp.py b/source_modelling/fsp.py index 0567f5a..c0acc0e 100644 --- a/source_modelling/fsp.py +++ b/source_modelling/fsp.py @@ -13,8 +13,8 @@ ---------- - 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/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/rupture_propagation.py b/source_modelling/rupture_propagation.py index 09b9cdc..e39ea12 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. @@ -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] p_tree = 1.0 for u, v in graph.edges: if tree.has_edge(u, v): @@ -72,6 +72,16 @@ def spanning_tree_with_probabilities( return trees, probabilities +@overload +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: # numpydoc ignore=GL08 + ... + + def sampled_spanning_tree( graph: nx.Graph, n_samples: int = 1 ) -> list[nx.Graph] | nx.Graph: @@ -209,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): + 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"] @@ -320,7 +330,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 +349,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 +400,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 3af5e8d..51b299b 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 @@ -21,6 +13,7 @@ import itertools import json import warnings +from collections.abc import Sequence from typing import NamedTuple, Self import networkx as nx @@ -110,8 +103,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 +431,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 +443,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 +516,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)) @@ -993,7 +988,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] + create_using=nx.DiGraph, ) try: self._validate_fault_plane_connectivity(points_into_graph) @@ -1171,12 +1167,12 @@ 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`. """ - return shapely.normalize( + return shapely.normalize( # ty: ignore[invalid-return-type] shapely.union_all([plane.geometry for plane in self.planes]) ) @@ -1247,8 +1243,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 +1330,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 4e7377e..68fe408 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. @@ -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 @@ -47,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 @@ -60,6 +57,46 @@ 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_EXTERNAL_FIELDS = {"VS", "DEN", "NT1", "SLIP2", "NT2", "SLIP3", "NT3"} + class Segments(Sequence): """A read-only view for SRF segments. @@ -85,7 +122,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 @@ -234,7 +273,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 ) @@ -262,7 +301,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 @@ -290,13 +329,65 @@ 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 | str, + ) -> None: + """Write the SRF file in SW4's SRF-HDF5 format. + + Parameters + ---------- + output_ffp : Path + The path to the output HDF5 file. + + 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 + ) + 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) + assert ( + SW4_POINTS_DTYPE.names is not None + ) + for field in SW4_POINTS_DTYPE.names: + if field in _SW4_POINTS_EXTERNAL_FIELDS: + continue + points_data[field] = self.points[ + "slip" if field == "SLIP1" else field.lower() + ].values.astype(SW4_POINTS_DTYPE[field].type) + + 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) + 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 +443,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) @@ -594,7 +687,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..a3ec59c 100644 --- a/tests/test_gsf.py +++ b/tests/test_gsf.py @@ -49,14 +49,14 @@ 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()) ) ) 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] def test_write_gsf(tmp_path: Path): @@ -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_magnitude_scaling.py b/tests/test_magnitude_scaling.py index 58d6041..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() + 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) + 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) + 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 ca1321f..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( +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): + 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, + 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,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] + 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..6dc78cd 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] 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] 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] min_size=2, max_size=10, ), diff --git a/tests/test_srf.py b/tests/test_srf.py index 847f539..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, # type: ignore + 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, # type: ignore + slipt1_array=None, # ty: ignore[invalid-argument-type] ) planes = mock_srf.planes @@ -500,6 +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, ( # type: ignore + 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(tmp_path: Path): + """Test that write_sw4_hdf5 preserves header, points, and slip data.""" + + original_srf = srf.read_srf(SRF_DIR / "3468575.srf") + + 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 + ) + + 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 + assert points["VS"] == pytest.approx(0.0) + assert points["DEN"] == pytest.approx(0.0) + + # 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 + 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 points[zero_field] == pytest.approx(0)