Write SRF in format SW4 needs

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"]
 

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

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

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

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(

source_modelling/sources.py

@@ -4,23 +4,16 @@
 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
 import dataclasses
 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