Fix minimum jumping depth

source_modelling/rupture_propagation.py

@@ -497,6 +497,7 @@ def sample_rupture_propagation(
 def jump_points_from_rupture_tree(
     source_map: dict[str, sources.IsSource],
     rupture_causality_tree: Tree,
+    min_depth: float | None = None,
 ) -> dict[str, JumpPair]:
     """
     Extract jump points between faults from a rupture causality tree.
@@ -510,6 +511,8 @@ def jump_points_from_rupture_tree(
         A mapping of fault names to their corresponding source objects.
     rupture_causality_tree : Tree
         A rupture causality tree.
+    min_depth : float | None, optional
+        The minimum depth to consider jumping between, or ``None`` to allow jumps at all depths.
 
     Returns
     -------
@@ -521,9 +524,28 @@ def jump_points_from_rupture_tree(
     for source, parent in rupture_causality_tree.items():
         if parent is None:
             continue
-        source_point, parent_point = sources.closest_point_between_sources(
-            source_map[source], source_map[parent]
-        )
+        elif min_depth is not None:
+            source_a = source_map[source]
+            source_b = source_map[parent]
+            depth = min(
+                min_depth,
+                # HACK: factor of 0.99 is used here because the closest points
+                # solver will not work if the minimum depth is precisely the
+                # bottom-edge of the fault. If the closest point is the bottom
+                # depth then this will still recover that, but a proper
+                # treatment of the degenerate case would require specialising
+                # the solver.
+                0.99 * source_a.bottom_m / 1000,
+                0.99 * source_b.bottom_m / 1000,
+            )
+            source_point, parent_point = sources.closest_points_beneath(
+                source_a, source_b, depth
+            )
+        else:
+            source_point, parent_point = sources.closest_point_between_sources(
+                source_map[source],
+                source_map[parent],
+            )
         jump_points[source] = JumpPair(parent_point, source_point)
     return jump_points
 

source_modelling/sources.py

@@ -59,6 +59,31 @@ class Point:
     dip: float
     dip_dir: float
 
+    @property
+    def top_m(self) -> float:  # numpydoc ignore=RT01
+        """float: The top of the point source pseudo-geometry"""
+        centroid_depth = self.bounds[-1]
+        #   -------------------------+--------
+        #     \-       /             |
+        #       \--   / dip          |
+        #          \-/               |
+        #            \--             |
+        #       fault   o------------+ centroid depth
+        #                 \--        |
+        #                    \-      |  sin(dip) / 2 * width
+        #                      \--   |
+        #                         \- |
+        #                           \+
+
+        return centroid_depth - self.width_m * np.sin(np.radians(self.dip)) / 2
+
+    @property
+    def bottom_m(self) -> float:  # numpydoc ignore=RT01
+        """float: The bottom of the point source pseudo-geometry"""
+        centroid_depth = self.bounds[-1]
+
+        return centroid_depth + self.width_m * np.sin(np.radians(self.dip)) / 2
+
     @classmethod
     def from_lat_lon_depth(cls, point_coordinates: np.ndarray, **kwargs) -> Self:
         """Construct a point source from a lat, lon, depth format.
@@ -1643,7 +1668,7 @@ def fault_coordinate_distance(
 
 
 def closest_points_beneath(
-    source_a: Fault | Plane, source_b: Fault | Plane, min_depth: float
+    source_a: IsSource, source_b: IsSource, min_depth: float
 ) -> tuple[np.ndarray, np.ndarray]:
     """Find the closest points between two sources beneath a minimum depth.
 

tests/test_rupture_propagation.py

@@ -502,7 +502,7 @@ def test_sample_rupture_propagation(
 
 
 @pytest.mark.parametrize(
-    "source_map, rupture_causality_tree, expected_jump_points",
+    "source_map, rupture_causality_tree, min_depth, expected_jump_points",
     [
         # Test case 1: Simple rupture causality tree
         (
@@ -518,6 +518,7 @@ def test_sample_rupture_propagation(
                 ),
             },
             {"A": None, "B": "A", "C": "B"},
+            None,
             {
                 "B": rupture_propagation.JumpPair(
                     np.array([0.5, 0.5]), np.array([0.5, 0.5])
@@ -541,6 +542,7 @@ def test_sample_rupture_propagation(
                 ),
             },
             {"A": None, "B": "A", "C": "A"},
+            None,
             {
                 "B": rupture_propagation.JumpPair(
                     np.array([0.5, 0.5]), np.array([0.5, 0.5])
@@ -555,11 +557,13 @@ def test_sample_rupture_propagation(
 def test_jump_points_from_rupture_tree(
     source_map: dict[str, sources.Point],
     rupture_causality_tree: dict[str, str | None],
+    min_depth: float | None,
     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]
         rupture_causality_tree,
+        min_depth,
     )
 
     # Check if the jump points match the expected values
@@ -568,3 +572,67 @@ def test_jump_points_from_rupture_tree(
             result_jump_points[fault].from_point, expected_jump.from_point
         )
         assert np.allclose(result_jump_points[fault].to_point, expected_jump.to_point)
+
+
+# The bottom depth (in metres) of the sources used in the minimum depth
+# tests below (planes 10km wide dipping at 45 degrees from the surface).
+_DIPPING_PLANE_BOTTOM_M = 10_000 * np.sin(np.radians(45))
+
+
+@pytest.mark.parametrize(
+    "min_depth, expected_jump_depth_m",
+    [
+        # No minimum depth: the sources dip away from each other, so the
+        # closest points are at the surface.
+        (None, 0),
+        # The minimum depth lies within both sources and so is respected exactly.
+        (5, 5000),
+        # The minimum depth is below the bottom of both sources and so is
+        # clamped to just above their bottom depth.
+        (100, 0.99 * _DIPPING_PLANE_BOTTOM_M),
+    ],
+)
+def test_jump_points_from_rupture_tree_min_depth(
+    min_depth: float | None, expected_jump_depth_m: float
+):
+    # Two planes dipping away from each other, so that the distance
+    # between them increases with depth.
+    source_map = {
+        "A": sources.Plane.from_centroid_strike_dip(
+            np.array([-41.2865, 174.7762]),
+            dip=45,
+            length=10,
+            width=10,
+            dtop=0,
+            strike_nztm=0,
+            dip_dir_nztm=270,
+        ),
+        "B": sources.Plane.from_centroid_strike_dip(
+            np.array([-41.2865, 174.8762]),
+            dip=45,
+            length=10,
+            width=10,
+            dtop=0,
+            strike_nztm=0,
+            dip_dir_nztm=90,
+        ),
+    }
+
+    result_jump_points = rupture_propagation.jump_points_from_rupture_tree(
+        source_map,  # ty: ignore[invalid-argument-type]
+        {"A": None, "B": "A"},
+        min_depth,
+    )
+
+    jump = result_jump_points["B"]
+    from_depth_m = source_map["A"].fault_coordinates_to_wgs_depth_coordinates(
+        jump.from_point
+    )[-1]
+    to_depth_m = source_map["B"].fault_coordinates_to_wgs_depth_coordinates(
+        jump.to_point
+    )[-1]
+
+    # The sources dip away from each other, so the jump is made as
+    # shallow as the minimum depth allows.
+    assert from_depth_m == pytest.approx(expected_jump_depth_m, abs=1)
+    assert to_depth_m == pytest.approx(expected_jump_depth_m, abs=1)

tests/test_sources.py

@@ -14,7 +14,7 @@
 
 from qcore import coordinates, geo
 from source_modelling import sources
-from source_modelling.sources import Fault, Plane, multi_fault_rx_ry_distance
+from source_modelling.sources import Fault, Plane, Point, multi_fault_rx_ry_distance
 
 DATA_PATH = Path("tests") / "data"
 np.random.seed(0)
@@ -74,6 +74,20 @@ def test_point_construction(
     assert np.allclose(point.centroid, point_coordinates)
 
 
+def test_top_bottom_point():
+    point = Point(
+        coordinates.nztm_to_wgs_depth(np.array([-43.0, 172.0, 1000.0])),
+        1000.0,
+        1000.0,
+        0,
+        60.0,
+        90.0,
+    )
+    sin_dip = np.sqrt(3) / 2
+    assert point.bottom_m == 1000.0 + sin_dip / 2 * 1000.0
+    assert point.top_m == 1000.0 - sin_dip / 2 * 1000.0
+
+
 @given(
     point_coordinates=st.builds(
         coordinate,