Perform basic sanity check on moment magnitude calculations

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(

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,35 +723,28 @@ 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(
     scaling_relation: ScalingRelation,
     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)

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.
 

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(

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.