Fix SRF writing logic

source_modelling/fsp.py

@@ -65,8 +65,6 @@ def _normalise_value(value: float) -> float | None:
 class FSPParseError(Exception):
     """Exception raised for errors in parsing FSP files."""
 
-    pass
-
 
 @dataclasses.dataclass
 class Segment:
@@ -278,11 +276,7 @@ def read_from_file(cls: Callable, fsp_ffp: Path) -> "FSPFile":
             for line in fsp_file_handle:
                 if line.startswith("% Data"):
                     break
-                if (
-                    line.startswith("% -")
-                    or line.strip() == "%"
-                    or line.startswith("% Event :")
-                ):
+                if line.startswith(("% -", "% Event :")) or line.strip() == "%":
                     continue
                 # Strip the leading "% ", and deduplicate the spaces in the line.
                 # This is required to normalise the string so that the parse

source_modelling/parse_utils.py

@@ -6,7 +6,6 @@
 class ParseError(Exception):
     """Error for parsing files in source modelling."""
 
-    pass
 
 
 def _is_seperator(char: str) -> bool:

source_modelling/rupture_propagation.py

@@ -59,7 +59,7 @@ def spanning_tree_with_probabilities(
     trees = []
     probabilities = []
 
-    for tree in mst.SpanningTreeIterator(graph):  # type: ignore
+    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):
@@ -223,7 +223,7 @@ def select_top_spanning_trees(
     cumulative_tree_weight = 0.0
     spanning_trees = []
 
-    for spanning_tree in mst.SpanningTreeIterator(weighted_graph, minimum=False):  # type: ignore
+    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"]

source_modelling/sources.py

@@ -1042,7 +1042,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,  # type: ignore
+            points_into_relation,  # ty: ignore[invalid-argument-type]
             create_using=nx.DiGraph,
         )
         try:
@@ -1221,12 +1221,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.Polygon | shapely.LineString:  # numpydoc ignore=RT01
+    def geometry(self) -> shapely.Geometry:  # 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(  # type: ignore
+        return shapely.normalize(
             shapely.union_all([plane.geometry for plane in self.planes])
         )
 
@@ -1376,8 +1376,8 @@ def rjb_distance(self, point: np.ndarray) -> float:
         float
             The Rjb distance (in metres) to the point.
         """
-        return self.geometry.distance(
-            shapely.Point(coordinates.wgs_depth_to_nztm(point))
+        return shapely.distance(
+            self.geometry, shapely.Point(coordinates.wgs_depth_to_nztm(point))
         )
 
     def rx_ry_distance(self, point: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
@@ -1512,15 +1512,14 @@ class CoordinateBounds(NamedTuple):
     """float: Maximum normalised dip coordinate, in the range of [0, 1]."""
 
 
+DEFAULT_BOUNDS = CoordinateBounds(min_strike=0, min_dip=0, max_strike=1, max_dip=1)
+
+
 def closest_point_between_sources(
     source_a: IsSource,
     source_b: IsSource,
-    source_a_coordinate_bounds: CoordinateBounds = CoordinateBounds(
-        min_strike=0, min_dip=0, max_strike=1, max_dip=1
-    ),
-    source_b_coordinate_bounds: CoordinateBounds = CoordinateBounds(
-        min_strike=0, min_dip=0, max_strike=1, max_dip=1
-    ),
+    source_a_coordinate_bounds: CoordinateBounds = DEFAULT_BOUNDS,
+    source_b_coordinate_bounds: CoordinateBounds = DEFAULT_BOUNDS,
 ) -> tuple[np.ndarray, np.ndarray]:
     """Find the closest point between two sources that have local coordinates.
 

source_modelling/stoch.py

@@ -14,7 +14,7 @@
 
 import dataclasses
 from pathlib import Path
-from typing import IO, NamedTuple, Self, TextIO, TypeAlias, cast
+from typing import IO, NamedTuple, Self, TextIO, cast
 
 import numpy as np
 from numpy.typing import NDArray
@@ -24,9 +24,9 @@
 from source_modelling.sources import Plane
 
 # Type aliases for NumPy arrays with specific shapes and dtypes
-FloatArray2D: TypeAlias = NDArray[np.float32]  # 2D array of float32
-LatLonArray: TypeAlias = NDArray[np.float64]  # Array for latitude/longitude
-CoordinateArray: TypeAlias = NDArray[np.float64]  # Array for coordinates
+type FloatArray2D = NDArray[np.float32]  # 2D array of float32
+type LatLonArray = NDArray[np.float64]  # Array for latitude/longitude
+type CoordinateArray = NDArray[np.float64]  # Array for coordinates
 
 
 class StochHeader(NamedTuple):

src/srf_parser.rs

@@ -53,6 +53,7 @@ fn read_srf_header(
         let ndip = scanner.next()?;
         let len = scanner.next()?;
         let wid = scanner.next()?;
+        scanner.expect_end_of_line()?;
         let stk = scanner.next()?;
         let dip = scanner.next()?;
         let dtop = scanner.next()?;
@@ -81,6 +82,7 @@ fn read_srf_header(
 const APPROX_BYTES_PER_SLIP_VALUE: usize = 12;
 
 // Tiny PointHeader used to construct actual points later.
+#[derive(Debug)]
 struct PointHeader {
     lon: f32,
     lat: f32,
@@ -206,11 +208,6 @@ fn read_srf_points_v2(
     let mut slipt1 = CsrMatrix::new(point_count, slipt1_capacity);
 
     for (i, plane) in planes.iter().enumerate() {
-        // In version 2.0 (and version 2.0 only), it is possible to construct SRFs
-        // with multiple POINT instantiations. Technically V1.0 SRFs could be
-        // constructed with multiple POINTS instantiations but we are deliberately
-        // parsing a stricter subset of the format.
-
         scanner.skip_token(b"POINTS")?;
         let plane_point_count = scanner.next()?;
         if plane.points() != plane_point_count {
@@ -293,7 +290,8 @@ mod tests {
 
     const SRF_V1: &[u8] = b"1.0\n\
 PLANE 1\n\
-0.0 0.0 2 1 4.0 2.0 90.0 45.0 0.0 0.0 1.0\n\
+0.0 0.0 2 1 4.0 2.0\n\
+90.0 45.0 0.0 0.0 1.0\n\
 POINTS 2\n\
 0.1 -43.0 5.0 90.0 45.0 1.0e10 0.5 0.1\n\
 30.0 1.5 3 0.0 0 0.0 0\n\
@@ -304,7 +302,8 @@ POINTS 2\n\
 
     const SRF_V2: &[u8] = b"2.0\n\
 PLANE 1\n\
-0.0 0.0 2 1 4.0 2.0 90.0 45.0 0.0 0.0 1.0\n\
+0.0 0.0 2 1 4.0 2.0\n\
+90.0 45.0 0.0 0.0 1.0\n\
 POINTS 2\n\
 0.1 -43.0 5.0 90.0 45.0 1.0e10 0.5 0.1 3.5 2.7\n\
 30.0 1.5 3 0.0 0 0.0 0\n\
@@ -347,8 +346,10 @@ POINTS 2\n\
     // Two 1x1 planes, each with its own POINTS block.
     const SRF_V2_TWO_PLANES: &[u8] = b"2.0\n\
 PLANE 2\n\
-0.0 0.0 1 1 4.0 2.0 90.0 45.0 0.0 0.0 1.0\n\
-0.5 0.5 1 1 4.0 2.0 90.0 45.0 0.0 0.0 1.0\n\
+0.0 0.0 1 1 4.0 2.0\n\
+90.0 45.0 0.0 0.0 1.0\n\
+0.5 0.5 1 1 4.0 2.0\n\
+90.0 45.0 0.0 0.0 1.0\n\
 POINTS 1\n\
 0.1 -43.0 5.0 90.0 45.0 1.0e10 0.5 0.1 3.5 2.7\n\
 30.0 1.5 2 0.0 0 0.0 0\n\
@@ -381,6 +382,34 @@ POINTS 1\n\
         assert_eq!(srf.slipt1.data, vec![0.1, 0.2, 0.4]);
     }
 
+    const SRF_V1_TWO_PLANES: &[u8] = b"1.0\n\
+PLANE 2\n\
+0.0 0.0 1 1 4.0 2.0\n\
+90.0 45.0 0.0 0.0 1.0\n\
+0.5 0.5 1 1 4.0 2.0\n\
+90.0 45.0 0.0 0.0 1.0\n\
+POINTS 2\n\
+0.1 -43.0 5.0 90.0 45.0 1.0e10 0.5 0.1\n\
+30.0 1.5 2 0.0 0 0.0 0\n\
+0.1 0.2\n\
+0.2 -43.1 5.5 90.0 45.0 1.0e10 0.6 0.1\n\
+45.0 2.0 1 0.0 0 0.0 0\n\
+0.4\n";
+
+    #[test]
+    fn parses_v1_with_multiple_planes() {
+        let mut scanner = scanner::Scanner::new(SRF_V1_TWO_PLANES);
+        let srf = read_srf_struct(&mut scanner).unwrap();
+        assert_eq!(srf.planes.len(), 2);
+        let metadata = match &srf.metadata {
+            SrfMetadataVersioned::V1(metadata) => metadata,
+            SrfMetadataVersioned::V2(_) => panic!("expected V1 metadata"),
+        };
+        assert_eq!(metadata.lon, vec![0.1, 0.2]);
+        assert_eq!(srf.slipt1.row_ptr, vec![0, 2, 3]);
+        assert_eq!(srf.slipt1.data, vec![0.1, 0.2, 0.4]);
+    }
+
     #[test]
     fn rejects_v1_point_count_mismatch() {
         // Plane declares 2x1 points but POINTS declares 3.
@@ -460,25 +489,16 @@ POINTS 1\n\
         assert!(read_srf_struct(&mut scanner).is_ok());
     }
 
+    const SRF_BAD_PLANE_HEADER: &[u8] = b"0.0 0.0 2 1 4.0 2.0 90.0 45.0 0.0 0.0 1.0";
+
     #[test]
-    fn mismatch_messages_read_correctly() {
-        let err = SrfParseError::PointCountMismatch {
-            declared: 3,
-            expected: 2,
-        };
-        assert_eq!(
-            err.to_string(),
-            "PLANE headers expect 2 total points but POINTS declares 3"
-        );
-        let err = SrfParseError::PlanePointCountMismatch {
-            plane: 1,
-            declared: 2,
-            expected: 1,
-        };
-        assert_eq!(
-            err.to_string(),
-            "plane 1 expects 1 points but its POINTS block declares 2"
-        );
+    fn plane_header_rejects_no_newline() {
+        let mut scanner = scanner::Scanner::new(SRF_BAD_PLANE_HEADER);
+        let err = read_srf_header(&mut scanner, 1).unwrap_err();
+        assert!(matches!(
+            err,
+            SrfParseError::Scanner(scanner::ScannerError::NoNewlineFound { .. })
+        ));
     }
 
     fn replace_once(data: &[u8], from: &[u8], to: &[u8]) -> Vec<u8> {

src/srf_writer.rs

@@ -113,18 +113,13 @@ fn write_plane_header<W: Write>(writer: &mut W, planes: &[SrfPlane]) -> Result<(
     for plane in planes {
         writeln!(
             writer,
-            "{} {} {} {} {} {} {} {} {} {} {}",
-            plane.elon,
-            plane.elat,
-            plane.nstk,
-            plane.ndip,
-            plane.len,
-            plane.wid,
-            plane.stk,
-            plane.dip,
-            plane.dtop,
-            plane.shyp,
-            plane.dhyp
+            "{} {} {} {} {} {}",
+            plane.elon, plane.elat, plane.nstk, plane.ndip, plane.len, plane.wid,
+        )?;
+        writeln!(
+            writer,
+            "{} {} {} {} {}",
+            plane.stk, plane.dip, plane.dtop, plane.shyp, plane.dhyp
         )?;
     }
     Ok(())
@@ -165,7 +160,8 @@ mod tests {
 
     const SRF_V1: &[u8] = b"1.0\n\
 PLANE 1\n\
-0.0 0.0 2 1 4.0 2.0 90.0 45.0 0.0 0.0 1.0\n\
+0.0 0.0 2 1 4.0 2.0\n\
+90.0 45.0 0.0 0.0 1.0\n\
 POINTS 2\n\
 0.1 -43.0 5.0 90.0 45.0 1.0e10 0.5 0.1\n\
 30.0 1.5 3 0.0 0 0.0 0\n\
@@ -176,8 +172,10 @@ POINTS 2\n\
 
     const SRF_V2_TWO_PLANES: &[u8] = b"2.0\n\
 PLANE 2\n\
-0.0 0.0 1 1 4.0 2.0 90.0 45.0 0.0 0.0 1.0\n\
-0.5 0.5 1 1 4.0 2.0 90.0 45.0 0.0 0.0 1.0\n\
+0.0 0.0 1 1 4.0 2.0\n\
+90.0 45.0 0.0 0.0 1.0\n\
+0.5 0.5 1 1 4.0 2.0\n\
+90.0 45.0 0.0 0.0 1.0\n\
 POINTS 1\n\
 0.1 -43.0 5.0 90.0 45.0 1.0e10 0.5 0.1 3.5 2.7\n\
 30.0 1.5 2 0.0 0 0.0 0\n\
@@ -275,7 +273,8 @@ POINTS 1\n\
     fn empty_slip_row_roundtrips() {
         let data = b"1.0\n\
 PLANE 1\n\
-0.0 0.0 1 1 4.0 2.0 90.0 45.0 0.0 0.0 1.0\n\
+0.0 0.0 1 1 4.0 2.0\n\
+90.0 45.0 0.0 0.0 1.0\n\
 POINTS 1\n\
 0.1 -43.0 5.0 90.0 45.0 1.0e10 0.5 0.1\n\
 30.0 1.5 0 0.0 0 0.0 0\n";

tests/test_magnitude_scaling.py

@@ -19,6 +19,7 @@
     magnitude_scaling.ScalingRelation.CONTRERAS_INTERFACE2017: magnitude_scaling.contreras_interface_magnitude_to_area,
     magnitude_scaling.ScalingRelation.CONTRERAS_SLAB2020: magnitude_scaling.strasser_slab_magnitude_to_area,
 }
+RELATIONS = tuple(MAGNITUDE_TO_AREA)
 
 AREA_TO_MAGNITUDE = {
     magnitude_scaling.ScalingRelation.LEONARD2014: magnitude_scaling.leonard_area_to_magnitude,
@@ -69,7 +70,7 @@ def test_rake_type(rake: float, expected: magnitude_scaling.RakeType):
 
 
 def relation_with_magnitude(
-    relations: list[magnitude_scaling.ScalingRelation] = list(MAGNITUDE_TO_AREA),
+    relations: tuple[magnitude_scaling.ScalingRelation, ...] = RELATIONS,
 ):
     @st.composite
     def sampler(
@@ -98,10 +99,10 @@ def sampler(
 # The coefficients are not invertible, so we cannot test the inversion of the area to magnitude function.
 @given(
     relation_with_magnitude(
-        [
+        (
             magnitude_scaling.ScalingRelation.LEONARD2014,
             magnitude_scaling.ScalingRelation.CONTRERAS_INTERFACE2017,
-        ]
+        )
     )
 )
 def test_inversion(
@@ -277,7 +278,7 @@ def test_normal_error_contreras_interface(area_to_mag: RandomFunction, area: flo
         sp.stats.norm,
         samples,
         statistic="ad",
-        known_params=dict(loc=area_to_mag(area), scale=0.73 / np.log(10)),
+        known_params={"loc": area_to_mag(area), "scale": 0.73 / np.log(10)},
     )
     assert result.pvalue > 0.05
 
@@ -306,7 +307,7 @@ def test_normal_error_contreras_interface_aspect_ratio(
         sp.stats.norm,
         np.log(samples),
         statistic="ad",
-        known_params=dict(loc=np.log(aspect_ratio(magnitude)), scale=sigma),
+        known_params={"loc": np.log(aspect_ratio(magnitude)), "scale": sigma},
     )
     assert result.pvalue > 0.05
 
@@ -358,10 +359,7 @@ def test_normal_error_strasser_slab(area_to_mag: RandomFunction, area: float):
     """Generate samples with random = True set on area_to_mag and check that it approximates the value with random = False."""
     samples = [area_to_mag(area, random=True) for _ in range(100)]
     result = sp.stats.goodness_of_fit(
-        sp.stats.norm,
-        samples,
-        statistic="ad",
-        known_params=dict(loc=area_to_mag(area)),
+        sp.stats.norm, samples, statistic="ad", known_params={"loc": area_to_mag(area)}
     )
     assert result.pvalue > 0.05
 
@@ -415,7 +413,7 @@ def test_normal_error_contreras_slab_aspect_ratio(
         sp.stats.norm,
         np.log(samples),
         statistic="ad",
-        known_params=dict(loc=np.log(aspect_ratio(magnitude)), scale=sigma),
+        known_params={"loc": np.log(aspect_ratio(magnitude)), "scale": sigma},
     )
     assert result.pvalue > 0.05