Add ry0_distance transform and fix multi-fault GC2 strike discordance

source_modelling/gc2_distances.py

@@ -389,29 +389,70 @@ 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 _trace_endpoints(
+    trace_points: np.ndarray, trace_indices: np.ndarray
 ) -> tuple[np.ndarray, np.ndarray]:
-    """Multi-trace rx and ry calculations.
+    """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
-    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.
+    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]:
+    """Calculate segment lengths and U-coordinate origins for multi-trace systems.
+
+    Parameters
+    ----------
+    trace_points : np.ndarray
+        Trace points, with discordant traces reversed by
+        `strike_corrected_traces`.
+    trace_indices : np.ndarray
+        Trace indices for end of segments
+
+    Returns
+    -------
+    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)
@@ -422,25 +463,23 @@ def multi_trace_rx_ry(
         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(
             "Invalid direction vectors calculated: likely indicates degenerate trace or invalid geometry."
         )
     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,
@@ -449,4 +488,39 @@ 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.
+
+    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
+    )
     return generalised_t_u_coordinates(segment_lengths, rx, ry, u_shift_origins)

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,15 +1580,67 @@ 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)
 
 
+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) or (2,).
+
+    Returns
+    -------
+    np.ndarray
+        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]
+    faults = list(faults)
+
+    if len(faults) == 1:
+        _, u = faults[0].rx_ry_distance(point)
+        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:
+        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, trace_indices = _multi_fault_traces(faults)
+        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()
+
+    ry0 = np.maximum(0.0, np.maximum(u_min - u, u - u_max))
+    return ry0.reshape(np.shape(point)[:-1])
+
+
 IsSource = Plane | Fault | Point