Add empirical calculations to intensity measures

pyproject.toml

@@ -14,11 +14,12 @@ dependencies = [
   "im-calculation>=2025.12.5",
   "velocity-modelling>=2026.2.1",
   "nshmdb>=2025.12.1",
-  "oq_wrapper>=2025.12.3",
+  "oq_wrapper>=2026.05.2",
   "qcore-utils>=2025.12.2",
   "source_modelling>=2026.07.2",
   # Data Formats
   "geopandas",
+  "netCDF4",               # Required for openquake shenanigans
   "pandas[parquet, hdf5]",
   "pyyaml",
   "xarray[io]",
@@ -30,10 +31,10 @@ dependencies = [
   "tqdm",
   "typer",
   # Misc.
-  "requests",  # For gcmt-to-realisation
-  "schema",    # For loading realisations
-  "structlog", # Logging.
-  "psutil",    # To get the CPU affinity for jobs
+  "requests",      # For gcmt-to-realisation
+  "schema",        # For loading realisations
+  "structlog",     # Logging.
+  "psutil",        # To get the CPU affinity for jobs
   "parse>=1.21.0",
   "rich>=14.3.2",
 ]

workflow/default_parameters/root/defaults.yaml

@@ -502,6 +502,9 @@ velocity_model_1d:
       rho: 3.33
       Qp: 394.80
       Qs: 197.40
+empirical:
+  tect_type: "active_shallow"
+  models: ["NSHM2022"]
 im:
   ims: ["PGA", "PGV", "PGD", "CAV", "AI", "Ds575", "Ds595", "pSA", "FAS"]
   valid_periods:

workflow/realisations.py

@@ -646,6 +646,23 @@ class Rakes(RealisationConfiguration):
     rakes: dict[str, float]
     """A map from faults to their rake angles."""
 
+    def as_vectors(self) -> dict[str, npt.NDArray[np.float64]]:
+        """Represent each rake angle as a unit vector.
+
+        Rakes are angles, so they cannot be averaged directly (the mean of
+        -179° and 179° is 0°, not 180°). Averaging the unit vectors and
+        recovering the angle with `arctan2` avoids this.
+
+        Returns
+        -------
+        dict
+            A map from faults to the unit vector of their rake angle.
+        """
+        return {
+            k: np.array([np.cos(np.radians(r)), np.sin(np.radians(r))])
+            for k, r in self.rakes.items()
+        }
+
     def __getitem__(self, key: str) -> float:
         """Get the rake for a fault name.
 
@@ -687,6 +704,38 @@ def __getitem__(self, key: str) -> BoldM:
         """
         return self.magnitudes[key]
 
+    @property
+    def moments(self) -> dict[str, float]:  # numpydoc ignore=RT01
+        """dict: a map from faults to their moment."""
+        return {
+            k: moment.magnitude_to_moment(mag, bold_m=True)
+            for k, mag in self.magnitudes.items()
+        }
+
+    def moment_averaged(self, values: dict[str, Any]) -> Any:
+        """Average per-fault quantities, weighted by fault moment.
+
+        Parameters
+        ----------
+        values : dict
+            A map from faults to the quantity to average. Every fault in
+            this realisation must be present. Values may be scalars or
+            arrays, provided they all share the same shape.
+
+        Returns
+        -------
+        Any
+            The moment-weighted average of `values`, with the same shape as
+            the individual values.
+        """
+        keys = list(self.magnitudes)
+        moments = self.moments
+        return np.average(
+            [values[key] for key in keys],
+            weights=[moments[key] for key in keys],
+            axis=0,
+        )
+
     @property
     def total_moment(self) -> float:  # numpydoc ignore=RT01
         """float: total moment of realisation"""
@@ -1223,6 +1272,21 @@ def to_dict(self) -> dict[str, Any]:
         return _dict
 
 
+@dataclasses.dataclass
+class EmpiricalParameters(RealisationConfiguration):
+    """Empirical (ground motion model) intensity measure parameters."""
+
+    _config_key: ClassVar[str] = "empirical"
+    _schema: ClassVar[Schema] = schemas.EMPIRICAL_PARAMETERS
+
+    # Types here are not explicitly declared so we do not pay the openquake tax
+    # importing this module.
+    tect_type: Any
+    """The tectonic type of the source (an `oq_wrapper.constants.TectType`)."""
+    models: list[Any]
+    """The ground motion models or logic trees to evaluate."""
+
+
 @dataclasses.dataclass
 class LogEntry:
     """Log entry for workflow utilities."""

workflow/schemas.py

@@ -1170,6 +1170,28 @@ def _corners_to_array(corners_spec: list[dict[str, float]]) -> np.ndarray:
 )
 
 
+EMPIRICAL_PARAMETERS = Schema(
+    {
+        Literal(
+            "tect_type",
+            description="Tectonic type of the source (one of oq_wrapper.constants.TectType)",
+        ): str,
+        Literal(
+            "models",
+            description=(
+                "Ground motion models or ground motion model logic trees to "
+                "evaluate (members of oq_wrapper.constants.GMM or "
+                "oq_wrapper.constants.GMMLogicTree)"
+            ),
+        ): [str],
+    }
+)
+# NOTE: The values of this schema are validated as plain strings rather than
+# `oq_wrapper.constants` enum members. Importing `oq_wrapper.constants` pulls in
+# OpenQuake, which is expensive (and must be precompiled), so the strings are
+# only resolved to enum members inside the IM calculation stage.
+
+
 LOG_ENTRY_SCHEMA = Schema(
     {
         Literal(

workflow/scripts/bb_sim.py

@@ -217,6 +217,7 @@ def combine_hf_and_lf(
             "y": ("station", lf.y.values),
             "latitude": ("station", lf.lat.values),
             "longitude": ("station", lf.lon.values),
+            "vs30": vs30_df["vsite"].to_xarray(),
         },
         attrs={
             "units": "g",