From 5bd305cbec33c84f197cce1ba15dfeae35e065be Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 16 Feb 2026 16:09:22 +1300 Subject: [PATCH 01/21] refactor(defaults): add PGD and im calc periods --- .../default_parameters/root/defaults.yaml | 84 ++++++++++++++++++- 1 file changed, 82 insertions(+), 2 deletions(-) diff --git a/workflow/default_parameters/root/defaults.yaml b/workflow/default_parameters/root/defaults.yaml index 223d049b..1dfcae6a 100644 --- a/workflow/default_parameters/root/defaults.yaml +++ b/workflow/default_parameters/root/defaults.yaml @@ -342,40 +342,120 @@ velocity_model_1d: Qp: 460.00 Qs: 230.00 im: - ims: ["PGA", "PGV", "CAV", "AI", "Ds575", "Ds595", "pSA", "FAS"] + ims: ["PGA", "PGV", "PGD", "CAV", "AI", "Ds575", "Ds595", "pSA", "FAS"] valid_periods: [ 0.01, 0.02, + 0.022, + 0.025, + 0.029, 0.03, + 0.032, + 0.035, + 0.036, 0.04, + 0.042, + 0.044, + 0.045, + 0.046, + 0.048, 0.05, + 0.055, + 0.06, + 0.065, + 0.067, + 0.07, 0.075, + 0.08, + 0.085, + 0.09, + 0.095, 0.1, + 0.11, 0.12, + 0.13, + 0.133, + 0.14, 0.15, + 0.16, 0.17, + 0.18, + 0.19, 0.2, + 0.22, + 0.24, 0.25, + 0.26, + 0.28, + 0.29, 0.3, + 0.32, + 0.34, + 0.35, + 0.36, + 0.38, 0.4, + 0.42, + 0.44, + 0.45, + 0.46, + 0.48, 0.5, + 0.55, 0.6, + 0.65, + 0.667, 0.7, 0.75, 0.8, + 0.85, 0.9, + 0.95, 1.0, - 1.25, + 1.1, + 1.2, + 1.3, + 1.4, 1.5, + 1.6, + 1.7, + 1.8, + 1.9, 2.0, + 2.2, + 2.4, 2.5, + 2.6, + 2.8, 3.0, + 3.2, + 3.4, + 3.5, + 3.6, + 3.8, 4.0, + 4.2, + 4.4, + 4.6, + 4.8, 5.0, + 5.5, 6.0, + 6.5, + 7.0, 7.5, + 8.0, + 8.5, + 9.0, + 9.5, 10.0, + 11.0, + 12.0, + 13.0, + 14.0, + 15.0, + 20.0, ] fas_frequencies: [ From f15a07fd57fd836c73e0adbabd761057805d34e9 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 20 Jul 2026 13:51:10 +1200 Subject: [PATCH 02/21] include rx/ry, source geometry, domain geometry, magnitudes --- workflow/realisations.py | 15 +++++- workflow/scripts/im_calc.py | 96 +++++++++++++++++++++++++++++++++++-- 2 files changed, 105 insertions(+), 6 deletions(-) diff --git a/workflow/realisations.py b/workflow/realisations.py index 62ee0af5..ba5f3fde 100644 --- a/workflow/realisations.py +++ b/workflow/realisations.py @@ -25,7 +25,7 @@ from schema import Schema from IM import im_calculation -from source_modelling import sources +from source_modelling import moment, sources from source_modelling.magnitude_scaling import BoldM from source_modelling.rupture_propagation import JumpPair from source_modelling.sources import IsSource @@ -687,6 +687,19 @@ def __getitem__(self, key: str) -> BoldM: """ return self.magnitudes[key] + @property + def total_moment(self) -> float: + """float: total moment of realisation""" + return sum( + moment.magnitude_to_moment(mag, bold_m=True) + for mag in self.magnitudes.values() + ) + + @property + def total_magnitude(self) -> BoldM: + """float: total magnitude of realisation""" + return moment.moment_to_magnitude(self.total_moment, bold_m=True) + @dataclasses.dataclass class RupturePropagationConfig(RealisationConfiguration): diff --git a/workflow/scripts/im_calc.py b/workflow/scripts/im_calc.py index 65645d89..bd341162 100644 --- a/workflow/scripts/im_calc.py +++ b/workflow/scripts/im_calc.py @@ -33,6 +33,7 @@ import numpy as np import pandas as pd +import shapely import tqdm import typer import xarray as xr @@ -40,9 +41,13 @@ from IM import im_reader, ims from IM.im_calculation import IM from qcore import cli, coordinates +from source_modelling import sources +from source_modelling.sources import IsSource from workflow import realisations, utils from workflow.realisations import ( + DomainParameters, IntensityMeasureCalculationParameters, + Magnitudes, RealisationMetadata, Resolution, RupturePropagationConfig, @@ -54,6 +59,57 @@ app = typer.Typer() +def _source_polygon(source_geometries: dict[str, IsSource]) -> shapely.Geometry: + """Extract source polygon in longitude, latitude format. + + Parameters + ---------- + source_geometries : dict[str, IsSource] + Realisation faults. + + Returns + ------- + Geometry + The union of all fault geometries. + """ + geometries = [] + for fault in source_geometries.values(): + geometry = fault.geometry + geometry = shapely.transform( + geometry, lambda c: coordinates.nztm_to_wgs_depth(c)[:, ::-1] + ) + + geometries.append(geometry) + return shapely.normalize(shapely.union_all(geometries)) + + +def _trace_polygon(source_geometries: dict[str, IsSource]) -> shapely.Geometry: + """Extract trace polygon in longitude, latitude format. + + Parameters + ---------- + source_geometries : dict[str, IsSource] + Realisation faults. + + Returns + ------- + Geometry + The union of all traces of geometries (for geometries that have traces). + """ + geometries = [] + for fault in source_geometries.values(): + if not hasattr(fault, "trace_geometry"): + continue + geometry = fault.trace_geometry + geometry = shapely.transform( + geometry, lambda c: coordinates.nztm_to_wgs_depth(c)[:, ::-1] + ) + + geometries.append(geometry) + + return shapely.normalize(shapely.union_all(geometries)) + + @cli.from_docstring(app) def calculate_instensity_measures( realisation_ffp: Annotated[ @@ -106,7 +162,9 @@ def calculate_instensity_measures( ) ) source_geometries = SourceConfig.read_from_realisation(realisation_ffp) + domain_parameters = DomainParameters.read_from_realisation(realisation_ffp) rup_prop_config = RupturePropagationConfig.read_from_realisation(realisation_ffp) + magnitudes = Magnitudes.read_from_realisation(realisation_ffp) broadband = xr.open_dataset(broadband_simulation_ffp) @@ -159,6 +217,7 @@ def calculate_instensity_measures( } latitude = broadband.latitude.values longitude = broadband.longitude.values + station_locations = np.stack((latitude, longitude), axis=-1) rrup = ( np.array( @@ -167,7 +226,7 @@ def calculate_instensity_measures( source.rrup_distance(np.append(station, 0)) for source in source_geometries.source_geometries.values() ) - for station in np.stack((latitude, longitude), axis=-1) + for station in station_locations ] ) / 1000 @@ -179,7 +238,7 @@ def calculate_instensity_measures( source.rjb_distance(np.append(station, 0)) for source in source_geometries.source_geometries.values() ) - for station in np.stack((latitude, longitude), axis=-1) + for station in station_locations ] ) / 1000 @@ -190,18 +249,27 @@ def calculate_instensity_measures( hyp = ( coordinates.distance_between_wgs_depth_coordinates( - np.stack((latitude, longitude, np.zeros_like(latitude)), axis=-1), + np.stack((station_locations, np.zeros_like(latitude)), axis=1), hypocentre, ) / 1000 ) epi = ( coordinates.distance_between_wgs_depth_coordinates( - np.stack((latitude, longitude), axis=-1), + station_locations, hypocentre[:2], ) / 1000 ) + all_faults_have_rx_ry = all( + isinstance(source, sources.Plane | sources.Fault) + for source in source_geometries.source_geometries.values() + ) + if all_faults_have_rx_ry: + rx, ry = sources.multi_fault_rx_ry_distance( + list(source_geometries.source_geometries.values()), # ty: ignore[invalid-argument-type] + station_locations, + ) dataset = xr.Dataset( coords={ @@ -210,12 +278,30 @@ def calculate_instensity_measures( "component", ["000", "090", "ver", "geom", "rotd0", "rotd50", "rotd100", "eas"], ), + "rx": ("station", rx), + "ry": ("station", ry), "rrup": ("station", rrup), "rjb": ("station", rjb), "hyp": ("station", hyp), "epi": ("station", epi), }, - attrs={"hypo_lat": hypocentre[0], "hypo_lon": hypocentre[1]}, + attrs={ + "hypo_lat": hypocentre[0], + "hypo_lon": hypocentre[1], + "source": shapely.to_wkt( + _source_polygon(source_geometries.source_geometries) + ), + "trace": shapely.to_wkt( + _trace_polygon(source_geometries.source_geometries) + ), + "domain": shapely.to_wkt( + shapely.transform( + domain_parameters.domain.polygon, lambda c: c[:, ::-1] + ) + ), + "magnitude": magnitudes.total_magnitude, + "event": metadata.name, + }, ) waveform = broadband.waveform.values.astype(np.float64) From 4ef0ffda06c8723911a58dacc7110a5259919300 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 20 Jul 2026 13:53:20 +1200 Subject: [PATCH 03/21] add PGD as default IM --- workflow/default_parameters/root/defaults.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/workflow/default_parameters/root/defaults.yaml b/workflow/default_parameters/root/defaults.yaml index 75868dd4..cd06c128 100644 --- a/workflow/default_parameters/root/defaults.yaml +++ b/workflow/default_parameters/root/defaults.yaml @@ -503,7 +503,7 @@ velocity_model_1d: Qp: 394.80 Qs: 197.40 im: - ims: ["PGA", "PGV", "CAV", "AI", "Ds575", "Ds595", "pSA", "FAS"] + ims: ["PGA", "PGV", "PGD", "CAV", "AI", "Ds575", "Ds595", "pSA", "FAS"] valid_periods: [ 0.01, From 2048160ad380b536625a1cd41f95ace6aebaac82 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 20 Jul 2026 14:13:19 +1200 Subject: [PATCH 04/21] fix hypocentre calculations --- workflow/scripts/im_calc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/workflow/scripts/im_calc.py b/workflow/scripts/im_calc.py index bd341162..36a50a3b 100644 --- a/workflow/scripts/im_calc.py +++ b/workflow/scripts/im_calc.py @@ -249,7 +249,7 @@ def calculate_instensity_measures( hyp = ( coordinates.distance_between_wgs_depth_coordinates( - np.stack((station_locations, np.zeros_like(latitude)), axis=1), + np.c_[station_locations, np.zeros_like(latitude)], hypocentre, ) / 1000 From ae2dfe34b0ad4657808e399bae4727246c3ef2f5 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 20 Jul 2026 14:22:46 +1200 Subject: [PATCH 05/21] fix rx, ry ordering --- workflow/scripts/im_calc.py | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/workflow/scripts/im_calc.py b/workflow/scripts/im_calc.py index 36a50a3b..666fb14b 100644 --- a/workflow/scripts/im_calc.py +++ b/workflow/scripts/im_calc.py @@ -75,6 +75,7 @@ def _source_polygon(source_geometries: dict[str, IsSource]) -> shapely.Geometry: geometries = [] for fault in source_geometries.values(): geometry = fault.geometry + geometry = shapely.transform( geometry, lambda c: coordinates.nztm_to_wgs_depth(c)[:, ::-1] ) @@ -261,25 +262,14 @@ def calculate_instensity_measures( ) / 1000 ) - all_faults_have_rx_ry = all( - isinstance(source, sources.Plane | sources.Fault) - for source in source_geometries.source_geometries.values() - ) - if all_faults_have_rx_ry: - rx, ry = sources.multi_fault_rx_ry_distance( - list(source_geometries.source_geometries.values()), # ty: ignore[invalid-argument-type] - station_locations, - ) - + stations = broadband.station.values dataset = xr.Dataset( coords={ - "station": ("station", broadband.station.values), + "station": ("station", stations), "component": ( "component", ["000", "090", "ver", "geom", "rotd0", "rotd50", "rotd100", "eas"], ), - "rx": ("station", rx), - "ry": ("station", ry), "rrup": ("station", rrup), "rjb": ("station", rjb), "hyp": ("station", hyp), @@ -304,6 +294,18 @@ def calculate_instensity_measures( }, ) + all_faults_have_rx_ry = all( + isinstance(source, sources.Plane | sources.Fault) + for source in source_geometries.source_geometries.values() + ) + if all_faults_have_rx_ry: + rx, ry = sources.multi_fault_rx_ry_distance( + list(source_geometries.source_geometries.values()), # ty: ignore[invalid-argument-type] + station_locations, + ) + dataset["rx"] = xr.DataArray(rx, dims="station", coords=dict(station=stations)) + dataset["ry"] = xr.DataArray(ry, dims="station", coords=dict(station=stations)) + waveform = broadband.waveform.values.astype(np.float64) for im_name in (pbar := tqdm.tqdm(intensity_measures)): From 3fe834b9a14dfc9bd968df580ddbb513cfd91dce Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 20 Jul 2026 14:23:24 +1200 Subject: [PATCH 06/21] Apply suggestions from code review Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- workflow/realisations.py | 3 +-- workflow/scripts/im_calc.py | 3 ++- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/workflow/realisations.py b/workflow/realisations.py index ba5f3fde..896a3537 100644 --- a/workflow/realisations.py +++ b/workflow/realisations.py @@ -697,8 +697,7 @@ def total_moment(self) -> float: @property def total_magnitude(self) -> BoldM: - """float: total magnitude of realisation""" - return moment.moment_to_magnitude(self.total_moment, bold_m=True) + """BoldM: total magnitude of realisation""" @dataclasses.dataclass diff --git a/workflow/scripts/im_calc.py b/workflow/scripts/im_calc.py index 666fb14b..d9111950 100644 --- a/workflow/scripts/im_calc.py +++ b/workflow/scripts/im_calc.py @@ -286,7 +286,8 @@ def calculate_instensity_measures( ), "domain": shapely.to_wkt( shapely.transform( - domain_parameters.domain.polygon, lambda c: c[:, ::-1] + domain_parameters.domain.polygon, + lambda c: coordinates.nztm_to_wgs_depth(c)[:, ::-1], ) ), "magnitude": magnitudes.total_magnitude, From f71466c718ef21e7e8df2b9bf3e2e8da7ee433e3 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 20 Jul 2026 14:26:07 +1200 Subject: [PATCH 07/21] fix name --- workflow/scripts/im_calc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/workflow/scripts/im_calc.py b/workflow/scripts/im_calc.py index d9111950..479da103 100644 --- a/workflow/scripts/im_calc.py +++ b/workflow/scripts/im_calc.py @@ -112,7 +112,7 @@ def _trace_polygon(source_geometries: dict[str, IsSource]) -> shapely.Geometry: @cli.from_docstring(app) -def calculate_instensity_measures( +def calculate_intensity_measures( realisation_ffp: Annotated[ Path, typer.Argument(exists=True, dir_okay=False, writable=True) ], From b5095a4308f3a728f6fd0ec7b439a5ad5ed71aa0 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 20 Jul 2026 14:35:06 +1200 Subject: [PATCH 08/21] Add in total magnitude --- workflow/realisations.py | 1 + 1 file changed, 1 insertion(+) diff --git a/workflow/realisations.py b/workflow/realisations.py index 896a3537..e71a8322 100644 --- a/workflow/realisations.py +++ b/workflow/realisations.py @@ -698,6 +698,7 @@ def total_moment(self) -> float: @property def total_magnitude(self) -> BoldM: """BoldM: total magnitude of realisation""" + return moment.moment_to_magnitude(self.total_moment, bold_m=True) @dataclasses.dataclass From 51657dcacb0a9c8b360c70c4e94d37bb62357297 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 20 Jul 2026 14:36:53 +1200 Subject: [PATCH 09/21] ignore transform type error --- workflow/scripts/im_calc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/workflow/scripts/im_calc.py b/workflow/scripts/im_calc.py index 479da103..a1e5b19f 100644 --- a/workflow/scripts/im_calc.py +++ b/workflow/scripts/im_calc.py @@ -104,7 +104,7 @@ def _trace_polygon(source_geometries: dict[str, IsSource]) -> shapely.Geometry: geometry = fault.trace_geometry geometry = shapely.transform( geometry, lambda c: coordinates.nztm_to_wgs_depth(c)[:, ::-1] - ) + ) # ty: ignore[no-matching-overload] geometries.append(geometry) From 17ee277b13b3a893cb22d7aefb7222b54128cd55 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 20 Jul 2026 14:37:31 +1200 Subject: [PATCH 10/21] ignore numpydoc error --- workflow/realisations.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/workflow/realisations.py b/workflow/realisations.py index e71a8322..0cb26307 100644 --- a/workflow/realisations.py +++ b/workflow/realisations.py @@ -688,7 +688,7 @@ def __getitem__(self, key: str) -> BoldM: return self.magnitudes[key] @property - def total_moment(self) -> float: + def total_moment(self) -> float: # numpydoc ignore=RT01 """float: total moment of realisation""" return sum( moment.magnitude_to_moment(mag, bold_m=True) @@ -696,7 +696,7 @@ def total_moment(self) -> float: ) @property - def total_magnitude(self) -> BoldM: + def total_magnitude(self) -> BoldM: # numpydoc ignore=RT01 """BoldM: total magnitude of realisation""" return moment.moment_to_magnitude(self.total_moment, bold_m=True) From e01b50229c1c5e1048223a595a402cd23de8efef Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 Jul 2026 09:51:16 +1200 Subject: [PATCH 11/21] add empirical class --- workflow/realisations.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/workflow/realisations.py b/workflow/realisations.py index 0cb26307..a5183e45 100644 --- a/workflow/realisations.py +++ b/workflow/realisations.py @@ -1223,6 +1223,17 @@ def to_dict(self) -> dict[str, Any]: return _dict +@dataclasses.dataclass +class EmpiricalParameters: + _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 + models: list[Any] + + @dataclasses.dataclass class LogEntry: """Log entry for workflow utilities.""" From f59df67760fcdc642cdf01c512888e9673734cfa Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 Jul 2026 11:18:46 +1200 Subject: [PATCH 12/21] extract out datatree logic --- workflow/scripts/im_calc.py | 200 ++++++++++++++++++++++++++++-------- 1 file changed, 157 insertions(+), 43 deletions(-) diff --git a/workflow/scripts/im_calc.py b/workflow/scripts/im_calc.py index a1e5b19f..b825da74 100644 --- a/workflow/scripts/im_calc.py +++ b/workflow/scripts/im_calc.py @@ -29,17 +29,22 @@ import functools from pathlib import Path -from typing import Annotated +from typing import Annotated, Any import numpy as np import pandas as pd + +# Importing pint_xarray registers pint units with xarray, allowing for +# unit-aware operations. It is not explicitly used, so we ignore the flake8 +# F401 error. +import pint_xarray # noqa: F401 import shapely import tqdm import typer import xarray as xr -from IM import im_reader, ims -from IM.im_calculation import IM +from IM import ims +from IM.ims import IM from qcore import cli, coordinates from source_modelling import sources from source_modelling.sources import IsSource @@ -59,6 +64,114 @@ app = typer.Typer() +COORDINATE_METADATA = { + "station": {"description": "Station identifiers"}, + "period": {"description": "Oscillation period", "units": "s"}, + "vs30": { + "description": "Average shear-wave velocity to 30m depth", + "units": "m/s", + }, + "epi": {"description": "Epicentral distance", "units": "km"}, + "hyp": {"description": "Hypocentral distance", "units": "km"}, + "rrup": {"description": "Rupture distance", "units": "km"}, + "rjb": {"description": "Joyner-Boore distance", "units": "km"}, + "rx": {"description": "Generalised strike-parallel distance", "units": "km"}, + "ry": {"description": "Generalised strike-normal distance", "units": "km"}, + "latitude": {"description": "Station latitude", "units": "degrees"}, + "longitude": {"description": "Station longitude", "units": "degrees"}, + "frequency": {"description": "Frequency of motion", "units": "Hz"}, + "period": {"description": "Period of motion", "units": "s"}, +} + +IM_METADATA = { + IM.PGA: "Peak ground acceleration", + IM.PGV: "Peak ground velocity", + IM.CAV: "Cumulative absolute velocity", + IM.CAV5: "Cumulative absolute velocity (above 5 cm/s)", + IM.AI: "Arias intensity", + IM.Ds575: "Significant duration (5-75%)", + IM.Ds595: "Significant duration (5-95%)", + IM.pSA: "Pseudo-spectral acceleration", + IM.FAS: "Fourier amplitude spectrum", +} + + +# The 'g0' unit is used for acceleration and is equivalent to 9.81 m/s^2. The +# reason for this is that 'g' is reserved for 'grams'. This is a decision +# made by the `pint` library, which is used to handle the units. +IM_UNITS = { + IM.PGA: "g0", + IM.PGV: "cm/s", + IM.CAV: "m/s", + IM.CAV5: "m/s", + IM.AI: "m/s", + IM.Ds575: "s", + IM.Ds595: "s", + IM.FAS: "g0 * s", + IM.pSA: "g0", +} + + +def add_distances( + dtree: xr.DataTree, distances: dict[str, xr.DataArray] +) -> xr.DataTree: + """Write intensity measures to a file, updating coordinate and variable metadata. + + Parameters + ---------- + dataset : xr.Dataset + The xarray dataset containing intensity measures to be written. + """ + + def distancify(dataset: xr.Dataset) -> xr.Dataset: + if "name" not in dataset.attrs: + return dataset + + dataset = dataset.copy(deep=False) + dataset.coords.update(distances) + return dataset + + dtree = dtree.map_over_datasets(distancify) + + return dtree + + +def add_units(dtree: xr.DataTree) -> xr.DataTree: + """Write intensity measures to a file, updating coordinate and variable metadata. + + Parameters + ---------- + dataset : xr.Dataset + The xarray dataset containing intensity measures to be written. + output_ffp : str or Path + The file path where the output dataset should be saved. + """ + + def unitify(dataset: xr.Dataset) -> xr.Dataset: + if "name" not in dataset.attrs: + return dataset + + dataset = dataset.copy(deep=False) + + for name, description in COORDINATE_METADATA.items(): + if name not in dataset.coords: + continue + dataset.coords[name].attrs.update(description) + + name = dataset.attrs["name"] + + for data_var in dataset.data_vars.values(): + data_var.attrs["units"] = IM_UNITS[name] + + description = IM_METADATA[name] + dataset.attrs["description"] = description + return dataset + + dtree = dtree.map_over_datasets(unitify) + + return dtree + + def _source_polygon(source_geometries: dict[str, IsSource]) -> shapely.Geometry: """Extract source polygon in longitude, latitude format. @@ -263,38 +376,12 @@ def calculate_intensity_measures( / 1000 ) stations = broadband.station.values - dataset = xr.Dataset( - coords={ - "station": ("station", stations), - "component": ( - "component", - ["000", "090", "ver", "geom", "rotd0", "rotd50", "rotd100", "eas"], - ), - "rrup": ("station", rrup), - "rjb": ("station", rjb), - "hyp": ("station", hyp), - "epi": ("station", epi), - }, - attrs={ - "hypo_lat": hypocentre[0], - "hypo_lon": hypocentre[1], - "source": shapely.to_wkt( - _source_polygon(source_geometries.source_geometries) - ), - "trace": shapely.to_wkt( - _trace_polygon(source_geometries.source_geometries) - ), - "domain": shapely.to_wkt( - shapely.transform( - domain_parameters.domain.polygon, - lambda c: coordinates.nztm_to_wgs_depth(c)[:, ::-1], - ) - ), - "magnitude": magnitudes.total_magnitude, - "event": metadata.name, - }, - ) - + distances: dict[str, Any] = { + "rrup": ("station", rrup), + "rjb": ("station", rjb), + "hyp": ("station", hyp), + "epi": ("station", epi), + } all_faults_have_rx_ry = all( isinstance(source, sources.Plane | sources.Fault) for source in source_geometries.source_geometries.values() @@ -304,11 +391,17 @@ def calculate_intensity_measures( list(source_geometries.source_geometries.values()), # ty: ignore[invalid-argument-type] station_locations, ) - dataset["rx"] = xr.DataArray(rx, dims="station", coords=dict(station=stations)) - dataset["ry"] = xr.DataArray(ry, dims="station", coords=dict(station=stations)) + rx /= 1000.0 + ry /= 1000.0 + distances["rx"] = xr.DataArray( + rx, dims="station", coords=dict(station=stations) + ) + distances["ry"] = xr.DataArray( + ry, dims="station", coords=dict(station=stations) + ) waveform = broadband.waveform.values.astype(np.float64) - + im_results: dict[str, xr.Dataset] = dict() for im_name in (pbar := tqdm.tqdm(intensity_measures)): pbar.set_description(im_name) im_fn = im_function_map[im_name] @@ -317,10 +410,31 @@ def calculate_intensity_measures( if isinstance(result, pd.DataFrame): result["station"] = broadband.station.values - result = result.set_index("station").to_xarray().to_array(dim="component") + result = result.set_index("station").to_xarray() elif isinstance(result, xr.DataArray): - result = result.assign_coords(station=broadband.station) - dataset[im_name] = result - im_reader.write_intensity_measures(dataset, output_path) - + result = result.assign_coords(station=broadband.station).to_dataset( + "component" + ) + result.attrs["name"] = im_name + im_results[im_name] = result + + dtree = xr.DataTree.from_dict(im_results, nested=True) + + dtree.attrs = { + "hypo_lat": hypocentre[0], + "hypo_lon": hypocentre[1], + "source": shapely.to_wkt(_source_polygon(source_geometries.source_geometries)), + "trace": shapely.to_wkt(_trace_polygon(source_geometries.source_geometries)), + "domain": shapely.to_wkt( + shapely.transform( + domain_parameters.domain.polygon, + lambda c: coordinates.nztm_to_wgs_depth(c)[:, ::-1], + ) + ), + "magnitude": magnitudes.total_magnitude, + "event": metadata.name, + } + dtree = add_distances(dtree, distances) + dtree = add_units(dtree) + dtree.to_netcdf(output_path) realisations.append_log_entry(realisation_ffp) From 7f2e90b4bf7bd67e32ada949f5b5ca23892c61f6 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 Jul 2026 12:48:25 +1200 Subject: [PATCH 13/21] add empricial measure calculations to ims --- pyproject.toml | 9 +- uv.lock | 2 +- .../default_parameters/root/defaults.yaml | 3 + workflow/realisations.py | 55 +- workflow/schemas.py | 22 + workflow/scripts/bb_sim.py | 1 + workflow/scripts/im_calc.py | 613 +++++++++++++++--- 7 files changed, 609 insertions(+), 96 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 8a4674bd..0abf4c95 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,6 +19,7 @@ dependencies = [ "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", ] diff --git a/uv.lock b/uv.lock index 48637ded..640c9aa6 100644 --- a/uv.lock +++ b/uv.lock @@ -225,7 +225,7 @@ name = "cffi" version = "2.0.0" source = { registry = "https://pypi.org/simple" } dependencies = [ - { name = "pycparser", marker = "implementation_name != 'PyPy'" }, + { name = "pycparser" }, ] sdist = { url = "https://files.pythonhosted.org/packages/eb/56/b1ba7935a17738ae8453301356628e8147c79dbb825bcbc73dc7401f9846/cffi-2.0.0.tar.gz", hash = "sha256:44d1b5909021139fe36001ae048dbdde8214afa20200eda0f64c068cac5d5529", size = 523588, upload-time = "2025-09-08T23:24:04.541Z" } wheels = [ diff --git a/workflow/default_parameters/root/defaults.yaml b/workflow/default_parameters/root/defaults.yaml index e02a5a76..f28a99c2 100644 --- a/workflow/default_parameters/root/defaults.yaml +++ b/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: diff --git a/workflow/realisations.py b/workflow/realisations.py index a5183e45..8f08ded8 100644 --- a/workflow/realisations.py +++ b/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""" @@ -1224,14 +1273,18 @@ def to_dict(self) -> dict[str, Any]: @dataclasses.dataclass -class EmpiricalParameters: +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 diff --git a/workflow/schemas.py b/workflow/schemas.py index 5df560cf..f22dee13 100644 --- a/workflow/schemas.py +++ b/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( diff --git a/workflow/scripts/bb_sim.py b/workflow/scripts/bb_sim.py index 9b03ae3e..66b572f1 100644 --- a/workflow/scripts/bb_sim.py +++ b/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", diff --git a/workflow/scripts/im_calc.py b/workflow/scripts/im_calc.py index b825da74..fe6842d6 100644 --- a/workflow/scripts/im_calc.py +++ b/workflow/scripts/im_calc.py @@ -27,21 +27,27 @@ See the output of `im-calc --help`. """ +import dataclasses import functools +import warnings from pathlib import Path -from typing import Annotated, Any +from typing import Annotated +# NOTE: netCDF4 must be imported before OpenQuake (which oq_wrapper imports). +# OpenQuake pulls in h5py, which is linked against a different build of HDF5 +# than netCDF4 is. Whichever of the two loads second fails to open our +# waveform files, so the import order here is load bearing. +import netCDF4 # noqa: F401 import numpy as np +import numpy.typing as npt +import oq_wrapper as oqw +import oq_wrapper.xarray as oqwx import pandas as pd - -# Importing pint_xarray registers pint units with xarray, allowing for -# unit-aware operations. It is not explicitly used, so we ignore the flake8 -# F401 error. -import pint_xarray # noqa: F401 import shapely import tqdm import typer import xarray as xr +from oq_wrapper.estimations import chiou_young_08_calc_z1p0, chiou_young_08_calc_z2p5 from IM import ims from IM.ims import IM @@ -51,8 +57,10 @@ from workflow import realisations, utils from workflow.realisations import ( DomainParameters, + EmpiricalParameters, IntensityMeasureCalculationParameters, Magnitudes, + Rakes, RealisationMetadata, Resolution, RupturePropagationConfig, @@ -66,11 +74,18 @@ COORDINATE_METADATA = { "station": {"description": "Station identifiers"}, - "period": {"description": "Oscillation period", "units": "s"}, "vs30": { "description": "Average shear-wave velocity to 30m depth", "units": "m/s", }, + "z1pt0": { + "description": "Depth to the 1.0 km/s shear-wave velocity horizon", + "units": "km", + }, + "z2pt5": { + "description": "Depth to the 2.5 km/s shear-wave velocity horizon", + "units": "km", + }, "epi": {"description": "Epicentral distance", "units": "km"}, "hyp": {"description": "Hypocentral distance", "units": "km"}, "rrup": {"description": "Rupture distance", "units": "km"}, @@ -112,43 +127,78 @@ } -def add_distances( - dtree: xr.DataTree, distances: dict[str, xr.DataArray] +EMPIRICAL_IM_NAMES = { + IM.PGA: "PGA", + IM.PGV: "PGV", + IM.CAV: "CAV", + IM.AI: "AI", + IM.Ds575: "Ds575", + IM.Ds595: "Ds595", + IM.pSA: "pSA", +} + + +EMPIRICAL_STATISTIC_METADATA = { + "mean": "Mean of the natural logarithm of {description}", + "std_Total": "Total standard deviation of the natural logarithm of {description}", + "std_Inter": "Between-event standard deviation of the natural logarithm of {description}", + "std_Intra": "Within-event standard deviation of the natural logarithm of {description}", +} + + +def add_station_parameters( + dtree: xr.DataTree, station_parameters: dict[str, xr.DataArray] ) -> xr.DataTree: - """Write intensity measures to a file, updating coordinate and variable metadata. + """Attach per-station parameters as coordinates on every leaf of the tree. Parameters ---------- - dataset : xr.Dataset - The xarray dataset containing intensity measures to be written. + dtree : xr.DataTree + The tree of intensity measure datasets. + station_parameters : dict + A map from parameter name (distance and site measures) to the + per-station values. + + Returns + ------- + xr.DataTree + The tree, with the parameters attached to every dataset containing + data. """ - def distancify(dataset: xr.Dataset) -> xr.Dataset: - if "name" not in dataset.attrs: + def parameterise(dataset: xr.Dataset) -> xr.Dataset: + if not dataset.data_vars: return dataset dataset = dataset.copy(deep=False) - dataset.coords.update(distances) + dataset.coords.update(station_parameters) return dataset - dtree = dtree.map_over_datasets(distancify) + dtree = dtree.map_over_datasets(parameterise) return dtree def add_units(dtree: xr.DataTree) -> xr.DataTree: - """Write intensity measures to a file, updating coordinate and variable metadata. + """Annotate coordinates and intensity measures with units and descriptions. + + Empirical datasets are left alone, because they are annotated as they + are calculated (their values are in log-space, so they do not share the + units of the simulated intensity measures). Parameters ---------- - dataset : xr.Dataset - The xarray dataset containing intensity measures to be written. - output_ffp : str or Path - The file path where the output dataset should be saved. + dtree : xr.DataTree + The tree of intensity measure datasets. + + Returns + ------- + xr.DataTree + The tree, with unit and description metadata attached. """ def unitify(dataset: xr.Dataset) -> xr.Dataset: - if "name" not in dataset.attrs: + if not dataset.data_vars: return dataset dataset = dataset.copy(deep=False) @@ -158,6 +208,9 @@ def unitify(dataset: xr.Dataset) -> xr.Dataset: continue dataset.coords[name].attrs.update(description) + if "name" not in dataset.attrs: + return dataset + name = dataset.attrs["name"] for data_var in dataset.data_vars.values(): @@ -224,6 +277,423 @@ def _trace_polygon(source_geometries: dict[str, IsSource]) -> shapely.Geometry: return shapely.normalize(shapely.union_all(geometries)) +@dataclasses.dataclass +class Distances: + """Source-to-site distance measures, in kilometres.""" + + rrup: xr.DataArray + """Shortest distance to the rupture plane.""" + rjb: xr.DataArray + """Shortest distance to the surface projection of the rupture.""" + hyp: xr.DataArray + """Distance to the hypocentre.""" + epi: xr.DataArray + """Distance to the epicentre.""" + rx: xr.DataArray | None = None + """Strike-parallel distance. Only defined for planar sources.""" + ry: xr.DataArray | None = None + """Strike-normal distance. Only defined for planar sources.""" + + def as_dict(self) -> dict[str, xr.DataArray]: + """Map distance measure name to distances, omitting undefined measures. + + Returns + ------- + dict + A map from distance measure name to per-station distances. + """ + return { + field.name: value + for field in dataclasses.fields(self) + if (value := getattr(self, field.name)) is not None + } + + +def calculate_distances( + source_geometries: SourceConfig, + hypocentre: np.ndarray, + broadband: xr.Dataset, +) -> Distances: + """Calculate source-to-site distances for every station in the broadband. + + Parameters + ---------- + source_geometries : SourceConfig + The source geometries of the realisation. + hypocentre : np.ndarray + The hypocentre, in latitude, longitude, depth format. + broadband : xr.Dataset + The broadband waveform dataset, supplying the station locations. + + Returns + ------- + Distances + The distance measures for each station, in kilometres. `rx` and `ry` + are only calculated if every source in the realisation is planar. + """ + latitude = broadband.latitude.values + longitude = broadband.longitude.values + station_locations = np.stack((latitude, longitude), axis=-1) + + rrup = xr.DataArray( + np.array( + [ + min( + source.rrup_distance(np.append(station, 0)) + for source in source_geometries.source_geometries.values() + ) + for station in station_locations + ] + ) + / 1000, + dims=["station"], + coords=dict(station=broadband.station), + ) + rjb = xr.DataArray( + np.array( + [ + min( + source.rjb_distance(np.append(station, 0)) + for source in source_geometries.source_geometries.values() + ) + for station in station_locations + ] + ) + / 1000, + dims=["station"], + coords=dict(station=broadband.station), + ) + + hyp = xr.DataArray( + coordinates.distance_between_wgs_depth_coordinates( + np.c_[station_locations, np.zeros_like(latitude)], + hypocentre, + ) + / 1000, + dims=["station"], + coords=dict(station=broadband.station), + ) + epi = xr.DataArray( + coordinates.distance_between_wgs_depth_coordinates( + station_locations, + hypocentre[:2], + ) + / 1000, + dims=["station"], + coords=dict(station=broadband.station), + ) + + distances = Distances(rrup=rrup, rjb=rjb, hyp=hyp, epi=epi) + all_faults_have_rx_ry = all( + isinstance(source, sources.Plane | sources.Fault) + for source in source_geometries.source_geometries.values() + ) + if all_faults_have_rx_ry: + rx, ry = sources.multi_fault_rx_ry_distance( + list(source_geometries.source_geometries.values()), # ty: ignore[invalid-argument-type] + station_locations, + ) + rx /= 1000.0 + ry /= 1000.0 + distances.rx = xr.DataArray( + rx, dims="station", coords=dict(station=broadband.station) + ) + distances.ry = xr.DataArray( + ry, dims="station", coords=dict(station=broadband.station) + ) + return distances + + +@dataclasses.dataclass +class SourceParameters: + """Rupture parameters describing the realisation as a single source.""" + + mag: float + """The total moment magnitude of the rupture.""" + avg_rake: float + """The moment-averaged rake angle (degrees).""" + avg_dip: float + """The moment-averaged dip angle (degrees).""" + avg_ztor: float + """The moment-averaged depth to the top of the rupture (km).""" + avg_zbot: float + """The moment-averaged depth to the bottom of the rupture (km).""" + hypo_depth: float + """The depth of the hypocentre (km).""" + + +def calculate_source_parameters( + source_config: SourceConfig, + magnitudes: Magnitudes, + rakes: Rakes, + hypocentre: np.ndarray, +) -> SourceParameters: + """Reduce a multi-fault realisation to a single set of rupture parameters. + + Ground motion models describe a rupture with a single magnitude, rake, + dip and depth. Multi-fault realisations are collapsed into these by + averaging each fault's contribution, weighted by its moment. + + Parameters + ---------- + source_config : SourceConfig + The source geometries of the realisation. + magnitudes : Magnitudes + The per-fault magnitudes, used for the moment weighting. + rakes : Rakes + The per-fault rake angles. + hypocentre : np.ndarray + The hypocentre, in latitude, longitude, depth (metres) format. + + Returns + ------- + SourceParameters + The rupture parameters of the realisation as a whole. + """ + mag = magnitudes.total_magnitude + + avg_rake_vector = magnitudes.moment_averaged(rakes.as_vectors()) + avg_rake = np.degrees(np.arctan2(avg_rake_vector[1], avg_rake_vector[0])) + + avg_dip_vector = magnitudes.moment_averaged( + { + k: np.array([np.cos(np.radians(f.dip)), np.sin(np.radians(f.dip))]) + for k, f in source_config.source_geometries.items() + } + ) + avg_dip = np.degrees(np.arctan2(avg_dip_vector[1], avg_dip_vector[0])) + + if all( + hasattr(f, "top_m") and hasattr(f, "bottom_m") + for f in source_config.source_geometries.values() + ): + avg_ztor = magnitudes.moment_averaged( + {k: f.top_m / 1000.0 for k, f in source_config.source_geometries.items()} # ty: ignore[unresolved-attribute] + ) + avg_zbot = magnitudes.moment_averaged( + {k: f.bottom_m / 1000.0 for k, f in source_config.source_geometries.items()} # ty: ignore[unresolved-attribute] + ) + else: + avg_ztor = magnitudes.moment_averaged( + { + k: f.centroid[-1] / 1000.0 + for k, f in source_config.source_geometries.items() + } + ) + avg_zbot = avg_ztor + + return SourceParameters( + mag=mag, + avg_rake=float(avg_rake), + avg_dip=float(avg_dip), + avg_ztor=avg_ztor, + avg_zbot=avg_zbot, + hypo_depth=float(hypocentre[2]) / 1000.0, + ) + + +@dataclasses.dataclass +class SiteParameters: + """Per-station site parameters.""" + + vs30: xr.DataArray + """Average shear-wave velocity to 30m depth (m/s).""" + z1pt0: xr.DataArray + """Depth to the 1.0 km/s shear-wave velocity horizon (km).""" + z2pt5: xr.DataArray + """Depth to the 2.5 km/s shear-wave velocity horizon (km).""" + + def as_dict(self) -> dict[str, xr.DataArray]: + """Map site parameter name to per-station values. + + Returns + ------- + dict + A map from site parameter name to per-station values. + """ + return { + field.name: getattr(self, field.name) for field in dataclasses.fields(self) + } + + +def calculate_site_parameters(vs30: xr.DataArray) -> SiteParameters: + """Estimate site parameters from vs30. + + Parameters + ---------- + vs30 : xr.DataArray + The per-station vs30 values (m/s). + + Returns + ------- + SiteParameters + The site parameters, with basin depths estimated using the Chiou + and Youngs (2008) relations. + """ + z1pt0 = chiou_young_08_calc_z1p0(vs30) + z2pt5 = chiou_young_08_calc_z2p5(z1pt0) + return SiteParameters(vs30=vs30, z1pt0=z1pt0, z2pt5=z2pt5) + + +def empirical_inputs( + source_parameters: SourceParameters, + site_parameters: SiteParameters, + distances: Distances, +) -> xr.Dataset: + """Assemble the rupture context ground motion models are evaluated against. + + Parameters + ---------- + source_parameters : SourceParameters + The rupture parameters of the realisation. + site_parameters : SiteParameters + The per-station site parameters. + distances : Distances + The per-station source-to-site distances. + + Returns + ------- + xr.Dataset + A dataset of OpenQuake rupture context variables. Site and distance + variables vary over the station dimension, rupture variables are + scalars. + """ + return xr.Dataset( + dict( + mag=source_parameters.mag, + dip=source_parameters.avg_dip, + rake=source_parameters.avg_rake, + ztor=source_parameters.avg_ztor, + zbot=source_parameters.avg_zbot, + hypo_depth=source_parameters.hypo_depth, + vs30=site_parameters.vs30, + z1pt0=site_parameters.z1pt0, + z2pt5=site_parameters.z2pt5, + vs30measured=False, + # TODO: Calculate backarc! + backarc=False, + ) + | distances.as_dict() + ) + + +def annotate_empirical(dataset: xr.Dataset, im_name: IM, model_name: str) -> xr.Dataset: + """Attach units and descriptions to an empirical intensity measure dataset. + + Ground motion models predict the distribution of the natural logarithm + of an intensity measure, so the values are dimensionless. The units of + the intensity measure itself are recorded in the `log_units` attribute. + + Parameters + ---------- + dataset : xr.Dataset + The dataset of statistics returned by `oq_wrapper`. + im_name : IM + The intensity measure the dataset describes. + model_name : str + The ground motion model (or logic tree) that produced the dataset. + + Returns + ------- + xr.Dataset + The dataset, with metadata attached. + """ + dataset = dataset.copy(deep=False) + description = IM_METADATA[im_name] + + for statistic, data_var in dataset.data_vars.items(): + data_var.attrs["units"] = "dimensionless" + data_var.attrs["log_units"] = IM_UNITS[im_name] + if statistic in EMPIRICAL_STATISTIC_METADATA: + data_var.attrs["description"] = EMPIRICAL_STATISTIC_METADATA[ # ty: ignore[invalid-argument-type] + statistic + ].format(description=description) + + dataset.attrs["intensity_measure"] = str(im_name) + dataset.attrs["model"] = model_name + dataset.attrs["description"] = ( + f"{description} predicted by the {model_name} ground motion model" + ) + return dataset + + +def calculate_empirical( + empirical_config: EmpiricalParameters, + source_parameters: SourceParameters, + site_parameters: SiteParameters, + distances: Distances, + intensity_measures: list[IM], + periods: npt.NDArray[np.float64], +) -> dict[str, xr.Dataset]: + """Calculate empirical intensity measures from ground motion models. + + Each model in the empirical configuration is evaluated for each + intensity measure a ground motion model can predict. Combinations that + a model does not support are skipped with a warning. + + Parameters + ---------- + empirical_config : EmpiricalParameters + The tectonic type and models to evaluate. + source_parameters : SourceParameters + The rupture parameters of the realisation. + site_parameters : SiteParameters + The per-station site parameters. + distances : Distances + The per-station source-to-site distances. + intensity_measures : list of IM + The intensity measures to calculate. + periods : np.ndarray + The periods to calculate pSA at. + + Returns + ------- + dict + A map from data tree path (`{im}/empirical/{model}`) to the log-mean + and log-standard deviation of that intensity measure. The paths are + chosen so this map can be merged with the simulated intensity + measures before building the output data tree. + """ + inputs = empirical_inputs(source_parameters, site_parameters, distances) + tect_type = oqw.constants.TectType(empirical_config.tect_type) + + empirical_results: dict[str, xr.Dataset] = {} + model_ims = [im for im in intensity_measures if im in EMPIRICAL_IM_NAMES] + + for model_name in empirical_config.models: + for im_name in (pbar := tqdm.tqdm(model_ims)): + pbar.set_description(f"{model_name} {im_name}") + try: + if model_name in oqw.constants.GMMLogicTree.__members__: + dataset = oqwx.run_gmm_logic_tree_xarray( + oqw.constants.GMMLogicTree[model_name], + tect_type, + inputs, + EMPIRICAL_IM_NAMES[im_name], + periods=periods.tolist(), + ) + else: + dataset = oqwx.run_gmm_xarray( + oqw.constants.GMM[model_name], + tect_type, + inputs, + EMPIRICAL_IM_NAMES[im_name], + periods=periods.tolist(), + ) + except (ValueError, KeyError, AttributeError) as e: + warnings.warn( + f"Skipping empirical {im_name} for {model_name}: {e}", + stacklevel=1, + ) + continue + + empirical_results[f"{im_name}/empirical/{model_name}"] = annotate_empirical( + dataset, im_name, model_name + ) + + return empirical_results + + @cli.from_docstring(app) def calculate_intensity_measures( realisation_ffp: Annotated[ @@ -240,6 +710,7 @@ def calculate_intensity_measures( ] = None, override_ims: Annotated[list[IM] | None, typer.Option("-i", "--im")] = None, cores: Annotated[int | None, typer.Option(min=1)] = None, + empirical: Annotated[bool, typer.Option()] = True, ) -> None: """Calculate intensity measures for simulation data. @@ -263,6 +734,10 @@ def calculate_intensity_measures( Set the number of cores for parallel processing of IMs. If set to `None`, will default to the available cores from `utils.get_available_cores`. + empirical : bool, default True + If passed, additionally estimate intensity measures from the ground + motion models in the realisation file. Requires the broadband + waveforms to carry a `vs30` coordinate. """ cores = cores or utils.get_available_cores() @@ -329,77 +804,15 @@ def calculate_intensity_measures( cores=cores, ), } - latitude = broadband.latitude.values - longitude = broadband.longitude.values - station_locations = np.stack((latitude, longitude), axis=-1) - - rrup = ( - np.array( - [ - min( - source.rrup_distance(np.append(station, 0)) - for source in source_geometries.source_geometries.values() - ) - for station in station_locations - ] - ) - / 1000 - ) - rjb = ( - np.array( - [ - min( - source.rjb_distance(np.append(station, 0)) - for source in source_geometries.source_geometries.values() - ) - for station in station_locations - ] - ) - / 1000 - ) hypocentre = source_geometries.source_geometries[ rup_prop_config.initial_fault ].fault_coordinates_to_wgs_depth_coordinates(rup_prop_config.hypocentre) - - hyp = ( - coordinates.distance_between_wgs_depth_coordinates( - np.c_[station_locations, np.zeros_like(latitude)], - hypocentre, - ) - / 1000 + distances = calculate_distances(source_geometries, hypocentre, broadband) + rakes = Rakes.read_from_realisation(realisation_ffp) + source_parameters = calculate_source_parameters( + source_geometries, magnitudes, rakes, hypocentre ) - epi = ( - coordinates.distance_between_wgs_depth_coordinates( - station_locations, - hypocentre[:2], - ) - / 1000 - ) - stations = broadband.station.values - distances: dict[str, Any] = { - "rrup": ("station", rrup), - "rjb": ("station", rjb), - "hyp": ("station", hyp), - "epi": ("station", epi), - } - all_faults_have_rx_ry = all( - isinstance(source, sources.Plane | sources.Fault) - for source in source_geometries.source_geometries.values() - ) - if all_faults_have_rx_ry: - rx, ry = sources.multi_fault_rx_ry_distance( - list(source_geometries.source_geometries.values()), # ty: ignore[invalid-argument-type] - station_locations, - ) - rx /= 1000.0 - ry /= 1000.0 - distances["rx"] = xr.DataArray( - rx, dims="station", coords=dict(station=stations) - ) - distances["ry"] = xr.DataArray( - ry, dims="station", coords=dict(station=stations) - ) - + site_parameters = calculate_site_parameters(broadband.vs30.astype(np.float64)) waveform = broadband.waveform.values.astype(np.float64) im_results: dict[str, xr.Dataset] = dict() for im_name in (pbar := tqdm.tqdm(intensity_measures)): @@ -418,6 +831,19 @@ def calculate_intensity_measures( result.attrs["name"] = im_name im_results[im_name] = result + if empirical: + empirical_parameters = EmpiricalParameters.read_from_realisation_or_defaults( + realisation_ffp, metadata.defaults_version + ) + im_results |= calculate_empirical( + empirical_parameters, + source_parameters, + site_parameters, + distances, + intensity_measures, + np.array(intensity_measure_parameters.valid_periods, dtype=np.float64), + ) + dtree = xr.DataTree.from_dict(im_results, nested=True) dtree.attrs = { @@ -433,8 +859,15 @@ def calculate_intensity_measures( ), "magnitude": magnitudes.total_magnitude, "event": metadata.name, + "rake": source_parameters.avg_rake, + "dip": source_parameters.avg_dip, + "ztor": source_parameters.avg_ztor, + "zbot": source_parameters.avg_zbot, + "hypo_depth": source_parameters.hypo_depth, } - dtree = add_distances(dtree, distances) + dtree = add_station_parameters( + dtree, distances.as_dict() | site_parameters.as_dict() + ) dtree = add_units(dtree) dtree.to_netcdf(output_path) realisations.append_log_entry(realisation_ffp) From 33ab6d2dec23cc21586ff3907a4b868a60a3ef8f Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 Jul 2026 16:29:57 +1200 Subject: [PATCH 14/21] include empirical parameter tectonic type --- workflow/scripts/im_calc.py | 1 + 1 file changed, 1 insertion(+) diff --git a/workflow/scripts/im_calc.py b/workflow/scripts/im_calc.py index fe6842d6..115cfb4d 100644 --- a/workflow/scripts/im_calc.py +++ b/workflow/scripts/im_calc.py @@ -864,6 +864,7 @@ def calculate_intensity_measures( "ztor": source_parameters.avg_ztor, "zbot": source_parameters.avg_zbot, "hypo_depth": source_parameters.hypo_depth, + "tect_type": str(empirical_parameters.tect_type), } dtree = add_station_parameters( dtree, distances.as_dict() | site_parameters.as_dict() From bc2fbe417c752b82776a027323b0ebf24652fb1c Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 Jul 2026 16:31:41 +1200 Subject: [PATCH 15/21] fix empirical tect type attributes --- workflow/scripts/im_calc.py | 35 ++++++++++++++++++----------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/workflow/scripts/im_calc.py b/workflow/scripts/im_calc.py index 115cfb4d..09f67ae5 100644 --- a/workflow/scripts/im_calc.py +++ b/workflow/scripts/im_calc.py @@ -831,22 +831,7 @@ def calculate_intensity_measures( result.attrs["name"] = im_name im_results[im_name] = result - if empirical: - empirical_parameters = EmpiricalParameters.read_from_realisation_or_defaults( - realisation_ffp, metadata.defaults_version - ) - im_results |= calculate_empirical( - empirical_parameters, - source_parameters, - site_parameters, - distances, - intensity_measures, - np.array(intensity_measure_parameters.valid_periods, dtype=np.float64), - ) - - dtree = xr.DataTree.from_dict(im_results, nested=True) - - dtree.attrs = { + attributes = { "hypo_lat": hypocentre[0], "hypo_lon": hypocentre[1], "source": shapely.to_wkt(_source_polygon(source_geometries.source_geometries)), @@ -864,8 +849,24 @@ def calculate_intensity_measures( "ztor": source_parameters.avg_ztor, "zbot": source_parameters.avg_zbot, "hypo_depth": source_parameters.hypo_depth, - "tect_type": str(empirical_parameters.tect_type), } + if empirical: + empirical_parameters = EmpiricalParameters.read_from_realisation_or_defaults( + realisation_ffp, metadata.defaults_version + ) + im_results |= calculate_empirical( + empirical_parameters, + source_parameters, + site_parameters, + distances, + intensity_measures, + np.array(intensity_measure_parameters.valid_periods, dtype=np.float64), + ) + attributes["tect_type"] = str(empirical_parameters.tect_type) + + dtree = xr.DataTree.from_dict(im_results, nested=True) + + dtree.attrs = attributes dtree = add_station_parameters( dtree, distances.as_dict() | site_parameters.as_dict() ) From aced5d1cf5759431bfb018ebadbe423c689c45b3 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 Jul 2026 16:35:53 +1200 Subject: [PATCH 16/21] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- workflow/scripts/im_calc.py | 1 + 1 file changed, 1 insertion(+) diff --git a/workflow/scripts/im_calc.py b/workflow/scripts/im_calc.py index 09f67ae5..f471f55a 100644 --- a/workflow/scripts/im_calc.py +++ b/workflow/scripts/im_calc.py @@ -101,6 +101,7 @@ IM_METADATA = { IM.PGA: "Peak ground acceleration", IM.PGV: "Peak ground velocity", + IM.PGD: "Peak ground displacement", IM.CAV: "Cumulative absolute velocity", IM.CAV5: "Cumulative absolute velocity (above 5 cm/s)", IM.AI: "Arias intensity", From 01772c5ba64539d073b325573c3449e311f15665 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 Jul 2026 16:36:03 +1200 Subject: [PATCH 17/21] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- workflow/scripts/im_calc.py | 1 + 1 file changed, 1 insertion(+) diff --git a/workflow/scripts/im_calc.py b/workflow/scripts/im_calc.py index f471f55a..464b7c98 100644 --- a/workflow/scripts/im_calc.py +++ b/workflow/scripts/im_calc.py @@ -118,6 +118,7 @@ IM_UNITS = { IM.PGA: "g0", IM.PGV: "cm/s", + IM.PGD: "cm", IM.CAV: "m/s", IM.CAV5: "m/s", IM.AI: "m/s", From 06617f3b372e98a3eae896b185c82bab9ab53264 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 Jul 2026 16:37:29 +1200 Subject: [PATCH 18/21] fix ci checks --- pyproject.toml | 2 +- workflow/scripts/im_calc.py | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 0abf4c95..45f4de1e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -14,7 +14,7 @@ dependencies = [ "im-calculation>=2025.12.5", "velocity-modelling>=2026.2.1", "nshmdb>=2025.12.1", - "oq_wrapper>=2025.12.3", + "oq_wrapper>=2025.05.2", "qcore-utils>=2025.12.2", "source_modelling>=2026.07.2", # Data Formats diff --git a/workflow/scripts/im_calc.py b/workflow/scripts/im_calc.py index 09f67ae5..c3e29a7b 100644 --- a/workflow/scripts/im_calc.py +++ b/workflow/scripts/im_calc.py @@ -166,7 +166,7 @@ def add_station_parameters( data. """ - def parameterise(dataset: xr.Dataset) -> xr.Dataset: + def parameterise(dataset: xr.Dataset) -> xr.Dataset: # numpydoc ignore=GL08 if not dataset.data_vars: return dataset @@ -197,7 +197,7 @@ def add_units(dtree: xr.DataTree) -> xr.DataTree: The tree, with unit and description metadata attached. """ - def unitify(dataset: xr.Dataset) -> xr.Dataset: + def unitify(dataset: xr.Dataset) -> xr.Dataset: # numpydoc ignore=GL08 if not dataset.data_vars: return dataset @@ -530,7 +530,7 @@ def calculate_site_parameters(vs30: xr.DataArray) -> SiteParameters: The site parameters, with basin depths estimated using the Chiou and Youngs (2008) relations. """ - z1pt0 = chiou_young_08_calc_z1p0(vs30) + z1pt0 = chiou_young_08_calc_z1p0(vs30) # ty: ignore[invalid-argument-type] z2pt5 = chiou_young_08_calc_z2p5(z1pt0) return SiteParameters(vs30=vs30, z1pt0=z1pt0, z2pt5=z2pt5) @@ -605,8 +605,8 @@ def annotate_empirical(dataset: xr.Dataset, im_name: IM, model_name: str) -> xr. data_var.attrs["units"] = "dimensionless" data_var.attrs["log_units"] = IM_UNITS[im_name] if statistic in EMPIRICAL_STATISTIC_METADATA: - data_var.attrs["description"] = EMPIRICAL_STATISTIC_METADATA[ # ty: ignore[invalid-argument-type] - statistic + data_var.attrs["description"] = EMPIRICAL_STATISTIC_METADATA[ + str(statistic) ].format(description=description) dataset.attrs["intensity_measure"] = str(im_name) From 11c87274b68fa97dac5e800edc6456963a77e9c9 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 Jul 2026 16:38:22 +1200 Subject: [PATCH 19/21] bump lock file --- uv.lock | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/uv.lock b/uv.lock index 640c9aa6..dc4cc779 100644 --- a/uv.lock +++ b/uv.lock @@ -3128,6 +3128,7 @@ source = { editable = "." } dependencies = [ { name = "geopandas" }, { name = "im-calculation" }, + { name = "netcdf4" }, { name = "nshmdb" }, { name = "numpy" }, { name = "oq-wrapper" }, @@ -3173,10 +3174,11 @@ requires-dist = [ { name = "geopandas" }, { name = "hypothesis", extras = ["numpy"], marker = "extra == 'test'", specifier = ">=6.0.0" }, { name = "im-calculation", specifier = ">=2025.12.5" }, + { name = "netcdf4" }, { name = "nshmdb", specifier = ">=2025.12.1" }, { name = "numpy" }, { name = "numpydoc", marker = "extra == 'dev'" }, - { name = "oq-wrapper", specifier = ">=2025.12.3" }, + { name = "oq-wrapper", specifier = ">=2025.5.2" }, { name = "pandas", extras = ["hdf5", "parquet"] }, { name = "pandas-stubs", marker = "extra == 'types'" }, { name = "parse", specifier = ">=1.21.0" }, From a514a227d9c7a1d87bac949593d40bd0dbae0127 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 Jul 2026 16:39:50 +1200 Subject: [PATCH 20/21] bump oq wrapper properly --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 45f4de1e..9ef98b3c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -14,7 +14,7 @@ dependencies = [ "im-calculation>=2025.12.5", "velocity-modelling>=2026.2.1", "nshmdb>=2025.12.1", - "oq_wrapper>=2025.05.2", + "oq_wrapper>=2026.05.2", "qcore-utils>=2025.12.2", "source_modelling>=2026.07.2", # Data Formats From 55088036e423c8577bfd4f4faa99cb5225d5dc1a Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 Jul 2026 16:40:33 +1200 Subject: [PATCH 21/21] bump lock file --- uv.lock | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/uv.lock b/uv.lock index dc4cc779..37100f88 100644 --- a/uv.lock +++ b/uv.lock @@ -1712,7 +1712,7 @@ numpy = [ [[package]] name = "oq-wrapper" -version = "2025.12.5" +version = "2026.5.2" source = { registry = "https://pypi.org/simple" } dependencies = [ { name = "numpy" }, @@ -1721,10 +1721,11 @@ dependencies = [ { name = "pyyaml" }, { name = "scipy" }, { name = "source-modelling" }, + { name = "xarray" }, ] -sdist = { url = "https://files.pythonhosted.org/packages/d4/8d/748f1fc683ff0d545ea1dfafefef8dda514f7bb14a6186bb4e5bbd753f6f/oq_wrapper-2025.12.5.tar.gz", hash = "sha256:5fa493cd532c7e0383dc07afbfdd823769f068fcc7357dd86d7751076b065f19", size = 179994, upload-time = "2026-01-23T02:23:35.49Z" } +sdist = { url = "https://files.pythonhosted.org/packages/54/55/d8ceb33b167b06f7413eaa5e98f35352b52459fd6247a7823a8e914440a5/oq_wrapper-2026.5.2.tar.gz", hash = "sha256:445a724adc6a7dbbeb05d919d898d08618c043889bb5272855355a9449f174e2", size = 207754, upload-time = "2026-05-14T03:26:39.69Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/96/82/9d2c938a44991f12fdcbca903f4b80fcc233af779dfd028a08c4eacdd552/oq_wrapper-2025.12.5-py3-none-any.whl", hash = "sha256:6fe11df916fe27c786e8690cbd8c90829a38eac5de2f31d86cc298b1fcee4da4", size = 29150, upload-time = "2026-01-23T02:23:34.12Z" }, + { url = "https://files.pythonhosted.org/packages/5c/6e/4a94bbe11b58629b0e1f08a163446b3eba92fbafaf6027d952c920d945d6/oq_wrapper-2026.5.2-py3-none-any.whl", hash = "sha256:ddade827c8eb668d8365aadb912765297c6ddfd5401133658b3b61621f25d99e", size = 32239, upload-time = "2026-05-14T03:26:38.307Z" }, ] [[package]] @@ -3178,7 +3179,7 @@ requires-dist = [ { name = "nshmdb", specifier = ">=2025.12.1" }, { name = "numpy" }, { name = "numpydoc", marker = "extra == 'dev'" }, - { name = "oq-wrapper", specifier = ">=2025.5.2" }, + { name = "oq-wrapper", specifier = ">=2026.5.2" }, { name = "pandas", extras = ["hdf5", "parquet"] }, { name = "pandas-stubs", marker = "extra == 'types'" }, { name = "parse", specifier = ">=1.21.0" },