Skip to content

Add empirical calculations to intensity measures - #127

Closed
lispandfound wants to merge 10 commits into
nzvm/bb-matched-filterfrom
empirical_ims
Closed

lispandfound wants to merge 10 commits into
nzvm/bb-matched-filterfrom
empirical_ims

Conversation

@lispandfound

@lispandfound lispandfound commented Jul 21, 2026 •

Copy link
Copy Markdown
Contributor

Adds empirical calculations to intensity measure outputs where empirical models support a given tectonic type and intensity measure.

To support this structure, I have also refactored im calc so that it finally outputs in xarray data tree format instead of dataset format. The upshot of doing this is that we can uncouple the components so that each intensity measure carries only components it actually computes. Especially for models like NSHM2022 that only support pSA it's kind of silly to have every intensity measure carry all-NaN EAS and empirical model components.

Copilot AI review requested due to automatic review settings July 21, 2026 04:30

@gemini-code-assist gemini-code-assist Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Code Review

This pull request introduces the capability to calculate empirical intensity measures from ground motion models using OpenQuake wrappers. It updates dependencies, default parameters, and schemas, and adds helper methods to average multi-fault parameters (rakes, magnitudes) weighted by fault moment. The im_calc.py script is refactored to calculate source-to-site distances, site parameters, and empirical IMs, outputting the results as an xr.DataTree in NetCDF format. The review feedback highlights critical issues where running the script with empirical=False or with broadband data lacking vs30 coordinates will result in NameError or AttributeError exceptions.

Important

The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.

Comment thread workflow/scripts/im_calc.py
Comment thread workflow/scripts/im_calc.py Outdated
Comment thread workflow/scripts/im_calc.py Outdated

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

Refactors intensity-measure (IM) calculation output to an xarray.DataTree netCDF structure and adds optional empirical (GMM-based) IM calculations driven by per-realisation configuration (tectonic type + model list), carrying through per-station parameters like Vs30 and derived basin depths.

Changes:

  • Refactor im-calc output from a single xarray.Dataset to a structured xarray.DataTree, with per-IM datasets and shared station/source metadata.
  • Add empirical IM calculation via oq_wrapper for supported IMs/tectonic types and store results under {im}/empirical/{model}.
  • Add vs30 to broadband waveform outputs and introduce realisation/schema/default support for empirical configuration.

Reviewed changes

Copilot reviewed 6 out of 7 changed files in this pull request and generated 4 comments.

Show a summary per file
File Description
workflow/scripts/im_calc.py Major refactor to DataTree output; adds empirical GMM evaluation and metadata/unit annotation.
workflow/scripts/bb_sim.py Adds vs30 coordinate to broadband netCDF output for downstream site parameter usage.
workflow/schemas.py Adds EMPIRICAL_PARAMETERS schema for validating empirical config in realisations/defaults.
workflow/realisations.py Adds moment-weighted averaging helpers and EmpiricalParameters realisation configuration.
workflow/default_parameters/root/defaults.yaml Introduces default empirical configuration (active_shallow + NSHM2022).
uv.lock Lockfile change related to dependency resolution (cffi/pycparser marker).
pyproject.toml Adds netCDF4 dependency to enforce safe import order vs OpenQuake/HDF5 stack.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Comment thread workflow/scripts/im_calc.py
Comment thread workflow/scripts/im_calc.py
Comment thread workflow/scripts/im_calc.py Outdated
Comment thread workflow/scripts/im_calc.py
Comment thread workflow/scripts/im_calc.py
The site parameters, with basin depths estimated using the Chiou
and Youngs (2008) relations.
"""
z1pt0 = chiou_young_08_calc_z1p0(vs30) # ty: ignore[invalid-argument-type]

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we want to estimate for all of the sites? Wondering if we have any "real sites" in this where we have a better measured / estimate value of z1.0 / 2.5 etc to use those instead from the site database.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah this is a stand-in at the moment, I intend to replace this with actual site database values at a later date.

@lispandfound
lispandfound changed the base branch from pegasus to nzvm/bb-matched-filter August 31, 2026 09:09

@github-actions github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please change the following dependencies for consistency.

Comment thread pyproject.toml Outdated
@@ -14,12 +14,13 @@ dependencies = [
"im-calculation @ git+https://github.com/ucgmsim/IM_calculation@no_parallel",

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
"im-calculation @ git+https://github.com/ucgmsim/IM_calculation@no_parallel",
"im-calculation @ git+https://github.com/ucgmsim/IM_calculation.git@no_parallel",

@github-actions github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please change the following dependencies for consistency.

Comment thread pyproject.toml Outdated
@@ -14,12 +14,13 @@ dependencies = [
"im-calculation @ git+https://github.com/ucgmsim/IM_calculation@no_parallel",

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
"im-calculation @ git+https://github.com/ucgmsim/IM_calculation@no_parallel",
"im-calculation @ git+https://github.com/ucgmsim/IM_calculation.git@no_parallel",

@lispandfound
lispandfound force-pushed the nzvm/bb-matched-filter branch from 427f137 to 85f21f1 Compare September 2, 2026 02:23
@lispandfound
lispandfound force-pushed the nzvm/bb-matched-filter branch from 85f21f1 to dc8667f Compare September 2, 2026 02:52
@lispandfound
lispandfound force-pushed the empirical_ims branch 2 times, most recently from 6e72aee to ed8468f Compare September 4, 2026 01:44
@lispandfound
lispandfound force-pushed the nzvm/bb-matched-filter branch from 09437b1 to 310ebe9 Compare September 8, 2026 22:56
@lispandfound
lispandfound force-pushed the empirical_ims branch 2 times, most recently from 786f419 to d4a7ec7 Compare September 9, 2026 00:09
@lispandfound
lispandfound force-pushed the nzvm/bb-matched-filter branch 2 times, most recently from 3b51442 to f410e8c Compare September 9, 2026 05:04
lispandfound and others added 10 commits September 23, 2026 12:43
Introduces the realisation sections an SW4 run is configured by, the
geometry derived from them, and the first shipped set of values. They
are one change because the geometry exists only to read the config
classes, and the defaults are the shipped instance of both -- the tests
that pin the geometry assert against those numbers, not synthetic ones.

**Configuration.** `Refinements` describes the stack of vertical grids
SW4 solves on, in the abstract -- the layers a domain of any depth would
be given. `refinements_for_depth` resolves that against a particular
domain: it truncates at the domain floor, extends with
`unbounded_refinement_resolution` if the listed layers do not reach that
far, and guarantees the last layer is at least two cells deep so a
domain ending just past a refinement boundary does not degenerate.

`SW4Parameters` holds SW4's input file as a list of `SW4Command`s rather
than a fixed set of typed fields. SW4 has a large and growing command
vocabulary and we do not want a field per command; anything SW4 accepts
can be written as `{"name": ..., "parameters": {...}}` without touching
this code. `render()` handles the two places SW4 differs from Python's
`str()`: booleans go out as 0/1, and `None` parameters are omitted
rather than rendered as "None".

**Geometry.** `workflow.sw4` is the single definition of the supergrid
sponge, which three later stages need and must not each re-derive.
`minimum_fault_buffer_m` is `sponge + 5h`, additive rather than a
multiple: the two terms have different origins -- the sponge is where
the equation changes, the five grid points are the source stencil plus
the dissipation operator applied to its outermost point -- and a
multiplicative margin would collapse as the grid refines while the
stencil still spans five points. `supergrid_width` follows SW4's own
precedence: `width=` beats `gp=`, and with neither, SW4's built-in 30
gridpoint default.

`absorbed_period` is advisory and nothing enforces it. Clearing the
sponge is not the same as the sponge working: the layer absorbs
adiabatically only while `W cos(theta) / lambda >> 0.431`, so a 12 km
sponge absorbs below ~8 s at normal incidence and ~4 s at 60 degrees.
That is a property of the period band being asked for, not of the
buffer.

**Values.** The load-bearing one is `fault_buffer: 14.0 km`, against the
root default of 2.0 km. Two criteria bear on it and only the first is
enforced (by a later commit):

  HARD FLOOR, 14 km. A source must sit outside the sponge plus the
  stencil margin: (30 + 5) * 400 m, at the widest sponge any supported
  domain produces. Below this the run is not a ground motion prediction
  -- it is the standing-wave failure seen in validation_results_24-08.

  LONG-PERIOD CRITERION, ~80 km. Sound runs in that campaign sat 83-113
  km from the domain edge, broken ones 4-5 km. Deliberately NOT
  enforced: enforcing it would reject every affordable domain.

The 2.0 km root default stays. It is the EMOD3D value, kept for
CyberShake reproducibility; v24.2.2.x have no absorbing layer to clear.
The sponge is stated as `width: 12000.0` rather than `gp: 30`, because
`gp` is measured on the coarsest grid and would vary with domain depth.

`SKIP_PAIRS` records the two genuinely version-specific sections:
`refinements` describes a grid the EMOD3D versions do not have, and
`resolution` is a single uniform spacing an SW4 run does not have.

No pipeline stage reads any of this yet.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
`create-sw4-input` renders the SW4 input file for a realisation: the
grid, the refinement stack, the source, the stations, and whatever
commands the `sw4` section carries.

The requested domain becomes the grid's *interior*. The grid is padded
laterally by one supergrid width on each side, so the sponge sits
outside the region anyone asked to simulate rather than eating into it.

`adjust_for_topography` handles the interaction between topography and
the refinement stack. Topography lifts the top of the grid above z=0,
which thins the first layer; SW4 requires every grid in the stack to
hold at least `nz_min` cells, so any layer left too thin is pushed down
until it does. The topography surface participates in that count as a
boundary but is not itself returned as a refinement -- the caller gets
back the real layers plus the depth used.

Commands come from the realisation rather than being hard-coded, so a
run can add or override any SW4 command without a code change. Output
commands that take a time window are merged with the simulation
duration on the way out.

`test_bottom_refinement_holds_the_sponge` asserts against the shipped
v26.7.1Hz configuration rather than a synthetic one: the invariant that
matters is that the bottom refinement of a *real* domain is thick enough
to hold the bottom sponge, at every depth we support.

Requires nzcvm (for `nzcvm.formats.sfile`, used to read the model's
lateral footprint), which in turn forces `requires-python >= 3.13`.
nzcvm is currently a git dependency and this cannot merge until it is
released.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
`create-nzvm-input` writes the config `nzcvm generate` consumes. One
realisation can produce several distinct velocity models, along two
independent axes:

  --format sw4 | emod3d       the grid the model is sampled onto: a
                              mesh-refined SW4 grid written as an sfile,
                              or a uniform EMOD3D binary grid
  --layers full | tomography  the layer stack queried to fill it: the
                              whole thing (basins, coastline, offshore,
                              Ely GTL taper) or just the background
                              tomography and the numerical clamps, which
                              is the reference no-bells-and-whistles model

Splitting those two axes is the point: comparing a full model against a
tomography-only one on the *same* grid is how a basin's contribution is
isolated, and it should not require hand-editing a config.

For SW4 the model is padded by more than the grid is -- one supergrid
width, as `create-sw4-input` uses, plus a few grid points of slack --
so SW4 never queries outside the sfile.
`test_model_padding_contains_the_padded_grid` asserts that ordering
directly, because if it ever inverts the symptom is a solver reading off
the end of the model rather than an error.

The `nzcvm` realisation section carries the layer stack as nzcvm's own
`LayerConfig` objects and the DEM surface as a `Path`, which is what the
`path_serialiser` added at the bottom of this stack exists to write.

The container image installs nzcvm explicitly. It is a bare requirement
resolved through `[tool.uv.sources]`, which pip does not read, so
`pip install workflow` alone cannot find it. The other git-sourced
dependencies are PEP 508 direct references in `dependencies` and pip
resolves those itself.

Blocked on an nzcvm release: this makes nzcvm a hard import of
`workflow.schemas`, so every workflow install pulls it in, not just the
SW4 path.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
nzcvm is a hard import of workflow.schemas but was never declared, so
deptry and pytest failed and `pip install ucgmsim-workflow` could not
find it. The container papered over this by installing nzcvm and
workflow from git, which undid #150's PyPI-based container build.

nzcvm 2026.9.1 is now on PyPI and requires Python 3.13, so declare it,
raise requires-python to match (the container already runs 3.13), and
restore the Dockerfile that installs ucgmsim-workflow==$WORKFLOW_VERSION.
Also declare pyproj, which nzvm_input_template imports directly.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
`lf-to-xarray --format sw4` reads SW4's HDF5 station file into the same
dataset shape the EMOD3D path produces, so everything downstream is
solver-agnostic. Stations are read in batches sized to a target chunk,
lazily via `map_blocks`, because a real run's station file does not fit
in memory.

Two things worth stating explicitly:

The datasets carry SW4's displacement-mode names (EW/NS/UP), but for an
SRF rupture source the time function SW4 receives is the slip *rate*, so
the nominal displacement output is physically velocity. Components are
reordered to EMOD3D's convention (x east, y north, z down).

SW4 reports, per station, how far into the supergrid sponge that station
sits. That number decides whether a trace is a ground motion prediction
at all, so it is carried all the way to the intensity measure file as a
station-dimension *coordinate*. That is load-bearing rather than
incidental: coordinates ride through `bb-sim` and `im-calc` untouched,
whereas data variables are dropped when chunks are recombined.

A station that reported nothing gets NaN, never 0.0. Zero means "checked
and found in the interior" -- a real claim about the run -- and an old
station file, or a solver with no absorbing layer, is not entitled to
make it. Stored as float32 for the same reason: readers open these files
with `mask_and_scale=False`, so an integer sentinel would read back raw
and become a plausible penetration depth.

`im-calc` attaches the coordinate for every solver, all-NaN where
nothing was reported, so the column always exists and its absence is
never mistaken for a clean run.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
A fault buffer too small to clear the supergrid absorbing layer puts the
source inside the region where SW4 solves a damped, coordinate-stretched
equation rather than the wave equation. The run completes and produces
waveforms; they are just not ground motion predictions. That failure has
already cost a campaign once (validation_results_24-08), and it is
cheapest to catch before a domain is written rather than after a
simulation has run.

`generate-domain` now checks the buffer against the sponge, using the
same `workflow.sw4` geometry `create-sw4-input` pads with, so the two
cannot disagree.

The depth is computed before the lateral domain so the refinements can
be resolved against it -- the sponge width depends on the coarsest grid,
which depends on how deep the domain goes. This is a reordering only;
`domain_max_depth` does not depend on the lateral extent.

Both `sw4` and `refinements` are optional and their absence is not an
error: the v24.2.2.x defaults are EMOD3D-only and have neither section,
so those runs skip the check entirely and are unaffected.

This is a behaviour change for SW4 realisations: a configuration that
previously produced a domain may now be rejected. That is the point.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Replaces the `srf2stoch` binary with an in-process conversion, and drops
it from the container build.

The conversion is a box average expressed as a sparse fractional-overlap
matrix: row j of the kernel holds the overlap weights between coarse bin
j and the fine cells it spans. The special case where the two grids span
exactly the same length is what srf2stoch.c implements; the sparse form
generalises it without materialising the empty overlaps, and handles the
padded case where they do not.

That padded case is why this was worth porting. The HF code demands a
single dx/dy across every SRF segment, and `nx = ceil(len / dx)` means
the coarse grid is generally *longer* than the plane it covers. The
overhang is split evenly between the two ends, so the grids stay
concentric -- the stoch format records only a centre point and an
`nx * dx` extent, so off-centre padding would sit the slip distribution
in the wrong place on the plane HF reconstructs. Edge bins then carry
weights summing to the covered fraction rather than to 1, which is what
conserves total slip*area rather than cell value.

Quantities that are not spread over a cell are handled separately.
Rupture time is a time, so it is divided by the covered fraction; every
cell is partially covered because nx and ny round up, so this never
divides by zero. Rise time is slip-averaged, and there the coverage
factor cancels, so it must *not* be divided out again.

Rake is averaged as a circular mean weighted by slip -- the arithmetic
mean of -179 and 179 degrees is 0, which is the opposite of the right
answer. A plane with no slip anywhere has no slip-weighted mean, so it
falls back to an unweighted average.

This changes every stoch file, and therefore every high-frequency
result. The previous binary's outputs were not moment-conserving for
segments whose length was not an exact multiple of the stoch dx.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Two changes to the same function body, which is why they are one commit.

**Larger than memory.** `bb-sim` read both waveform files entirely into
memory and looped over the three components in Python. Now LF and HF are
opened lazily, the common stations are selected on the *backend* arrays
before chunking -- selecting after chunking makes station reordering an
all-to-all dask shuffle, which materialises the whole array -- and the
recombination runs per chunk under `map_blocks`. Chunking is over
stations only, so each chunk holds complete traces for resampling,
alignment and filtering.

`resample_signal` and `align_datasets` replace the old pad-and-align:
the two legs no longer have to share a timestep, so an SW4 LF run at one
dt can be combined with HF at another. `relabel_hf_components` maps HF's
090/000/ver onto LF's x/y/z once, up front.

**Site amplification.** Replaces `qcore.siteamp_models.cb_amp_multi`
with the `site-calculation` models, selected by `bb.site_amp_version`:
CB2014 as before, or BA2018. The amplification is now constrained to an
explicit [fmin, fmax] band with logarithmic tapers at both ends, which
is what the two new `fhightop` and `fmax` parameters are for -- the old
code had a lowpass taper only.

`site_amp_version` becomes an enum rather than a free string. It was
previously declared, defaulted to "2014", and read by nothing at all.

The shared broadband parameters move from each defaults version into
root; the versions now carry only `flo`, which is the one value that
genuinely differs between them.

Station-dimension coordinates -- `supergrid_depth` among them -- ride
through `map_blocks` untouched, which is why the supergrid penetration
is carried as a coordinate rather than a data variable.

This changes broadband results: different site amplification model, and
a highpass taper where there was none.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Both solvers low-pass their *source* time functions: SW4 through its
`prefilter` command, EMOD3D through `bfilt` with a corner derived from
`min_vs` and the grid spacing. The wave equation is linear, so filtering
the source is equivalent to filtering every trace -- which means the LF
leg arrives at `bb-sim` already low-passed, and the matched pair then
low-passes it a second time.

The result is a hole around the merge frequency: the matched high-pass
and low-pass are power-complementary by construction, but only if each
leg is filtered once. `tests/test_bb_filters.py` measures the hole and
then measures that each correction closes it.

Because the filter is applied to the source rather than the output,
`bb-sim` can divide it back out here instead of the solver having to be
re-run. `--solver` reads the filter the realisation actually specifies
(rather than assuming one), and refuses to guess: an SW4 prefilter that
is not a lowpass, or an EMOD3D configuration that also high-passes, is
an error rather than a correction for the wrong thing.

`--filter` chooses which leg absorbs it:

  lf    restores the low-frequency leg exactly. Needs the largest boost.
  hf    fills the missing power from the high-frequency side instead.
        Cannot restore below the merge frequency, where the HF leg
        carries nothing to scale up.
  both  one factor on both legs. Best conditioned -- `current_power` is
        bounded below by the high-pass leg -- and restores the power sum
        exactly, at the cost of touching both legs.

`MAX_BOOST` caps any correction at 10x. The two real configurations need
about 2.4x; anything near the cap means the source filter rolls off
faster than the target and the leg is being reconstructed from content
that is not there. `warn_if_ill_conditioned` says so up front.

Without `--solver` the recombination is bit-for-bit unchanged: the
correction is applied after the matched pair rather than in place of it,
and the gains are exactly 1. Verified against the previous commit's
output.

What was applied is written into the broadband file's attributes.
Without that the correction is invisible downstream and two files that
differ by it look identical.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
`im-calc` now also evaluates the ground motion models named in the
realisation, writing them alongside the simulated measures at
`{im}/empirical/{model}` in the same data tree. Having both in one file
is the point: comparing a simulation against a GMM should not require
joining two files on station and period.

A multi-fault realisation has to be reduced to the single rupture a GMM
describes -- one magnitude, rake, dip, ztor, zbot -- which
`calculate_source_parameters` does by moment-weighted averaging, with
rakes and dips averaged as unit vectors rather than as angles. Site
parameters come from vs30 via the Chiou & Young z1pt0/z2pt5
estimators, and the distances are the ones already computed for the
simulated measures, so both sides see identical geometry.

Combinations a model does not support are skipped with a warning rather
than failing the run -- the NSHM2022 logic tree configures pSA but not
PGA or PGV, and that is not an error.

Two deliberate awkwardnesses:

`netCDF4` is imported before anything that pulls in OpenQuake. OpenQuake
brings h5py, which is linked against a different HDF5 build than netCDF4
is, and whichever loads second cannot open our waveform files. The
import order is load bearing, not stylistic.

`EmpiricalParameters` types its fields as `Any` and its schema validates
plain strings rather than `oq_wrapper.constants` enum members. Importing
that module pulls in OpenQuake, which is expensive and must be
precompiled; the strings are resolved to enum members inside the
calculation instead, so `workflow.realisations` stays cheap to import.

vs30 is loaded eagerly rather than left as a dask array -- it is one
float per station, and letting it stay lazy propagates dask into every
station coordinate.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@lispandfound
lispandfound force-pushed the nzvm/bb-matched-filter branch 2 times, most recently from 813ab2b to 81450af Compare September 23, 2026 01:35
@lispandfound

Copy link
Copy Markdown
Contributor Author

Folded into #139 (the empirical commit now sits at the tip of nzvm/im-calc-metadata); #145 is retargeted onto nzvm/bb-matched-filter.

@lispandfound
lispandfound removed this pull request from stack #146 September 23, 2026 01:38
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants