diff --git a/source_modelling/rupture_propagation.py b/source_modelling/rupture_propagation.py index d74b3b3..d0f8853 100644 --- a/source_modelling/rupture_propagation.py +++ b/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 diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 5bffe43..162971c 100644 --- a/source_modelling/sources.py +++ b/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. diff --git a/tests/test_rupture_propagation.py b/tests/test_rupture_propagation.py index 3dfb722..7ab6cf5 100644 --- a/tests/test_rupture_propagation.py +++ b/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) diff --git a/tests/test_sources.py b/tests/test_sources.py index 758bad9..72b304a 100644 --- a/tests/test_sources.py +++ b/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,