From 1fba3ed7a4873df920dc692a2dd608c1601cb941 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 19 Oct 2026 13:32:07 +1300 Subject: [PATCH 1/3] Add ry0_distance transform to sources Implements Ry0 = max(0, -U, U - L) as a standalone function. For multi-fault ruptures the rupture extent is taken as the span of the fault traces along the GC2 U axis, and a RuntimeWarning is raised since the transform is not well-defined there. Factors the segment length / U-origin calculation out of multi_trace_rx_ry into multi_trace_u_origins so the rupture extent can be computed without evaluating GC2 on the trace itself. Co-Authored-By: Claude Opus 5.5 --- source_modelling/gc2_distances.py | 49 ++++++++++++----- source_modelling/sources.py | 58 +++++++++++++++++++++ tests/test_rx_ry.py | 87 +++++++++++++++++++++++++++++++ 3 files changed, 182 insertions(+), 12 deletions(-) diff --git a/source_modelling/gc2_distances.py b/source_modelling/gc2_distances.py index 14ca6ad8..cc03f4b7 100644 --- a/source_modelling/gc2_distances.py +++ b/source_modelling/gc2_distances.py @@ -389,13 +389,10 @@ def strike_corrected_directions( return trace_directions -def multi_trace_rx_ry( - trace_points: np.ndarray, - trace_indices: np.ndarray, - rx: np.ndarray, - ry: np.ndarray, +def multi_trace_u_origins( + trace_points: np.ndarray, trace_indices: np.ndarray ) -> tuple[np.ndarray, np.ndarray]: - """Multi-trace rx and ry calculations. + """Calculate segment lengths and U-coordinate origins for multi-trace systems. Parameters ---------- @@ -403,15 +400,13 @@ def multi_trace_rx_ry( Trace points. trace_indices : np.ndarray Trace indices for end of segments - rx : np.ndarray - The rx values calculated to points. - ry : np.ndarray - The ry values calculated to points. Returns ------- - tuple[np.ndarray, np.ndarray] - The (T, U) generalised coordinates for Rx and Ry calculations. + segment_lengths : np.ndarray + The length of each segment, shape (m,). + u_shift_origins : np.ndarray + The U-coordinate origin of each segment, shape (m,). """ direction_vectors = diff_reduction(trace_points, trace_indices) segment_lengths = np.linalg.norm(direction_vectors, axis=1) @@ -449,4 +444,34 @@ def multi_trace_rx_ry( b_hat, ) + return segment_lengths, u_shift_origins + + +def multi_trace_rx_ry( + trace_points: np.ndarray, + trace_indices: np.ndarray, + rx: np.ndarray, + ry: np.ndarray, +) -> tuple[np.ndarray, np.ndarray]: + """Multi-trace rx and ry calculations. + + Parameters + ---------- + trace_points : np.ndarray + Trace points. + trace_indices : np.ndarray + Trace indices for end of segments + rx : np.ndarray + The rx values calculated to points. + ry : np.ndarray + The ry values calculated to points. + + Returns + ------- + tuple[np.ndarray, np.ndarray] + The (T, U) generalised coordinates for Rx and Ry calculations. + """ + segment_lengths, u_shift_origins = multi_trace_u_origins( + trace_points, trace_indices + ) return generalised_t_u_coordinates(segment_lengths, rx, ry, u_shift_origins) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 5f34c824..ef7cea2e 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -1561,6 +1561,64 @@ def multi_fault_rx_ry_distance( return gc2_distances.multi_trace_rx_ry(trace_points, trace_indices, rx, ry) +def ry0_distance( + faults: Fault | Plane | Sequence[Fault | Plane], point: np.ndarray +) -> np.ndarray: + """Calculate the Ry0 distance between a rupture and a set of points. + + Ry0 is the distance from the point to the ends of the rupture + measured parallel to strike, i.e. Ry0 = max(0, -U, U - L) where U + is the generalised ry distance and L is the length of the rupture + in the U coordinate system. Points with a projection along strike + falling within the rupture have Ry0 = 0. + + For multi-fault ruptures, the rupture is taken to span the + range of U values covered by the fault traces. The transformation + is not well-defined for these ruptures, so a RuntimeWarning is + raised. + + Parameters + ---------- + faults : Fault | Plane | Sequence[Fault | Plane] + The fault, or faults, of the rupture. + point : np.ndarray + Points to calculate distance to, has shape (n, 2). + + Returns + ------- + np.ndarray + The Ry0 distance (in metres) between the rupture and the points. Has shape (n,) + """ + if isinstance(faults, (Fault, Plane)): + faults = [faults] + faults = list(faults) + + traces = [fault.trace[:, :2] for fault in faults] + if len(faults) == 1: + _, u = faults[0].rx_ry_distance(point) + segments = traces[0].reshape((-1, 2, 2)) + u_min = 0.0 + u_max = np.linalg.norm(segments[:, 1] - segments[:, 0], axis=-1).sum() + else: + warnings.warn( + "Ry0 is not well-defined for multi-fault ruptures." + " Using the extent of the fault traces along the GC2 U axis as the rupture length.", + RuntimeWarning, + ) + _, u = multi_fault_rx_ry_distance(faults, point) + trace_points = np.concatenate(traces, axis=0) + trace_indices = np.cumulative_sum( + [len(trace) for trace in traces], include_initial=True + ) + segment_lengths, u_origins = gc2_distances.multi_trace_u_origins( + trace_points, trace_indices + ) + u_min = u_origins.min() + u_max = (u_origins + segment_lengths).max() + + return np.maximum(0.0, np.maximum(u_min - u, u - u_max)) + + IsSource = Plane | Fault | Point diff --git a/tests/test_rx_ry.py b/tests/test_rx_ry.py index bb3388f0..aad362c3 100644 --- a/tests/test_rx_ry.py +++ b/tests/test_rx_ry.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np import pytest import scipy as sp @@ -426,3 +428,88 @@ def test_multi_trace_rx_ry_local_shift_accumulation() -> None: # Seg 1 Global U: 10 + 0 = 10 # Seg 2 Global U: 5 + 5 = 10 assert u.item() == pytest.approx(10.0) + + +def test_ry0_plane() -> None: + """Test that Ry0 = max(0, -U, U - L) for a single plane.""" + plane = sources.Plane.from_centroid_strike_dip( + np.array([-43.538, 172.6474, 5.0]), 90.0, 10.0, width=5.0, strike=0.0 + ) + points = np.array( + [ + [-43.538, 172.6174], # beside the plane + [-43.60, 172.6574], # south of the plane + [-43.48, 172.6574], # north of the plane + ] + ) + _, ry = plane.rx_ry_distance(points) + length = plane.length_m + + ry0 = sources.ry0_distance(plane, points) + + assert ry0 == pytest.approx(np.maximum(0.0, np.maximum(-ry, ry - length)), abs=1e-3) + assert ry0[0] == pytest.approx(0.0, abs=1e-3) + assert ry0[1] > 0 + assert ry0[2] > 0 + + +def test_ry0_fault() -> None: + """Test that Ry0 = max(0, -U, U - L) for a multi-plane fault.""" + fault = sources.Fault.from_trace_points( + np.array( + [ + [-43.58302058, 172.64739978], + [-43.49297907, 172.6474], + [-43.44832311516278, 172.7160675056574], + ] + ), + dtop=0.0, + dbottom=10.0, + dip=90.0, + dip_dir=90.0, + ) + points = np.array( + [ + [-43.538, 172.6174], + [-43.60, 172.6574], + [-43.40, 172.75], + ] + ) + _, ry = fault.rx_ry_distance(points) + trace = fault.trace[:, :2] + length = np.linalg.norm(trace[1::2] - trace[::2], axis=1).sum() + + ry0 = sources.ry0_distance(fault, points) + + assert ry0 == pytest.approx(np.maximum(0.0, np.maximum(-ry, ry - length)), abs=1e-3) + assert ry0[0] == pytest.approx(0.0, abs=1e-3) + assert ry0[1] > 0 + assert ry0[2] > 0 + + +def test_ry0_single_fault_does_not_warn() -> None: + """Test that Ry0 does not warn for single-fault ruptures.""" + plane = sources.Plane.from_centroid_strike_dip( + np.array([-43.538, 172.6474, 5.0]), 90.0, 10.0, width=5.0, strike=0.0 + ) + with warnings.catch_warnings(): + warnings.simplefilter("error") + sources.ry0_distance([plane], np.array([[-43.538, 172.6174]])) + + +def test_ry0_multi_fault_warns() -> None: + """Test that Ry0 raises a RuntimeWarning for multi-fault ruptures.""" + plane_a = sources.Plane.from_centroid_strike_dip( + np.array([-43.538, 172.6474, 5.0]), 90.0, 10.0, width=5.0, strike=0.0 + ) + plane_b = sources.Plane.from_centroid_strike_dip( + np.array([-43.40, 172.70, 5.0]), 90.0, 10.0, width=5.0, strike=0.0 + ) + points = np.array([[-43.538, 172.6174], [-43.70, 172.65], [-43.20, 172.70]]) + + with pytest.warns(RuntimeWarning, match="multi-fault"): + ry0 = sources.ry0_distance([plane_a, plane_b], points) + + assert ry0.shape == (3,) + assert ry0[0] == pytest.approx(0.0, abs=1e-3) + assert np.all(ry0[1:] > 0) From a5ba73f40c608a5da7dec2ec6409d2cce3832258 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 19 Oct 2026 13:33:42 +1300 Subject: [PATCH 2/3] Make ry0_distance output shape match input points and test vectorised use Co-Authored-By: Claude Opus 5.5 --- source_modelling/sources.py | 8 +++++--- tests/test_rx_ry.py | 41 +++++++++++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+), 3 deletions(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index ef7cea2e..109cf1e5 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -1582,12 +1582,13 @@ def ry0_distance( faults : Fault | Plane | Sequence[Fault | Plane] The fault, or faults, of the rupture. point : np.ndarray - Points to calculate distance to, has shape (n, 2). + Points to calculate distance to, has shape (n, 2) or (2,). Returns ------- np.ndarray - The Ry0 distance (in metres) between the rupture and the points. Has shape (n,) + The Ry0 distance (in metres) between the rupture and the + points. Has shape (n,), or () for a single point of shape (2,). """ if isinstance(faults, (Fault, Plane)): faults = [faults] @@ -1616,7 +1617,8 @@ def ry0_distance( u_min = u_origins.min() u_max = (u_origins + segment_lengths).max() - return np.maximum(0.0, np.maximum(u_min - u, u - u_max)) + ry0 = np.maximum(0.0, np.maximum(u_min - u, u - u_max)) + return ry0.reshape(np.shape(point)[:-1]) IsSource = Plane | Fault | Point diff --git a/tests/test_rx_ry.py b/tests/test_rx_ry.py index aad362c3..b2818b87 100644 --- a/tests/test_rx_ry.py +++ b/tests/test_rx_ry.py @@ -513,3 +513,44 @@ def test_ry0_multi_fault_warns() -> None: assert ry0.shape == (3,) assert ry0[0] == pytest.approx(0.0, abs=1e-3) assert np.all(ry0[1:] > 0) + + +@pytest.mark.parametrize("multi_fault", [False, True]) +def test_ry0_vectorised(multi_fault: bool) -> None: + """Test that Ry0 over many points matches a per-point calculation.""" + fault = sources.Fault.from_trace_points( + np.array( + [ + [-43.58302058, 172.64739978], + [-43.49297907, 172.6474], + [-43.44832311516278, 172.7160675056574], + ] + ), + dtop=0.0, + dbottom=10.0, + dip=90.0, + dip_dir=90.0, + ) + faults: list[sources.Fault | sources.Plane] = [fault] + if multi_fault: + faults.append( + sources.Plane.from_centroid_strike_dip( + np.array([-43.30, 172.80, 5.0]), 90.0, 10.0, width=5.0, strike=20.0 + ) + ) + rng = np.random.default_rng(0) + points = np.column_stack( + [rng.uniform(-43.8, -43.2, 200), rng.uniform(172.4, 173.0, 200)] + ) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + ry0 = sources.ry0_distance(faults, points) + ry0_loop = np.array([sources.ry0_distance(faults, p) for p in points]) + ry0_single_row = sources.ry0_distance(faults, points[:1]) + + assert ry0.shape == (200,) + assert ry0_loop.shape == (200,) + assert ry0_single_row.shape == (1,) + assert ry0 == pytest.approx(ry0_loop) + assert np.any(ry0 == 0) and np.any(ry0 > 0) From e11fcdff255f8533ca1ad09dc97d00f9d485f294 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 19 Oct 2026 13:42:01 +1300 Subject: [PATCH 3/3] Reverse discordant traces in multi-fault GC2 Spudich & Chiou (2015, p. 5-6) reverse the point order of each trace whose end-minus-start vector disagrees with the system's general strike before calculating T and U. Previously only segment direction vectors were flipped (and only to compute the nominal strike), so discordant traces contributed rx/ry with the wrong orientation and U was not monotonic along strike (e.g. the Darfield rupture, where Greendale NW-SE and the fault near Hororata are discordant). - Add gc2_distances.strike_corrected_traces, deciding discordance per trace and reversing the trace points. - Use corrected traces in multi_fault_rx_ry_distance and ry0_distance. - Take the nominal strike as the sum of corrected trace vectors (Eq. 11) and place the U origin at the back end of it. Co-Authored-By: Claude Opus 5.5 --- source_modelling/gc2_distances.py | 73 ++++++++++++++++++++++++++----- source_modelling/sources.py | 42 +++++++++++++----- tests/test_rx_ry.py | 58 ++++++++++++++++++++++++ 3 files changed, 150 insertions(+), 23 deletions(-) diff --git a/source_modelling/gc2_distances.py b/source_modelling/gc2_distances.py index cc03f4b7..158c8cbd 100644 --- a/source_modelling/gc2_distances.py +++ b/source_modelling/gc2_distances.py @@ -389,6 +389,51 @@ def strike_corrected_directions( return trace_directions +def _trace_endpoints( + trace_points: np.ndarray, trace_indices: np.ndarray +) -> tuple[np.ndarray, np.ndarray]: + """Return the first and last point of each trace.""" + return trace_points[trace_indices[:-1]], trace_points[trace_indices[1:] - 1] + + +def strike_corrected_traces( + trace_points: np.ndarray, trace_indices: np.ndarray +) -> np.ndarray: + """Reverse the point order of discordant traces. + + Implements the algorithm described in page 6. A trace is discordant + if its end-minus-start vector disagrees with the general strike + direction of the entire fault system. Discordant traces must be + reversed before calculating rx and ry for each segment. + + Parameters + ---------- + trace_points : np.ndarray + Trace points. + trace_indices : np.ndarray + Trace indices for end of segments + + Returns + ------- + np.ndarray + Trace points with the points of each discordant trace in reverse order. + """ + trace_starts, trace_ends = _trace_endpoints(trace_points, trace_indices) + p_origin, p_end = trial_strike_vector(np.concatenate((trace_starts, trace_ends))) + trial_unit_vector = (p_end - p_origin) / np.linalg.norm(p_end - p_origin) + trace_directions = trace_ends - trace_starts + corrected_directions = strike_corrected_directions( + trace_directions, trial_unit_vector + ) + discordant = np.any(corrected_directions != trace_directions, axis=1) + + trace_points = trace_points.copy() + for j, (start, end) in enumerate(itertools.pairwise(trace_indices)): + if discordant[j]: + trace_points[start:end] = trace_points[start:end][::-1] + return trace_points + + def multi_trace_u_origins( trace_points: np.ndarray, trace_indices: np.ndarray ) -> tuple[np.ndarray, np.ndarray]: @@ -397,7 +442,8 @@ def multi_trace_u_origins( Parameters ---------- trace_points : np.ndarray - Trace points. + Trace points, with discordant traces reversed by + `strike_corrected_traces`. trace_indices : np.ndarray Trace indices for end of segments @@ -417,18 +463,11 @@ def multi_trace_u_origins( directions_per_trace, include_initial=True ) - trace_starts = trace_points[trace_indices[:-1]] - end_indices = trace_indices[1:] - 1 - trace_ends = trace_points[end_indices] - trace_endpoints = np.concatenate((trace_starts, trace_ends)) + trace_starts, trace_ends = _trace_endpoints(trace_points, trace_indices) - p_origin, p_end = trial_strike_vector(trace_endpoints) - trial_unit_vector = p_end - p_origin - trial_unit_vector /= np.linalg.norm(trial_unit_vector) - direction_vectors = strike_corrected_directions( - direction_vectors, trial_unit_vector - ) - b_hat = np.sum(direction_vectors, axis=0) + # Equation 11: the nominal strike is the sum of the (corrected) + # end-minus-start vectors of every trace. + b_hat = np.sum(trace_ends - trace_starts, axis=0) b_hat_norm = np.linalg.norm(b_hat) if np.isclose(b_hat_norm, 0.0): raise ValueError( @@ -436,6 +475,11 @@ def multi_trace_u_origins( ) b_hat /= b_hat_norm + # The origin is the farthest endpoint at the back end of the nominal strike. + p_origin, p_end = trial_strike_vector(np.concatenate((trace_starts, trace_ends))) + if np.dot(p_end - p_origin, b_hat) < 0: + p_origin = p_end + u_shift_origins = calculate_gc2_u_origins( segment_lengths, trace_direction_start_indices, @@ -470,6 +514,11 @@ def multi_trace_rx_ry( ------- tuple[np.ndarray, np.ndarray] The (T, U) generalised coordinates for Rx and Ry calculations. + + Notes + ----- + The ``trace_points`` (and the rx and ry values calculated from them) + must have discordant traces reversed by `strike_corrected_traces`. """ segment_lengths, u_shift_origins = multi_trace_u_origins( trace_points, trace_indices diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 109cf1e5..89faad6d 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -1532,6 +1532,34 @@ def ry_distance(self, point: np.ndarray) -> np.ndarray: return self.rx_ry_distance(point)[1] +def _multi_fault_traces( + faults: Sequence[Fault | Plane], +) -> tuple[np.ndarray, np.ndarray]: + """Combine the fault traces of a multi-fault rupture for GC2 calculations. + + Parameters + ---------- + faults : Sequence[Fault | Plane] + Faults to combine. + + Returns + ------- + trace_points : np.ndarray + The NZTM trace points of every fault, with discordant traces reversed. + trace_indices : np.ndarray + The index of the start of each trace in ``trace_points``, plus the total number of points. + """ + traces = [fault.trace[:, :2] for fault in faults] + trace_points = np.concatenate(traces, axis=0) + trace_indices = np.cumulative_sum( + [len(trace) for trace in traces], include_initial=True + ) + return ( + gc2_distances.strike_corrected_traces(trace_points, trace_indices), + trace_indices, + ) + + def multi_fault_rx_ry_distance( faults: list[Fault | Plane], point: np.ndarray ) -> tuple[np.ndarray, np.ndarray]: @@ -1552,11 +1580,7 @@ def multi_fault_rx_ry_distance( The generalised ry distance (in metres) between the faults and the points. Has shape (n,) """ point = coordinates.wgs_depth_to_nztm(point[..., :2]) - traces = [fault.trace[:, :2] for fault in faults] - trace_points = np.concatenate(traces, axis=0) - trace_indices = np.cumulative_sum( - [len(trace) for trace in traces], include_initial=True - ) + trace_points, trace_indices = _multi_fault_traces(faults) rx, ry = gc2_distances.segment_rx_ry(trace_points, point) return gc2_distances.multi_trace_rx_ry(trace_points, trace_indices, rx, ry) @@ -1594,10 +1618,9 @@ def ry0_distance( faults = [faults] faults = list(faults) - traces = [fault.trace[:, :2] for fault in faults] if len(faults) == 1: _, u = faults[0].rx_ry_distance(point) - segments = traces[0].reshape((-1, 2, 2)) + segments = faults[0].trace[:, :2].reshape((-1, 2, 2)) u_min = 0.0 u_max = np.linalg.norm(segments[:, 1] - segments[:, 0], axis=-1).sum() else: @@ -1607,10 +1630,7 @@ def ry0_distance( RuntimeWarning, ) _, u = multi_fault_rx_ry_distance(faults, point) - trace_points = np.concatenate(traces, axis=0) - trace_indices = np.cumulative_sum( - [len(trace) for trace in traces], include_initial=True - ) + trace_points, trace_indices = _multi_fault_traces(faults) segment_lengths, u_origins = gc2_distances.multi_trace_u_origins( trace_points, trace_indices ) diff --git a/tests/test_rx_ry.py b/tests/test_rx_ry.py index b2818b87..c8d0d301 100644 --- a/tests/test_rx_ry.py +++ b/tests/test_rx_ry.py @@ -554,3 +554,61 @@ def test_ry0_vectorised(multi_fault: bool) -> None: assert ry0_single_row.shape == (1,) assert ry0 == pytest.approx(ry0_loop) assert np.any(ry0 == 0) and np.any(ry0 > 0) + + +def _gc2_t_u( + trace_points: np.ndarray, points: np.ndarray +) -> tuple[np.ndarray, np.ndarray]: + """Calculate GC2 (T, U) for two-point traces in NZTM (northing, easting).""" + trace_indices = np.arange(0, len(trace_points) + 1, 2) + trace_points = gc2_distances.strike_corrected_traces(trace_points, trace_indices) + rx, ry = gc2_distances.segment_rx_ry(trace_points, points) + return gc2_distances.multi_trace_rx_ry(trace_points, trace_indices, rx, ry) + + +def test_strike_corrected_traces_reverses_discordant_trace() -> None: + """Test that discordant traces are reversed as a whole (page 6).""" + # Trace 1 strikes east with a short westward kink, so only a + # per-trace (not per-segment) check leaves it alone. Trace 2 + # strikes west and is discordant. + trace_points = np.array( + [ + [0.0, 0.0], + [0.0, 20.0], + [0.0, 20.0], + [1.0, 19.0], + [0.0, 40.0], + [0.0, 30.0], + ] + ) + trace_indices = np.array([0, 4, 6]) + + corrected = gc2_distances.strike_corrected_traces(trace_points, trace_indices) + + assert corrected[:4] == pytest.approx(trace_points[:4]) + assert corrected[4:] == pytest.approx(trace_points[4:][::-1]) + + +def test_multi_trace_t_u_invariant_to_discordant_trace_order() -> None: + """Test that reversing a trace's point order does not change T and U.""" + concordant = np.array([[0.0, 0.0], [0.0, 15.0], [0.0, 20.0], [5.0, 30.0]]) + discordant = np.array([[0.0, 0.0], [0.0, 15.0], [5.0, 30.0], [0.0, 20.0]]) + points = np.array([[3.0, -10.0], [-4.0, 5.0], [2.0, 15.0], [-6.0, 40.0]]) + + t_concordant, u_concordant = _gc2_t_u(concordant, points) + t_discordant, u_discordant = _gc2_t_u(discordant, points) + + assert t_discordant == pytest.approx(t_concordant) + assert u_discordant == pytest.approx(u_concordant) + # U increases along strike across both traces. + assert np.all(np.diff(u_concordant) > 0) + + +def test_multi_trace_u_origin_at_back_of_nominal_strike() -> None: + """Test that U = 0 at the start of a westward striking system.""" + trace_points = np.array([[0.0, 30.0], [0.0, 20.0], [0.0, 10.0], [0.0, 0.0]]) + points = np.array([[0.0, 30.0], [0.0, 25.0], [0.0, 5.0]]) + + _, u = _gc2_t_u(trace_points, points) + + assert u == pytest.approx([0.0, 5.0, 25.0])