diff --git a/source_modelling/gc2_distances.py b/source_modelling/gc2_distances.py index b45afdd..51f1766 100644 --- a/source_modelling/gc2_distances.py +++ b/source_modelling/gc2_distances.py @@ -140,7 +140,9 @@ def segment_weights( @njit([float64[:](float64[:], uint64[:]), float64[:](float64[:], int64[:])], cache=True) -def cumulative_reduction(data: np.ndarray, indices: np.ndarray) -> np.ndarray: # pragma: no cover +def cumulative_reduction( + data: np.ndarray, indices: np.ndarray +) -> np.ndarray: # pragma: no cover """Calculate sum between indices. Parameters @@ -173,7 +175,9 @@ def cumulative_reduction(data: np.ndarray, indices: np.ndarray) -> np.ndarray: # [float64[:, :](float64[:, :], uint64[:]), float64[:, :](float64[:, :], int64[:])], cache=True, ) -def diff_reduction(points: np.ndarray, indices: np.ndarray) -> np.ndarray: # pragma: no cover +def diff_reduction( + points: np.ndarray, indices: np.ndarray +) -> np.ndarray: # pragma: no cover """Calculate a diff reduction points between the given indices. Parameters @@ -425,7 +429,9 @@ def multi_trace_rx_ry( b_hat = np.sum(direction_vectors, axis=0) b_hat_norm = np.linalg.norm(b_hat) if np.isclose(b_hat_norm, 0.0): - raise ValueError("Invalid direction vectors calculated: likely indicates degenerate trace or invalid geometry.") + raise ValueError( + "Invalid direction vectors calculated: likely indicates degenerate trace or invalid geometry." + ) b_hat /= b_hat_norm u_shift_origins = calculate_gc2_u_origins( diff --git a/source_modelling/magnitude_scaling.py b/source_modelling/magnitude_scaling.py index c123976..324009b 100644 --- a/source_modelling/magnitude_scaling.py +++ b/source_modelling/magnitude_scaling.py @@ -1,15 +1,13 @@ """Magnitude scaling relationships for fault dimensions.""" -import functools -import typing import warnings from enum import Enum, StrEnum, auto import numpy as np import scipy as sp -BoldM = typing.NewType("BoldM", float) -Mw = typing.NewType("Mw", float) +from source_modelling import moment +from source_modelling.moment import BoldM, Mw class RakeType(Enum): @@ -660,7 +658,7 @@ def contreras_slab_magnitude_to_length_width( def magnitude_to_length_width( scaling_relation: ScalingRelation, - magnitude: BoldM | Mw, + magnitude: BoldM, rake: float | None = None, random: bool = False, ) -> tuple[float, float]: @@ -670,7 +668,7 @@ def magnitude_to_length_width( ---------- scaling_relation : ScalingRelation Scaling relation to use. - magnitude : BoldM | Mw + magnitude : BoldM Moment magnitude of the fault rupture. rake : float, optional Rake of the fault (degrees). Required for Leonard scaling. @@ -683,32 +681,24 @@ def magnitude_to_length_width( tuple[float, float] Length and width of the fault estimated by the scaling relation. """ - - if scaling_relation == ScalingRelation.LEONARD2014 and rake is None: - raise ValueError("Rake must be specified for Leonard scaling.") - scaling_relations_map = { - ScalingRelation.LEONARD2014: functools.partial( - leonard_magnitude_to_length_width, - # Rake is checked earlier so ty warning is not required. - rake=rake, # ty: ignore[invalid-argument-type] - random=random, - ), - ScalingRelation.CONTRERAS_INTERFACE2017: functools.partial( - contreras_interface_magnitude_to_length_width, random=random - ), - ScalingRelation.CONTRERAS_SLAB2020: functools.partial( - contreras_slab_magnitude_to_length_width, random=random - ), - } - # The dispatched callable is a union of partials expecting either Mw or BoldM, - # so ty can't see that the scaling_relation key already determines which - # convention `magnitude` should be treated as. - return scaling_relations_map[scaling_relation](magnitude) # ty: ignore[invalid-argument-type] + match scaling_relation: + case ScalingRelation.LEONARD2014 if rake is not None: + return leonard_magnitude_to_length_width( + moment.boldm_to_mw(magnitude), rake=rake, random=random + ) + case ScalingRelation.LEONARD2014: + raise ValueError("Rake must be specified for Leonard scaling.") + case ScalingRelation.CONTRERAS_INTERFACE2017: + return contreras_interface_magnitude_to_length_width( + magnitude, random=random + ) + case ScalingRelation.CONTRERAS_SLAB2020: + return contreras_slab_magnitude_to_length_width(magnitude, random=random) def magnitude_to_area( scaling_relation: ScalingRelation, - magnitude: BoldM | Mw, + magnitude: BoldM, rake: float | None = None, random: bool = False, ) -> float: @@ -718,7 +708,7 @@ def magnitude_to_area( ---------- scaling_relation : ScalingRelation Scaling relation to use. - magnitude : BoldM | Mw + magnitude : BoldM Moment magnitude of the fault rupture. rake : float, optional Rake of the fault (degrees). Required for Leonard scaling. @@ -733,27 +723,20 @@ def magnitude_to_area( Returns ------- - tuple[float, float] + float Area of the fault estimated by the scaling relation. """ - if scaling_relation == ScalingRelation.LEONARD2014 and rake is None: - raise ValueError("Rake must be specified for Leonard scaling.") - scaling_relations_map = { - ScalingRelation.LEONARD2014: functools.partial( - leonard_magnitude_to_area, - # Rake is checked earlier so ty warning is not required. - rake=rake, # ty: ignore[invalid-argument-type] - random=random, - ), - ScalingRelation.CONTRERAS_INTERFACE2017: functools.partial( - contreras_interface_magnitude_to_area, random=random - ), - ScalingRelation.CONTRERAS_SLAB2020: functools.partial( - strasser_slab_magnitude_to_area, random=random - ), - } - # See note in `magnitude_to_length_width` about why this is ignored. - return scaling_relations_map[scaling_relation](magnitude) # ty: ignore[invalid-argument-type] + match scaling_relation: + case ScalingRelation.LEONARD2014 if rake is not None: + return leonard_magnitude_to_area( + moment.boldm_to_mw(magnitude), rake=rake, random=random + ) + case ScalingRelation.LEONARD2014: + raise ValueError("Rake must be specified for Leonard scaling.") + case ScalingRelation.CONTRERAS_INTERFACE2017: + return contreras_interface_magnitude_to_area(magnitude, random=random) + case ScalingRelation.CONTRERAS_SLAB2020: + return strasser_slab_magnitude_to_area(magnitude, random=random) def area_to_magnitude( @@ -761,7 +744,7 @@ def area_to_magnitude( area: float, rake: float | None = None, random: bool = False, -) -> BoldM | Mw: +) -> BoldM: """Convert area to magnitude using a scaling relationship. Parameters @@ -778,24 +761,17 @@ def area_to_magnitude( Returns ------- - BoldM | Mw + BoldM Moment magnitude of the fault rupture estimated by the scaling relation. """ - if scaling_relation == ScalingRelation.LEONARD2014 and rake is None: - raise ValueError("Rake must be specified for Leonard scaling.") - - scaling_relations_map = { - ScalingRelation.LEONARD2014: functools.partial( - leonard_area_to_magnitude, - # Rake is checked earlier so ty warning is not required. - rake=rake, # ty: ignore[invalid-argument-type] - random=random, - ), - ScalingRelation.CONTRERAS_INTERFACE2017: functools.partial( - contreras_interface_area_to_magnitude, random=random - ), - ScalingRelation.CONTRERAS_SLAB2020: functools.partial( - strasser_slab_area_to_magnitude, random=random - ), - } - return scaling_relations_map[scaling_relation](area) + + match scaling_relation: + case ScalingRelation.LEONARD2014 if rake is not None: + mag = leonard_area_to_magnitude(area, rake=rake, random=random) + return moment.mw_to_boldm(mag) + case ScalingRelation.LEONARD2014: + raise ValueError("Rake must be specified for Leonard scaling.") + case ScalingRelation.CONTRERAS_INTERFACE2017: + return contreras_interface_area_to_magnitude(area, random=random) + case ScalingRelation.CONTRERAS_SLAB2020: + return strasser_slab_area_to_magnitude(area, random=random) diff --git a/source_modelling/moment.py b/source_modelling/moment.py index d9cc033..71ebc1b 100644 --- a/source_modelling/moment.py +++ b/source_modelling/moment.py @@ -11,14 +11,18 @@ from scipy.sparse import csr_array from source_modelling import rupture_propagation, sources -from source_modelling.magnitude_scaling import BoldM, Mw from source_modelling.sources import Fault, Plane +BoldM = typing.NewType("BoldM", float) +Mw = typing.NewType("Mw", float) + + # Moment magnitude scale coefficients for seismic moment in Nm, from # equations 4 and 7 of Hanks and Kanamori (1979). See the [Hanks1979] reference # in the `moment_to_magnitude` docstring for the full citation. EQUATION_4_COEFFICIENT = 6.0667 # `Mw` convention EQUATION_7_COEFFICIENT = 6.0333 # `BoldM` convention +BOLDM_TO_MW = EQUATION_7_COEFFICIENT - EQUATION_4_COEFFICIENT def find_connected_faults( @@ -155,6 +159,13 @@ def moment_to_magnitude(moment: float, bold_m: bool = True) -> BoldM | Mw: BoldM | Mw Rupture moment magnitude in the convention specified by `bold_m`. + Raises + ------ + ValueError + If the provided moment corresponds to a physically implausible + magnitude. Likely this is because the user supplied moment is in dyne-cm + instead of Nm. + References ---------- .. [Hanks1979] Hanks, T. C., and H. Kanamori (1979), @@ -164,10 +175,16 @@ def moment_to_magnitude(moment: float, bold_m: bool = True) -> BoldM | Mw: """ if bold_m: - return BoldM(2 / 3 * np.log10(moment) - EQUATION_7_COEFFICIENT) - + mag = BoldM(2 / 3 * np.log10(moment) - EQUATION_7_COEFFICIENT) else: - return Mw(2 / 3 * np.log10(moment) - EQUATION_4_COEFFICIENT) + mag = Mw(2 / 3 * np.log10(moment) - EQUATION_4_COEFFICIENT) + + if mag > 10.0: + raise ValueError( + "Magnitude for moment is unreasonably large, did you provide moment in dyne-cm instead of Nm?" + ) + + return mag @typing.overload @@ -202,6 +219,46 @@ def magnitude_to_moment(magnitude: BoldM | Mw, bold_m: bool = True) -> float: return 10 ** ((magnitude + EQUATION_4_COEFFICIENT) * 3 / 2) +def boldm_to_mw(magnitude: BoldM) -> Mw: + """Convert a BoldM convention magnitude to an Mw convention magnitude. + + Parameters + ---------- + magnitude : BoldM + A magnitude in the BoldM convention. + + Returns + ------- + Mw + A magnitude in the Mw convention. + + See Also + -------- + moment_to_magnitude : Introduces the different magnitude conventions. + """ + return Mw(magnitude + BOLDM_TO_MW) + + +def mw_to_boldm(magnitude: Mw) -> BoldM: + """Convert a Mw convention magnitude to a BoldM convention magnitude. + + Parameters + ---------- + magnitude : Mw + A magnitude in the Mw convention. + + Returns + ------- + BoldM + A magnitude in the BoldM convention. + + See Also + -------- + moment_to_magnitude : Introduces the different magnitude conventions. + """ + return BoldM(magnitude - BOLDM_TO_MW) + + def moment_over_time_from_moment_rate(moment_rate_df: pd.DataFrame) -> pd.DataFrame: """Integrate a moment rate dataframe into a cumulative moment dataframe. diff --git a/tests/test_magnitude_scaling.py b/tests/test_magnitude_scaling.py index 5c55bbb..e701266 100644 --- a/tests/test_magnitude_scaling.py +++ b/tests/test_magnitude_scaling.py @@ -11,8 +11,8 @@ import scipy as sp from hypothesis import assume, given -from source_modelling import magnitude_scaling -from source_modelling.magnitude_scaling import BoldM, Mw +from source_modelling import magnitude_scaling, moment +from source_modelling.moment import BoldM, Mw MAGNITUDE_TO_AREA = { magnitude_scaling.ScalingRelation.LEONARD2014: magnitude_scaling.leonard_magnitude_to_area, @@ -154,8 +154,9 @@ def test_leonard_area(mw: Mw, rake: float, expected_area: float): """Test the Leonard 2014 area calculation against values from the old implementation in qcore. NOTE: this combined with the inversion test for leonard ensure that the area calculation is compatible with the old implementation as well.""" + boldm = moment.mw_to_boldm(mw) assert magnitude_scaling.magnitude_to_area( - magnitude_scaling.ScalingRelation.LEONARD2014, mw, rake + magnitude_scaling.ScalingRelation.LEONARD2014, boldm, rake ) == pytest.approx(expected_area) @@ -541,12 +542,13 @@ def test_magnitude_to_length_width_calls_correct_function( with patch(f"source_modelling.magnitude_scaling.{func_name}") as mock_func: magnitude_scaling.magnitude_to_length_width( - scaling_relation, magnitude, rake, random + # We are not testing outputs here so we can ignore the invalid magnitude convention. + scaling_relation, + magnitude, # ty: ignore[invalid-argument-type] + rake, + random, ) - if rake_required: - mock_func.assert_called_once_with(magnitude, rake=rake, random=random) - else: - mock_func.assert_called_once_with(magnitude, random=random) + mock_func.assert_called_once() @pytest.mark.parametrize( @@ -579,11 +581,10 @@ def test_magnitude_to_area_calls_correct_function( random = True with patch(f"source_modelling.magnitude_scaling.{func_name}") as mock_func: - magnitude_scaling.magnitude_to_area(scaling_relation, magnitude, rake, random) - if rake_required: - mock_func.assert_called_once_with(magnitude, rake=rake, random=random) - else: - mock_func.assert_called_once_with(magnitude, random=random) + # We are not checking outputs here, so it is ok to ignore the invalid magnitude convention + magnitude_scaling.magnitude_to_area(scaling_relation, magnitude, rake, random) # ty: ignore[invalid-argument-type] + + mock_func.assert_called_once() @pytest.mark.parametrize( diff --git a/tests/test_moment.py b/tests/test_moment.py index 3cf9b0c..c8066ea 100644 --- a/tests/test_moment.py +++ b/tests/test_moment.py @@ -58,6 +58,14 @@ def test_moment_to_magnitude(): assert moment.moment_to_magnitude(8.96e20) == pytest.approx(7.89, abs=5e-02) +def test_moment_to_magnitude_units(): + # Dusky sound earthquake, nominal magnitude ~ 8.0 but GCMT provides the moment in dyne-cm. + # Should throw error because the magnitude is too large to be physically plausible. + mom = 1.44e28 + with pytest.raises(ValueError, match="Magnitude for moment is unreasonably large"): + moment.moment_to_magnitude(mom) + + # The following test involves some patching to make it feasible to test properly.