ADD: Adding changes for original organization repo to this repo - #1
Merged
Conversation
Turn the repository into an installable Python package so the advective
interpolation and spectral (FFT) gridding work has somewhere to land. This
commit is packaging only -- the subpackages are documented but empty, so the
API surface can be settled before ~2000 lines of research code moves in.
Build and layout:
- pyproject.toml with a setuptools + setuptools-scm backend. Versions come
from git tags; fallback_version keeps builds working from a tag-less clone
or an sdist with no git metadata.
- src/ layout, so tests run against the installed package rather than the
working directory.
- Subpackages radar_palette.{advection,gridding,util,testing}, each carrying
the scope, conventions and status of the capability that will fill it.
py.typed is shipped.
- Optional-dependency extras keep the core install to numpy/scipy/arm_pyart:
[advection] pulls scikit-image (optical flow), [spectral] pulls finufft
(non-uniform FFT). Plus [all], [test], [docs], [dev].
Quality gates:
- pytest config with slow/network markers; packaging tests assert every name
a subpackage advertises in __all__ is actually importable.
- tests/conftest.py provides requires_skimage / requires_finufft skip markers
so a minimal install still collects a green suite.
- ruff for lint, import sort, pyupgrade, numpy rules and numpydoc docstring
style; pre-commit runs it alongside the standard hygiene hooks.
- Sphinx docs stub (autosummary + numpydoc + myst-nb). Builds with -W;
RADAR_PALETTE_DOCS_OFFLINE=1 drops intersphinx for offline builds.
The GitHub Actions workflow is deliberately not in this commit: pushing files
under .github/workflows/ requires a separate token permission, so CI arrives
as a follow-up commit rather than blocking the packaging work.
Documentation records the two conventions that are easy to get wrong and
expensive to rediscover: reflectivity is interpolated in dBZ rather than
linear Z, and displacement is the physical echo displacement from the first
volume to the second, so velocity = displacement / dt with no sign flip.
Verified locally: editable install imports, 10 tests pass, ruff clean, sdist
and wheel build and pass twine check --strict, docs build warning-free.
Co-authored-by: Scott Collis <scollis.acrf@gmail.com>
Co-authored-by: Claude <noreply@anthropic.com>
Scaffold the package: setuptools build, src layout and docs
Claude created CI
An advection-interpolated volume represents an instant that was never observed,
but the interpolation built its output by deep-copying the earlier bracketing
volume -- so the reconstructed volume inherited that volume's clock. A volume
interpolated halfway across a 7-minute gap reported the time of the scan 3.5
minutes before it.
The failure is silent. The volume is structurally valid, plots without
complaint, and survives a cfradial round-trip, so nothing surfaces the error at
the point it is introduced. It corrupts any time-ordered downstream use --
rainfall accumulation, cell tracking, matching against another instrument --
by up to the full inter-volume interval, and the cause is hard to trace back
from the symptom.
radar_palette.advection.timing:
- volume_reference_time reduces a volume to one representative instant (the
mean ray time, the centroid of the acquisition window, so it does not drift
with scan strategy the way a first- or last-ray convention would).
- interpolate_ray_times builds the per-ray time dictionary of the interpolated
volume.
- apply_interpolated_time stamps it onto the output and refreshes the derived
time metadata.
Two details the implementation is careful about:
Intra-volume structure is preserved. A volume is not an instant; its rays are
acquired in sequence over minutes. That per-ray pattern is a property of the
scan strategy, not of the interpolation, so the ray pattern is carried through
and the volume as a whole is shifted, rather than collapsing every ray to a
single stamp.
Arithmetic is done on absolute times decoded from each volume's own CF units.
Two volumes read from separate files are not guaranteed to share an epoch, and
differencing raw offsets across mismatched epochs would silently produce a
wrong interval.
Also adds radar_palette.testing.{make_empty_ppi_volume,assign_scan_times},
which build bare PPI geometries with explicitly controlled acquisition times so
timing is tested independently of any field data -- a failing test localises to
the time handling rather than to optical flow or sampling.
Test driven: the 35 unit tests were written first and run red (the module did
not exist), then against the implementation. One test failed for an instructive
reason -- it compared raw offsets, which are meaningless across a re-based
units epoch -- and is now split into a comparison on absolute times plus an
explicit guard that the units epoch itself moves, which is the hazard a
consumer reading time["data"] while ignoring time["units"] would hit.
38 tests in the new timing module: parameterised alpha sweeps, endpoint recovery
at alpha 0 and 1, linearity, extrapolation beyond the bracket, monotonicity,
scan-pattern preservation, mismatched epochs, masked ray times, input
immutability, array aliasing, CF key completeness, error paths, metadata refresh
without fabrication, and three cfradial round-trip integration tests -- because
an in-memory fix that does not survive serialisation would not be a fix. The
repository suite is 48 tests, the other 10 being the pre-existing packaging
checks.
Verified: 48 passed (38 timing + 10 packaging), ruff check and format clean,
docs build warning-free,
sdist and wheel pass twine check --strict, and advection.timing imports with
scikit-image absent (it does not depend on the optional extra).
Co-authored-by: Scott Collis <scollis.acrf@gmail.com>
Co-authored-by: Claude <noreply@anthropic.com>
Users arrive with one of two object families describing the same radar data:
Py-ART's Radar/Grid, or xradar's DataTree and xarray Dataset. Requiring callers
to convert before and after every operator is friction that also invites
mistakes, so the entry points now accept either family and return the family the
caller expects.
Default policy, per the request:
- advection: the output family mirrors the input. Radar in gives Radar out,
DataTree in gives DataTree out.
- gridding: Py-ART input gives a pyart.core.Grid, xradar input gives an xarray
Dataset.
Either default is overridden per call with output_flavor ('pyart', 'xradar',
'xarray'; case-insensitive, and 'xradar' is accepted as an alias for 'xarray'
when naming a grid flavour, since a caller working in that ecosystem names it
that way).
Conversions delegate to Py-ART's own interoperability layer rather than
reimplementing the mapping between the two data models: DataTree.pyart.to_radar()
for the Py-ART surface over a tree, Grid.to_xarray() and pyart.xradar.Xgrid for
grids, and xradar's cfradial readers for the reverse direction. That layer is
maintained alongside the formats and handles metadata this package has no
business duplicating.
Two sharp edges were found by measuring the installed libraries rather than
assuming, and both are absorbed here with regression tests:
Assigning Xradar.time does NOT write through to the underlying DataTree, even
though add_field does. Correcting an interpolated volume's time through the
Py-ART wrapper alone would leave an xradar caller's tree on its original clock --
the same silent staleness the timing module exists to prevent. Time is now
written back explicitly, per sweep, by write_ray_times_to_datatree.
pyart.xradar.Xgrid rejects a dataset whose time has been decoded to datetimes,
but Grid.to_xarray produces exactly that, so a direct Xgrid(dataset) raises. The
conversion re-encodes time to CF offsets before wrapping.
A third was caught by a test rather than by inspection: xradar registers an
accessor named 'xradar' on DataTree, so a hasattr(obj, 'xradar') probe intended
to unwrap an Xradar wrapper is also true for a bare tree and yields the accessor
instead of the data. The check is now an explicit isinstance, with a comment
recording why.
Test driven: 71 tests across the new io layer and entry points were written
first and run red (the modules did not exist), then against the implementation.
Coverage includes flavour detection for every input type, the default-mirroring
policy in both directions, explicit overrides, string aliases and case handling,
round-trip geometry/field/reference-time preservation, per-sweep time write-back,
intra-sweep ray-pattern preservation, mixed-flavour bracketing volumes, wrong-
family and unknown-flavour error paths, and an equivalence check that flavour
changes the container but not the numbers.
The repository suite is 121 tests (10 packaging, 38 timing, 50 io flavours,
21 entry points, 2 new subpackage-export cases).
radar_palette.gridding.grid_volume is currently backed by
pyart.map.grid_from_radars. The spectral operator is not yet ported; settling the
flavour contract first means the port replaces the backing without changing the
public signature or its tests.
Verified: 121 passed, ruff check and format clean, docs build warning-free, sdist
and wheel pass twine check --strict, and all entry points import with
scikit-image and finufft absent.
Co-authored-by: Scott Collis <scollis.acrf@gmail.com>
Co-authored-by: Claude <noreply@anthropic.com>
Advection timing
… floor to 3.11 Three related fixes to the packaging of the merged dual-flavour work, squashed because the second and third exist only to correct the first. 1. Declare xarray, xradar and netCDF4 as direct dependencies. radar_palette.io imports xradar and xarray directly, and both it and advection.timing import netCDF4, but none was declared -- all three arrived transitively through arm_pyart. That works wherever arm_pyart was installed first and breaks for a user whose resolver chose differently, failing inside a conversion rather than at install. Two packaging tests guard it: one walks the AST of every module under src/radar_palette and asserts each non-stdlib import is declared; the other asserts xarray.DataTree exists rather than trusting the version string. The guard was checked by deliberately removing a declaration and confirming it fails, and it was the guard rather than inspection that found the netCDF4 gap. 2. Raise requires-python to >=3.11. CI then failed on py3.10: pip cannot satisfy xradar>=0.12.0 there, because that release is the first to require 3.11. Reproduced locally in a real 3.10 environment with the identical error, so the floor added in (1) was a genuine regression. Keeping 3.10 was viable and was measured, not guessed: a 3.10 environment resolves arm_pyart 2.2.0 / xradar 0.11.1 / xarray 2025.6.1, that stack has everything radar_palette.io needs (xarray.DataTree, Xradar, Xgrid, the cfradial readers), and 122 of 123 tests passed on it. It was rejected because Py-ART made the same move -- arm_pyart 2.2.1 dropped 3.10 -- so supporting it means testing against a stack two minor versions behind what users resolve. The 3.10 run exposed a second defect in (1), independent of the resolver: the new import-guard test imports tomllib, stdlib only from 3.11, so it raised ModuleNotFoundError rather than skipping. Now an explicit pytest.importorskip. Raising ruff's target-version to py311 unlocked UP017; 22 datetime.timezone.utc usages are now datetime.UTC. The matching CI matrix change (dropping "3.10") is NOT in this commit: GitHub refuses a personal access token pushing .github/workflows without the classic workflow scope. It must be applied separately, and the CHANGELOG says so. 3. Document the version fallback as observed rather than assumed. The README claimed a tag-less clone reports 0.1.dev<N>+g<sha>. That was never observed -- every build here reports 0.0.0, and setuptools-scm raises LookupError in this environment even with a tag present, tested by creating a throwaway tag. The README now documents only what happens: the build takes fallback_version, and a 0.0.0 wheel is not a release. Verified on a real Python 3.11 environment (xarray 2026.7.0, xradar 0.12.0, arm_pyart 2.2.4): 123 passed. On 3.12: 123 passed, ruff check and format clean, docs build warning-free, wheel metadata reports Requires-Python: >=3.11 with 3.11/3.12/3.13 classifiers, twine check --strict passes. Co-authored-by: Scott Collis <scollis.acrf@gmail.com> Co-authored-by: Claude <noreply@anthropic.com>
Declare xarray, xradar and netCDF4 as direct dependencies
First step of the spectral gridding port (see SPECTRAL_PORT_API_PROPOSAL.md). This
is self-contained and useful on its own: the census answers "is this volume
spectrally griddable, and if not why" without gridding anything.
radar_palette.gridding.geometry
-------------------------------
antenna_to_cartesian_43 reproduces pyart.core.antenna_to_cartesian exactly, with
one interface difference: metres in and out rather than kilometres in.
cartesian_to_antenna_43 exists because pyart.core.cartesian_to_antenna is NOT the
inverse of Py-ART's own forward transform. The forward returns great-circle arc
length on the 4/3 earth; the Py-ART inverse computes a flat straight-line
sqrt(x^2+y^2+z^2) with no curvature term, so any product inverting a Cartesian grid
back to antenna coordinates inherits a silent systematic bias.
This is the one place the project's "use Py-ART's conversion layers" rule does not
apply, because the layer is wrong for the purpose.
The discrepancy was re-measured here rather than copied from the earlier findings
report, and the measurement corrected a misreading of it. The reported figure of
-286 m at 118 km is a maximum OVER ELEVATION, not a value at typical scan
elevations: the error peaks near 35 degrees at every range, and at 1 degree the
same range shows only -19 m. Measured against our forward transform:
range err at 1 deg worst err elevation err
10 km -0.1 m -2.3 m +34 mdeg
118 km -19 m -313 m +398 mdeg
250 km -109 m -1392 m +843 mdeg
460 km -496 m -4649 m +1552 mdeg
Azimuth agrees to ~1e-14 degrees; only range and elevation are affected. Our
inverse closes the round trip to ~1e-10 m.
Two Py-ART interface traps are pinned by tests, both found while writing them
rather than by inspection: cartesian_to_antenna assigns into its azimuth result, so
it raises TypeError on scalar input and must be handed arrays; and it returns range
in METRES while the forward transform takes KILOMETRES, so a naive round trip is
wrong by 1000x in the opposite direction to the usual mistake.
radar_palette.gridding.census
-----------------------------
census_sweep and census_radar measure each sweep and classify it as
EXACT_UNIFORM_PERIODIC (full 360 on a closing lattice, so a plain FFT is valid),
UNIFORM_PARTIAL_SECTOR (uniform but not periodic) or NON_UNIFORM. census_radar also
tags split-cut groups: NEXRAD VCPs repeat the low tilts, and de-duplicating them is
mandatory before vertical interpolation or a tilt is interpolated against itself.
Two deliberate departures from the research code:
- Returns a list of dataclasses, not a pandas DataFrame. pandas is not a declared
dependency of this package and adding one for a return type is not warranted;
as_row() feeds a DataFrame in one line for anyone who wants a table. This settles
question 3 of the port proposal.
- Accepts either object flavour by routing through radar_palette.io, so the census
works on an xradar DataTree as well as a Py-ART Radar.
SweepClass is an enum.StrEnum, which ruff's UP042 correctly prefers now the floor
is 3.11, and which keeps as_row() CSV-serialisable without enum handling.
Test driven: 90 collected cases (81 test functions; the round-trip closure test is
parametrized over 5 ranges x 5 elevations, and the Py-ART range-error bound over 3
ranges) written first and run red. They cover analytic properties of both
transforms, agreement with Py-ART's forward transform, the documented inverse
discrepancy including its elevation dependence, and classifier decisions placed
deliberately on both sides of every tolerance boundary -- jitter just inside and
just outside AZ_DEV_TOL_FRAC, exact and jittered lattices, partial sectors, wrap
seams, multiple revolutions, non-monotonic azimuth, and duplicate fixed angles at
and beyond SPLIT_CUT_TOL_DEG.
Repository suite is 213 (14 packaging, 38 timing, 50 io flavours, 21 entry points,
90 geometry).
Verified on Python 3.11 (xarray 2026.7.0, xradar 0.12.0, arm_pyart 2.2.4) and 3.12:
213 passed on both. ruff check and format clean, docs build warning-free, twine
check --strict passes.
Co-authored-by: Scott Collis <scollis.acrf@gmail.com>
Co-authored-by: Claude <noreply@anthropic.com>
Port 4/3-earth geometry and sweep geometry census
Step 2 of the spectral gridding port. SweepSpectralEvaluator turns a sampled sweep
into a continuous band-limited function, with three azimuth paths selected by the
sweep class the census assigns, and range mirror-extended rather than treated as
periodic.
Two defaults differ from the research code, both because measurement said so.
1. Conjugate-gradient refinement is ON by default (n_cg=12).
Convolution gridding computes the ADJOINT of the sampling operator, not its inverse.
That is not an academic distinction. Measured on a constant field sampled with
+/-0.3-spacing azimuth jitter, unrefined gridding reproduces the constant to only 60%
relative deviation, and the error scales with jitter: 3.6% at +/-0.02, 9.1% at
+/-0.05, 19% at +/-0.1, 39% at +/-0.2.
CG on the density-weighted normal equations converges fast, and the curve was
measured rather than assumed:
n_cg 0 2 4 6 8 10 12 16 20
err 60.3% 32.6% 6.06% 1.84% 0.79% 0.16% 0.07% 0.05% 0.05%
It plateaus at 12, costing ~50% more build time (11 ms -> 17 ms on a 120x34 sweep).
A default returning 60% error on a constant field is not defensible, so 12 is the
default; n_cg=0 recovers the unrefined operator and is still tested.
This only works because adjoint() is the exact transpose of forward(). An earlier
version of the operator carried a stray factor of M or N0 and its adjoint test
returned ~1 instead of ~1e-15 -- invisible to any test that merely checked the output
looked smooth. The test is now asserted at < 1e-10 and recorded in the report.
2. The linear-Z guard tests MAGNITUDE, not dynamic range.
Interpolation must run in dBZ. Reproduced here on a 55 dBZ core against a -10 dBZ
background -- a hard edge, which is the geometry that provokes it -- interpolating in
linear Z put 42.6% of evaluated samples below zero, worst excursion -0.16x the field
maximum, every one of which clips on conversion and becomes indistinguishable from
real weak echo.
The first guard I wrote used max/min dynamic range, and it was wrong in both
directions: a dBZ field spanning 0 to 50 has a ratio of 1e5 and would be falsely
flagged, while linear Z from a 6.8-31.8 dBZ scene has a ratio of only 316 and would
be missed. Magnitude separates them cleanly -- no meteorological dBZ value approaches
200, while linear Z passes 200 above about 23 dBZ. Known blind spot, documented: a
linear field peaking below ~23 dBZ is not caught.
Also documented rather than glossed: dBZ does NOT remove Gibbs ringing, and the
ringing is NOT symmetric. On that same hard edge (360 rays, 100 gates, 600x600 probe)
the dBZ interpolant reaches 74.8 dBZ against a true maximum of 55.0 and -20.6 against
a true minimum of -10.0 -- that is +19.8 dB above the data on the high side but only
-10.6 dB below it on the low side, so the high-side excursion is roughly twice the
low-side one. The magnitudes also depend on sampling: at 180 rays by 60 gates the same
edge gives +17.5 and -8.7 dB, so every quoted figure now states its configuration. A
regression test asserts the asymmetry, because an earlier draft of these docstrings
summarised both directions as "about 20 dB", overstating the low side by ~2x.
What dBZ buys is that an overshooting dBZ value is still interpretable and clippable,
where a negative reflectivity factor is not. The guarantee is narrow and the
docstrings now say so; band_frac trades resolution against overshoot, which is tested.
A smooth field shows ZERO negative samples under linear-Z interpolation, so the
sharp-edge fixture is load-bearing: an earlier version of that check used a smooth
field and would have suggested the dBZ requirement was unnecessary.
Test driven: 87 collected cases written first and run red. They cover exact recovery
at and off the lattice for the exact path, periodicity across the branch cut,
linearity, ray-order independence, mirror-extension behaviour on a range ramp, fast-
vs-direct-path convergence, both sector gap-fill modes, the adjoint identity, CG
convergence, density compensation, the guard and its message, all three fill methods
preserving measured samples, and the band-limited lattice round trip.
Repository suite is 300 (14 packaging, 38 timing, 50 io flavours, 21 entry points,
90 geometry, 87 evaluator).
Verified on Python 3.11.15 (xarray 2026.7.0, xradar 0.12.0, arm_pyart 2.2.4, scipy
1.17.1) and 3.12.13 (arm_pyart 2.2.5, scipy 1.18.0): 300 passed on both. ruff check
and format clean, docs build warning-free, twine check --strict passes, and the
gridding API imports with scikit-image and finufft absent.
Co-authored-by: Scott Collis <scollis.acrf@gmail.com>
Co-authored-by: Claude <noreply@anthropic.com>
Port the sweep spectral evaluator, NUFFT operators and dBZ guard
Step 1 of the advection port, and the piece that makes the capability real. Until
now radar_palette.advection could only CORRECT the clock on a volume something else
had interpolated; it could not interpolate. This adds the two operators that do.
radar_palette.advection.flow
----------------------------
grid_optical_flow estimates a dense TV-L1 motion field between two gridded volumes,
independently per vertical level. Per-level rather than per-volume because storm
motion is height dependent, and dense rather than rigid because in an MCS the broad
stratiform region and the embedded convective cells move at different speeds -- a
single displacement is dominated by the larger, slower area and can read near zero
while cells move at 20 m/s.
radar_palette.advection.morph
-----------------------------
advection_interpolate reconstructs a volume at a time between two observed volumes.
The morph runs PER GATE in native (azimuth, slant range, elevation) coordinates, so
the result is an ordinary volume, not a Cartesian grid; only the motion estimation
passes through a grid, because optical flow needs a regular raster.
THE SIGN IS ESTABLISHED BY MEASUREMENT, NOT INSPECTION.
This is the one thing in the module that is easy to get backwards, and getting it
backwards produces output that looks entirely plausible while being worse than a
naive time average. So it was pinned end-to-end before the implementation was
written: warp the earlier volume by the FULL displacement and check it reproduces
the later volume.
warp direction RMSE vs the later volume
position - 1.0 * D 0.46 dBZ
position + 1.0 * D 6.09 dBZ
A factor of 13. The implementation follows that measurement. The flow sign was
verified the same way, on a synthetic scene translating due east at a known
20 m/s over 300 s: recovered +5966 m against a true +6000 m, with the cross-track
component at -10 m against a true 0.
Accuracy, measured against an ANALYTIC truth rather than another approximation. The
test scene is a Gaussian blob in rigid translation, so the field at any fractional
time is known exactly:
alpha=0.5 advected 0.047 dBZ RMSE
alpha=0.5 naive average 0.171 dBZ RMSE (73% worse)
alpha=0.0 vs earlier 0.015 dBZ
alpha=1.0 vs later 0.016 dBZ
That 73% is the EASY case and the docstrings say so. A rigid translation is the most
favourable possible scene; real convection grows, decays and changes shape between
volumes, and none of that is motion an advection scheme can capture. On real storm
scenes the improvement over a time average is a few percent, not seventy. Overstating
this would be the most tempting error in the module, so the honest figure is in the
module docstring rather than only in a commit message.
Timing is wired in at the source. The output is built by deepcopy of the earlier
volume, so apply_interpolated_time is applied before returning -- otherwise the
reconstruction would report having been observed up to a full volume interval before
the instant it depicts. Verified exact at every alpha on a 420 s bracket: 0/105/210/
315/420 s offsets for alpha 0/0.25/0.5/0.75/1.0, with intra-volume ray spacing
preserved and the reference epoch rebased.
Two interface facts found by testing, both now carrying comments so they are not
re-broken:
- to_pyart_radar returns (radar, flavour), so the separate detect_radar_flavor
call it replaced was redundant.
- pyart.map.grid_from_radars calls len() on its first argument, which the Xradar
wrapper does not support as a bare object -- it must be passed as a 1-tuple.
Without this, xradar input fails while Py-ART input works.
Also fixes a latent flaw in the packaging guard. test_third_party_imports_are_
declared_dependencies rejected any correctly-gated OPTIONAL import, so the new
scikit-image import failed it. Accepting extras is the right fix, but on its own it
would let a contributor promote an optional dependency to required unnoticed, so a
companion test asserts scikit-image and finufft stay out of project.dependencies.
Both guards were verified non-vacuous by deleting a declaration and confirming each
fires with the offending name.
Test driven: 28 cases written first and run red (ImportError on the absent
operators), covering flow sign and magnitude, argument-order reversal, return order,
shape and mismatch validation, warp direction end-to-end, endpoint recovery at both
ends, beating the naive baseline, output geometry and value bounds, field renaming
and metadata, target_time-to-alpha conversion, and both object flavours.
Repository suite is 329 (15 packaging, 38 timing, 28 advection interpolation, 50 io
flavours, 21 entry points, 90 geometry, 87 evaluator).
Verified on Python 3.11.15 (xarray 2026.7.0, xradar 0.12.0, arm_pyart 2.2.4, scipy
1.17.1, scikit-image 0.26.0) and 3.12.13: 329 passed on both. ruff check and format
clean, docs build warning-free, twine check --strict passes, and the advection API
imports with scikit-image absent.
Co-authored-by: Scott Collis <scollis.acrf@gmail.com>
Co-authored-by: Claude <noreply@anthropic.com>
Port the advection interpolation operator
Step 3 of the spectral gridding port. The evaluator turns one sweep into a
continuous surface; this turns a stack of surfaces into a 3-D field, which is the
last component missing before grid_volume can be wired to the spectral path.
radar_palette.gridding.cones
----------------------------
A sweep at fixed elevation does not sample a plane. On a curved earth the beam climbs
faster than a straight line, so the sample locus is a CONE, and every consequence in
these modules follows from that. cone_range_height is the exact inverse of
antenna_to_cartesian_43 at fixed elevation; target_lattice builds the shared
horizontal lattice each tilt is gridded onto, carrying that tilt's own height surface
and validity mask; build_cones grids every sweep onto it.
Two details that are load-bearing rather than cosmetic, both now commented:
- Lattice axis order is row=y, column=x, matching pyart's (nz, ny, nx). Transposing
it does not announce itself on a roughly isotropic field -- it produces an error
equal to the field's own standard deviation, not a visible flip.
- The validity mask excludes the near-radar hole because the range axis is
mirror-extended for the FFT and the mirror is EVEN about the first gate, so it
reflects the echo into the hole and fabricates a plausible near-radar return.
dedup_sweeps collapses split cuts by MEASURED range rather than by cut naming, so it
needs no knowledge of any particular radar's scan-strategy vocabulary: within a
split-cut group it keeps the member reaching furthest, ties falling back to the lower
sweep index, and records the group size so dropped sweeps are visible.
radar_palette.gridding.vertical
-------------------------------
interp_column_stack interpolates the cone stack onto target heights, per column,
because cone curvature means both the tilt heights and the gaps between them vary
from column to column.
Extrapolation policy is NONE, and the invariant is exact: every finite value carries
VerticalFlag.INTERPOLATED, and every non-interpolated cell is NaN. A test asserts the
two sets are equal rather than merely consistent -- a NaN wearing an interpolated flag
is a silent hole, and a finite value without one is a silent extrapolation.
INTERIOR_GAP handles the case that is easy to get wrong. A column's valid tilt run is
usually contiguous, since at fixed arc length slant range grows with elevation. But
maximum valid range is a property of the ECHO, not the geometry, and can be
non-monotonic in elevation -- a higher tilt reaching further than the one below it --
leaving an annulus where the lower tilt is absent and the higher present. Such targets
are flagged rather than interpolated across two non-adjacent surfaces.
MEASURED: SAMPLING DOMINATES SCHEME CHOICE.
The default is pchip_z, but the more useful result is how little that matters. On a
synthetic bright band at 2500 m with a 400 m standard deviation -- a feature sharper
than the tilt spacing -- reconstructed onto five height levels:
scheme 6 tilts 23 tilts
linear_z 4.086 dB 0.295 dB
linear_elev 4.085 dB --
pchip_z 4.024 dB 0.202 dB
nearest_tilt 5.194 dB --
Refining the scan strategy (median inter-cone thickness 680 m -> 213 m) improves the
result by a factor of 14. Choosing the best scheme over the worst improves it by 1.2x.
Within one strategy the error tracks the local layer thickness directly: median
absolute error 0.42 dB where adjacent cones are 400-700 m apart, 3.81 dB where they
are 1200-2500 m apart.
So the honest reading of a gridded value is alongside its gap_m, and no scheme
recovers structure the sampling did not capture -- asserted as a test rather than left
as advice. pchip_z is the default because it measured best in both configurations and,
being monotone, cannot invent an extremum between two cones.
A field LINEAR in height cannot discriminate the schemes at all: all three smooth
schemes reproduce it to 0.001 dB. The bright-band fixture is therefore load-bearing,
and there is a test asserting that a linear field shows no difference -- otherwise a
future comparison built on the easy field would suggest the choice is free.
Test driven: 69 collected cases written first and run red (ImportError on the absent
module). Coverage includes cone-surface inversion against the forward transform,
curvature (a planar surface is explicitly rejected), the 0-degree rise, lattice axis
order and both validity cuts, split-cut selection and tie-breaking, cone ordering and
per-tilt diagnostics, the flag/value equivalence, no extrapolation below the lowest or
above the highest cone, interior gaps versus above-coverage, layer thickness growing
with range, monotone non-overshoot, nearest-tilt returning only observed values, and
the sampling-versus-scheme comparison above.
Repository suite is 398 (15 packaging, 38 timing, 28 advection interpolation, 50 io
flavours, 21 entry points, 90 geometry, 87 evaluator, 69 vertical).
Verified on Python 3.11.15 (xradar 0.12.0, arm_pyart 2.2.4, scipy 1.17.1) and
3.12.13: 398 passed on both. ruff check and format clean, docs build warning-free,
twine check --strict passes, and the vertical API imports with scikit-image and
finufft absent.
Co-authored-by: Scott Collis <scollis.acrf@gmail.com>
Co-authored-by: Claude <noreply@anthropic.com>
Port cone geometry and vertical assembly
Final step of the spectral gridding port. Every component was upstream but nothing
connected them: gridder.py mentioned grid_from_radars three times and
SweepSpectralEvaluator zero times. grid_volume now takes method={"pyart","spectral"},
and the FFT gridding capability works end to end from either object family.
The spectral path runs the census, evaluates each sweep as a band-limited interpolant
on its own cone, resamples the cone stack onto the caller's requested axes, and
interpolates vertically per column.
WHY THE DEFAULT STAYS "pyart"
The two operators are different, not ranked, so switching the default would change
every existing caller's numbers in exchange for a trade-off they did not ask for.
Measured on a hard azimuthal echo edge (a 10 -> 50 dBZ wedge):
method output range overshoot
spectral [ 4.67, 55.52] +5.52 / -5.33 dB
pyart [10.00, 50.00] +0.00 / -0.00 dB
Distance weighting is an average, so it cannot exceed its inputs. The spectral
interpolant is exact where the sampling supports it but rings at a discontinuity.
On a smooth field the two agree to 0.35 dB median where both produce data.
The gridder docstring originally quoted the per-sweep evaluator's +19.8 / -10.6 dB
excursion, which is a different measurement on a different scene; it now states this
path's own figure and cross-references the evaluator's, noting the smaller number
reflects vertical interpolation and horizontal resampling averaging some of it away
rather than a gentler operator.
Two tests pin this rather than leaving it as a docstring claim: one asserts the
default reproduces a direct pyart.map.grid_from_radars call byte for byte, and one
asserts the default is the non-ringing operator -- so flipping the constant fails
with a reason. Reversing the recommendation is a one-line change.
HONEST COVERAGE
The spectral path emits a coverage_flag field carrying VerticalFlag per cell. A radar
volume does not observe a box, and an unexplained NaN leaves a user unable to tell a
data gap from a geometric coverage limit; the flag distinguishes below-lowest-tilt,
above-highest-tilt, hole-in-the-tilt-run and too-few-tilts. A test asserts the
finite-value set equals the INTERPOLATED-flag set end to end, and that the flag
survives conversion to xarray.
REFUSALS RATHER THAN PLAUSIBLE OUTPUT
Vertical interpolation needs two cones to bracket a target, so a single-tilt volume
raises rather than returning an all-NaN grid. Merging several volumes spectrally
would need a cross-radar combination rule this operator does not define, so that
raises NotImplementedError pointing at method="pyart" -- inventing a merge rule here
would be a research decision disguised as plumbing.
Test driven: 22 cases written first and run red (ImportError on the absent
constant), then 9 more added after measuring the edge behaviour and the flavour
composition. Coverage includes method dispatch and rejection, default-equals-pyart
byte equality, spectral shape and axis fidelity, the coverage flag and its
invariant, both refusals, that the two backings genuinely differ (so a mis-wired
dispatch running Py-ART for both cannot pass) while broadly agreeing on smooth
fields, all four flavour x method combinations, and that both object families give
the same spectral answer.
Also corrects the README, which still said the modules were "documented but not yet
implemented", and adds a usage section. The example was wrong on first writing --
grid_optical_flow requires a field argument and advection_interpolate takes volumes
plus alpha rather than a precomputed flow -- and is now executed verbatim as written.
Repository suite is 429 (15 packaging, 38 timing, 28 advection interpolation, 50 io
flavours, 21 entry points, 90 geometry, 87 evaluator, 69 vertical, 31 gridder).
Verified on Python 3.11.15 (xradar 0.12.0, arm_pyart 2.2.4, scipy 1.17.1) and
3.12.13: 429 passed on both. ruff check and format clean, docs build warning-free,
twine check --strict passes, and the gridding API imports with scikit-image and
finufft absent.
Co-authored-by: Scott Collis <scollis.acrf@gmail.com>
Co-authored-by: Claude <noreply@anthropic.com>
Wire the spectral operator behind grid_volume
Found by running the advection operator on real data for the first time. interpolate_ray_times raised ValueError on three consecutive ARM C-SAPR2 volumes because they have 13955, 13951 and 13954 rays: a sweep occasionally records one ray more or fewer, so consecutive volumes of the SAME 15-sweep strategy do not repeat their ray count. The guard therefore rejected every real pair, and advection_interpolate was unusable on anything but synthetic input. The arithmetic never needed the restriction. Nothing in the function pairs rays between the two volumes: the output carries radar_early's acquisition pattern shifted by a single scalar derived from the two volumes' reference times, so radar_late contributes exactly one number and there is no correspondence to violate. The guard encoded an assumption the implementation did not rely on. Removed it, documented in the Notes section why ray counts are deliberately not required to match, and cited the measured counts so the next reader does not reinstate it. The pre-existing test asserting the guard is replaced rather than deleted, and renamed to record that the requirement was superseded rather than forgotten. Four new tests in TestUnequalRayCounts cover the case real data presents: a one-ray difference, output length always following the earlier volume across several later-volume lengths, a shift that varies only through the reference times and not the counts, and apply_interpolated_time not reimposing the restriction at the public boundary. Verified end to end afterwards: advection_interpolate now runs on the real volumes and reconstructs the withheld middle time to 0.000 s. Scored against that withheld volume on a common grid, echo above 15 dBZ inside 105 km (n = 369), RMSE is 6.23 dBZ against 7.86 for a naive time average and 13.00 for persistence. 433 passed on 3.12, ruff check and format clean, docs build warning-free. Co-authored-by: Scott Collis <scollis.acrf@gmail.com> Co-authored-by: Claude <noreply@anthropic.com>
Second bug found by running the advection operator on a real sequence. Sampling raised on the 13-volume ARM C-SAPR2 run: ValueError: The points in dimension 0 must be strictly ascending or descending _sample_sweep tiles the azimuth axis at -360, 0 and +360 degrees so a query that lands on the 0/360 seam interpolates between real neighbouring rays. That tiled axis goes to RegularGridInterpolator, which requires it to be STRICTLY ascending. _sweep_sampler sorted the azimuths but did not deduplicate them, so a sweep containing the same azimuth twice produced equal consecutive entries and the whole reconstruction raised. Real antennas do record a repeated azimuth: measured across these 13 volumes, 3 of 195 sweeps contain exactly one duplicate (17:07:36 sweep 13, 17:23:36 sweeps 1 and 4), which is what a brief dwell looks like in the data. A synthetic volume built from arange never does, which is why the synthetic suite could not reach this. Note the azimuths are otherwise well behaved -- zero sweeps are non-monotonic once the 360 wrap is accounted for -- so deduplication, not resorting, is the fix. Duplicates are now collapsed to one ray by averaging the rays that share an angle, rather than discarding one: they are independent samples of the same beam position, so their mean is the better estimate of it. The average is a nanmean per group, because a duplicated ray can be masked in one copy and not the other, and an all-nan group must stay nan rather than collapsing to zero. Four tests in TestRepeatedAzimuths, written first and confirmed red on the exact ValueError above: one duplicate does not raise; the reconstruction still carries usable echo rather than merely not raising; a ten-ray dwell also works; and strictly-ascending input is bit-identical, so the well-behaved path is untouched. The class carries @requires_skimage, like every other class in this file that reaches the optical-flow path. Omitting it was my error and CI caught it: the four tests failed on the py3.12-minimal job with ImportError, because advection_interpolate needs scikit-image and that job installs without the [advection] extra deliberately. Verified on a purpose-built minimal environment (arm_pyart and xradar present, scikit-image and finufft absent) rather than by reasoning about the marker: 405 passed and 32 skipped there, with these four among the skips and a reason attached. On the full install the same four RUN and pass, so the gate does not hide them from the job that can actually exercise them. 437 passed on 3.12, ruff check and format clean, docs build warning-free. Co-authored-by: Scott Collis <scollis.acrf@gmail.com> Co-authored-by: Claude <noreply@anthropic.com>
Adds carry_fields to advection_interpolate, so a single optical-flow estimate can be applied to several variables in one call. Motivation is physical, not ergonomic. The operator estimated flow from the field it was interpolating, so reconstructing four variables meant four calls -- and three of them would have derived motion from a variable that is not a tracer. Differential reflectivity is a shape measure, correlation coefficient is closer to a quality flag, and Doppler velocity is signed and folded at the Nyquist interval; tracking any of them tracks the wrong thing. Motion is a property of the scene, so it should be estimated once from reflectivity and applied to everything. The warp loop was already field-agnostic -- displacement, resampling coordinates and blend weights depend only on geometry and alpha -- so the change is to loop the sampler over the requested fields inside the existing per-sweep loop, rather than to add a second code path. Four fields therefore cost one flow estimate and one gridding pass, not four. Each carried field keeps its own metadata dict, so units are not overwritten with the tracer's: a carried ZDR field stays dB rather than becoming dBZ. interp_field_name renames only the tracer, since renaming several outputs from one string is not well defined. Documented caveat for folded velocity: the blend is a weighted mean, so averaging across the Nyquist interval returns a value near zero where the truth wraps. Measured on the ARM C-SAPR2 pair used to exercise this (17:55:37 and 18:03:37 on 8 Aug 2026, Nyquist 16.30 m/s), gates changing by more than one Nyquist interval between scans are 7.14% within echo above 15 dBZ and 20.38% across all valid gates, the remainder being noise whose velocity is essentially uniform. The docstring says to dealias first if the interpolated velocity is to be read quantitatively. Seven tests in TestCarryFields, written first and confirmed red. Beyond presence and error handling they pin the properties that a plausible-looking wiring mistake would break: the carried field is warped rather than copied from the earlier volume (the fixture blob translates, so a copy sits at the start position); it keeps its own units; adding it does not perturb the tracer's own reconstruction; and because the fixture makes it an exact affine function of the tracer, the output must satisfy that same relation -- which holds only if motion, resampling and blend weights are identical. Default behaviour is unchanged: no carry_fields gives one field out. The class carries @requires_skimage, verified on a real minimal environment rather than assumed: 405 passed and 39 skipped with scikit-image absent, and 444 passed on the full install where the seven run. Co-authored-by: Scott Collis <scollis.acrf@gmail.com> Co-authored-by: Claude <noreply@anthropic.com>
The three preceding commits changed public behaviour but left CHANGELOG.md untouched, which breaks the convention every prior feature commit in this repository follows (6225994, 20c3779, 2dcffbf all update it in the same commit as the change). Added: carry_fields on advection_interpolate, with the reason it exists -- the polarimetric variables are not tracers, so the flow must come from reflectivity and be applied to them -- and the folded-velocity caveat with its measured figures. Fixed: the equal-ray-count guard and the repeated-azimuth failure, each with the measurement that motivated it and the note that both were found by running the operator on a real volume sequence, unreachable from synthetic volumes built on arange. No source change. 444 passed on 3.12, ruff clean, docs build warning-free with the changelog included. Co-authored-by: Scott Collis <scollis.acrf@gmail.com> Co-authored-by: Claude <noreply@anthropic.com>
Third real-data failure in the advection operator, and the one a user hits first.
volume_reference_time returns a UTC-aware instant, but pyart.util.datetime_from_radar
returns a naive cftime.DatetimeGregorian -- so the most natural call a Py-ART user
can make,
target = pyart.util.datetime_from_radar(earlier) + half_the_gap
advection_interpolate(earlier, later, target_time=target)
raised TypeError: can't subtract offset-naive and offset-aware datetimes. Only alpha
was usable on real volumes.
Every existing target_time test built its target FROM volume_reference_time, so all
of them were aware-on-aware and none could reach this. That is why the gap survived
the port.
A naive target is now read as UTC. CF times are UTC by construction, and this is the
rule _decode_ray_times already applies to decoded ray times, so the two agree rather
than each inventing a convention.
The normalisation rebuilds the instant field by field instead of calling
replace(tzinfo=UTC). A first version did use replace(), passed a naive-datetime test,
and still raised on real Py-ART input: cftime.DatetimeGregorian is not a
datetime.datetime subclass and its replace() rejects tzinfo outright. The cftime case
is therefore pinned as its own test rather than assumed to follow from the naive one.
Three tests: a naive datetime target matches an aware one; they agree off the
midpoint too, where a sign or offset slip would not cancel; and a genuine
cftime.DatetimeGregorian target matches the equivalent alpha. Verified on real
volumes -- the call above now returns a volume whose reference time equals the
requested target exactly.
447 passed on 3.12, ruff clean, docs build warning-free.
Co-authored-by: Scott Collis <scollis.acrf@gmail.com>
Co-authored-by: Claude <noreply@anthropic.com>
Ray count guard
The gridding module docstring justified the opt-in spectral path by contrasting
two costs: the spectral operator rings at a discontinuity, while Cressman
weighting "cannot overshoot, and it degrades gracefully where sampling is poor".
The first half is right. The second is wrong, and it flatters the default.
A radius of influence is a fixed length. The vertical separation between
adjacent tilts grows with range. Beyond the range where the separation exceeds
the radius, no gate is within reach of a cell and it is left empty -- an
interior hole, with valid data both above and below it in the same column. That
is not graceful degradation; it is a hole where the operator could not reach.
Measured on a C-SAPR2 volume (15 tilts, 1.5 to 42 degrees) at 500 m over a
60 km box, against an observed maximum of 68.0 dBZ:
roi_func filled interior holes peak
dist_beam (pyart default) 70.7% 570 63.1 dBZ
dist, z_factor=0.1 95.1% 6 57.0 dBZ
dist, z_factor=0.2 96.3% 0 53.1 dBZ
constant_roi=4000.0 98.3% 0 43.0 dBZ
The trade-off runs in both directions, which is why the docstring now states it
as a trade-off rather than as graceful behaviour: closing the gaps entirely
costs 15 dB of peak. The spectral path reaches zero interior holes at a cost of
4.7 dB, because vertical assembly interpolates between whichever cones bracket
the target height instead of reaching a fixed distance.
Worth being explicit that this is not a mis-set argument. Py-ART's own defaults
for roi_func="dist_beam" produce the identical result -- 570 holes either way --
so the honest statement is about the operator, and a caller who wants filled
columns has to choose a point on the trade-off deliberately.
grid_volume's **gridding_kwargs entry now names roi_func as the consequential
argument for the "pyart" backing and points at the table.
Five tests in TestDistanceWeightingLeavesInterConeGaps, written before the
docstring change and verified non-vacuous:
- interior holes appear with the beam-width radius
- the Py-ART defaults reproduce them exactly, so it is not the arguments
- a wider radius closes them
- closing them lowers the peak (both directions asserted)
- the spectral path has zero interior holes on the same volume
Non-vacuity probes: a closely-spaced-tilt volume (0.5 to 2.0 degrees) yields 0
interior holes against 228 for the widely-spaced fixture, so the fixture's tilt
spacing is load-bearing rather than the assertion being trivially true; and the
spectral section is 439 of 984 cells finite, so its zero hole count is a real
measurement rather than an empty array.
452 passed on Python 3.11.15 and 3.12.13, ruff check and format clean, docs
build warning-free with the new list-table, twine check --strict passes.
Co-authored-by: Scott Collis <scollis.acrf@gmail.com>
Co-authored-by: Claude <noreply@anthropic.com>
Correct the claim that distance weighting degrades gracefully
nufft.py hand-rolls a Kaiser-Bessel kernel, its Fourier transform by quadrature, and a CG solve on the density-weighted normal equations. It is correct and is left untouched as the reference. This factors the same operator into two independent axes so faster machinery can be used without changing the mathematics: an *engine* supplies forward/adjoint, a *solver* turns those products into a lattice. Engines: reference, scipy (default), dense, finufft, ducc0, torch. Availability is checked at construction, never at import, so a minimal install keeps working and an absent engine raises a message naming the extra rather than an ImportError. The default engine needs no new dependency. It is the reference algorithm with two exact substitutions: the spreading stencil becomes one CSR matrix (np.add.at is an unbuffered scatter and was 38% of solve time) and the spread-to-spectrum transform becomes rfft (the spread field is real, so the reference computed a conjugate-symmetric half it discarded). Agreement with nufft.py is 1.5e-15 -- round-off, not approximation -- for 1.9-2.3x on sweeps of 360 rays and up. The dense engine has no interpolation kernel at all: it forms the DFT matrix. Affordable here and only here, because on the azimuth axis that matrix is n_rays x n_modes and both are ray counts (8 MB at 720 rays). So it is exact to round-off, needs no dependency, and is the fastest engine below ~720 rays; finufft overtakes it past ~1440, the O(n log n) vs O(n^2) crossover. solver='direct' factorises the normal equations instead of iterating on them, reaching the least-squares solution CG approaches (checked against CG run to 1e-14). Worth 10-30x, against ~2x for the fastest engine change -- the solver is the larger of the two effects, because n_lattice here is a ray count and a Cholesky at that size is microseconds. The accuracy gain belongs to the engine/solver PAIR, not the solver. On a field exactly band-limited on the lattice, so truth is known in closed form: direct on dense or finufft recovers it to ~1e-10 against ~3e-5 for 12 CG iterations, but direct on a Kaiser-Bessel engine converges to that engine's own ~3e-4 kernel error and gains nothing at low jitter. Tests assert each half; a reader who took the headline figure for a property of solver='direct' alone would be wrong. What direct improves on every engine is degraded sampling, where the conditioning fails rather than the kernel: at +/-0.45-spacing jitter, 12 iterations reach 1.7e-2 against 4.5e-4 on the same transforms. Its ridge is likewise load-bearing -- a sector's normal matrix is singular by construction (a 30-120 deg sector at 1 deg spacing fits 91 rays onto a 366-point lattice) and a bare Cholesky fails on it. Regularised, the direct solve fits the measured rays to 8.9e-11 where 12 iterations reach 3.7e-4. Note the engines are not all round-off equivalent to nufft.py: dense, finufft, ducc0 and torch each shift the answer by ~5e-4, replacing the reference's Kaiser-Bessel kernel (own error 2.5e-4) with a different or exact transform. Towards exact arithmetic, but a shift -- which is why the default is the engine that reproduces the reference bit-for-bit. torchkbnufft's phase convention cost a 69% error that its adjoint test passed cleanly, a consistently wrong pair of operators being still a consistent pair. Hence every engine is tested against nufft.py and, where it claims to beat that kernel, against exact trig-polynomial arithmetic computed in closed form. Test suite mutation-checked: 7 of 8 seeded faults fail a test, and the 8th (dropping the Nyquist split) turned out to be a no-op for real fields, which is now recorded as a test with the measurement behind it. benchmarks/bench_nufft_engines.py regenerates every figure quoted; timings are from an arm64 laptop and hardware-dependent.
The docstring said 2-4x; the benchmark spans 1.54x (1440x1832) to 4.12x (360x1832), so both ends were outside the stated bounds. The narrowing at the largest shape is the informative part -- it is where the transform starts to outweigh the per-call setup -- so it is now stated rather than rounded away.
The docstrings and CHANGELOG said dense is fastest below ~720 rays with finufft overtaking past ~1440. The saved benchmark says otherwise: finufft leads at BOTH ends of the tested range (22.5x vs 18.8x at 120x34, 32.7x vs 21.5x at 1440x1832) and dense leads only the middle band (35.3x at 360x500, 44.0x at 360x1832, 32.2x at 720x1200). So the ranking is not monotone in ray count and the single-threshold framing was wrong. Two different mechanisms bound the dense engine at the two ends -- O(n_rays * n_modes) losing to O(n log n) at the top, and the per-solve BLAS work failing to amortise the built matrix at the bottom -- and the full table is now in the module docstring so the reader picks by measurement instead of extrapolating from one number. Operational 1 deg and 0.5 deg volumes sit in the band where dense wins, which is the useful part of the finding and survives the correction. Also corrects the mutation-testing record in the previous commit message, which misnamed the surviving fault. Two rounds were run. In the final round of 8 (torch omega sign, scipy rfft conj, remove ridge, dense sign flip, no symmetrise, scipy skip deapod, wrong lattice scale, finufft isign flip) the survivor was "no symmetrise" -- and investigating it showed the symmetrise line is defensive rather than load-bearing: the transforms leave only ~3e-16 relative asymmetry and cho_factor reads one triangle anyway, which is now stated in the comment there rather than the overstatement it replaced. The Nyquist-split survivor came from the earlier round; it is a genuine no-op for real fields (adjoint identical to the bit) and is recorded as a test with the measurement behind it.
Everything so far was synthetic: sums of azimuthal harmonics, band-limited by construction, errors relative to np.ptp. Validating by held-out ray prediction on three storm-filled sweeps (SPOL S-band convective line, SWX C-band widespread precipitation, CSAPR2 C-band convective, all 0.5 deg and all inspected as PPIs first) changes two conclusions. The ridge default of 1e-10 was wrong by orders of magnitude, and not conservatively wrong -- it is the WORST value tested on all three sweeps. Held-out median error, best ridge vs 1e-10: 6.70 vs 7.80 dB (SPOL), 4.79 vs 5.08 (SWX), 8.52 vs 10.19 (CSAPR2). So ridge now defaults to 'auto', which reads the eigenvalue spectrum and caps the condition number, resolving to the 1e-10 floor on a full-rank synthetic geometry and ~1.6e-2 on a real sweep. No constant serves both. The driver is NOT rank deficiency, which was my first hypothesis and is wrong: two of the three real sweeps are full-rank and still want a large ridge. It is that real reflectivity is not band-limited. Speckle, clutter and echo edges put energy in modes the sweep cannot determine, a small ridge fits that energy at the measured rays, and the interpolant between them is noise. A synthetic field built from a few harmonics has nothing in those modes and shows none of this, which is exactly why synthetic tuning gave an answer eight orders too small. A parametrised test adds broadband noise to a synthetic field and watches the optimum leave the floor, which isolates band-limitedness as the variable rather than rank or conditioning. Second correction, and the more important one: the direct solver's headline accuracy gain does not transfer. With a well-chosen ridge it is COMPARABLE to 12 CG iterations on real data (6.70 vs 7.80, 4.79 vs 4.92, 8.52 vs 8.47) -- it wins two and loses one, by tenths of a dB. The ~1e-10 recovery stands but is a statement about representing a field this operator can represent exactly, and real reflectivity is not such a field. The module docstring now carries a warning to that effect, and the honest remaining case for the direct solver on real data is speed (10-30x) and having no iteration count to tune, not accuracy. Rank deficiency is still real, still geometric, and now reported: normal_null_dim counts the undetermined modes, since a caller cannot distinguish an under-determined answer from a converged one by looking at the values. Operational azimuth sampling leaves gaps of ~2x nominal, giving 3 undetermined modes of 414 on SPOL and 23 of 993 on a GUC X-band sweep. Adds benchmarks/bench_nufft_engines.py --real to reproduce all of it, and five tests covering the ridge contract. Three existing tests asserted the old small-ridge behaviour and are updated rather than deleted; one of them now says explicitly that data fit is the right yardstick only because its field is band-limited, since on real data a better data fit is the failure mode.
The DEFAULT_RIDGE comment collapsed both wins to '~0.15 dB' three lines under a table giving SPOL as 6.70 against CG's 7.80. Only SWX is ~0.15 (0.13); SPOL is 1.10. The conclusion is unchanged -- comparable, not orders better, with the direct solver's real case being speed rather than accuracy -- but understating the one case where it does clearly win is no better than overstating it.
Redone against the ARM BNF C-SAPR2 case from the performance notebook: a 15-tilt volume, all 15 cones NON_UNIFORM, ~800 rays each, so five tilts sharing a scan strategy and differing only in what they observe. Three findings, and two of them are corrections. Stratification is not optional. Over half the held-out gates in a real sweep sit below 0 dBZ -- noise and the gap-fill floor -- and an unstratified median is dominated by them. It also ranks the solvers differently: on BNF sweep 11 a ridge of 1e-1 is the BEST value on all gates (4.14 dB) and the WORST on gates above 10 dBZ (8.09 dB). Every dB figure in this work is now reported both ways, and the echo column is the one to read. 'auto' does nothing on a well-conditioned sweep, which is most of them. At 4.5-42 deg these tilts have condition numbers of 15-32, so the condition-number term never binds and the rule returns its floor -- yet the held-out optimum on echo gates is 1e-3..1e-2 at every tilt. Conditioning is a correlate of out-of-band energy, not that quantity, and on a clean geometry the correlation fails in the unsafe direction. Two fixes were tried and both reverted, with the reasoning kept at MIN_RIDGE because the trade is the kind a later change would silently break. Raising the floor to 1e-2 wins on this volume (5.330 dB against 5.446, mean of per-tilt echo medians, versus 5.272 for a per-sweep oracle and 5.389 for CG) but costs eight orders of magnitude on a band-limited field -- a bad trade for a default, since it hard-codes an assumption about the data into a geometry-only rule. Generalised cross-validation recovers 1e-10 exactly on the synthetic case and then selects 1e-10 on four BNF tilts whose optimum is 1e-2, because it scores on training residual and overfitting is what minimises that. Any rule fitted on the measured rays has the same defect. So the honest position is that no automatic rule tested is safe as a silent default: pass ridge=1e-2 for real reflectivity, and use --real to measure it. Eight real sweeps now, graded on echo gates with a per-sweep oracle ridge: direct beats 12 CG iterations on 6 of 8 by a mean of 0.13 dB. A wash, consistent with the earlier single-sweep result. Speed is unchanged and remains the reason to adopt this. On BNF sweep 0 (798 rays x 1050 gates) against nufft.py: scipy+cg 1.6x, dense+direct 27x, finufft+direct 32x. Adds evaluate_lattice(), which evaluates the recovered interpolant at arbitrary azimuths -- forward() only evaluates at the measured rays, which cannot distinguish fitting from overfitting. It is a module function rather than an engine method because the interpolant belongs to the lattice. This also exposed that the reference engine lacked the mode_low attribute every other engine sets; it alone failed on that path. --real gains the BNF volume via RADAR_PALETTE_BNF_VOLUME (the file is not downloadable, so it is skipped when unset).
The engine layer stopped at SweepSpectralEvaluator, so grid_volume() -- the entry point every caller and notebook actually uses -- could not reach it. build_cones and _grid_spectrally now forward the three arguments, left unset by default so the evaluator's own defaults remain the single source of truth for an unconfigured call. This is the consequential place for them. Spectral gridding builds one evaluator per sweep and that build dominates, so the total is nearly flat in output resolution: on this volume the published ramp reads ~41.5 s from 2000 m down to 250 m, a 3000x change in output cells. The cost is set by the azimuth solve, which until now was not selectable from the entry point. Also corrects the GCV note in CHANGELOG.md and nufft_engines.py: it returns the floor on four of five BNF tilts, but those tilts' measured optima are 1e-2 on three and 1e-3 on one, not 1e-2 on four. The conclusion is unchanged -- a training-residual criterion cannot see overfitting -- but the two counts are different numbers and were conflated.
build_cones ran its 15 tilts serially while the machine had 14 cores and BLAS alone was using about 5 of them. The tilts are independent -- each builds its own evaluator, writes its own cone, shares nothing -- so this is the one place in the spectral path where hardware helps without changing the algorithm. Threads rather than processes: the per-tilt cost is NumPy/SciPy linear algebra and FFTs, which release the GIL, and a process pool would have to pickle the radar object to every worker (hundreds of megabytes for a research volume, and 727 MB for the one this was measured on). The cap is memory, not cores, and that distinction is the substance of resolve_tilt_workers. One in-flight tilt holds an upsampled lattice that reaches ~1 GB on a 171 km domain at 250 m, so a core-count default would try to hold fifteen of those and swap -- slower than serial, and hard to diagnose. The resolver sizes the pool from free memory and the per-tilt lattice, and warns when it reduces a request rather than silently running fewer workers than asked, which would look like the parallelism not working. On this machine n_jobs=-1 gives 14 workers at 2 km and 8 at 250 m. Default stays serial (n_jobs=None), so the memory profile of previous releases is unchanged unless a caller opts in. The loop body moved into a closure to make it schedulable, which is a large diff for no behaviour change; the 550 pre-existing tests passing unaltered is the evidence for that. Six new tests cover what could genuinely go wrong: parallel and serial grids must be BIT-identical (any difference at all would be a race rather than a tolerable reordering, so a tolerance here would hide the bug it is meant to catch), the cone stack must stay in ascending-elevation order regardless of completion order, the memory cap must bind and warn on a fine lattice, and n_jobs=0 must be rejected rather than clamped to serial.
The n_jobs work in d156439 was correct but incomplete, and the gap was bad enough to take a workstation down. ThreadPoolExecutor started 14 workers, and each worker's BLAS and FFT calls started their own pool sized to the whole machine: ~200 threads on 14 cores, load average above 500, and the processes outlived SIGTERM because native thread pools do not respond to it. Found by watching the load during a benchmark rather than by any test -- the bit-identical test passed throughout, because the output was never wrong, only the machine. resolve_tilt_workers now returns (workers, blas_threads) and guarantees their product does not exceed the core count. build_cones applies the thread limit with threadpoolctl inside the pool. The memory cap it had before was the right instinct aimed at the wrong resource; both bind, and the function now says so in its name and its docstring. Measured, 15 tilts over a 60 km domain at 1 km, dense engine, direct solver, 14 cores: n_jobs=1 1 worker x 14 threads 17.47 s 1.00x n_jobs=2 2 workers x 7 threads 10.70 s 1.63x n_jobs=4 4 workers x 3 threads 8.44 s 2.07x n_jobs=8 8 workers x 1 thread 5.65 s 3.09x Worth recording that the fastest setting gives each worker a SINGLE thread. The azimuth solve is a few hundred rows, too small for a wide GEMM to pay for its own synchronisation, so spreading whole tilts across cores beats splitting each tilt's linear algebra. Speedup stays short of linear because tilts share memory bandwidth and the vertical interpolation afterwards is serial. threadpoolctl is optional (extra: parallel). The serial default never touches the native pools so it must not become a hard requirement; without it the parallel path is still correct, just over-subscribed, and warns loudly -- a silent fallback would look exactly like "parallelism does not help on this machine", which is the failure it exists to prevent. Tests assert the property that was missing rather than the symptom: workers * threads <= cores across the plausible request range, since either factor alone looked fine and it was the product that was wrong.
Asked which settings produce nonsensical grids, and the answer is worse than expected: several, and one of them needs no new arguments at all. Measured on the BNF C-SAPR2 volume, 468k cells, spectral path: default, n_cg=12 -61 .. 66 dBZ 0 cells |Z|>80 n_cg=25 -1907 .. 2471 dBZ 1294 n_cg=200 -1907 .. 2471 dBZ 1294 az_ridge=1e-10 -1903 .. 2466 dBZ 1292 az_ridge=1e-6 -98 .. 138 dBZ 26 az_ridge=1e-2 -61 .. 58 dBZ 0 az_ridge='auto' -61 .. 59 dBZ 0 n_cg is the trap. It predates this branch, its docstring invites tuning, and raising it looks like asking for a better answer -- but the iteration count is acting as regularisation rather than as a convergence budget, so past about a dozen iterations the solve amplifies out-of-band content instead of converging. Twelve to twenty-five iterations is a cliff, not a dial. The default is not exempt either. Per-cone, one tilt of fifteen overshoots by 107 dB and reaches 169 dBZ at n_cg=12; it is invisible in the final grid only because the vertical interpolation happens to drop those cells. az_solver='direct' with az_ridge=1e-2 removes it and leaves the other fourteen tilts within 1 dB, which is a stronger argument for the direct solver than the speed was. The check warns rather than raising or clipping. Raising would abort a fifteen-tilt volume for one bad tilt; clipping would hide a settings problem behind a field that looks plausible -- a scientist who sees 2471 dBZ knows the run is wrong, one who sees a silently clipped 100 does not. Its own warning category so a production pipeline can escalate it to an error while interactive work still gets the field and the message. It reuses the per-tilt overshoot_db that cones.py already computed and then discarded, and it runs after gridding rather than before, because the failure happens BETWEEN measured rays: every ray position can be fine and the interpolant between them still diverge, so no check on inputs can see it. Also plumbs cg_tol and field_units through grid_volume, which were reachable on the evaluator but not from the public entry point -- cg_tol matters here because it and n_cg jointly decide whether the solve diverges. On the tests: they exercise the predicate on constructed reports rather than through a synthetic volume, and that is deliberate. Several attempts to provoke divergence synthetically -- jittered azimuths, a sharp azimuthal wedge, deliberately tightened ray pairs at 0.003x to 0.09x nominal spacing -- all stayed bounded. The real volume's ill-conditioning comes from ray pairs 11x to 300x closer than nominal, and a fixture reproducing that faithfully enough to diverge would be fitting the test to one file. Testing the rule directly is honest about what is verified; a monkeypatched spy covers the wiring, which was the one seeded fault the first version of these tests missed -- the synthetic fixture never trips the guard, so a guard that was never called looked identical to a happy one.
Two changes the real-data work argued for. DEFAULT_ENGINE scipy -> dense, DEFAULT_SOLVER cg -> direct. This moves numbers for existing callers, which is why it was not done when the solver landed. It is done now because the numbers it moves were partly wrong. Per cone on the 15-tilt BNF C-SAPR2 volume, the old default at n_cg=12 had tilt 1 of 15 overshooting its input range by 107 dB and reaching 169 dBZ -- invisible in the gridded volume only because the vertical interpolation happened to drop those cells -- and across the resolution ramp the whole-volume range reached -299..297 dBZ. dense/direct holds -99..66 dBZ on the same grids with the other fourteen tilts within 1 dB. The deeper reason is that n_cg was doing regularisation's job. CG reaches well-conditioned directions first, so the iteration count is implicit early stopping; 12 lands near the optimum by luck, and 25 is off a cliff (-1907..2471 dBZ). A default whose safety depends on nobody raising an argument named "number of iterations" is not a good default. direct takes an explicit ridge, so the regularisation is a parameter rather than a side effect. dense rather than scipy because it needs no new dependency and has no kernel error: on this axis the DFT matrix is n_rays x n_modes, both ray counts, so forming it exactly costs 8 MB at 720 rays. az_engine='scipy', az_solver='cg' reproduces the old behaviour bit-for-bit, diagnostics included, and a test asserts it still does. EngineConcurrencyWarning covers finufft with n_jobs>1. Diagnosed rather than assumed: serially finufft is steady (0.29 s per tilt over six runs, 17.9 s per volume over three, 1.0x spread), but at n_jobs=4 three consecutive runs took 88.2, 6.5 and 16.3 s -- 13.6x -- while dense and torch held 1.0-1.1x on the identical path. So it is contention in finufft's own global planner state, not a slow transform. Kept and warned rather than dropped or silently serialised: the results are correct, a user may have measured their own case, and forcing n_jobs=1 for one engine would make an engine comparison quietly unfair. Nine tests failed on the default change and each was fixed by intent rather than by making it pass. Seven were about the CG iteration and now pin az_solver='cg' explicitly, with a comment saying why -- the default ignores n_cg, so without the pin both arms of those comparisons would have been the same direct solve and the tests would have passed vacuously. Two encoded the old default as the contract: "the default engine is scipy or reference" became a check against the dependency map (so changing the default cannot silently break minimal installs, which a name list permitted), and "the default matches the reference to round-off" became two tests -- the reference path is still exactly reproducible, and the new default differs from the old by a bounded amount, asserted in both directions so it proves the change happened AND that it did not move the field. Also corrects the GPU documentation, which claimed device='cuda' as though CUDA were the only option: device is passed to torch, so 'mps' works on Apple silicon too. Both GPU paths are explicitly marked untested -- every figure in that module is CPU -- and the docstring now notes torchkbnufft defaults its tables to single precision, so a GPU run is a precision change as well as a device change.
Both engines had claims attached to them that measurement did not support. finufft. Ran both engines at n_jobs=1 and 4, three repeats, eight sweep shapes. Serially finufft and dense are the same speed -- medians 0.97-1.00, 1.00-1.03x spread -- not 1.9x apart as this branch previously said. That earlier figure compared serial finufft against PARALLEL dense; reading across the wrong columns of my own table nearly cost the engine its place. Under n_jobs=4 it is bimodal: six of eight shapes make it the fastest CPU configuration measured (~8.3 s vs dense's ~10.1 s), and the rest take 12-19x longer. Best and worst single runs were 7.9 s and 177.8 s at fixed workload, against dense's 10.0-11.9 s over the same 24 runs, and slow runs appear in later repeats as well as first ones -- so a race in finufft's global planner state, not warm-up. The warning is reworded accordingly: it previously implied finufft was a poor choice, when it is often the fastest and only ever unpredictable. GPU. The torch engine's device= argument shipped documented but never executed. On an NVIDIA A10 it is slower than CPU on 15 of 15 tilts: 93.3 ms (float64) and 91.3 ms (float32) against 83.9 ms, with cuda.synchronize() around every timing so this measures execution rather than launch. The reason is the problem shape -- a 931 x 1050 solve is tens of milliseconds of arithmetic, so transfer and launch dominate. That is the same argument that makes dense competitive here: the azimuth axis is a ray count, not an image dimension. Docstrings now say measured-and-slower instead of untested, which is a more useful claim than either the old silence or a guess. Accuracy on the device is fine if anyone wants it anyway: float64 agrees with CPU to 1e-15 and float32 is 1.5x worse at 3.3e-4, because the Kaiser-Bessel kernel error dominates round-off. Setting that up exposed that dtype was hardcoded to complex128, so the single-precision path a GPU is actually fast at could not be expressed; it is now a parameter. Not committed, deliberately: an upper bound on torch. torchkbnufft segfaulted at import inside torch.jit.script on the CUDA container I first used, and I briefly pinned torch<2.6 with that as evidence. The cause was my own base image -- nvidia/cuda:12.4.1-devel fighting the host CUDA layer -- not the torch version. On a clean debian-slim image with plain pip torch, torch 2.13.0 imports torchkbnufft and runs a GPU transform without complaint. The pin was reverted; a plausible-looking traceback is not a substitute for varying the variable you chose arbitrarily.
census_sweep's valid_field fallback was next(iter(radar.fields)). Py-ART builds that dict from a set, so under hash randomisation the first key differs between processes for a byte-identical file. On a NEXRAD split cut this is decisive: reflectivity spans the full surveillance sweep (~460 km) while velocity stops at ~2 km, so range_max_valid_m flipped and dedup_sweeps kept the wrong member of the same elevation. Measured on KLOT 27 Jul 2026 VCP-212, which has six SAILS repeats at 0.48 deg whose peaks range from 66.5 to 85.0 dBZ: the gridded volume maximum alternated between 67.60 and 86.86 dBZ across identical runs, 9 of 20 fresh processes taking the wrong sweep. Serial and parallel were equally affected, which ruled out the tilt pool early. Nine consecutive fresh processes now agree exactly. Resolved via VALID_FIELD_PREFERENCE -- reflectivity-like names in order, then a sorted fallback -- so no path depends on dict iteration order. tests/gridding/test_census.py is new: 8 tests covering order invariance, the preference list's own contents (a Doppler moment there would silently reinstate the bug), and split-cut selection either way round. Mutation-checked: reverting to dict order fails 3, putting velocity first fails 5, an unsorted fallback fails 1, flipping the dedup tiebreak fails 2. The C-SAPR2 results already in the changelog are unaffected -- that strategy has no split cuts, so each elevation group has a single member.
Spectral nufft engines
|
You are seeing this message because GitHub Code Scanning has recently been set up for this repository, or this pull request contains the workflow file for the Code Scanning tool. What Enabling Code Scanning Means:
For more information about GitHub Code Scanning, check out the documentation. |
scollis
approved these changes
Aug 25, 2026
AdamTheisen
approved these changes
Aug 25, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
4 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.