ADD: NEXRAD Level 3 (NIDS) radial-product reader

[mgrover1]

Adds an xarray backend (engine="nexradlevel3") and open_nexradlevel3_datatree for NEXRAD Level 3 (NIDS) radial products. Level 3 is available on AWS in real time and archived back to ~2020 (s3://unidata-nexrad-level3), so this opens a large, previously unreadable dataset to the xradar stack.

Coverage

27 message codes across the three radial packet formats:

• Digital radial (packet 16): super-res reflectivity/velocity/spectrum width (153/154/155), legacy digital 94/99, dual-pol ZDR/CC/KDP (159/161/163), hydro class (165), digital accumulations (170/172-175)
• Run-length encoded (packet AF1F): legacy reflectivity/velocity/spectrum width (19/20/25/27/28/30), storm-relative velocity (56), legacy accumulations (78-80)
• Generic data packet (packet 28, XDR): DPR precip rate (176), hybrid hydro class (177) — via a small inline XDR unpacker, no new dependency

Deferred products (VIL/EET special encodings, TDWR, legacy hybrid variants) raise a clear NotImplementedError naming the product; follow-up tracking issue to come. A multi-tilt volume of the same product (e.g. N0B-N3B) assembles into a DataTree ordered by fixed angle, following the IMD multi-file precedent.

Deliberate divergences from Py-ART

Each verified against the ICD (2620001) and real data:

What	This PR	Py-ART
azimuth coord	ray centers (start + width/2)	ray start angles
range coord	bin centers at the true ICD bin size	edges at packet "range scale" — that field is actually floor(1000*cos(elevation)), so ranges drift up to ~580 m by the last bin
units	uniform m/s (VRADH/WRADH/SRMV) and mm (ACCUM/RATE)	legacy velocity in knots, precip in inches
legacy 16-level thresholds	per-flag sign/scale incl. the x0.01 flag	flags[0] applied to all levels
below-threshold vs range-folded	flag counts read from PDB halfwords 36-38; RF exposed as a <moment>_range_folded mask	fixed raw < 1/raw < 2 masks — on a real DPR file this marks ~92% of bins (genuine zero rain rate) as missing

model.py gains SRMV/ACCUM/HCLASS (and RATE) in the canonical moment names; the two test_io.py expectation updates reflect real moments in existing sample files that the canonical set previously missed.

Verification

• 37 synthetic unit tests against byte-built files with exact ground truth (every decode scheme, packet format, mask semantics, metadata field, bz2 round-trip, error paths)
• 11 real-file integration tests with pinned statistics (skip-guarded until the open-radar-data samples land — companion data PR to follow: 12 LOT products, ~0.9 MB, one coherent scene)
• All 12 real LOT products cross-checked against Py-ART 2.2.5 and MetPy 1.7.1: decoded values bit-identical after the documented unit conversions
• Example notebook opens live AWS data, plots a georeferenced PPI, and assembles a 4-tilt volume

Performance (M-series laptop, warm, mean of 20)

Case	xradar	Py-ART	MetPy
N0B super-res (720x1840)	27 ms	31	32
N0S legacy RLE	4.7 ms	6.0	3.4
DPR (XDR)	14 ms	76	24
60-file mixed batch, per file	14 ms	20	13

Known limitation

Level 3 carries no per-ray times, so every ray holds the volume scan start time: multi-tilt volumes export to CfRadial2 but not CfRadial1 (documented in the docstring).

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 94.63%. Comparing base (79ba495) to head (9c8826c).

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #392      +/-   ##
==========================================
+ Coverage   94.20%   94.63%   +0.43%     
==========================================
  Files          29       30       +1     
  Lines        6417     6934     +517     
==========================================
+ Hits         6045     6562     +517     
  Misses        372      372              

Flag	Coverage Δ
notebooktests	0.00% <0.00%> (ø)
unittests	94.63% <100.00%> (+0.43%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[aladinor]

Hey @mgrover1 — really nice work on this one. I did a deep review pass: re-ran the suite in a clean env (all green), checked the coverage claim (it's real — 99% branch, 0 missed statements), and verified the decode math against the ICD (linear halfword scale/offset, float scale/offset, legacy-16 per-flag sign/scale, the unit conversions, and the hw36–38 flag-count handling all check out). The structured-dtype packet-16 path and the vectorized RLE expansion also hold up nicely.

That said, I went adversarial on the edge cases and found a handful of things I think need fixing before merge. Everything in the first section I reproduced with the PR's own test helpers.

Bugs (reproduced)

1. Multi-component generic packets crash with a raw AttributeError
_unpack_components() returns a bare component when num == 1 but a list otherwise, and _read_generic_packet only handles the bare shape (component.radials). A packet-28 file with two components dies with:

packet = _generic_packet28(2, 2, data, ncomponents=2)
NEXRADLevel3File(io.BytesIO(build_level3_file(msg_code=176, packet=packet)))
# AttributeError: 'list' object has no attribute 'radials'

test_xdr_multiple_components covers this input but drives _Level3XDRParser directly, so it never hits the crash. I'd drop the num == 1 collapse entirely (always return a list, caller takes components[0] with an explicit check) — that also removes the test_xdr_single_parameter_collapses quirk-pinning test.

2. mask_and_scale=False resurrects the Py-ART zero-rain bug on DPR
_build_l3_sweep sets data_attrs["_FillValue"] = 0 unconditionally in the raw path. For DPR, raw 0 is a genuine zero rain rate (leading flag count 0) — the PR description's own headline divergence from Py-ART. But run the raw output through xr.decode_cf and:

get_data()[0]:   [ 0.    2.54 10.16]
decode_cf()[0]:  [ nan   2.54 10.16]

So the decoded path and the raw+CF path disagree, and the raw path re-masks exactly the bins the PR says it rescues. _FillValue should follow the leading-flag count (omit it when leading == 0), same as the valid_min logic right below.

3. open_nexradlevel3_datatree combines files it should reject
Only the message code is checked, so both of these "work":

• KLOT N0B + KMKX N0B → merges silently; the KMKX sweep inherits KLOT root coordinates, so georeferencing is quietly wrong for that sweep.
• The same tilt from two different volume scans → a "volume" with sweep_fixed_angle = [0.5, 0.5].

Both come up naturally when someone globs LOT_N?B_* from the bucket. I'd validate site (lat/lon or the radar id from the text header) and volume scan time, and error on duplicate fixed angles.

4. Truncated bz2 escapes the except OSError
bz2.decompress raises ValueError for truncation ("Compressed data ended before the end-of-stream marker...") and OSError only for corrupt streams (checked on CPython 3.12). So the "file is corrupt or truncated" message never fires for the truncated case — the one it names. except (OSError, ValueError, EOFError) fixes it; test_corrupt_bz2_raises only covers the corrupt-tail case.

5. Files truncated inside the headers leak raw struct.error
A file cut inside the MHB/PDB/symbology header raises struct.error: unpack_from requires a buffer of at least ... from _unpack_from_buf — and struct.error is not a ValueError subclass, so callers catching ValueError per the class's contract miss it. A length check up front (or wrapping the header unpacks) would keep the error contract consistent.

6. The nbins-override rule turns ICD halfword padding into a phantom bin
The ICD requires an even byte count per packet-16 radial, so an odd-nbins product would carry nbytes = nbins + 1 — and the "byte count wins" rule then decodes a 5-bin radial as 6 bins with a garbage pad column (reproduced), growing range by one gate. All currently supported products have even bin counts, so it's latent — but the override also only peeks at radial 0, and the ragged fallback silently truncates any later radial wider than that (also reproduced: a 10-byte radial 1 behind an 8-byte radial 0 loses its last two bins with zero warnings). Suggest never widening by exactly +1 and warning when the counts disagree.

7. Raw-mode attrs vs decoded-mode mask disagree for linear_hw
In the raw path _build_l3_sweep hardcodes leading, valid_max = 2, None for linear_hw, while get_data calls get_flag_counts() for the same products and will honor plausible hw36–38 values including valid_max. Reproduced on a synthetic 153 file with max_val=255, leading=2, trailing=1: get_data NaNs raw 255, the raw-mode attrs carry only valid_min=2, and the two paths disagree on 2 of 5 bins. One source of truth (get_flag_counts) feeding both would fix it.

While verifying this I hit a bigger problem with the raw-mode strategy itself: xarray's decode_cf does not apply valid_min/valid_max at all (only _FillValue/missing_value and scale/offset), so the inline comment "valid_min/max are in packed (raw) units so CF-aware decoders mask flag levels (below-threshold, range-folded) too" doesn't hold for xarray itself. Concretely, the range-folded level (raw 1) survives xr.decode_cf as a real-looking −32.5 dBZ for every linear_hw product. The raw path probably needs to rely on _FillValue (derived from the leading-flag count, per bug 2) rather than valid-range attrs, or document that flag levels must be masked manually via range_folded_raw_value.

Design / consolidation (should-fix, not blocking)

• ProductSpec.post_scale is dead for the float family. get_scale_offset overwrites factor = spec.post_scale with hardcoded 0.01 * IN_TO_MM / IN_TO_MM for "precip"/"rate", while those table rows carry post_scale=1.0 — contradicting the field's own docstring. Put the factors in the table and "precip"/"rate" stop being decode schemes at all; the enum drops to 4 and several if-chains shrink.
• SUPPORTED_VERSION_NUMBERS is a fully parallel 26-key dict of PRODUCT_TABLE — a max_version field on ProductSpec (with defaults=) removes the drift risk. Two related nits: the ProductSpec docstring omits the "rate" scheme, and entry 177's bin_size=250.0 is never read (packet-28 geometry comes from the component) — clearer as None like 176.
• packet_header shape-punning: _read_generic_packet grafts first_bin/range_scale onto the dict with different units/types than the radial-packet fields (meters-float vs index-int), plus the getattr(self, "_range_from_centers", False) sentinel — and packet_header["nbins"] goes stale when the nbins override fires, so the documented attribute lies. Uniform parse-time attrs (self._first_center_m, self._gate_spacing_m) would collapse get_range to one line.
• Re-raise keyed on message text: if isinstance(err, ValueError) and "expand" in str(err): raise breaks silently if the message is ever reworded — a distinct exception type (or moving the expansion check out of the try) is sturdier.
• Reuse of the shared helpers: _build_l3_root duplicates common._assign_root almost line for line (I checked — reuse works fine with L3's uniform times), _sweep_for_datatree duplicates common._attach_sweep_groups, and the (name, fmt) tuple tables + _unpack_from_buf re-implement the iris._unpack_dictionary/_get_fmt_string pattern that both nexrad_level2.py and uf.py already import.
• Dead code: the _apply_site_as_coords call in _sweep_for_datatree is provably a no-op (station vars are dropped on the line before), which makes the datatree's site_as_coords kwarg do nothing while its docstring promises sweep-level coords. L2's kwarg is equally inert so this is parity — but the docstring should be fixed at minimum.
• File-handle machinery: the whole file is read eagerly in __init__, so keeping the handle open (one fd per sweep in a datatree), __del__ = close, and the context-manager plumbing serve nothing — and ds.close() currently closes user-provided file-likes. Closing path-opened handles right after read() deletes most of it.

Performance / memory

• The packet-28 path pins the entire XDR payload until ds.close(): unpack_int_array returns np.frombuffer views and the _RadialData tuples in gen_data_pack keep them alive alongside the raw_data copy (~2.5× necessary memory per DPR/HHC open). Dropping the radial references after raw_data is built fixes it.
• (raw * scale + offset) with Python-float scale promotes u1 → float64 (~10 MB temp on 720×1840, a few ms of the advertised 27). np.float32(scale)/np.float32(offset) keeps it f32 end-to-end.

model.py side effects worth a look

• imd.py's comment "HCLASS does not [have a canonical]" is now false, and its attrs-overlay loop will now overwrite IMD's file-derived HCLASS attrs with the NEXRAD-flavored canonical ones — untested behavior change in an untouched backend.
• util.get_sweep_dataset_vars(..., non_standard=False) now returns RATE for Furuno and HCLASS for CfRadial1/IMD files — probably worth its own changelog line.
• ACCUM/SRMV are NEXRAD-invented names entering the canonical set; ODIM's accumulation quantity is ACRR. Might be worth an explicit openradar decision before downstream code standardizes on them.

Questions / scope

• The reader requires an SDUS header in the first 80 bytes and a 30-byte text header. Checking the reference implementations' source: Py-ART has the identical buf.find(b"SDUS") requirement (so this is parity with Py-ART — I originally claimed otherwise here, corrected), but MetPy is more permissive: its wmo_finder regex accepts (?:NX|SD|NO)US and the header is optional (if match:), so headerless and NXUS/NOUS files that MetPy reads get "Not a valid NEXRAD Level 3 file" here. Fine to scope to the AWS bucket, but the docstring should say so — and note the TDWR NotImplementedError may be unreachable if TDWR files ship under non-SDUS headings, since the header check fires first.
• Legacy-16 coded values (flag 0x80) all become NaN, so range folding in legacy velocity/SRMV products is merged into "missing" while digital products get a dedicated <moment>_range_folded mask — worth a docstring note or a follow-up issue.

Test suggestions (they line up with the bugs)

• Route the multi-component XDR test through NEXRADLevel3File (catches Next initialization steps #1).
• Add an xr.decode_cf round-trip on mask_and_scale=False output (catches ADD: build sdist/wheel upload to pypi #2).
• Datatree negative tests: mixed sites, duplicate tilts, mixed volume times (catches Update env to include pydata sphinx #3).
• Header-level truncation tests (catches FIX: Move to workflow dispatch #5), and a synthetic test for product 177 through the file-level generic-packet path — only 176 is covered until the real HHC sample lands.

Happy to open follow-up issues for the deferred items (legacy RF, WMO header flexibility) or help with any of the fixes. The core decode work here is really solid — most of the above is hardening around it.

---

Edited: corrected the Py-ART header claim in "Questions / scope" (Py-ART requires SDUS exactly like this PR — only MetPy is more permissive; verified in both sources), and added executed-repro evidence to items 6 and 7, including that xr.decode_cf ignores valid_min/valid_max.