From 0d8c9ae6c8b85df920241db8a722e2ecfe3dd812 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 27 Oct 2026 10:54:21 +1300 Subject: [PATCH 1/6] Add Haynie et al. (2026) interface/intraslab tectonic classification Within the Slab2 footprint, the probability of a subduction interface event is the product of an on-interface probability, exp(-dz^2 / 2 sigma^2) with sigma = 10 km for the vertical distance to the interface, a subduction slip probability, exp(-K^2 / 2 sigma^2) with sigma = 30 degrees for the minimum Kagan angle to slip on the interface in pure reverse or along the Pacific-Australia relative plate motion (NUVEL-1A), and whether the centroid is no deeper than 70 km. The remainder is upper plate above the interface and intraslab below it. Fixed depths and missing mechanisms leave the corresponding factor at 1. Adds kagan_angle, plate_motion_azimuth and subduction_slip_planes, and CMTClassifier.haynie_tectonic_type(_probabilities). Co-Authored-By: Claude Opus 5.5 --- source_modelling/focal_mechanism.py | 285 +++++++++++++++++++++++++++- tests/test_focal_mechanism.py | 136 +++++++++++++ 2 files changed, 420 insertions(+), 1 deletion(-) diff --git a/source_modelling/focal_mechanism.py b/source_modelling/focal_mechanism.py index 03e6f9b..c663a32 100644 --- a/source_modelling/focal_mechanism.py +++ b/source_modelling/focal_mechanism.py @@ -18,7 +18,12 @@ in the depth below the interface. Alternatively, the modified NGA-SUB (2020) rule used by the NZGMDB classifies from the location and depth alone, by the event's position relative to the up-dip, seismogenic and -down-dip zones of the interface. +down-dip zones of the interface. A third rule (Haynie et al. 2026) scores +interface events by the product of Gaussian likelihoods in the vertical +distance to the interface and the Kagan angle between the mechanism and +interface slip (pure reverse, or along the relative plate motion), +splitting the remainder into upper plate and intraslab by the side of the +interface. Examples -------- @@ -34,6 +39,8 @@ >>> classifier.nga_sub_tectonic_type(centroid) +>>> classifier.haynie_tectonic_type(centroid, nodal_plane_1) # 9 km below Slab2 + """ from __future__ import annotations @@ -99,6 +106,18 @@ CENTROID_DEPTH_SIGMA_KM = 7.0 """Standard deviation of the centroid depth for the NGA-SUB rule.""" +ON_INTERFACE_SIGMA_KM = 10.0 +"""Width of the Gaussian on-interface probability in the vertical distance to the interface (Haynie et al. 2026).""" + +SUBDUCTION_SLIP_SIGMA = 30.0 +"""Width in degrees of the Gaussian subduction slip probability in the Kagan angle (Haynie et al. 2026).""" + +SUBDUCTION_MAX_DEPTH_KM = 70.0 +"""Maximum centroid depth of subduction interface events (Haynie et al. 2026).""" + +PACIFIC_AUSTRALIA_EULER_POLE = (-60.1, -178.3, 1.0744) +"""Latitude, longitude (degrees) and rate (degrees/Myr) of Pacific rotation relative to Australia (NUVEL-1A, DeMets et al. 1994).""" + FAULT_NEIGHBOURS = 8 """Number of nearest CFM fault segments used for the strike, dip and rake misfits.""" @@ -171,6 +190,91 @@ def _plane_normal(strike: float, dip: float) -> npt.NDArray[np.float64]: ) +def _slip_vector(strike: float, dip: float, rake: float) -> npt.NDArray[np.float64]: + """Unit slip vector of the hanging wall in (north, east, down) coordinates (Aki & Richards eq. 4.88).""" + strike, dip, rake = np.radians(strike), np.radians(dip), np.radians(rake) + return np.array( + [ + np.cos(rake) * np.cos(strike) + np.sin(rake) * np.cos(dip) * np.sin(strike), + np.cos(rake) * np.sin(strike) - np.sin(rake) * np.cos(dip) * np.cos(strike), + -np.sin(rake) * np.sin(dip), + ] + ) + + +def _principal_axes(plane: NodalPlane) -> npt.NDArray[np.float64]: + """Rows T, P and B = T x P of a double couple in (north, east, down) coordinates.""" + normal = _plane_normal(plane.strike, plane.dip) + slip = _slip_vector(plane.strike, plane.dip, plane.rake) + t_axis = (normal + slip) / np.sqrt(2.0) + p_axis = (normal - slip) / np.sqrt(2.0) + return np.stack([t_axis, p_axis, np.cross(t_axis, p_axis)]) + + +def kagan_angle(plane_a: NodalPlane, plane_b: NodalPlane) -> float: + """Kagan angle between two double-couple focal mechanisms. + + The Kagan angle (Kagan 1991) is the smallest rotation taking one + mechanism's principal axes onto the other's, allowing for the + symmetries of the double couple. Either nodal plane of a mechanism gives + the same angle. + + Parameters + ---------- + plane_a : NodalPlane + A nodal plane of the first mechanism. + plane_b : NodalPlane + A nodal plane of the second mechanism. + + Returns + ------- + float + The Kagan angle in degrees, between 0 and 120. + """ + # Diagonal of the rotation between the two principal axis frames; the + # double-couple symmetries flip the signs of two axes at a time. + m = np.diag(_principal_axes(plane_a) @ _principal_axes(plane_b).T) + trace = max( + m[0] + m[1] + m[2], m[0] - m[1] - m[2], -m[0] + m[1] - m[2], -m[0] - m[1] + m[2] + ) + return float(np.degrees(np.arccos(np.clip((trace - 1.0) / 2.0, -1.0, 1.0)))) + + +def plate_motion_azimuth( + lat: float, + lon: float, + euler_pole: tuple[float, float, float] = PACIFIC_AUSTRALIA_EULER_POLE, +) -> float: + """Azimuth of relative plate motion at a location. + + Parameters + ---------- + lat : float + Latitude in degrees. + lon : float + Longitude in degrees. + euler_pole : tuple[float, float, float], optional + Latitude, longitude (degrees) and rate of rotation of one plate + relative to the other. Defaults to Pacific relative to Australia. + + Returns + ------- + float + Azimuth in degrees of the rotating plate's motion relative to the + fixed plate. + """ + pole_lat, pole_lon, rate = euler_pole + velocity = np.cross( + rate * _unit_vectors(pole_lat, pole_lon), _unit_vectors(lat, lon) + ) + lat, lon = np.radians(lat), np.radians(lon) + north = np.array( + [-np.sin(lat) * np.cos(lon), -np.sin(lat) * np.sin(lon), np.cos(lat)] + ) + east = np.array([-np.sin(lon), np.cos(lon), 0.0]) + return float(np.degrees(np.arctan2(velocity @ east, velocity @ north)) % 360.0) + + class SlabQuery(NamedTuple): """Slab surface properties at a location.""" @@ -653,6 +757,119 @@ def nga_sub_tectonic_type_probabilities( return probabilities +def subduction_slip_planes( + lat: float, lon: float, slab: SlabQuery +) -> tuple[NodalPlane, NodalPlane]: + """Mechanisms of slip on the subduction interface at a location. + + Parameters + ---------- + lat : float + Latitude in degrees. + lon : float + Longitude in degrees. + slab : SlabQuery + Slab surface beneath the location. + + Returns + ------- + tuple[NodalPlane, NodalPlane] + Slip on the slab surface (a) in pure reverse (rake 90) and (b) in + oblique reverse with the horizontal slip parallel to the relative + Pacific-Australia plate motion. + """ + # Components of the plate motion along strike and down the dip direction. + relative_azimuth = np.radians(plate_motion_azimuth(lat, lon) - slab.strike) + along_strike, down_dip = np.cos(relative_azimuth), np.sin(relative_azimuth) + # Take the sense of the motion that is reverse slip, whichever plate + # is subducting. + if down_dip > 0: + along_strike, down_dip = -along_strike, -down_dip + # The horizontal slip of a plane is cos(rake) along strike and + # -sin(rake) cos(dip) down the dip direction. + rake = np.degrees( + np.arctan2(-down_dip / np.cos(np.radians(slab.dip)), along_strike) + ) + return ( + NodalPlane(slab.strike, slab.dip, 90.0), + NodalPlane(slab.strike, slab.dip, float(rake)), + ) + + +def haynie_tectonic_type_probabilities( + depth_above_slab: float, + depth: float, + subduction_kagan_angle: float | None = None, + fixed_depth: bool = False, + on_interface_sigma: float = ON_INTERFACE_SIGMA_KM, + subduction_slip_sigma: float = SUBDUCTION_SLIP_SIGMA, +) -> dict[TectonicType, float]: + """Tectonic type probabilities by the rule of Haynie et al. (2026). + + Events outside the slab footprint are crustal (upper plate). Inside it, + the subduction interface probability is the product of + + - the on-interface probability, ``exp(-dz**2 / (2 * on_interface_sigma**2))`` + for vertical distance ``dz`` to the interface (1 for a fixed depth), + - the subduction slip probability, + ``exp(-kagan**2 / (2 * subduction_slip_sigma**2))`` for the Kagan + angle to interface slip (1 without a focal mechanism), and + - whether the event is no deeper than `SUBDUCTION_MAX_DEPTH_KM`. + + The remaining probability is upper plate for events above the + interface and intraslab for events below it. + + Parameters + ---------- + depth_above_slab : float + Vertical distance in km from the centroid up to the slab surface + (negative below it), NaN outside the slab footprint. + depth : float + Centroid depth in km. + subduction_kagan_angle : float | None, optional + Minimum Kagan angle in degrees between the event's mechanism and + slip on the interface (see `subduction_slip_planes`), None without + a focal mechanism. + fixed_depth : bool, optional + Whether the centroid depth was fixed rather than solved for, so says + nothing about the distance to the interface. + on_interface_sigma : float, optional + Width in km of the on-interface probability. + subduction_slip_sigma : float, optional + Width in degrees of the subduction slip probability. + + Returns + ------- + dict[TectonicType, float] + Probability of each tectonic type. + """ + if np.isnan(depth_above_slab): + return { + TectonicType.CRUSTAL: 1.0, + TectonicType.INTERFACE: 0.0, + TectonicType.SLAB: 0.0, + } + on_interface = ( + 1.0 + if fixed_depth + else np.exp(-(depth_above_slab**2) / (2.0 * on_interface_sigma**2)) + ) + subduction_slip = ( + 1.0 + if subduction_kagan_angle is None + else np.exp(-(subduction_kagan_angle**2) / (2.0 * subduction_slip_sigma**2)) + ) + interface = float( + on_interface * subduction_slip * (depth <= SUBDUCTION_MAX_DEPTH_KM) + ) + above_slab = depth_above_slab >= 0 + return { + TectonicType.CRUSTAL: (1.0 - interface) * above_slab, + TectonicType.INTERFACE: interface, + TectonicType.SLAB: (1.0 - interface) * (not above_slab), + } + + def depth_temperature( depth: npt.ArrayLike, scale: float, depth_rate: float ) -> npt.NDArray[np.float64]: @@ -963,3 +1180,69 @@ def nga_sub_tectonic_type(self, centroid: npt.ArrayLike) -> TectonicType: """ probabilities = self.nga_sub_tectonic_type_probabilities(centroid, sigma=0.0) return max(probabilities, key=probabilities.__getitem__) + + def haynie_tectonic_type_probabilities( + self, + centroid: npt.ArrayLike, + nodal_plane: NodalPlane | None = None, + fixed_depth: bool = False, + ) -> dict[TectonicType, float]: + """Tectonic type probabilities by the Haynie et al. (2026) rule (see the module-level `haynie_tectonic_type_probabilities`). + + Parameters + ---------- + centroid : npt.ArrayLike + Centroid as (lat, lon, depth_km). If depth is omitted a default + is assumed with a warning and treated as fixed. + nodal_plane : NodalPlane | None, optional + Either nodal plane of the event's focal mechanism, None if it + has none. + fixed_depth : bool, optional + Whether the centroid depth was fixed rather than solved for. + + Returns + ------- + dict[TectonicType, float] + Probability of each tectonic type. + """ + fixed_depth = fixed_depth or np.size(centroid) == 2 + lat, lon, depth = _parse_centroid(centroid) + slab = self.slab_model.query(lat, lon) + subduction_kagan_angle = None + if nodal_plane is not None and slab.present: + subduction_kagan_angle = min( + kagan_angle(nodal_plane, subduction_plane) + for subduction_plane in subduction_slip_planes(lat, lon, slab) + ) + return haynie_tectonic_type_probabilities( + slab.depth - depth, depth, subduction_kagan_angle, fixed_depth + ) + + def haynie_tectonic_type( + self, + centroid: npt.ArrayLike, + nodal_plane: NodalPlane | None = None, + fixed_depth: bool = False, + ) -> TectonicType: + """Most probable tectonic type by the Haynie et al. (2026) rule (see `haynie_tectonic_type_probabilities`). + + Parameters + ---------- + centroid : npt.ArrayLike + Centroid as (lat, lon, depth_km). If depth is omitted a default + is assumed with a warning and treated as fixed. + nodal_plane : NodalPlane | None, optional + Either nodal plane of the event's focal mechanism, None if it + has none. + fixed_depth : bool, optional + Whether the centroid depth was fixed rather than solved for. + + Returns + ------- + TectonicType + The tectonic type. + """ + probabilities = self.haynie_tectonic_type_probabilities( + centroid, nodal_plane, fixed_depth + ) + return max(probabilities, key=probabilities.__getitem__) diff --git a/tests/test_focal_mechanism.py b/tests/test_focal_mechanism.py index 554b168..86653c7 100644 --- a/tests/test_focal_mechanism.py +++ b/tests/test_focal_mechanism.py @@ -166,3 +166,139 @@ def test_nga_sub_tectonic_type_probabilities( ) assert sum(probabilities.values()) == pytest.approx(1.0) assert max(probabilities, key=probabilities.__getitem__) == expected + + +@pytest.mark.parametrize( + "plane_a, plane_b, expected", + [ + (PLANE_1, PLANE_1, 0.0), + (NodalPlane(0, 45, 90), NodalPlane(180, 45, 90), 0.0), # Auxiliary plane + (NodalPlane(0, 45, 90), NodalPlane(0, 45, -90), 90.0), # Thrust vs normal + (NodalPlane(0, 45, 90), NodalPlane(0, 45, 60), 30.0), # Rake rotation + (NodalPlane(0, 90, 0), NodalPlane(90, 90, 0), 90.0), # Strike rotation + ], +) +def test_kagan_angle(plane_a: NodalPlane, plane_b: NodalPlane, expected: float): + assert focal_mechanism.kagan_angle(plane_a, plane_b) == pytest.approx( + expected, abs=1e-4 + ) + assert focal_mechanism.kagan_angle(plane_b, plane_a) == pytest.approx( + expected, abs=1e-4 + ) + + +def test_kagan_angle_is_at_most_120(): + rng = np.random.default_rng(0) + planes = rng.uniform([0, 0, -180], [360, 90, 180], size=(500, 3)) + angles = [ + focal_mechanism.kagan_angle(NodalPlane(*a), NodalPlane(*b)) + for a, b in zip(planes[::2], planes[1::2]) + ] + assert 0.0 <= min(angles) and max(angles) <= 120.0 + + +def test_plate_motion_azimuth(): + # Pacific moves west relative to Australia at the Hikurangi margin. + assert 250 < focal_mechanism.plate_motion_azimuth(-38.7, 178.5) < 280 + assert 230 < focal_mechanism.plate_motion_azimuth(-47.0, 165.5) < 250 + + +def test_subduction_slip_planes(): + slab_model = focal_mechanism.SlabModel.load() + for lat, lon in [(-41.3, 174.8), (-38.7, 177.5), (-46.0, 166.0)]: + slab = slab_model.query(lat, lon) + pure, oblique = focal_mechanism.subduction_slip_planes(lat, lon, slab) + assert pure == NodalPlane(slab.strike, slab.dip, 90.0) + assert (oblique.strike, oblique.dip) == (slab.strike, slab.dip) + assert 0 < oblique.rake < 180 # Reverse sense + # The horizontal slip is parallel to the plate motion. + slip = focal_mechanism._slip_vector(*oblique) + slip_azimuth = np.degrees(np.arctan2(slip[1], slip[0])) + difference = focal_mechanism._angular_difference( + slip_azimuth, focal_mechanism.plate_motion_azimuth(lat, lon), 180.0 + ) + assert difference == pytest.approx(0.0, abs=1e-6) + + +@pytest.mark.parametrize( + "depth_above_slab, depth, kagan, fixed_depth, expected", + [ + (float("nan"), 10.0, 0.0, False, {TectonicType.CRUSTAL: 1.0}), + (0.0, 20.0, 0.0, False, {TectonicType.INTERFACE: 1.0}), + (0.0, 20.0, None, False, {TectonicType.INTERFACE: 1.0}), + (-1.0, 80.0, 0.0, False, {TectonicType.SLAB: 1.0}), # Too deep + ( + 10.0, + 20.0, + 30.0, + False, + { + TectonicType.CRUSTAL: 1.0 - np.exp(-1.0), + TectonicType.INTERFACE: np.exp(-1.0), + }, + ), + ( + -10.0, + 40.0, + None, + False, + { + TectonicType.SLAB: 1.0 - np.exp(-0.5), + TectonicType.INTERFACE: np.exp(-0.5), + }, + ), + # A fixed depth says nothing about the distance to the interface. + (-30.0, 50.0, 0.0, True, {TectonicType.INTERFACE: 1.0}), + (-30.0, 50.0, 90.0, True, {TectonicType.SLAB: 1.0 - np.exp(-4.5)}), + ], +) +def test_haynie_tectonic_type_probabilities( + depth_above_slab: float, + depth: float, + kagan: float | None, + fixed_depth: bool, + expected: dict[TectonicType, float], +): + probabilities = focal_mechanism.haynie_tectonic_type_probabilities( + depth_above_slab, depth, kagan, fixed_depth + ) + assert sum(probabilities.values()) == pytest.approx(1.0) + for tectonic_type, probability in expected.items(): + assert probabilities[tectonic_type] == pytest.approx(probability) + + +def test_haynie_tectonic_type(classifier: focal_mechanism.CMTClassifier): + # Thrust on the Hikurangi interface beneath Wellington. + slab = classifier.slab_model.query(-41.3, 174.8) + on_interface = np.array([-41.3, 174.8, slab.depth]) + thrust = NodalPlane(slab.strike, slab.dip, 100.0) + assert classifier.haynie_tectonic_type(on_interface, thrust) == ( + TectonicType.INTERFACE + ) + # Either nodal plane gives the same probabilities. + np.testing.assert_allclose( + list( + classifier.haynie_tectonic_type_probabilities( + on_interface, PLANE_1 + ).values() + ), + list( + classifier.haynie_tectonic_type_probabilities( + on_interface, PLANE_2 + ).values() + ), + atol=0.01, + ) + # Strike-slip near the interface is not interface slip. + strike_slip = NodalPlane(40.0, 85.0, 175.0) + shallower = on_interface - [0, 0, 2.0] + assert classifier.haynie_tectonic_type(shallower, strike_slip) == ( + TectonicType.CRUSTAL + ) + deeper = on_interface + [0, 0, 2.0] + assert classifier.haynie_tectonic_type(deeper, strike_slip) == TectonicType.SLAB + # Outside the slab footprint (Christchurch). + assert ( + classifier.haynie_tectonic_type([-43.5, 172.6, 10.0], thrust) + == TectonicType.CRUSTAL + ) From 8de75e21530608a961848cc478a43bb72c2ed3c8 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 27 Oct 2026 14:16:01 +1300 Subject: [PATCH 2/6] Replace NSHM 2022 and NGA-SUB tectonic rules with Haynie et al. (2026) CMTClassifier.tectonic_type(centroid, nodal_plane_1, nodal_plane_2), the call gcmt-to-realisation makes, now classifies by the Haynie et al. rule. The planes are optional (no focal mechanism) and a fixed_depth keyword drops the on-interface factor. The Kagan angle is taken over both planes so the result does not depend on their order. Removes the NSHM 2022 rule (tectonic_type_probabilities' old signature, _is_interface_like and its constants), the NGA-SUB rule (nga_sub_tectonic_type(_probabilities), SlabZone, SlabZoneQuery, SlabModel.zone) and the haynie_* names. Co-Authored-By: Claude Opus 5.5 --- source_modelling/focal_mechanism.py | 473 ++-------------------------- tests/test_focal_mechanism.py | 156 ++------- 2 files changed, 69 insertions(+), 560 deletions(-) diff --git a/source_modelling/focal_mechanism.py b/source_modelling/focal_mechanism.py index c663a32..7b066fe 100644 --- a/source_modelling/focal_mechanism.py +++ b/source_modelling/focal_mechanism.py @@ -13,17 +13,13 @@ ``scripts/train_cmt_classifiers.py``. The tectonic type (crustal, subduction interface or intraslab) follows the -NZ NSHM 2022 rule (Rollins et al. 2022) against the Slab2 interface -geometry (Hayes et al. 2018), with probabilities from Gaussian uncertainty -in the depth below the interface. Alternatively, the modified NGA-SUB -(2020) rule used by the NZGMDB classifies from the location and depth -alone, by the event's position relative to the up-dip, seismogenic and -down-dip zones of the interface. A third rule (Haynie et al. 2026) scores -interface events by the product of Gaussian likelihoods in the vertical -distance to the interface and the Kagan angle between the mechanism and -interface slip (pure reverse, or along the relative plate motion), -splitting the remainder into upper plate and intraslab by the side of the -interface. +probabilistic rule of Haynie et al. (2026) against the Slab2 interface +geometry (Hayes et al. 2018). The interface probability is the product of +Gaussian likelihoods in the vertical distance to the interface and in the +Kagan angle between the mechanism and interface slip (pure reverse, or +along the relative plate motion), for events no deeper than 70 km. The +remainder is crustal (upper plate) above the interface and intraslab below +it. Examples -------- @@ -36,10 +32,6 @@ >>> classifier.most_likely_nodal_plane(centroid, nodal_plane_1, nodal_plane_2) NodalPlane(strike=20, dip=35, rake=79) >>> classifier.tectonic_type(centroid, nodal_plane_1, nodal_plane_2) - ->>> classifier.nga_sub_tectonic_type(centroid) - ->>> classifier.haynie_tectonic_type(centroid, nodal_plane_1) # 9 km below Slab2 """ @@ -65,47 +57,10 @@ NodalPlane, get_community_fault_model, ) -from source_modelling.magnitude_scaling import RakeType, rake_type DEFAULT_CENTROID_DEPTH_KM = 10.0 """Depth assumed when a centroid is given without a depth.""" -INTERFACE_DEPTH_TOLERANCE_KM = 10.0 -"""Vertical distance from the interface within which events may be interface events (Rollins et al. 2022).""" - -INTERFACE_RAKE_TOLERANCE = 60.0 -"""Maximum deviation of rake from pure reverse (90 degrees) for an interface-like plane.""" - -INTERFACE_MAX_DIP = 75.0 -"""Maximum dip of an interface-like plane.""" - -INTERFACE_MAX_SLAB_ANGLE = 35.0 -"""Maximum angle between a plane and the local slab surface for an interface-like plane.""" - -DEPTH_BELOW_SLAB_SIGMA_KM = 7.0 -"""Standard deviation of the depth below the interface (combined centroid and Slab2 depth uncertainty).""" - -SEISMOGENIC_ZONE_DEPTHS_KM = {"ker": (10.0, 47.0), "puy": (11.0, 30.0)} -"""Interface depths bounding the seismogenic zone of each Slab2 region (Hayes et al. 2018).""" - -SLAB_ZONE_SEARCH_RADIUS_KM = 10.0 -"""Horizontal distance from a slab zone within which an event is assigned to it.""" - -SEISMOGENIC_CRUSTAL_MAX_DEPTH_KM = 20.0 -"""Maximum depth of crustal events above the seismogenic zone (NGA-SUB rule).""" - -DOWNDIP_CRUSTAL_MAX_DEPTH_KM = 30.0 -"""Depth above which events down-dip of the seismogenic zone are always crustal (NGA-SUB rule).""" - -UNDETERMINED_CRUSTAL_MAX_DEPTH_KM = 50.0 -"""Depth splitting NGA-SUB "undetermined" events into crustal and intraslab (as in the NZGMDB).""" - -INTERFACE_MAX_DEPTH_KM = 60.0 -"""Maximum depth of interface events (NGA-SUB rule).""" - -CENTROID_DEPTH_SIGMA_KM = 7.0 -"""Standard deviation of the centroid depth for the NGA-SUB rule.""" - ON_INTERFACE_SIGMA_KM = 10.0 """Width of the Gaussian on-interface probability in the vertical distance to the interface (Haynie et al. 2026).""" @@ -153,19 +108,6 @@ class TectonicType(Enum): """Intraslab (including outer-rise) earthquake.""" -class SlabZone(Enum): - """Part of a subduction interface in the NGA-SUB (2020) classification.""" - - UPDIP = "updip" - """Interface shallower than the seismogenic zone, near the trench (region A).""" - - SEISMOGENIC = "seismogenic" - """Seismogenic zone of the interface (region B).""" - - DOWNDIP = "downdip" - """Interface deeper than the seismogenic zone (region C).""" - - def _unit_vectors(lat: npt.ArrayLike, lon: npt.ArrayLike) -> npt.NDArray[np.float64]: """Unit vectors from the Earth's centre to points on a sphere.""" lat, lon = np.radians(lat), np.radians(lon) @@ -293,24 +235,12 @@ def present(self) -> bool: return bool(np.isfinite(self.depth)) -class SlabZoneQuery(NamedTuple): - """Slab zone near a location.""" - - zone: SlabZone | None - """The slab zone, None away from all zones.""" - - depth: float - """Depth of the slab surface in km at the nearest point of the zone, NaN away from all zones.""" - - class SlabModel: """Slab2 subduction interface geometry cropped to New Zealand. Contains the Kermadec-Hikurangi ("ker") and Puysegur ("puy") regions of Slab2 (Hayes et al. 2018), interpolated bilinearly. Strike is - interpolated as a unit vector to avoid wraparound. The grid nodes are - also split into the NGA-SUB slab zones by their depth relative to the - region's seismogenic zone. + interpolated as a unit vector to avoid wraparound. Parameters ---------- @@ -318,17 +248,9 @@ class SlabModel: Mapping with keys ``_lat``, ``_lon``, ``_depth``, ``_dip`` and ``_strike`` for each region. - seismogenic_zone_depths : dict[str, tuple[float, float]], optional - Interface depths in km bounding the seismogenic zone of each region. """ - def __init__( - self, - grids: dict[str, npt.NDArray], - seismogenic_zone_depths: dict[ - str, tuple[float, float] - ] = SEISMOGENIC_ZONE_DEPTHS_KM, - ): + def __init__(self, grids: dict[str, npt.NDArray]): """Create a slab model from named grids. Parameters @@ -337,12 +259,8 @@ def __init__( Mapping with keys ``_lat``, ``_lon``, ``_depth``, ``_dip`` and ``_strike`` for each region. - seismogenic_zone_depths : dict[str, tuple[float, float]], optional - Interface depths in km bounding the seismogenic zone of each region. """ self._regions = [] - zone_points = {zone: [] for zone in SlabZone} - zone_depths = {zone: [] for zone in SlabZone} for region in sorted({key.split("_")[0] for key in grids}): lats = grids[f"{region}_lat"] lons = grids[f"{region}_lon"] @@ -357,24 +275,6 @@ def __init__( ) self._regions.append((bool(lons[-1] > 180.0), interpolator)) - points = _unit_vectors(*np.meshgrid(lats, lons, indexing="ij")) - top, bottom = seismogenic_zone_depths[region] - # NaN depths (outside the footprint) fall in no zone. - for zone, in_zone in [ - (SlabZone.UPDIP, depth < top), - (SlabZone.SEISMOGENIC, (depth >= top) & (depth <= bottom)), - (SlabZone.DOWNDIP, depth > bottom), - ]: - zone_points[zone].append(points[in_zone]) - zone_depths[zone].append(depth[in_zone]) - self._zones = { - zone: ( - sp.spatial.KDTree(np.concatenate(zone_points[zone])), - np.concatenate(zone_depths[zone]), - ) - for zone in SlabZone - } - @classmethod @functools.cache def load(cls) -> SlabModel: @@ -414,35 +314,6 @@ def query(self, lat: float, lon: float) -> SlabQuery: return SlabQuery(float(depth), float(dip), float(strike)) return SlabQuery(float("nan"), float("nan"), float("nan")) - def zone(self, lat: float, lon: float) -> SlabZoneQuery: - """NGA-SUB slab zone near a location. - - A location is in a zone if it is within `SLAB_ZONE_SEARCH_RADIUS_KM` - of one of its grid nodes. Near zone boundaries the seismogenic zone - takes precedence, then the down-dip zone. - - Parameters - ---------- - lat : float - Latitude in degrees. - lon : float - Longitude in degrees (either -180..180 or 0..360). - - Returns - ------- - SlabZoneQuery - The zone and the slab depth at its nearest node. - """ - point = _unit_vectors(lat, lon) - # Chord length on the unit sphere, equal to arc length at this scale. - radius = SLAB_ZONE_SEARCH_RADIUS_KM / geo.R_EARTH - for zone in (SlabZone.SEISMOGENIC, SlabZone.DOWNDIP, SlabZone.UPDIP): - tree, depths = self._zones[zone] - distance, index = tree.query(point, distance_upper_bound=radius) - if np.isfinite(distance): - return SlabZoneQuery(zone, float(depths[index])) - return SlabZoneQuery(None, float("nan")) - class FaultSegmentIndex: """Straight-line segments of the CFM fault traces with their attributes. @@ -581,182 +452,6 @@ def _parse_centroid(centroid: npt.ArrayLike) -> tuple[float, float, float]: return float(centroid[0]), float(centroid[1]), float(centroid[2]) -def _is_interface_like(plane: NodalPlane, slab: SlabQuery) -> bool: - """Whether a plane is thrust-like, not too steep and near parallel to the slab surface.""" - if not slab.present: - return False - # Angle between the plane and the slab surface (between their normals). - cosine = abs( - float( - _plane_normal(plane.strike, plane.dip) - @ _plane_normal(slab.strike, slab.dip) - ) - ) - slab_angle = np.degrees(np.arccos(min(cosine, 1.0))) - return bool( - _angular_difference(plane.rake, 90.0) <= INTERFACE_RAKE_TOLERANCE - and plane.dip <= INTERFACE_MAX_DIP - and slab_angle <= INTERFACE_MAX_SLAB_ANGLE - ) - - -def tectonic_type_probabilities( - depth_below_slab: float, - interface_like: bool, - normal_like: bool, - sigma: float = DEPTH_BELOW_SLAB_SIGMA_KM, -) -> dict[TectonicType, float]: - """Tectonic type probabilities by the NZ NSHM 2022 rule (Rollins et al. 2022). - - Events outside the slab footprint are crustal. Inside it, events within - `INTERFACE_DEPTH_TOLERANCE_KM` of the slab surface are interface if a - plane is interface-like and intraslab if the mechanism is normal - faulting (bending of the subducting plate, whose CMT centroids are often - too shallow). Other events are crustal above the slab surface and - intraslab below it. The rule is applied with the depth below the slab - normally distributed. - - Parameters - ---------- - depth_below_slab : float - Centroid depth below the slab surface in km (negative above it), - NaN outside the slab footprint. - interface_like : bool - Whether either nodal plane is interface-like. - normal_like : bool - Whether the mechanism is normal faulting. - sigma : float, optional - Standard deviation of the depth below the slab in km. With - ``sigma=0`` the deterministic rule is applied. - - Returns - ------- - dict[TectonicType, float] - Probability of each tectonic type. - """ - if np.isnan(depth_below_slab): - return { - TectonicType.CRUSTAL: 1.0, - TectonicType.INTERFACE: 0.0, - TectonicType.SLAB: 0.0, - } - tolerance = INTERFACE_DEPTH_TOLERANCE_KM - # Probability that the depth below the slab is above the tolerance band, - # above the slab surface, and above the bottom of the band. - if sigma > 0: - above_band, above_slab, above_band_bottom = np.asarray( - sp.stats.norm.cdf([-tolerance, 0.0, tolerance], depth_below_slab, sigma) - ).tolist() - else: - above_band = float(depth_below_slab < -tolerance) - above_slab = float(depth_below_slab <= 0.0) - above_band_bottom = float(depth_below_slab <= tolerance) - band_crustal = above_slab - above_band - band_slab = above_band_bottom - above_slab - interface = 0.0 - if interface_like: - interface, band_crustal, band_slab = band_crustal + band_slab, 0.0, 0.0 - elif normal_like: - band_crustal, band_slab = 0.0, band_crustal + band_slab - return { - TectonicType.CRUSTAL: float(above_band + band_crustal), - TectonicType.INTERFACE: float(interface), - TectonicType.SLAB: float(1.0 - above_band_bottom + band_slab), - } - - -def _nga_sub_depth_intervals( - zone: SlabZone | None, slab_depth: float -) -> list[tuple[float, TectonicType]]: - """Maximum depth of each tectonic type under the NGA-SUB rule, shallowest first.""" - tolerance = INTERFACE_DEPTH_TOLERANCE_KM - match zone: - case None: - return [ - (UNDETERMINED_CRUSTAL_MAX_DEPTH_KM, TectonicType.CRUSTAL), - (np.inf, TectonicType.SLAB), - ] - case SlabZone.UPDIP: - # Outer-rise events above 60 km, intraslab below. - return [(np.inf, TectonicType.SLAB)] - case SlabZone.SEISMOGENIC: - return [ - ( - min(slab_depth - tolerance, SEISMOGENIC_CRUSTAL_MAX_DEPTH_KM), - TectonicType.CRUSTAL, - ), - ( - min(slab_depth + tolerance, INTERFACE_MAX_DEPTH_KM), - TectonicType.INTERFACE, - ), - (np.inf, TectonicType.SLAB), - ] - case SlabZone.DOWNDIP: - # Undetermined events between the crustal and slab depths are - # split at UNDETERMINED_CRUSTAL_MAX_DEPTH_KM. - crustal_max_depth = max( - DOWNDIP_CRUSTAL_MAX_DEPTH_KM, - min(slab_depth - tolerance, UNDETERMINED_CRUSTAL_MAX_DEPTH_KM), - ) - return [ - (crustal_max_depth, TectonicType.CRUSTAL), - (np.inf, TectonicType.SLAB), - ] - - -def nga_sub_tectonic_type_probabilities( - depth: float, - zone: SlabZone | None, - slab_depth: float, - sigma: float = CENTROID_DEPTH_SIGMA_KM, -) -> dict[TectonicType, float]: - """Tectonic type probabilities by the modified NGA-SUB (2020) rule of the NZGMDB. - - The rule depends only on the location and depth, not the mechanism: - - - Up-dip zone: intraslab (outer-rise above 60 km). - - Seismogenic zone: crustal above both 20 km and 10 km above the slab; - otherwise interface above both 60 km and 10 km below the slab; - otherwise intraslab. - - Down-dip zone: crustal above 30 km, intraslab within 10 km above the - slab or deeper, and undetermined otherwise. - - Away from all zones: crustal above 30 km, intraslab below 60 km, and - undetermined otherwise. - - Undetermined events are crustal above `UNDETERMINED_CRUSTAL_MAX_DEPTH_KM` - and intraslab below, as in the NZGMDB. The rule is applied with the - centroid depth normally distributed. - - Parameters - ---------- - depth : float - Centroid depth in km. - zone : SlabZone | None - Slab zone near the centroid, None away from all zones. - slab_depth : float - Slab surface depth in km at the nearest point of the zone (unused - when `zone` is None). - sigma : float, optional - Standard deviation of the centroid depth in km. With ``sigma=0`` - the deterministic rule is applied. - - Returns - ------- - dict[TectonicType, float] - Probability of each tectonic type. - """ - intervals = _nga_sub_depth_intervals(zone, slab_depth) - bounds = np.array([-np.inf] + [max_depth for max_depth, _ in intervals]) - if sigma > 0: - cdf = sp.stats.norm.cdf(bounds, depth, sigma) - else: - cdf = (depth <= bounds).astype(float) - probabilities = dict.fromkeys(TectonicType, 0.0) - for (_, tectonic_type), probability in zip(intervals, np.diff(cdf)): - probabilities[tectonic_type] += float(probability) - return probabilities - - def subduction_slip_planes( lat: float, lon: float, slab: SlabQuery ) -> tuple[NodalPlane, NodalPlane]: @@ -796,7 +491,7 @@ def subduction_slip_planes( ) -def haynie_tectonic_type_probabilities( +def tectonic_type_probabilities( depth_above_slab: float, depth: float, subduction_kagan_angle: float | None = None, @@ -1076,127 +771,21 @@ def most_likely_nodal_plane( def tectonic_type_probabilities( self, centroid: npt.ArrayLike, - nodal_plane_1: NodalPlane, - nodal_plane_2: NodalPlane, - sigma: float = DEPTH_BELOW_SLAB_SIGMA_KM, - ) -> dict[TectonicType, float]: - """Tectonic type probabilities (see the module-level `tectonic_type_probabilities`). - - Parameters - ---------- - centroid : npt.ArrayLike - Centroid as (lat, lon, depth_km). If depth is omitted a default - is assumed with a warning. - nodal_plane_1 : NodalPlane - First nodal plane. - nodal_plane_2 : NodalPlane - Second nodal plane. - sigma : float, optional - Standard deviation of the depth below the slab in km. - - Returns - ------- - dict[TectonicType, float] - Probability of each tectonic type. - """ - lat, lon, depth = _parse_centroid(centroid) - slab = self.slab_model.query(lat, lon) - # The planes share the slip sense; the plane with more dip-slip - # decides whether the mechanism is normal or strike-slip. - dominant = max( - nodal_plane_1, nodal_plane_2, key=lambda p: abs(np.sin(np.radians(p.rake))) - ) - return tectonic_type_probabilities( - depth - slab.depth, - _is_interface_like(nodal_plane_1, slab) - or _is_interface_like(nodal_plane_2, slab), - rake_type(dominant.rake) in (RakeType.NORMAL, RakeType.NORMAL_OBLIQUE), - sigma, - ) - - def tectonic_type( - self, - centroid: npt.ArrayLike, - nodal_plane_1: NodalPlane, - nodal_plane_2: NodalPlane, - ) -> TectonicType: - """Tectonic type by the NZ NSHM 2022 rule (see `tectonic_type_probabilities`). - - Parameters - ---------- - centroid : npt.ArrayLike - Centroid as (lat, lon, depth_km). If depth is omitted a default - is assumed with a warning. - nodal_plane_1 : NodalPlane - First nodal plane. - nodal_plane_2 : NodalPlane - Second nodal plane. - - Returns - ------- - TectonicType - The tectonic type. - """ - probabilities = self.tectonic_type_probabilities( - centroid, nodal_plane_1, nodal_plane_2, sigma=0.0 - ) - return max(probabilities, key=probabilities.__getitem__) - - def nga_sub_tectonic_type_probabilities( - self, centroid: npt.ArrayLike, sigma: float = CENTROID_DEPTH_SIGMA_KM - ) -> dict[TectonicType, float]: - """Tectonic type probabilities by the NGA-SUB rule (see the module-level `nga_sub_tectonic_type_probabilities`). - - Parameters - ---------- - centroid : npt.ArrayLike - Centroid as (lat, lon, depth_km). If depth is omitted a default - is assumed with a warning. - sigma : float, optional - Standard deviation of the centroid depth in km. - - Returns - ------- - dict[TectonicType, float] - Probability of each tectonic type. - """ - lat, lon, depth = _parse_centroid(centroid) - zone, slab_depth = self.slab_model.zone(lat, lon) - return nga_sub_tectonic_type_probabilities(depth, zone, slab_depth, sigma) - - def nga_sub_tectonic_type(self, centroid: npt.ArrayLike) -> TectonicType: - """Tectonic type by the NGA-SUB rule (see `nga_sub_tectonic_type_probabilities`). - - Parameters - ---------- - centroid : npt.ArrayLike - Centroid as (lat, lon, depth_km). If depth is omitted a default - is assumed with a warning. - - Returns - ------- - TectonicType - The tectonic type. - """ - probabilities = self.nga_sub_tectonic_type_probabilities(centroid, sigma=0.0) - return max(probabilities, key=probabilities.__getitem__) - - def haynie_tectonic_type_probabilities( - self, - centroid: npt.ArrayLike, - nodal_plane: NodalPlane | None = None, + nodal_plane_1: NodalPlane | None = None, + nodal_plane_2: NodalPlane | None = None, fixed_depth: bool = False, ) -> dict[TectonicType, float]: - """Tectonic type probabilities by the Haynie et al. (2026) rule (see the module-level `haynie_tectonic_type_probabilities`). + """Tectonic type probabilities (see the module-level `tectonic_type_probabilities`). Parameters ---------- centroid : npt.ArrayLike Centroid as (lat, lon, depth_km). If depth is omitted a default is assumed with a warning and treated as fixed. - nodal_plane : NodalPlane | None, optional - Either nodal plane of the event's focal mechanism, None if it - has none. + nodal_plane_1 : NodalPlane | None, optional + First nodal plane, None without a focal mechanism. + nodal_plane_2 : NodalPlane | None, optional + Second nodal plane, None without a focal mechanism. fixed_depth : bool, optional Whether the centroid depth was fixed rather than solved for. @@ -1208,32 +797,38 @@ def haynie_tectonic_type_probabilities( fixed_depth = fixed_depth or np.size(centroid) == 2 lat, lon, depth = _parse_centroid(centroid) slab = self.slab_model.query(lat, lon) + # Both planes describe the same mechanism; using both keeps the + # result symmetric when catalogue rounding makes them disagree. + planes = [p for p in (nodal_plane_1, nodal_plane_2) if p is not None] subduction_kagan_angle = None - if nodal_plane is not None and slab.present: + if planes and slab.present: subduction_kagan_angle = min( - kagan_angle(nodal_plane, subduction_plane) + kagan_angle(plane, subduction_plane) + for plane in planes for subduction_plane in subduction_slip_planes(lat, lon, slab) ) - return haynie_tectonic_type_probabilities( + return tectonic_type_probabilities( slab.depth - depth, depth, subduction_kagan_angle, fixed_depth ) - def haynie_tectonic_type( + def tectonic_type( self, centroid: npt.ArrayLike, - nodal_plane: NodalPlane | None = None, + nodal_plane_1: NodalPlane | None = None, + nodal_plane_2: NodalPlane | None = None, fixed_depth: bool = False, ) -> TectonicType: - """Most probable tectonic type by the Haynie et al. (2026) rule (see `haynie_tectonic_type_probabilities`). + """Most probable tectonic type (see `tectonic_type_probabilities`). Parameters ---------- centroid : npt.ArrayLike Centroid as (lat, lon, depth_km). If depth is omitted a default is assumed with a warning and treated as fixed. - nodal_plane : NodalPlane | None, optional - Either nodal plane of the event's focal mechanism, None if it - has none. + nodal_plane_1 : NodalPlane | None, optional + First nodal plane, None without a focal mechanism. + nodal_plane_2 : NodalPlane | None, optional + Second nodal plane, None without a focal mechanism. fixed_depth : bool, optional Whether the centroid depth was fixed rather than solved for. @@ -1242,7 +837,7 @@ def haynie_tectonic_type( TectonicType The tectonic type. """ - probabilities = self.haynie_tectonic_type_probabilities( - centroid, nodal_plane, fixed_depth + probabilities = self.tectonic_type_probabilities( + centroid, nodal_plane_1, nodal_plane_2, fixed_depth ) return max(probabilities, key=probabilities.__getitem__) diff --git a/tests/test_focal_mechanism.py b/tests/test_focal_mechanism.py index 86653c7..b99acdb 100644 --- a/tests/test_focal_mechanism.py +++ b/tests/test_focal_mechanism.py @@ -3,7 +3,7 @@ from source_modelling import focal_mechanism from source_modelling.community_fault_model import NodalPlane -from source_modelling.focal_mechanism import SlabZone, TectonicType +from source_modelling.focal_mechanism import TectonicType # 2003 Fiordland earthquake (GeoNet 2103645), a Puysegur interface event. CENTROID = np.array([-45.1929, 166.83, 22.0]) @@ -70,104 +70,6 @@ def test_slab_query(): assert not slab_model.query(-43.5, 172.6).present # Christchurch -def test_tectonic_type(classifier: focal_mechanism.CMTClassifier): - assert ( - classifier.tectonic_type(CENTROID, PLANE_1, PLANE_2) == TectonicType.INTERFACE - ) - # Same mechanism well above the interface is crustal. - shallow = np.array([-45.1929, 166.83, 2.0]) - assert classifier.tectonic_type(shallow, PLANE_1, PLANE_2) == TectonicType.CRUSTAL - - -@pytest.mark.parametrize( - "depth_below_slab, interface_like, normal_like, expected", - [ - (float("nan"), True, False, TectonicType.CRUSTAL), - (-20.0, True, False, TectonicType.CRUSTAL), - (-5.0, True, False, TectonicType.INTERFACE), - (5.0, True, False, TectonicType.INTERFACE), - (-5.0, False, False, TectonicType.CRUSTAL), - (5.0, False, False, TectonicType.SLAB), - (-5.0, False, True, TectonicType.SLAB), - (20.0, True, False, TectonicType.SLAB), - ], -) -def test_tectonic_type_probabilities( - depth_below_slab: float, - interface_like: bool, - normal_like: bool, - expected: TectonicType, -): - rule = focal_mechanism.tectonic_type_probabilities( - depth_below_slab, interface_like, normal_like, sigma=0.0 - ) - assert rule[expected] == 1.0 - probabilities = focal_mechanism.tectonic_type_probabilities( - depth_below_slab, interface_like, normal_like - ) - assert sum(probabilities.values()) == pytest.approx(1.0) - assert max(probabilities, key=probabilities.__getitem__) == expected - - -@pytest.mark.parametrize( - "lat, lon, expected", - [ - (-41.3, 174.8, SlabZone.SEISMOGENIC), # Wellington - (-38.7, 176.1, SlabZone.DOWNDIP), # Taupo - (-41.5, 177.5, SlabZone.UPDIP), # Near the Hikurangi trough - (-43.5, 172.6, None), # Christchurch - ], -) -def test_slab_zone(lat: float, lon: float, expected: SlabZone | None): - slab_model = focal_mechanism.SlabModel.load() - query = slab_model.zone(lat, lon) - assert query.zone == expected - assert np.isfinite(query.depth) == (expected is not None) - wrapped = slab_model.zone(lat, lon - 360) - assert wrapped.zone == query.zone - np.testing.assert_equal(wrapped.depth, query.depth) - - -def test_nga_sub_tectonic_type(classifier: focal_mechanism.CMTClassifier): - assert classifier.nga_sub_tectonic_type(CENTROID) == TectonicType.INTERFACE - shallow = np.array([-41.3, 174.8, 5.0]) # Above the slab beneath Wellington - assert classifier.nga_sub_tectonic_type(shallow) == TectonicType.CRUSTAL - - -@pytest.mark.parametrize( - "depth, zone, slab_depth, expected", - [ - (5.0, None, float("nan"), TectonicType.CRUSTAL), - (45.0, None, float("nan"), TectonicType.CRUSTAL), - (55.0, None, float("nan"), TectonicType.SLAB), - (5.0, SlabZone.UPDIP, 8.0, TectonicType.SLAB), - (5.0, SlabZone.SEISMOGENIC, 25.0, TectonicType.CRUSTAL), - (18.0, SlabZone.SEISMOGENIC, 25.0, TectonicType.INTERFACE), - (22.0, SlabZone.SEISMOGENIC, 40.0, TectonicType.INTERFACE), - (40.0, SlabZone.SEISMOGENIC, 25.0, TectonicType.SLAB), - (25.0, SlabZone.DOWNDIP, 100.0, TectonicType.CRUSTAL), - (45.0, SlabZone.DOWNDIP, 100.0, TectonicType.CRUSTAL), - (55.0, SlabZone.DOWNDIP, 100.0, TectonicType.SLAB), - (40.0, SlabZone.DOWNDIP, 48.0, TectonicType.SLAB), - ], -) -def test_nga_sub_tectonic_type_probabilities( - depth: float, - zone: SlabZone | None, - slab_depth: float, - expected: TectonicType, -): - rule = focal_mechanism.nga_sub_tectonic_type_probabilities( - depth, zone, slab_depth, sigma=0.0 - ) - assert rule[expected] == 1.0 - probabilities = focal_mechanism.nga_sub_tectonic_type_probabilities( - depth, zone, slab_depth - ) - assert sum(probabilities.values()) == pytest.approx(1.0) - assert max(probabilities, key=probabilities.__getitem__) == expected - - @pytest.mark.parametrize( "plane_a, plane_b, expected", [ @@ -252,14 +154,14 @@ def test_subduction_slip_planes(): (-30.0, 50.0, 90.0, True, {TectonicType.SLAB: 1.0 - np.exp(-4.5)}), ], ) -def test_haynie_tectonic_type_probabilities( +def test_tectonic_type_probabilities( depth_above_slab: float, depth: float, kagan: float | None, fixed_depth: bool, expected: dict[TectonicType, float], ): - probabilities = focal_mechanism.haynie_tectonic_type_probabilities( + probabilities = focal_mechanism.tectonic_type_probabilities( depth_above_slab, depth, kagan, fixed_depth ) assert sum(probabilities.values()) == pytest.approx(1.0) @@ -267,38 +169,50 @@ def test_haynie_tectonic_type_probabilities( assert probabilities[tectonic_type] == pytest.approx(probability) -def test_haynie_tectonic_type(classifier: focal_mechanism.CMTClassifier): +def test_tectonic_type(classifier: focal_mechanism.CMTClassifier): # Thrust on the Hikurangi interface beneath Wellington. slab = classifier.slab_model.query(-41.3, 174.8) on_interface = np.array([-41.3, 174.8, slab.depth]) thrust = NodalPlane(slab.strike, slab.dip, 100.0) - assert classifier.haynie_tectonic_type(on_interface, thrust) == ( - TectonicType.INTERFACE - ) + assert classifier.tectonic_type(on_interface, thrust) == (TectonicType.INTERFACE) # Either nodal plane gives the same probabilities. np.testing.assert_allclose( - list( - classifier.haynie_tectonic_type_probabilities( - on_interface, PLANE_1 - ).values() - ), - list( - classifier.haynie_tectonic_type_probabilities( - on_interface, PLANE_2 - ).values() - ), + list(classifier.tectonic_type_probabilities(on_interface, PLANE_1).values()), + list(classifier.tectonic_type_probabilities(on_interface, PLANE_2).values()), atol=0.01, ) # Strike-slip near the interface is not interface slip. strike_slip = NodalPlane(40.0, 85.0, 175.0) shallower = on_interface - [0, 0, 2.0] - assert classifier.haynie_tectonic_type(shallower, strike_slip) == ( - TectonicType.CRUSTAL - ) + assert classifier.tectonic_type(shallower, strike_slip) == (TectonicType.CRUSTAL) deeper = on_interface + [0, 0, 2.0] - assert classifier.haynie_tectonic_type(deeper, strike_slip) == TectonicType.SLAB + assert classifier.tectonic_type(deeper, strike_slip) == TectonicType.SLAB # Outside the slab footprint (Christchurch). assert ( - classifier.haynie_tectonic_type([-43.5, 172.6, 10.0], thrust) - == TectonicType.CRUSTAL + classifier.tectonic_type([-43.5, 172.6, 10.0], thrust) == TectonicType.CRUSTAL ) + + +def test_tectonic_type_with_both_planes(classifier: focal_mechanism.CMTClassifier): + # The call gcmt-to-realisation makes, which must not depend on plane order. + forward = classifier.tectonic_type_probabilities(CENTROID, PLANE_1, PLANE_2) + backward = classifier.tectonic_type_probabilities(CENTROID, PLANE_2, PLANE_1) + assert forward == backward + assert sum(forward.values()) == pytest.approx(1.0) + # 9 km below the Slab2 surface, so mostly intraslab. + assert classifier.tectonic_type(CENTROID, PLANE_1, PLANE_2) == TectonicType.SLAB + + +def test_tectonic_type_without_mechanism_or_depth( + classifier: focal_mechanism.CMTClassifier, +): + slab = classifier.slab_model.query(-41.3, 174.8) + on_interface = np.array([-41.3, 174.8, slab.depth]) + assert classifier.tectonic_type(on_interface) == TectonicType.INTERFACE + # A fixed depth leaves the interface probability to the mechanism. + deep = on_interface + [0, 0, 30.0] + assert classifier.tectonic_type(deep) == TectonicType.SLAB + assert classifier.tectonic_type(deep, fixed_depth=True) == TectonicType.INTERFACE + with pytest.warns(UserWarning, match="without depth"): + probabilities = classifier.tectonic_type_probabilities(on_interface[:2]) + assert probabilities[TectonicType.INTERFACE] == 1.0 From c3750bcff7ff76cccb197a7db028864386ad6470 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 27 Oct 2026 20:09:47 +1300 Subject: [PATCH 3/6] Build random nodal planes from Python floats so ty can type them Co-Authored-By: Claude Opus 5.5 --- tests/test_focal_mechanism.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/tests/test_focal_mechanism.py b/tests/test_focal_mechanism.py index b99acdb..9c7af14 100644 --- a/tests/test_focal_mechanism.py +++ b/tests/test_focal_mechanism.py @@ -91,10 +91,14 @@ def test_kagan_angle(plane_a: NodalPlane, plane_b: NodalPlane, expected: float): def test_kagan_angle_is_at_most_120(): rng = np.random.default_rng(0) - planes = rng.uniform([0, 0, -180], [360, 90, 180], size=(500, 3)) + planes = [ + NodalPlane(strike, dip, rake) + for strike, dip, rake in rng.uniform( + [0, 0, -180], [360, 90, 180], size=(500, 3) + ).tolist() + ] angles = [ - focal_mechanism.kagan_angle(NodalPlane(*a), NodalPlane(*b)) - for a, b in zip(planes[::2], planes[1::2]) + focal_mechanism.kagan_angle(a, b) for a, b in zip(planes[::2], planes[1::2]) ] assert 0.0 <= min(angles) and max(angles) <= 120.0 From 7ae0e2aab85a339d01aef3d83377ba54db1e86f2 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 27 Oct 2026 20:28:54 +1300 Subject: [PATCH 4/6] Follow CatSep (Haynie et al. 2026) as published, without the Moho grid Replace the Gaussian interface likelihoods and plate-motion reference slip with CatSep v1.0.0's ramp functions: the interface score is the flex-zone score (within 15 km of Slab2), times a ramp from 1 at the seismogenic zone depth to 0 15 km below, times a Kagan ramp (1 at 30 degrees to 0.25 at 53) against slip at the region's seismogenic rake. Crustal is the remainder times a ramp from 20 km above to 20 km below the slab, and intraslab takes the rest. The CRUST1.0 Moho is not used: moho_depth defaults to infinity, so the crust extends to the slab and the mantle-wedge branch never applies. SlabQuery now carries its Slab2 region to look up the seismogenic zone (Hayes et al. 2018; Puysegur 30 km, 90 degrees). Co-Authored-By: Claude Opus 5.5 --- source_modelling/focal_mechanism.py | 298 +++++++++++++--------------- tests/test_focal_mechanism.py | 156 ++++++--------- 2 files changed, 198 insertions(+), 256 deletions(-) diff --git a/source_modelling/focal_mechanism.py b/source_modelling/focal_mechanism.py index 7b066fe..2af0e1c 100644 --- a/source_modelling/focal_mechanism.py +++ b/source_modelling/focal_mechanism.py @@ -13,13 +13,14 @@ ``scripts/train_cmt_classifiers.py``. The tectonic type (crustal, subduction interface or intraslab) follows the -probabilistic rule of Haynie et al. (2026) against the Slab2 interface -geometry (Hayes et al. 2018). The interface probability is the product of -Gaussian likelihoods in the vertical distance to the interface and in the -Kagan angle between the mechanism and interface slip (pure reverse, or -along the relative plate motion), for events no deeper than 70 km. The -remainder is crustal (upper plate) above the interface and intraslab below -it. +CatSep plausibility scores of Haynie et al. (2026) against the Slab2 +interface geometry and seismogenic zone depths (Hayes et al. 2018). The +interface score is the product of ramp functions in the vertical distance +to the interface, the depth below the seismogenic zone and the Kagan angle +between the mechanism and interface slip. The crustal score ramps down +from 20 km above to 20 km below the interface, and intraslab takes the +remainder. CatSep's CRUST1.0 Moho is not used: the crust is taken to +extend down to the slab. Examples -------- @@ -32,7 +33,7 @@ >>> classifier.most_likely_nodal_plane(centroid, nodal_plane_1, nodal_plane_2) NodalPlane(strike=20, dip=35, rake=79) >>> classifier.tectonic_type(centroid, nodal_plane_1, nodal_plane_2) - + """ from __future__ import annotations @@ -61,17 +62,62 @@ DEFAULT_CENTROID_DEPTH_KM = 10.0 """Depth assumed when a centroid is given without a depth.""" -ON_INTERFACE_SIGMA_KM = 10.0 -"""Width of the Gaussian on-interface probability in the vertical distance to the interface (Haynie et al. 2026).""" +FLEX_ZONE_KM = 15.0 +"""Half-width of the interface flex zone about the slab surface, and the buffer below the seismogenic zone (Haynie et al. 2026).""" -SUBDUCTION_SLIP_SIGMA = 30.0 -"""Width in degrees of the Gaussian subduction slip probability in the Kagan angle (Haynie et al. 2026).""" -SUBDUCTION_MAX_DEPTH_KM = 70.0 -"""Maximum centroid depth of subduction interface events (Haynie et al. 2026).""" +class Ramp(NamedTuple): + """Ramp function: `p1` at or below `x1`, `p2` at or above `x2`, linear between (Haynie et al. 2026, eq. 1).""" -PACIFIC_AUSTRALIA_EULER_POLE = (-60.1, -178.3, 1.0744) -"""Latitude, longitude (degrees) and rate (degrees/Myr) of Pacific rotation relative to Australia (NUVEL-1A, DeMets et al. 1994).""" + x1: float + """Start of the ramp.""" + + p1: float + """Score at or below `x1`.""" + + x2: float + """End of the ramp.""" + + p2: float + """Score at or above `x2`.""" + + def __call__(self, x: float) -> float: + """Score at `x`.""" + if x <= self.x1: + return self.p1 + if x >= self.x2: + return self.p2 + return self.p1 + (self.p2 - self.p1) * (x - self.x1) / (self.x2 - self.x1) + + +KAGAN_RAMP = Ramp(30.0, 1.0, 53.0, 0.25) +"""Interface score in the Kagan angle (degrees) to interface slip (Haynie et al. 2026, Table 1).""" + +CRUSTAL_DEPTH_DIFFERENCE_RAMP = Ramp(-20.0, 1.0, 20.0, 0.0) +"""Crustal score in the depth (km) below the slab surface (Haynie et al. 2026, Table 1).""" + +MANTLE_WEDGE_FLEX_SCORES = (0.0, 0.4) +"""Interface flex scores at the Moho and the top of the flex zone for events in the mantle wedge (Haynie et al. 2026, Table 1).""" + + +class SeismogenicZone(NamedTuple): + """Seismogenic zone of a Slab2 region.""" + + depth: float + """Depth in km of the deep limit of the seismogenic zone.""" + + rake: float + """Typical rake of interface slip in degrees.""" + + +SEISMOGENIC_ZONES = { + "ker": SeismogenicZone(47.0, 93.0), + "puy": SeismogenicZone(30.0, 90.0), +} +"""Seismogenic zone of each Slab2 region (Hayes et al. 2018). Puysegur has no published rake, so it takes CatSep's default.""" + +DEFAULT_SEISMOGENIC_ZONE = SeismogenicZone(40.0, 90.0) +"""Seismogenic zone of Slab2 regions without one in `SEISMOGENIC_ZONES` (CatSep's default).""" FAULT_NEIGHBOURS = 8 """Number of nearest CFM fault segments used for the strike, dip and rake misfits.""" @@ -108,14 +154,6 @@ class TectonicType(Enum): """Intraslab (including outer-rise) earthquake.""" -def _unit_vectors(lat: npt.ArrayLike, lon: npt.ArrayLike) -> npt.NDArray[np.float64]: - """Unit vectors from the Earth's centre to points on a sphere.""" - lat, lon = np.radians(lat), np.radians(lon) - return np.stack( - [np.cos(lat) * np.cos(lon), np.cos(lat) * np.sin(lon), np.sin(lat)], axis=-1 - ) - - def _angular_difference( a: npt.ArrayLike, b: npt.ArrayLike, period: float = 360.0 ) -> npt.NDArray[np.float64]: @@ -182,41 +220,6 @@ def kagan_angle(plane_a: NodalPlane, plane_b: NodalPlane) -> float: return float(np.degrees(np.arccos(np.clip((trace - 1.0) / 2.0, -1.0, 1.0)))) -def plate_motion_azimuth( - lat: float, - lon: float, - euler_pole: tuple[float, float, float] = PACIFIC_AUSTRALIA_EULER_POLE, -) -> float: - """Azimuth of relative plate motion at a location. - - Parameters - ---------- - lat : float - Latitude in degrees. - lon : float - Longitude in degrees. - euler_pole : tuple[float, float, float], optional - Latitude, longitude (degrees) and rate of rotation of one plate - relative to the other. Defaults to Pacific relative to Australia. - - Returns - ------- - float - Azimuth in degrees of the rotating plate's motion relative to the - fixed plate. - """ - pole_lat, pole_lon, rate = euler_pole - velocity = np.cross( - rate * _unit_vectors(pole_lat, pole_lon), _unit_vectors(lat, lon) - ) - lat, lon = np.radians(lat), np.radians(lon) - north = np.array( - [-np.sin(lat) * np.cos(lon), -np.sin(lat) * np.sin(lon), np.cos(lat)] - ) - east = np.array([-np.sin(lon), np.cos(lon), 0.0]) - return float(np.degrees(np.arctan2(velocity @ east, velocity @ north)) % 360.0) - - class SlabQuery(NamedTuple): """Slab surface properties at a location.""" @@ -229,6 +232,14 @@ class SlabQuery(NamedTuple): strike: float """Strike of the slab surface in degrees, NaN outside the footprint.""" + region: str | None = None + """Slab2 region code ("ker" or "puy"), None outside the footprint.""" + + @property + def seismogenic_zone(self) -> SeismogenicZone: + """Seismogenic zone of the slab region.""" + return SEISMOGENIC_ZONES.get(self.region or "", DEFAULT_SEISMOGENIC_ZONE) + @property def present(self) -> bool: """Whether the location is inside a slab footprint.""" @@ -273,7 +284,7 @@ def __init__(self, grids: dict[str, npt.NDArray]): interpolator = sp.interpolate.RegularGridInterpolator( (lats, lons), values, bounds_error=False, fill_value=np.nan ) - self._regions.append((bool(lons[-1] > 180.0), interpolator)) + self._regions.append((region, bool(lons[-1] > 180.0), interpolator)) @classmethod @functools.cache @@ -306,12 +317,12 @@ def query(self, lat: float, lon: float) -> SlabQuery: SlabQuery Slab properties, NaN-valued outside the slab footprints. """ - for positive_lon, interpolator in self._regions: + for region, positive_lon, interpolator in self._regions: query_lon = lon % 360.0 if positive_lon else ((lon + 180.0) % 360.0) - 180.0 depth, dip, sin_strike, cos_strike = interpolator([lat, query_lon])[0] if np.isfinite(depth): strike = np.degrees(np.arctan2(sin_strike, cos_strike)) % 360.0 - return SlabQuery(float(depth), float(dip), float(strike)) + return SlabQuery(float(depth), float(dip), float(strike), region) return SlabQuery(float("nan"), float("nan"), float("nan")) @@ -452,116 +463,88 @@ def _parse_centroid(centroid: npt.ArrayLike) -> tuple[float, float, float]: return float(centroid[0]), float(centroid[1]), float(centroid[2]) -def subduction_slip_planes( - lat: float, lon: float, slab: SlabQuery -) -> tuple[NodalPlane, NodalPlane]: - """Mechanisms of slip on the subduction interface at a location. - - Parameters - ---------- - lat : float - Latitude in degrees. - lon : float - Longitude in degrees. - slab : SlabQuery - Slab surface beneath the location. - - Returns - ------- - tuple[NodalPlane, NodalPlane] - Slip on the slab surface (a) in pure reverse (rake 90) and (b) in - oblique reverse with the horizontal slip parallel to the relative - Pacific-Australia plate motion. - """ - # Components of the plate motion along strike and down the dip direction. - relative_azimuth = np.radians(plate_motion_azimuth(lat, lon) - slab.strike) - along_strike, down_dip = np.cos(relative_azimuth), np.sin(relative_azimuth) - # Take the sense of the motion that is reverse slip, whichever plate - # is subducting. - if down_dip > 0: - along_strike, down_dip = -along_strike, -down_dip - # The horizontal slip of a plane is cos(rake) along strike and - # -sin(rake) cos(dip) down the dip direction. - rake = np.degrees( - np.arctan2(-down_dip / np.cos(np.radians(slab.dip)), along_strike) - ) - return ( - NodalPlane(slab.strike, slab.dip, 90.0), - NodalPlane(slab.strike, slab.dip, float(rake)), - ) - - def tectonic_type_probabilities( - depth_above_slab: float, depth: float, + slab_depth: float, + seismogenic_depth: float = DEFAULT_SEISMOGENIC_ZONE.depth, subduction_kagan_angle: float | None = None, - fixed_depth: bool = False, - on_interface_sigma: float = ON_INTERFACE_SIGMA_KM, - subduction_slip_sigma: float = SUBDUCTION_SLIP_SIGMA, + moho_depth: float = np.inf, ) -> dict[TectonicType, float]: - """Tectonic type probabilities by the rule of Haynie et al. (2026). + """Tectonic type plausibility scores by CatSep (Haynie et al. 2026). - Events outside the slab footprint are crustal (upper plate). Inside it, - the subduction interface probability is the product of + Events outside the slab footprint are crustal. Inside it, the interface + score is the product (eq. 3) of - - the on-interface probability, ``exp(-dz**2 / (2 * on_interface_sigma**2))`` - for vertical distance ``dz`` to the interface (1 for a fixed depth), - - the subduction slip probability, - ``exp(-kagan**2 / (2 * subduction_slip_sigma**2))`` for the Kagan - angle to interface slip (1 without a focal mechanism), and - - whether the event is no deeper than `SUBDUCTION_MAX_DEPTH_KM`. + - the flex score: 1 within `FLEX_ZONE_KM` of the slab surface, if no + deeper than `FLEX_ZONE_KM` below the seismogenic zone or above the + Moho; ramping from 0 at the Moho to 0.4 at the top of the flex zone + in the mantle wedge; 0 otherwise, + - the seismogenic score, ramping from 1 at the seismogenic zone depth + to 0 at `FLEX_ZONE_KM` below it, and + - the Kagan score (`KAGAN_RAMP`), when the mechanism is known. - The remaining probability is upper plate for events above the - interface and intraslab for events below it. + The crustal score (eq. 4) is the remainder scaled by + `CRUSTAL_DEPTH_DIFFERENCE_RAMP` and by a ramp from 1 at the Moho to 0 + at the slab surface. Intraslab takes what is left (eq. 5). Parameters ---------- - depth_above_slab : float - Vertical distance in km from the centroid up to the slab surface - (negative below it), NaN outside the slab footprint. depth : float Centroid depth in km. + slab_depth : float + Depth of the slab surface in km beneath the centroid, NaN outside + the slab footprint. + seismogenic_depth : float, optional + Depth in km of the deep limit of the seismogenic zone. subduction_kagan_angle : float | None, optional - Minimum Kagan angle in degrees between the event's mechanism and - slip on the interface (see `subduction_slip_planes`), None without - a focal mechanism. - fixed_depth : bool, optional - Whether the centroid depth was fixed rather than solved for, so says - nothing about the distance to the interface. - on_interface_sigma : float, optional - Width in km of the on-interface probability. - subduction_slip_sigma : float, optional - Width in degrees of the subduction slip probability. + Kagan angle in degrees between the event's mechanism and interface + slip, None without a focal mechanism. + moho_depth : float, optional + Moho depth in km. The default puts the Moho below the slab, so the + crust extends down to the slab surface and there is no mantle wedge. Returns ------- dict[TectonicType, float] - Probability of each tectonic type. + Score of each tectonic type, summing to 1. """ - if np.isnan(depth_above_slab): + if np.isnan(slab_depth): return { TectonicType.CRUSTAL: 1.0, TectonicType.INTERFACE: 0.0, TectonicType.SLAB: 0.0, } - on_interface = ( - 1.0 - if fixed_depth - else np.exp(-(depth_above_slab**2) / (2.0 * on_interface_sigma**2)) + flex_top = slab_depth - FLEX_ZONE_KM + in_flex_zone = flex_top <= depth <= slab_depth + FLEX_ZONE_KM + if in_flex_zone and ( + depth <= seismogenic_depth + FLEX_ZONE_KM or depth <= moho_depth + ): + flex = 1.0 + elif moho_depth <= depth < flex_top: + flex = Ramp( + moho_depth, + MANTLE_WEDGE_FLEX_SCORES[0], + flex_top, + MANTLE_WEDGE_FLEX_SCORES[1], + )(depth) + else: + flex = 0.0 + seismogenic = Ramp(seismogenic_depth, 1.0, seismogenic_depth + FLEX_ZONE_KM, 0.0)( + depth ) - subduction_slip = ( - 1.0 - if subduction_kagan_angle is None - else np.exp(-(subduction_kagan_angle**2) / (2.0 * subduction_slip_sigma**2)) + kagan = ( + 1.0 if subduction_kagan_angle is None else KAGAN_RAMP(subduction_kagan_angle) ) - interface = float( - on_interface * subduction_slip * (depth <= SUBDUCTION_MAX_DEPTH_KM) + interface = flex * seismogenic * kagan + crustal = ( + (1.0 - interface) + * CRUSTAL_DEPTH_DIFFERENCE_RAMP(depth - slab_depth) + * Ramp(moho_depth, 1.0, slab_depth, 0.0)(depth) ) - above_slab = depth_above_slab >= 0 return { - TectonicType.CRUSTAL: (1.0 - interface) * above_slab, + TectonicType.CRUSTAL: crustal, TectonicType.INTERFACE: interface, - TectonicType.SLAB: (1.0 - interface) * (not above_slab), + TectonicType.SLAB: 1.0 - interface - crustal, } @@ -773,42 +756,38 @@ def tectonic_type_probabilities( centroid: npt.ArrayLike, nodal_plane_1: NodalPlane | None = None, nodal_plane_2: NodalPlane | None = None, - fixed_depth: bool = False, ) -> dict[TectonicType, float]: - """Tectonic type probabilities (see the module-level `tectonic_type_probabilities`). + """Tectonic type scores (see the module-level `tectonic_type_probabilities`). Parameters ---------- centroid : npt.ArrayLike Centroid as (lat, lon, depth_km). If depth is omitted a default - is assumed with a warning and treated as fixed. + is assumed with a warning. nodal_plane_1 : NodalPlane | None, optional First nodal plane, None without a focal mechanism. nodal_plane_2 : NodalPlane | None, optional Second nodal plane, None without a focal mechanism. - fixed_depth : bool, optional - Whether the centroid depth was fixed rather than solved for. Returns ------- dict[TectonicType, float] - Probability of each tectonic type. + Score of each tectonic type. """ - fixed_depth = fixed_depth or np.size(centroid) == 2 lat, lon, depth = _parse_centroid(centroid) slab = self.slab_model.query(lat, lon) + zone = slab.seismogenic_zone # Both planes describe the same mechanism; using both keeps the # result symmetric when catalogue rounding makes them disagree. planes = [p for p in (nodal_plane_1, nodal_plane_2) if p is not None] subduction_kagan_angle = None if planes and slab.present: + interface_slip = NodalPlane(slab.strike, slab.dip, zone.rake) subduction_kagan_angle = min( - kagan_angle(plane, subduction_plane) - for plane in planes - for subduction_plane in subduction_slip_planes(lat, lon, slab) + kagan_angle(plane, interface_slip) for plane in planes ) return tectonic_type_probabilities( - slab.depth - depth, depth, subduction_kagan_angle, fixed_depth + depth, slab.depth, zone.depth, subduction_kagan_angle ) def tectonic_type( @@ -816,21 +795,18 @@ def tectonic_type( centroid: npt.ArrayLike, nodal_plane_1: NodalPlane | None = None, nodal_plane_2: NodalPlane | None = None, - fixed_depth: bool = False, ) -> TectonicType: - """Most probable tectonic type (see `tectonic_type_probabilities`). + """Most plausible tectonic type (see `tectonic_type_probabilities`). Parameters ---------- centroid : npt.ArrayLike Centroid as (lat, lon, depth_km). If depth is omitted a default - is assumed with a warning and treated as fixed. + is assumed with a warning. nodal_plane_1 : NodalPlane | None, optional First nodal plane, None without a focal mechanism. nodal_plane_2 : NodalPlane | None, optional Second nodal plane, None without a focal mechanism. - fixed_depth : bool, optional - Whether the centroid depth was fixed rather than solved for. Returns ------- @@ -838,6 +814,6 @@ def tectonic_type( The tectonic type. """ probabilities = self.tectonic_type_probabilities( - centroid, nodal_plane_1, nodal_plane_2, fixed_depth + centroid, nodal_plane_1, nodal_plane_2 ) return max(probabilities, key=probabilities.__getitem__) diff --git a/tests/test_focal_mechanism.py b/tests/test_focal_mechanism.py index 9c7af14..1e85588 100644 --- a/tests/test_focal_mechanism.py +++ b/tests/test_focal_mechanism.py @@ -103,120 +103,86 @@ def test_kagan_angle_is_at_most_120(): assert 0.0 <= min(angles) and max(angles) <= 120.0 -def test_plate_motion_azimuth(): - # Pacific moves west relative to Australia at the Hikurangi margin. - assert 250 < focal_mechanism.plate_motion_azimuth(-38.7, 178.5) < 280 - assert 230 < focal_mechanism.plate_motion_azimuth(-47.0, 165.5) < 250 - - -def test_subduction_slip_planes(): - slab_model = focal_mechanism.SlabModel.load() - for lat, lon in [(-41.3, 174.8), (-38.7, 177.5), (-46.0, 166.0)]: - slab = slab_model.query(lat, lon) - pure, oblique = focal_mechanism.subduction_slip_planes(lat, lon, slab) - assert pure == NodalPlane(slab.strike, slab.dip, 90.0) - assert (oblique.strike, oblique.dip) == (slab.strike, slab.dip) - assert 0 < oblique.rake < 180 # Reverse sense - # The horizontal slip is parallel to the plate motion. - slip = focal_mechanism._slip_vector(*oblique) - slip_azimuth = np.degrees(np.arctan2(slip[1], slip[0])) - difference = focal_mechanism._angular_difference( - slip_azimuth, focal_mechanism.plate_motion_azimuth(lat, lon), 180.0 - ) - assert difference == pytest.approx(0.0, abs=1e-6) +def test_ramp(): + decreasing = focal_mechanism.Ramp(30.0, 1.0, 53.0, 0.25) + assert decreasing(0.0) == 1.0 + assert decreasing(30.0) == 1.0 + assert decreasing(41.5) == pytest.approx(0.625) + assert decreasing(53.0) == 0.25 + assert decreasing(120.0) == 0.25 + increasing = focal_mechanism.Ramp(5.0, 0.0, 15.0, 0.4) + assert increasing(10.0) == pytest.approx(0.2) @pytest.mark.parametrize( - "depth_above_slab, depth, kagan, fixed_depth, expected", + "depth, slab_depth, kagan, moho_depth, expected", [ - (float("nan"), 10.0, 0.0, False, {TectonicType.CRUSTAL: 1.0}), - (0.0, 20.0, 0.0, False, {TectonicType.INTERFACE: 1.0}), - (0.0, 20.0, None, False, {TectonicType.INTERFACE: 1.0}), - (-1.0, 80.0, 0.0, False, {TectonicType.SLAB: 1.0}), # Too deep - ( - 10.0, - 20.0, - 30.0, - False, - { - TectonicType.CRUSTAL: 1.0 - np.exp(-1.0), - TectonicType.INTERFACE: np.exp(-1.0), - }, - ), - ( - -10.0, - 40.0, - None, - False, - { - TectonicType.SLAB: 1.0 - np.exp(-0.5), - TectonicType.INTERFACE: np.exp(-0.5), - }, - ), - # A fixed depth says nothing about the distance to the interface. - (-30.0, 50.0, 0.0, True, {TectonicType.INTERFACE: 1.0}), - (-30.0, 50.0, 90.0, True, {TectonicType.SLAB: 1.0 - np.exp(-4.5)}), + # Outside the slab footprint. + (10.0, float("nan"), None, np.inf, (1.0, 0.0, 0.0)), + # On the interface, above the seismogenic limit. + (20.0, 20.0, None, np.inf, (0.0, 1.0, 0.0)), + (20.0, 20.0, 41.5, np.inf, (0.1875, 0.625, 0.1875)), + (20.0, 20.0, 90.0, np.inf, (0.375, 0.25, 0.375)), + # Halfway down the ramp below the seismogenic zone. + (37.5, 30.0, None, np.inf, (0.15625, 0.5, 0.34375)), + # Above and below the flex zone. + (5.0, 30.0, None, np.inf, (1.0, 0.0, 0.0)), + (60.0, 30.0, None, np.inf, (0.0, 0.0, 1.0)), + # In the mantle wedge, between the Moho and the top of the flex zone. + (10.0, 30.0, None, 5.0, (0.64, 0.2, 0.16)), ], ) def test_tectonic_type_probabilities( - depth_above_slab: float, depth: float, + slab_depth: float, kagan: float | None, - fixed_depth: bool, - expected: dict[TectonicType, float], + moho_depth: float, + expected: tuple[float, float, float], ): probabilities = focal_mechanism.tectonic_type_probabilities( - depth_above_slab, depth, kagan, fixed_depth + depth, slab_depth, 30.0, kagan, moho_depth ) - assert sum(probabilities.values()) == pytest.approx(1.0) - for tectonic_type, probability in expected.items(): - assert probabilities[tectonic_type] == pytest.approx(probability) + assert [ + probabilities[TectonicType.CRUSTAL], + probabilities[TectonicType.INTERFACE], + probabilities[TectonicType.SLAB], + ] == pytest.approx(expected) + + +def test_seismogenic_zone(): + slab_model = focal_mechanism.SlabModel.load() + assert slab_model.query(-41.3, 174.8).region == "ker" # Wellington + assert slab_model.query(*CENTROID[:2]).region == "puy" + assert slab_model.query(-41.3, 174.8).seismogenic_zone.depth == 47.0 + assert slab_model.query(*CENTROID[:2]).seismogenic_zone.depth == 30.0 + outside = slab_model.query(-43.5, 172.6) # Christchurch + assert outside.seismogenic_zone == focal_mechanism.DEFAULT_SEISMOGENIC_ZONE def test_tectonic_type(classifier: focal_mechanism.CMTClassifier): - # Thrust on the Hikurangi interface beneath Wellington. + # The 2003 Fiordland earthquake ruptured the Puysegur interface, and the + # result does not depend on the order of the planes. + assert classifier.tectonic_type(CENTROID, PLANE_1, PLANE_2) == ( + TectonicType.INTERFACE + ) + assert classifier.tectonic_type_probabilities( + CENTROID, PLANE_1, PLANE_2 + ) == classifier.tectonic_type_probabilities(CENTROID, PLANE_2, PLANE_1) + + # Thrusting on the Hikurangi interface beneath Wellington. slab = classifier.slab_model.query(-41.3, 174.8) on_interface = np.array([-41.3, 174.8, slab.depth]) - thrust = NodalPlane(slab.strike, slab.dip, 100.0) - assert classifier.tectonic_type(on_interface, thrust) == (TectonicType.INTERFACE) - # Either nodal plane gives the same probabilities. - np.testing.assert_allclose( - list(classifier.tectonic_type_probabilities(on_interface, PLANE_1).values()), - list(classifier.tectonic_type_probabilities(on_interface, PLANE_2).values()), - atol=0.01, - ) - # Strike-slip near the interface is not interface slip. + thrust = NodalPlane(slab.strike, slab.dip, 95.0) + assert classifier.tectonic_type(on_interface, thrust) == TectonicType.INTERFACE + # Strike-slip near the interface is crustal above it and intraslab below. strike_slip = NodalPlane(40.0, 85.0, 175.0) - shallower = on_interface - [0, 0, 2.0] - assert classifier.tectonic_type(shallower, strike_slip) == (TectonicType.CRUSTAL) - deeper = on_interface + [0, 0, 2.0] - assert classifier.tectonic_type(deeper, strike_slip) == TectonicType.SLAB + above = on_interface - [0, 0, 10.0] + below = on_interface + [0, 0, 10.0] + assert classifier.tectonic_type(above, strike_slip) == TectonicType.CRUSTAL + assert classifier.tectonic_type(below, strike_slip) == TectonicType.SLAB + # Without a mechanism, location alone decides. + assert classifier.tectonic_type(on_interface) == TectonicType.INTERFACE # Outside the slab footprint (Christchurch). assert ( classifier.tectonic_type([-43.5, 172.6, 10.0], thrust) == TectonicType.CRUSTAL ) - - -def test_tectonic_type_with_both_planes(classifier: focal_mechanism.CMTClassifier): - # The call gcmt-to-realisation makes, which must not depend on plane order. - forward = classifier.tectonic_type_probabilities(CENTROID, PLANE_1, PLANE_2) - backward = classifier.tectonic_type_probabilities(CENTROID, PLANE_2, PLANE_1) - assert forward == backward - assert sum(forward.values()) == pytest.approx(1.0) - # 9 km below the Slab2 surface, so mostly intraslab. - assert classifier.tectonic_type(CENTROID, PLANE_1, PLANE_2) == TectonicType.SLAB - - -def test_tectonic_type_without_mechanism_or_depth( - classifier: focal_mechanism.CMTClassifier, -): - slab = classifier.slab_model.query(-41.3, 174.8) - on_interface = np.array([-41.3, 174.8, slab.depth]) - assert classifier.tectonic_type(on_interface) == TectonicType.INTERFACE - # A fixed depth leaves the interface probability to the mechanism. - deep = on_interface + [0, 0, 30.0] - assert classifier.tectonic_type(deep) == TectonicType.SLAB - assert classifier.tectonic_type(deep, fixed_depth=True) == TectonicType.INTERFACE - with pytest.warns(UserWarning, match="without depth"): - probabilities = classifier.tectonic_type_probabilities(on_interface[:2]) - assert probabilities[TectonicType.INTERFACE] == 1.0 From 9004ce93a6dc696dd5c3c7a9f248e01c97bd6599 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 27 Oct 2026 22:45:10 +1300 Subject: [PATCH 5/6] Document Ramp.__call__ parameters and return for npdlint Co-Authored-By: Claude Opus 5.5 --- source_modelling/focal_mechanism.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/source_modelling/focal_mechanism.py b/source_modelling/focal_mechanism.py index 2af0e1c..1de29e2 100644 --- a/source_modelling/focal_mechanism.py +++ b/source_modelling/focal_mechanism.py @@ -82,7 +82,18 @@ class Ramp(NamedTuple): """Score at or above `x2`.""" def __call__(self, x: float) -> float: - """Score at `x`.""" + """Score at a value. + + Parameters + ---------- + x : float + The value to score. + + Returns + ------- + float + The score at `x`. + """ if x <= self.x1: return self.p1 if x >= self.x2: From 582658e5e418bca6f44a37adadc02e527aae1cda Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 27 Oct 2026 22:54:24 +1300 Subject: [PATCH 6/6] Use the CRUST1.0 Moho and simplify the CatSep classifier Package CatSep's CRUST1.0 depth-to-Moho grid cropped to New Zealand (crust1_moho_nz.npz, 2.7 kB) behind a MohoModel, and pass the Moho depth at the centroid from CMTClassifier, completing CatSep's p_moho and mantle-wedge flex terms. Simplify: the ramps are np.interp calls on (xs, ps) constants instead of a Ramp class; the Kagan angle uses scipy's D2 rotation group instead of enumerating the double-couple symmetries by hand; SlabModel takes the seismogenic zones again and SlabQuery carries its zone instead of a region string looked up in a global; the redundant depth conditions on the flex score are dropped. Scores match the previous implementation to 1e-15 over 200,000 random inputs and Kagan angles to 1e-12 degrees. Co-Authored-By: Claude Opus 5.5 --- source_modelling/NZ_CFM/crust1_moho_nz.npz | Bin 0 -> 2717 bytes source_modelling/focal_mechanism.py | 251 ++++++++++++--------- tests/test_focal_mechanism.py | 34 ++- 3 files changed, 163 insertions(+), 122 deletions(-) create mode 100644 source_modelling/NZ_CFM/crust1_moho_nz.npz diff --git a/source_modelling/NZ_CFM/crust1_moho_nz.npz b/source_modelling/NZ_CFM/crust1_moho_nz.npz new file mode 100644 index 0000000000000000000000000000000000000000..b3e63080dbfeb405d3049b5db22755aa1a9c398f GIT binary patch literal 2717 zcmZ9Oc{CK-AHYXujAh8P49_w|rG&AR-55zPi9|7B7)B$Rtl6?=PogYY8YV@OiScZa ztW&5#DDld^j5L`U+n~XZSN+cEeeb#Fd(OSz@4e@9zW;ozEx3WI004k@XYT{Vm2;H0 ze+C#J2Ee;sQ}gx>jXtzTO|D@b=BI2CypKk&}LJuH9p zpx~a%d7I1}38G78I=@J;1@*8MC|SZ|KeTAiy433g{XBjC8vbc+HJ*egrVxlLQU*wg zMAAVNZzLT^2}0sQ6dV!{q+CamK$I&;5|HALQ~^=kkSah*IFbsY;E_}yB}zzBJ)CP& z2gv~@F#s%ac?37*FQGjE>)(i4?$S8F5&3xkBI4!xfs*=@D8EyWzYPh9pQv1cs7>kL zO!b*kq|l29d;H#Bb~3S=vt> zy&7(W)7EHiN`==K(r9bAwN#py{o8C#lo?CSY!ZDEbm7!{>%X`ytKTKWa7)NuE9^6s zNL33#SCkS0N>yFOQ-n2jm4$L_?HD7HYRtK0mg0YLYt>PeN z7rA3ilEo4VL$ECZ#p1(ew@fJ^kF<1@wXP#4*Ky60#Am$M|D?>`4Si43}UzS;N zK+YhVk>4U>i=D7yKA7fUT_?#Hk8I`Hm}HqtCBY^Yf3ALL#=Ntk)^M3Zf<*2_8rHQF z)0Y-Mb9ork%6Pwes1^ zgUtSTBfg2^oLeou2ttYM=UNL+roJ z(sr~dw)9l19FJTWWvhsB7&D|;ul40?T(gy7@Tr0`8xO}HgBYqw_GNSHi{0Z7&#k_) zFlgOMIgJ&Q0>-JXCWy6;gr zTGWqj<*;x%Pt-n%)HWkus`pxiXQk!jwCZZC;(@?v)rIWmnHjP%N5<)|>Rv8-T?_+B z=skV!xH)AIb^^DzOBMjO5 zT&`DEfOB5#&lZRroOC64An{NR+bkRUwk3iF)5W2W%SaT%)kj~+OpOVCG74F^d)>`j z-ER1c^F$dzMwnE)sX`XI@xpt}m^VLGx5JO+dNJ)p{)}wjd7SN0`sVw!=9Rs=gG*T& zI2V*)jl9IFxeQ@O4f8b{(ca)j?*!K?TS87nR6f25my_NDez|rJPxzRdu3!EfLnGcPY6C{K2GFZsQ0kw71(Dd2kZfy*TD13#F$NwBvHnGBFPKIk$ zr***R9r|@>q_$1lMP*&>H=jNQ-@G&uw*_8{5I{M1Y;_z^8s1KdeW|2Xfu%(1X2A9+ zDKR5Y`rb*SFH3vT@|Le2FE<1@M<$NID`^8yKR@xc_3e^*o)g`TRoI9C_edR=Gmnmi zwYG8gNu-GSWoArI$zLNTj@hR?b37T+q}jC1n=dEa+)HMdP7yytlBh*zxwa2eRwe0# z1fqQ*v2C~@yYYjDXQEX$##M-HN0@Q$65N2#3YmplHA5TiKWm|SBkk{FH$>iJn;2ll>^V$HBm zZ@a75h^6pj!)=vCY2ttnR*Ws`Ll1c5;-M z7u*8Ae~je#2`1LJi6t(dUez7c2osOryu}U(+AyPs4SeVrsjqaFD=SQ^qzSA$Gspa% zNnqF-(^QK!1n|E3kk4J%=`-|m^u1+;LZ#?YlOBK61~4F#tzOR?){Bh#9B=wv&9bP_ zIALjsysqXBZ6IRC*r;us<;RPW!mF?t^z4^I_&_Xti9w#;L_K%J0cVvM(u1l~?brLN z+6C@lLtwTGVcoSIuMX<1hLpE0T%a%`2w0)8W=&SVQ6QNqJXqU_M&H-J6#;*f($!E+ zk|$EhZUT3?L%J@jlV5dx7NOoK!kWR-BQ0lQsz{)&yj3iP8okf540(5^eB z&Ga*$&qUmV8UVrHGYf`A{6C_gEHn2}wQ8TwKx)iLam)4h-YueRY+#Ag4flNg%c93y%Q+Ilz5=um{;vGZ^( z=b#O8Ze1OJYP6_2`^z7>i8&>Nvt7mk3D4*12!@C2DS3me;9ij|J%wk3jB-rYos#Q+ zjMB1bD>a>Yva5$8h_9(BqS;O>(2CBNrVGmu{i-?O2_W5<|!14&wPZO9Wbh} zcez(m1Evp_vJ`C4xGY1jtM1s2%DL50tKDW@Iy3)#ZU_RcHP6pH=zv*jaWKA=HV`=A zzSeXu9y>T~=CtA0coZa)m=bvmd^c>|w?gnOp4#+yJ8H11x8ZiFps$m7INi9cum|F6!uRW&ogB@TyYtNdHOATk#It+Ly|W5;^6bC0 F{{wy71bF}e literal 0 HcmV?d00001 diff --git a/source_modelling/focal_mechanism.py b/source_modelling/focal_mechanism.py index 1de29e2..a155080 100644 --- a/source_modelling/focal_mechanism.py +++ b/source_modelling/focal_mechanism.py @@ -14,13 +14,13 @@ The tectonic type (crustal, subduction interface or intraslab) follows the CatSep plausibility scores of Haynie et al. (2026) against the Slab2 -interface geometry and seismogenic zone depths (Hayes et al. 2018). The -interface score is the product of ramp functions in the vertical distance -to the interface, the depth below the seismogenic zone and the Kagan angle -between the mechanism and interface slip. The crustal score ramps down -from 20 km above to 20 km below the interface, and intraslab takes the -remainder. CatSep's CRUST1.0 Moho is not used: the crust is taken to -extend down to the slab. +interface geometry and seismogenic zone depths (Hayes et al. 2018) and the +CRUST1.0 Moho (Laske et al. 2013). The interface score is the product of +ramp functions in the vertical distance to the interface, the depth below +the seismogenic zone and the Kagan angle between the mechanism and +interface slip. The crustal score ramps down from 20 km above to 20 km +below the interface and from the Moho to the slab, and intraslab takes the +remainder. Examples -------- @@ -66,49 +66,11 @@ """Half-width of the interface flex zone about the slab surface, and the buffer below the seismogenic zone (Haynie et al. 2026).""" -class Ramp(NamedTuple): - """Ramp function: `p1` at or below `x1`, `p2` at or above `x2`, linear between (Haynie et al. 2026, eq. 1).""" +KAGAN_RAMP = ((30.0, 53.0), (1.0, 0.25)) +"""Interface score ramp (`np.interp` points) in the Kagan angle in degrees to interface slip (Haynie et al. 2026, Table 1).""" - x1: float - """Start of the ramp.""" - - p1: float - """Score at or below `x1`.""" - - x2: float - """End of the ramp.""" - - p2: float - """Score at or above `x2`.""" - - def __call__(self, x: float) -> float: - """Score at a value. - - Parameters - ---------- - x : float - The value to score. - - Returns - ------- - float - The score at `x`. - """ - if x <= self.x1: - return self.p1 - if x >= self.x2: - return self.p2 - return self.p1 + (self.p2 - self.p1) * (x - self.x1) / (self.x2 - self.x1) - - -KAGAN_RAMP = Ramp(30.0, 1.0, 53.0, 0.25) -"""Interface score in the Kagan angle (degrees) to interface slip (Haynie et al. 2026, Table 1).""" - -CRUSTAL_DEPTH_DIFFERENCE_RAMP = Ramp(-20.0, 1.0, 20.0, 0.0) -"""Crustal score in the depth (km) below the slab surface (Haynie et al. 2026, Table 1).""" - -MANTLE_WEDGE_FLEX_SCORES = (0.0, 0.4) -"""Interface flex scores at the Moho and the top of the flex zone for events in the mantle wedge (Haynie et al. 2026, Table 1).""" +CRUSTAL_DEPTH_DIFFERENCE_RAMP = ((-20.0, 20.0), (1.0, 0.0)) +"""Crustal score ramp (`np.interp` points) in the depth in km below the slab surface (Haynie et al. 2026, Table 1).""" class SeismogenicZone(NamedTuple): @@ -222,13 +184,12 @@ def kagan_angle(plane_a: NodalPlane, plane_b: NodalPlane) -> float: float The Kagan angle in degrees, between 0 and 120. """ - # Diagonal of the rotation between the two principal axis frames; the - # double-couple symmetries flip the signs of two axes at a time. - m = np.diag(_principal_axes(plane_a) @ _principal_axes(plane_b).T) - trace = max( - m[0] + m[1] + m[2], m[0] - m[1] - m[2], -m[0] + m[1] - m[2], -m[0] - m[1] + m[2] + rotation = sp.spatial.transform.Rotation.from_matrix( + _principal_axes(plane_a) @ _principal_axes(plane_b).T ) - return float(np.degrees(np.arccos(np.clip((trace - 1.0) / 2.0, -1.0, 1.0)))) + # The double couple is symmetric under half turns about its principal axes. + symmetries = sp.spatial.transform.Rotation.create_group("D2") + return float(np.degrees((symmetries * rotation).magnitude().min())) class SlabQuery(NamedTuple): @@ -243,13 +204,8 @@ class SlabQuery(NamedTuple): strike: float """Strike of the slab surface in degrees, NaN outside the footprint.""" - region: str | None = None - """Slab2 region code ("ker" or "puy"), None outside the footprint.""" - - @property - def seismogenic_zone(self) -> SeismogenicZone: - """Seismogenic zone of the slab region.""" - return SEISMOGENIC_ZONES.get(self.region or "", DEFAULT_SEISMOGENIC_ZONE) + seismogenic_zone: SeismogenicZone = DEFAULT_SEISMOGENIC_ZONE + """Seismogenic zone of the slab region.""" @property def present(self) -> bool: @@ -270,9 +226,16 @@ class SlabModel: Mapping with keys ``_lat``, ``_lon``, ``_depth``, ``_dip`` and ``_strike`` for each region. + seismogenic_zones : dict[str, SeismogenicZone], optional + Seismogenic zone of each region; others take + `DEFAULT_SEISMOGENIC_ZONE`. """ - def __init__(self, grids: dict[str, npt.NDArray]): + def __init__( + self, + grids: dict[str, npt.NDArray], + seismogenic_zones: dict[str, SeismogenicZone] = SEISMOGENIC_ZONES, + ): """Create a slab model from named grids. Parameters @@ -281,6 +244,9 @@ def __init__(self, grids: dict[str, npt.NDArray]): Mapping with keys ``_lat``, ``_lon``, ``_depth``, ``_dip`` and ``_strike`` for each region. + seismogenic_zones : dict[str, SeismogenicZone], optional + Seismogenic zone of each region; others take + `DEFAULT_SEISMOGENIC_ZONE`. """ self._regions = [] for region in sorted({key.split("_")[0] for key in grids}): @@ -295,7 +261,8 @@ def __init__(self, grids: dict[str, npt.NDArray]): interpolator = sp.interpolate.RegularGridInterpolator( (lats, lons), values, bounds_error=False, fill_value=np.nan ) - self._regions.append((region, bool(lons[-1] > 180.0), interpolator)) + zone = seismogenic_zones.get(region, DEFAULT_SEISMOGENIC_ZONE) + self._regions.append((zone, bool(lons[-1] > 180.0), interpolator)) @classmethod @functools.cache @@ -328,15 +295,83 @@ def query(self, lat: float, lon: float) -> SlabQuery: SlabQuery Slab properties, NaN-valued outside the slab footprints. """ - for region, positive_lon, interpolator in self._regions: + for zone, positive_lon, interpolator in self._regions: query_lon = lon % 360.0 if positive_lon else ((lon + 180.0) % 360.0) - 180.0 depth, dip, sin_strike, cos_strike = interpolator([lat, query_lon])[0] if np.isfinite(depth): strike = np.degrees(np.arctan2(sin_strike, cos_strike)) % 360.0 - return SlabQuery(float(depth), float(dip), float(strike), region) + return SlabQuery(float(depth), float(dip), float(strike), zone) return SlabQuery(float("nan"), float("nan"), float("nan")) +class MohoModel: + """Moho depth, interpolated bilinearly from a grid. + + The packaged grid is the CRUST1.0 depth to Moho (Laske et al. 2013), as + shipped with CatSep, cropped to 1 degree cell centres from 54.5S to + 24.5S and 156.5E to 189.5E. + + Parameters + ---------- + lat : npt.NDArray + Increasing grid latitudes in degrees. + lon : npt.NDArray + Increasing grid longitudes in degrees (0..360). + depth : npt.NDArray + Moho depth in km, shape ``(len(lat), len(lon))``. + """ + + def __init__(self, lat: npt.NDArray, lon: npt.NDArray, depth: npt.NDArray): + """Create a Moho model from a grid. + + Parameters + ---------- + lat : npt.NDArray + Increasing grid latitudes in degrees. + lon : npt.NDArray + Increasing grid longitudes in degrees (0..360). + depth : npt.NDArray + Moho depth in km, shape ``(len(lat), len(lon))``. + """ + # Outside the grid there is no Moho above the slab. + self._interpolator = sp.interpolate.RegularGridInterpolator( + (lat, lon), depth, bounds_error=False, fill_value=np.inf + ) + + @classmethod + @functools.cache + def load(cls) -> MohoModel: + """Load the packaged CRUST1.0 Moho grid. + + Returns + ------- + MohoModel + The loaded Moho model (cached). + """ + with ( + (_DATA_DIR / "crust1_moho_nz.npz").open("rb") as handle, + np.load(handle) as data, + ): + return cls(data["lat"], data["lon"], data["depth"]) + + def depth(self, lat: float, lon: float) -> float: + """Moho depth beneath a location. + + Parameters + ---------- + lat : float + Latitude in degrees. + lon : float + Longitude in degrees (either -180..180 or 0..360). + + Returns + ------- + float + Moho depth in km, infinite outside the grid. + """ + return float(self._interpolator([lat, lon % 360.0])[0]) + + class FaultSegmentIndex: """Straight-line segments of the CFM fault traces with their attributes. @@ -481,22 +516,15 @@ def tectonic_type_probabilities( subduction_kagan_angle: float | None = None, moho_depth: float = np.inf, ) -> dict[TectonicType, float]: - """Tectonic type plausibility scores by CatSep (Haynie et al. 2026). + """Tectonic type plausibility scores by CatSep (Haynie et al. 2026, eqs. 3-5). Events outside the slab footprint are crustal. Inside it, the interface - score is the product (eq. 3) of - - - the flex score: 1 within `FLEX_ZONE_KM` of the slab surface, if no - deeper than `FLEX_ZONE_KM` below the seismogenic zone or above the - Moho; ramping from 0 at the Moho to 0.4 at the top of the flex zone - in the mantle wedge; 0 otherwise, - - the seismogenic score, ramping from 1 at the seismogenic zone depth - to 0 at `FLEX_ZONE_KM` below it, and - - the Kagan score (`KAGAN_RAMP`), when the mechanism is known. - - The crustal score (eq. 4) is the remainder scaled by - `CRUSTAL_DEPTH_DIFFERENCE_RAMP` and by a ramp from 1 at the Moho to 0 - at the slab surface. Intraslab takes what is left (eq. 5). + score is the product of a flex score (1 within `FLEX_ZONE_KM` of the + slab surface; in the mantle wedge, ramping from 0 at the Moho to 0.4 at + the top of the flex zone), a ramp from 1 at the seismogenic zone depth + to 0 `FLEX_ZONE_KM` below it, and `KAGAN_RAMP`. The crustal score is the + remainder scaled by `CRUSTAL_DEPTH_DIFFERENCE_RAMP` and a ramp from 1 at + the Moho to 0 at the slab surface. Intraslab takes what is left. Parameters ---------- @@ -511,8 +539,8 @@ def tectonic_type_probabilities( Kagan angle in degrees between the event's mechanism and interface slip, None without a focal mechanism. moho_depth : float, optional - Moho depth in km. The default puts the Moho below the slab, so the - crust extends down to the slab surface and there is no mantle wedge. + Moho depth in km. The default (no Moho) extends the crust down to + the slab surface, leaving no mantle wedge. Returns ------- @@ -526,31 +554,33 @@ def tectonic_type_probabilities( TectonicType.SLAB: 0.0, } flex_top = slab_depth - FLEX_ZONE_KM - in_flex_zone = flex_top <= depth <= slab_depth + FLEX_ZONE_KM - if in_flex_zone and ( - depth <= seismogenic_depth + FLEX_ZONE_KM or depth <= moho_depth - ): + # Below the seismogenic zone the seismogenic ramp is already zero, so + # CatSep's extra depth conditions on the flex zone change nothing. + if flex_top <= depth <= slab_depth + FLEX_ZONE_KM: flex = 1.0 elif moho_depth <= depth < flex_top: - flex = Ramp( - moho_depth, - MANTLE_WEDGE_FLEX_SCORES[0], - flex_top, - MANTLE_WEDGE_FLEX_SCORES[1], - )(depth) + flex = np.interp(depth, [moho_depth, flex_top], [0.0, 0.4]) else: flex = 0.0 - seismogenic = Ramp(seismogenic_depth, 1.0, seismogenic_depth + FLEX_ZONE_KM, 0.0)( - depth + seismogenic = np.interp( + depth, [seismogenic_depth, seismogenic_depth + FLEX_ZONE_KM], [1.0, 0.0] ) kagan = ( - 1.0 if subduction_kagan_angle is None else KAGAN_RAMP(subduction_kagan_angle) + 1.0 + if subduction_kagan_angle is None + else np.interp(subduction_kagan_angle, *KAGAN_RAMP) ) - interface = flex * seismogenic * kagan - crustal = ( + interface = float(flex * seismogenic * kagan) + # Below a Moho that is itself below the slab, the score is 0. + above_moho = ( + 1.0 + if depth <= moho_depth + else np.interp(depth, [moho_depth, max(moho_depth, slab_depth)], [1.0, 0.0]) + ) + crustal = float( (1.0 - interface) - * CRUSTAL_DEPTH_DIFFERENCE_RAMP(depth - slab_depth) - * Ramp(moho_depth, 1.0, slab_depth, 0.0)(depth) + * np.interp(depth - slab_depth, *CRUSTAL_DEPTH_DIFFERENCE_RAMP) + * above_moho ) return { TectonicType.CRUSTAL: crustal, @@ -598,6 +628,8 @@ class CMTClassifier: Log-odds multiplier at the surface (see `depth_temperature`). temperature_depth_rate : float, optional Decay rate of the log-odds multiplier per km of centroid depth. + moho_model : MohoModel | None, optional + Moho depth for the tectonic type, None for no Moho. """ def __init__( @@ -607,6 +639,7 @@ def __init__( weights: npt.NDArray[np.float64], temperature_scale: float = 1.0, temperature_depth_rate: float = 0.0, + moho_model: MohoModel | None = None, ): """Create a classifier. @@ -623,12 +656,15 @@ def __init__( Log-odds multiplier at the surface (see `depth_temperature`). temperature_depth_rate : float, optional Decay rate of the log-odds multiplier per km of centroid depth. + moho_model : MohoModel | None, optional + Moho depth for the tectonic type, None for no Moho. """ self.segments = segments self.slab_model = slab_model self.weights = weights self.temperature_scale = temperature_scale self.temperature_depth_rate = temperature_depth_rate + self.moho_model = moho_model @classmethod def from_faults(cls, faults: list[CommunityFault]) -> CMTClassifier: @@ -654,12 +690,13 @@ def from_faults(cls, faults: list[CommunityFault]) -> CMTClassifier: np.asarray(model["weights"], dtype=float), model["temperature"]["scale"], model["temperature"]["depth_rate"], + MohoModel.load(), ) @classmethod @functools.cache def load(cls) -> CMTClassifier: - """Load the classifier with the packaged CFM, slab grids and model. + """Load the classifier with the packaged CFM, slab and Moho grids and model. Returns ------- @@ -787,18 +824,26 @@ def tectonic_type_probabilities( """ lat, lon, depth = _parse_centroid(centroid) slab = self.slab_model.query(lat, lon) - zone = slab.seismogenic_zone # Both planes describe the same mechanism; using both keeps the # result symmetric when catalogue rounding makes them disagree. planes = [p for p in (nodal_plane_1, nodal_plane_2) if p is not None] subduction_kagan_angle = None if planes and slab.present: - interface_slip = NodalPlane(slab.strike, slab.dip, zone.rake) + interface_slip = NodalPlane( + slab.strike, slab.dip, slab.seismogenic_zone.rake + ) subduction_kagan_angle = min( kagan_angle(plane, interface_slip) for plane in planes ) + moho_depth = ( + np.inf if self.moho_model is None else self.moho_model.depth(lat, lon) + ) return tectonic_type_probabilities( - depth, slab.depth, zone.depth, subduction_kagan_angle + depth, + slab.depth, + slab.seismogenic_zone.depth, + subduction_kagan_angle, + moho_depth, ) def tectonic_type( diff --git a/tests/test_focal_mechanism.py b/tests/test_focal_mechanism.py index 1e85588..abc342e 100644 --- a/tests/test_focal_mechanism.py +++ b/tests/test_focal_mechanism.py @@ -103,17 +103,6 @@ def test_kagan_angle_is_at_most_120(): assert 0.0 <= min(angles) and max(angles) <= 120.0 -def test_ramp(): - decreasing = focal_mechanism.Ramp(30.0, 1.0, 53.0, 0.25) - assert decreasing(0.0) == 1.0 - assert decreasing(30.0) == 1.0 - assert decreasing(41.5) == pytest.approx(0.625) - assert decreasing(53.0) == 0.25 - assert decreasing(120.0) == 0.25 - increasing = focal_mechanism.Ramp(5.0, 0.0, 15.0, 0.4) - assert increasing(10.0) == pytest.approx(0.2) - - @pytest.mark.parametrize( "depth, slab_depth, kagan, moho_depth, expected", [ @@ -130,6 +119,9 @@ def test_ramp(): (60.0, 30.0, None, np.inf, (0.0, 0.0, 1.0)), # In the mantle wedge, between the Moho and the top of the flex zone. (10.0, 30.0, None, 5.0, (0.64, 0.2, 0.16)), + # Below the seismogenic zone and below a Moho that is below the slab. + (48.0, 40.0, None, np.inf, (0.3, 0.0, 0.7)), + (48.0, 40.0, None, 45.0, (0.0, 0.0, 1.0)), ], ) def test_tectonic_type_probabilities( @@ -149,14 +141,18 @@ def test_tectonic_type_probabilities( ] == pytest.approx(expected) -def test_seismogenic_zone(): - slab_model = focal_mechanism.SlabModel.load() - assert slab_model.query(-41.3, 174.8).region == "ker" # Wellington - assert slab_model.query(*CENTROID[:2]).region == "puy" - assert slab_model.query(-41.3, 174.8).seismogenic_zone.depth == 47.0 - assert slab_model.query(*CENTROID[:2]).seismogenic_zone.depth == 30.0 - outside = slab_model.query(-43.5, 172.6) # Christchurch - assert outside.seismogenic_zone == focal_mechanism.DEFAULT_SEISMOGENIC_ZONE +def test_seismogenic_zone(classifier: focal_mechanism.CMTClassifier): + slab_model = classifier.slab_model + assert slab_model.query(-41.3, 174.8).seismogenic_zone.depth == 47.0 # Hikurangi + assert slab_model.query(*CENTROID[:2]).seismogenic_zone.depth == 30.0 # Puysegur + + +def test_moho_model(): + moho_model = focal_mechanism.MohoModel.load() + wellington = moho_model.depth(-41.3, 174.8) + assert 15 < wellington < 45 + assert moho_model.depth(-37.0, 179.5 - 360) == moho_model.depth(-37.0, 179.5) + assert moho_model.depth(0.0, 0.0) == np.inf def test_tectonic_type(classifier: focal_mechanism.CMTClassifier):