diff --git a/AGENTS.md b/AGENTS.md new file mode 100644 index 0000000..27a9dd8 --- /dev/null +++ b/AGENTS.md @@ -0,0 +1,16 @@ +Keep code concise with minimal but useful documentation. + +Use uv for dependency management and `uv run` to run any python code, including small inline runs. + +Put any session notes, longer text analyses, purely temporary code, or any other session artifacts you want to keep but that aren't part of the main flows in `claude_notes/`. Put plans in `claude_plans/`, prefix them with the date in YYYYMMDD-title.md format, and update these when you have completed them for future reference. + +Do not edit anything outside of this directory. Put any outputs into an `outputs/` directory. + +When long tasks are required, use subagents to perform them. + +Review the `docs/` directory before starting on any task. `docs/architecture.md` describes the basic outline. + +Run the tests (using pytest) after making any changes. Add tests to verify new functionality you've been asked to add, but do not modify existing tests without asking. + +Useful references: +* Required surface SNR calculations: https://github.com/thomasteisberg/required_surface_snr \ No newline at end of file diff --git a/CLAUDE.md b/CLAUDE.md deleted file mode 100644 index f77d461..0000000 --- a/CLAUDE.md +++ /dev/null @@ -1,14 +0,0 @@ -Keep code concise with minimal but useful documentation. - -Use uv for dependency management and `uv run` to run any python code, including small inline runs. - -Put any session notes, longer text analyses, purely temporary code, or any other session artifacts you want to keep but that aren't part of the main flows in `claude_notes/`. Put plans in `claude_plans/`, prefix them with the date in YYYYMMDD-title.md format, and update these when you have completed them for future reference. - -Do not edit anything outside of this directory. Put any outputs into an `outputs/` directory. - -Review the `docs/` directory before starting on any task. `docs/architecture.md` describes the basic outline. - -Run the tests (using pytest) after making any changes. Add tests to verify new functionality you've been asked to add, but do not modify existing tests without asking. - -Useful references: -* Required surface SNR calculations: https://github.com/thomasteisberg/required_surface_snr \ No newline at end of file diff --git a/CLAUDE.md b/CLAUDE.md new file mode 120000 index 0000000..47dc3e3 --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1 @@ +AGENTS.md \ No newline at end of file diff --git a/claude_notes/calib_img_availability_probe.py b/claude_notes/calib_img_availability_probe.py new file mode 100644 index 0000000..1cc2eb7 --- /dev/null +++ b/claude_notes/calib_img_availability_probe.py @@ -0,0 +1,46 @@ +"""Probe img_01 file availability for every season in the production stores. + +Answers: how often would the img-overlap check be impossible (fallback needed)? +Uses HTTP HEAD against the public OPR site; no data downloads. +""" +import icechunk +import zarr +import numpy as np +import requests +from xopr import OPRConnection + +seasons = set() +for prefix in ["antarctica", "greenland", "ase", "utig", "crosssystem"]: + storage = icechunk.s3_storage( + bucket="opr-radar-metrics", prefix=f"icechunk/{prefix}", + region="us-west-2", anonymous=True, + ) + repo = icechunk.Repository.open(storage=storage) + root = zarr.open_group(repo.readonly_session(branch="main").store, mode="r") + colls = root.attrs.asdict().get("frame_collections") + if colls is None: + print(f"{prefix}: no frame_collections attr, skipping") + continue + seasons.update(colls) +print(f"{len(seasons)} unique seasons\n") + +opr = OPRConnection() +for coll in sorted(seasons): + try: + frames = opr.query_frames(collections=[coll], max_items=3) + except Exception as e: + print(f"{coll}: query failed ({type(e).__name__})") + continue + statuses = [] + for fid in frames.index[:2]: + # Data_20161014_03_001 -> segment 20161014_03 + parts = fid.replace("Data_", "").rsplit("_", 1) + seg = parts[0] + url = (f"https://data.cresis.ku.edu/data/rds/{coll}/CSARP_standard/" + f"{seg}/Data_img_01_{fid.replace('Data_', '')}.mat") + try: + r = requests.head(url, timeout=30, allow_redirects=True) + statuses.append(str(r.status_code)) + except Exception as e: + statuses.append(type(e).__name__) + print(f"{coll}: img_01 HEAD -> {statuses}") diff --git a/claude_notes/calib_weights_survey.py b/claude_notes/calib_weights_survey.py new file mode 100644 index 0000000..bf4e75a --- /dev/null +++ b/claude_notes/calib_weights_survey.py @@ -0,0 +1,56 @@ +"""Survey img_comb_weights / img_comb_weights_mode across all store seasons. + +Loads one combined frame per season and searches every param_* attr dict for +weights-related fields. +""" +import numpy as np +import icechunk +import zarr +from xopr import OPRConnection + +seasons = set() +for prefix in ["antarctica", "greenland", "ase", "utig"]: + storage = icechunk.s3_storage( + bucket="opr-radar-metrics", prefix=f"icechunk/{prefix}", + region="us-west-2", anonymous=True, + ) + repo = icechunk.Repository.open(storage=storage) + root = zarr.open_group(repo.readonly_session(branch="main").store, mode="r") + colls = root.attrs.asdict().get("frame_collections") + if colls: + seasons.update(colls) + + +def find_fields(d, needle, prefix="", depth=0, out=None): + if out is None: + out = [] + if depth > 5 or not isinstance(d, dict): + return out + for k, v in d.items(): + if isinstance(v, dict): + find_fields(v, needle, prefix + k + ".", depth + 1, out) + elif needle in k.lower(): + out.append((prefix + k, v)) + return out + + +opr = OPRConnection() +for coll in sorted(seasons): + try: + frames = opr.query_frames(collections=[coll], max_items=1) + ds = opr.load_frame(frames.iloc[0], data_product="CSARP_standard") + except Exception as e: + print(f"{coll}: load failed ({type(e).__name__}: {str(e)[:80]})") + continue + hits = [] + for pk in [k for k in ds.attrs if k.startswith("param")]: + for path, v in find_fields(ds.attrs[pk], "img_comb_weight"): + s = str(v).replace("\n", " ")[:80] + hits.append(f"{pk}.{path}={s}") + if hits: + # dedupe identical values + print(f"{coll} ({frames.index[0]}):") + for h in sorted(set(hits)): + print(f" {h}") + else: + print(f"{coll} ({frames.index[0]}): no img_comb_weights fields") diff --git a/claude_notes/calibration_baseline_results.md b/claude_notes/calibration_baseline_results.md new file mode 100644 index 0000000..0595883 --- /dev/null +++ b/claude_notes/calibration_baseline_results.md @@ -0,0 +1,152 @@ +# Calibration baseline results (2026-08-31) + +Steps 1–4 of `claude_plans/20260831-radiometric-calibration-checks.md`. +Code: `src/radar_return_statistics/calibration.py` + `calibration_cli.py`; +tests: `tests/unit/test_calibration.py` (16 tests; full unit suite 64 passing). +Outputs: `outputs/calibration/baseline/` and `baseline_extra/` +(parquet + figures + `manifest.json`). No store writes anywhere. + +## Check 1: image-combine residual seam offsets + +36 frames across 5 seasons (2 frames × 5 segments for 2018_Greenland_P3, +2016_Antarctica_DC8, 2018_Antarctica_DC8; 3 frames 2014_Greenland_P3; +3 frames 2022_Antarctica_BaslerMKB). Zero load failures; statuses: +27 `ok`, 8 `insufficient_overlap` (real geometry, see below), 1 more in +extra run. All four param layouts exercised (`param_combine.array_param`, +`param_array.array`, empty-`img_comb`-era, BaslerMKB full-gate). + +### Per-season residual seam offsets (worst pair, frame means, `ok` frames) + +| season | pair | n | mean (dB) | frame-to-frame std | weights | +|---|---|---|---|---|---| +| 2018_Greenland_P3 | 1 | 8 | **+1.07** | 0.21 | all 0 (raw stitch) | +| 2016_Antarctica_DC8 | 1 | 8 | **+1.45** | 0.65 | all 0 (raw stitch) | +| 2016_Antarctica_DC8 | 2 | 2 | +0.12 | 0.09 | | +| 2018_Antarctica_DC8 | 1 | 8 | **−0.57** | 0.99 | w2/w3 ≈ **106–111 dB, per-frame** | +| 2018_Antarctica_DC8 | 2 | 2 | −0.07 | 0.49 | | +| 2022_Antarctica_BaslerMKB | 1 | 3 | +0.55 | 0.04 | all 0 | + +- Weight recovery was rock solid everywhere: per-frame spread ≤ 0.05 dB on + every frame, no `weight_recovery_failed`. +- **2018_Antarctica_DC8's auto weights vary frame to frame (106–111 dB)** — + empirical recovery is mandatory, not optional, for that season; the + declared params could never supply these. Even after auto correction one + frame retains a **−2.92 dB** residual seam (Data_20181016_01_051) — auto + mode does not always fix the seam. +- Distribution vs the ~3 dB downstream threshold: max |frame mean| = 2.92 dB, + 2/27 frames above 2 dB, 0 above 3 dB. A 3 dB cut would currently reject + borderline frames only; the bulk sits at |offset| ≤ 1.5 dB. +- The healthy-frame ~+1.8 dB from exploration is confirmed as season-typical + for 2016 DC8 (+1.45 ± 0.65) and 2018 Greenland (+1.07 ± 0.21): real, stable + seam steps in the combined products. + +### Mean vs median sensitivity (Codex #16 check) + +Max |mean − median| over ok frames = **1.21 dB** (Data_20161031_06_015: mean +2.66 / median 1.45, MAD 1.04, drift +2.02 — a genuinely drifting frame). +Typical difference ≪ 0.3 dB. The mean headline stands (user decision); high +MAD/drift flags the frames where it matters, and both center measures ship in +the parquet/store. + +### insufficient_overlap is real geometry, not a bug + +2014_Greenland_P3 comes back `insufficient_overlap` on all frames at typical +AGL: its `img_comb` (T_comb=3 µs, T_guard=2.64 µs, found *non-empty* at +`param_combine.array_param` — the survey's "empty" was the `combine` path) +plus img1 ending at 7.85 µs leaves **no artifact-free overlap window** when +td_surf ≈ 3.2 µs (window [6.2, 5.2] µs is empty): the combine switched earlier +than the T_comb-safe point, so any measurement would include img2's saturated +surface response. The scattered insufficient_overlap frames in other seasons +(5/36) are the same geometry at high AGL. Honest NaN + status is the designed +behavior; weights still recover (2014: all ≈ 0 → raw stitch). + +## Check 2: surface saturation (all five stores, read-only) + +### Season fits + +| store | season | status | level (dB) | overall slope (dB/dec) | pileup | +|---|---|---|---|---|---| +| antarctica | 2008_Antarctica_BaslerJKB | fit_ok | −15.3 | +2.8 | 0.03 | +| antarctica | 2012_Antarctica_DC8 | **fit_ok** | **−33.1** | −9.8 | **0.17** | +| antarctica | 2013_Antarctica_Basler | no_plateau | | −16.6 | | +| antarctica | 2013_Antarctica_P3 | insufficient_support | | | | +| antarctica | 2014_Antarctica_DC8 | no_plateau | | −25.9 | | +| antarctica | 2016_Antarctica_DC8 | no_plateau | | −19.0 | | +| antarctica | 2017_Antarctica_Basler | fit_ok | −25.5 | −4.1 | 0.02 | +| antarctica | 2017_Antarctica_BaslerJKB | insufficient_support (n=178) | | | | +| antarctica | 2017_Antarctica_P3 | no_plateau | | −27.9 | | +| antarctica | 2018_Antarctica_DC8 | fit_ok | −21.8 | +3.9 | 0.02 | +| antarctica | 2019_Antarctica_GV | insufficient_support (span 0.42 dec) | | | | +| antarctica | 2022_Antarctica_BaslerMKB | no_plateau | | −20.2 | | +| antarctica | 2023_Antarctica_BaslerMKB | no_plateau | | −15.8 | | +| greenland | 2013_Greenland_P3 | **fit_ok** | **−50.4** | +6.2 | 0.01 | +| greenland | 2014_Greenland_P3 | no_plateau (two-regime, see below) | | −24.0 | | +| greenland | 2016_Greenland_P3 | insufficient_support | | | | +| greenland | 2017_Greenland_P3 | no_plateau | | −24.3 | | +| greenland | 2018_Greenland_P3 | no_plateau | | −21.1 | | +| greenland | 2019_Greenland_P3 | no_plateau | | −29.5 | | +| ase | 2012_Antarctica_DC8 | **fit_ok** | −33.2 | −9.0 | **0.11** | +| ase | 2014_Antarctica_DC8 | fit_ok | −30.3 | −10.5 | 0.02 | +| ase | 2016_Antarctica_DC8 | no_plateau | | −12.5 | | +| ase | 2018_Antarctica_DC8 | fit_ok | −25.1 | −6.7 | 0.08 | +| utig | 2008_Antarctica_BaslerJKB | fit_ok | −14.4 | −1.4 | 0.01 | +| utig | 2017_Antarctica_BaslerJKB | insufficient_support (n=276) | | | | + +Headline: **2012_Antarctica_DC8 is the most saturated season** — flat top at +−33 dB across 400 m–10 km with 11–17% of traces within 1 dB of the ceiling +(consistent with the ASE saturation discussion around Chu et al. 2021). +2013_Greenland_P3 confirms the exploration finding (razor-flat top, −50.4 dB). +Every real-data `fit_ok` came from the **flat-envelope rule** (overall +upper-quantile slope far shallower than −20 dB/decade); the piecewise +partial-saturation model never won on real data. + +### Fit hardening during the baseline + +The first fit version produced a false negative (2014_Greenland) and false +positives (2018_Greenland with a rising, then a −57 dB/dec plunging second +segment). Root cause: **two-regime seasons** — a flat upper envelope at one +level below ~1 km range, a ~15–20 dB step, then a second flat envelope +(altitude-dependent receiver gain/attenuation, not a clip-level crossing). +Three constraints were added to `SAT_DEFAULTS` and unit-tested against a +synthetic two-regime population: + +- `max_slope_beyond_db_per_decade = −8` — the post-breakpoint segment must fall; +- `min_slope_beyond_db_per_decade = −35` — but not plunge like a regime cliff; +- `continuity_tol_db = 3` — it must continue from the plateau level. + +### Two-regime seasons (key finding for review) + +`figures/saturation_greenland_2014_Greenland_P3.png`: two *flat-topped* bands +(~−30 dB below ~800 m range, ~−48 dB beyond). **Both regimes individually +look ceiling-limited.** The season fit correctly refuses a single ceiling +(`no_plateau`) but thereby under-reports real saturation. Per-segment fits +cannot rescue this: segments fly at near-constant altitude, so 34/35 segment +fits returned `insufficient_support` (no range lever arm within a segment) — +season (or altitude-regime) level is the only workable granularity. + +## Proposed thresholds + +- **Image combine (downstream guidance, not enforced)**: |frame mean residual| + > 3 dB, as the user anticipated — currently rejects 0/27 measured baseline + frames (max 2.92); 2 dB would flag 2/27. +- **Saturation defaults** (all recorded in the manifest): flat-envelope rule + at slope ≥ −12 dB/decade; piecewise acceptance needs SSE improvement ≥ 30%, + plateau internal slope ≤ 5, slope-beyond in [−35, −8], continuity ≤ 3 dB; + support minimums 50 traces/bin, 6 bins, 0.5 decades. + +## Open questions needing user decisions + +1. **Two-regime seasons**: stratify by altitude band (e.g. split at the + step) or by receiver attenuation settings from params, fitting a ceiling + per regime? Without this, 2014_Greenland_P3 (and possibly 2017/2019 GL) + get no margins despite visibly ceiling-limited populations. +2. **insufficient_overlap frames** (2014_Greenland at typical AGL; scattered + high-AGL frames elsewhere): accept honest NaN, or add a clearly-labeled + relaxed window (which would include part of img2's saturated surface + response)? +3. **Piecewise model**: it never fired on real data (all detections via the + flat rule). Keep it (tested, harmless) or simplify to flat-rule-only? +4. 2019_Antarctica_GV fails on range span (0.42 < 0.5 decades) — relax + `min_span_decades` for high-altitude-only platforms, or leave? +5. 2016 GL / 2013 Antarctica P3 / 2017 BaslerJKB fail support in their + stores but may pass in a merged view — worth a cross-store union fit? diff --git a/claude_notes/calibration_integration_notes.md b/claude_notes/calibration_integration_notes.md new file mode 100644 index 0000000..fa3343f --- /dev/null +++ b/claude_notes/calibration_integration_notes.md @@ -0,0 +1,183 @@ +# Calibration integration (steps 5–6) — 2026-08-31 + +Branch `calibration-checks`, uncommitted. No S3 writes; all runs local. + +## What changed + +- `config.py` — `normalize_config()` shared by `load_config` and + `process_frame`; `processing.calibration.{img_combine: true, + image_load_retries: 2}` defaults ON at every entry point. +- `calibration.py` — + - images now load **sequentially** and convert to numpy immediately + (dataset closed per image); + - `_align_traces`: images reindexed to the combined frame's (decimated) + slow_time by nearest within 1 s; unmatched traces → NaN. Needed because + the pipeline decimates the combined frame but images are full-resolution + (the CLI passed full-res frames so this never triggered before). + - `min_traces_per_pair` capped proportionally + (`min(30, max(5, 0.25·n_traces))`) — the absolute 30 was tuned against + full-res frames (~1% of traces) and was a far stricter relative bar for + ~40-trace decimated frames. + - new: `fit_season_by_source` (img1 / img2 populations fit separately, + `min_span_decades_source_restricted: 0.25`), `cross_cap_step` + (empirical img2 bias), `margins_by_source`. +- `processing.py` — `run_frame_calibration()` (never fails a frame; statuses + degrade); emits `img_comb_offset_dB` (float32), `img_comb_pair` (int8), + `surface_source_image_index` (int8), `surface_ceiling_margin_dB` (all-NaN + float32) every frame regardless of enablement; frame attrs + `frame_img_comb_status` / `frame_img_comb_weights_mode` / + `frame_img_comb_offset_dB` (only when measured). `decimate_frame()` factored + out for backfill reuse. Offsets stored only for `ok` / `partial_images` + (`weight_recovery_failed` residuals stay parquet-only, per plan). +- `store.py` — `_zarr_append` auto-migrates (missing variables created sized + to the store, fill NaN/−1/False/"", loud warning) instead of silently + skipping; `science_fingerprint` (n_traces + sha256 of sorted + processed_frames), `saturation_stale`, `write_saturation_results`, + `update_frame_calibration` (in-place backfill; contiguous slice or oindex). +- `runner.py` — `FRAME_ATTR_CASTS` typed per-frame attrs (strings ride along; + NaN floats → None so attrs stay JSON-valid); stale-second-pass warning after + the final commit. +- New commands: `saturation_pass.py` (`--dry-run`; fits per season by source + image, writes margins + `saturation` root attr with fingerprint), + `calibration_backfill.py` (`--statuses load_error,disabled`, `--limit`). +- Viewer: `toFloat64Array` byte-reinterpret bug fixed (int arrays now convert + element-wise); int8 sentinel (−1) mapped to NaN on load + (`INT_SENTINEL_VARIABLES`); `integer` formatting flag; variables added to + config + dropdown (`img_comb_pair` loadable but not a dropdown layer — + it's bookkeeping for the offset). +- Docs: architecture.md (variables, frame/season attrs, pipeline step, cost), + processing.md (calibration config + commands), data_access.md (downstream + filtering summary, user's "provenance not validity" framing). +- Tests: `tests/unit/test_calibration_integration.py` (9 tests: migration, + enabled/disabled transitions, placeholder emission, degrade-to-status, + in-place backfill, split-population fit + cross-cap recovery, second-pass + write, fingerprint staleness). Suite: **73 passed**. +- **Existing-test modification (flagged)**: the exact-variable-set assertion + in `test_process_frame_output_variables` gained the four new variable + names — unavoidable given the user-approved schema; nothing else touched. + +## Verification (local smoke store, outputs/icechunk_store_calib_smoke) + +`test_config.yaml`-equivalent run (David/Drygalski, 3 frames of +2013_Antarctica_P3, calibration on): statuses `['ok', 'insufficient_overlap', +'insufficient_overlap']`, frame offset +0.51 dB on the ok frame, per-trace +offsets finite exactly on that frame, `surface_source_image_index` all 1, +margins all-NaN. Second pass: `by_source_image` mode, +`insufficient_support` both populations (117 traces — expected), attrs + +fingerprint written, `saturation_stale` False after its own commit. Backfill +CLI: correctly reports nothing retryable. + +## Benchmark (3 real 2018_Greenland_P3 frames, sequential, streaming cache) + +| | wall time | downloaded | per-frame cache high-water | +|---|---|---|---| +| calibration OFF | 21.1 s | 200 MB | 118 MB | +| calibration ON | 65.3 s | 512 MB | 318 MB | + +≈ **3.1× wall time, 2.6× download volume, 2.7× temp disk** per frame with the +check on (this season has 3 images). Worker per-frame cache cleanup already +bounds disk (unchanged). For production planning: multiply current per-store +processing time/bandwidth by ~3. + +Benchmark statuses: Data_20180426_03_001 `ok` (+1.1 dB, matches baseline); +both Data_20180404_02_* frames `insufficient_overlap` with n_valid=0 at full +resolution too — real geometry, not a decimation artifact (verified +explicitly). + +## Reprocessing regression verification (2026-08-31, post-review) + +`claude_notes/reprocess_regression_check.py`: re-ran `process_frame` with the +branch code (calibration ON, production configs) on 5 frames — 2018/2013 +Greenland P3, 2022 BaslerMKB, 2016/2012 Antarctica DC8 — and compared all 25 +pre-existing per-trace variables trace-by-trace against the production store +values (NaN-pattern equality + float tolerance). **0 mismatches on all 5 +frames**; trace counts identical. Calibration statuses: ok (2016 DC8, +BaslerMKB), insufficient_overlap (2018 GL), params_missing (2013 GL — no +imgs; 2012 DC8 — see below). + +Follow-up worth one look pre-production: **2012_Antarctica_DC8 returns +`params_missing`** even though its img files exist (availability probe was +200/200). Its 2012-era attrs presumably carry neither a recognized img_comb +path nor usable Tpd. Since 2012 DC8 is the most saturated season, recovering +its seam check (extend param search or Tpd-default window) may be worthwhile. + +## 2012-era param support + noise-floor fix (2026-08-31, post-review) + +**2012_Antarctica_DC8 `params_missing` resolved** — three era differences, all +handled in `calibration.py` (+ unit tests, suite at 75): +1. Combine params live under `param_combine_wf_chan` (new search path). +2. Old 2-element-per-pair `img_comb` format `[T_comb, T_guard]` (no T_blank), + normalized to modern `[T_comb, -inf, T_guard]` — applied only under the + 2012-era attr so modern 6-element vectors are never misread. Verified + empirically: combined follows img1 until ≈ surf + 10 µs (T_comb = Tpd2), + weight recovery w2 = +0.025 dB (spread 0.004). +3. `radar.wfs` is a list of per-waveform dicts (Tpd extracted per dict); + `imgs` entries are (N,2) wf-adc rows (wf axis chosen as the constant one). + +**Noise-floor contamination fix**: the first rerun then hit +`insufficient_overlap` because img1's record ends mid-ice (18.8 µs), putting +the record-tail noise window (6.8–11.8 µs) in strong englacial signal — floor +inflated ~45 dB, SNR gate starved. `image_noise_floor_db` now takes the +elementwise min of the tail median and the global 10th percentile (a +slightly-low floor only admits extra bins into a median — safe). Not +2012-specific: any short-gate image whose record ends above the bed had the +same hazard. + +**Final regression check: all 5 frames PASS** (0/25 variable mismatches each); +2012 DC8 now status `ok`. Note the floor change makes the SNR gate slightly +more permissive than in the earlier baseline parquet — offsets are +medians-over-bins so ok-frame numbers move negligibly, and the production run +supersedes the baseline anyway. + +## ASE production rehearsal (2026-08-31, completed) + +Full local-then-sync flow validated end to end on the ASE store: +1. Scratch-prefix sync smoke test: local icechunk layout uploads to S3 intact + (32/32 arrays identical, anonymous read OK). +2. ASE synced down; baseline sha256 snapshot of all 28 arrays + attrs. +3. In-place parallel backfill of all 815 frames (statuses: ~678 ok, + ~120 insufficient_overlap, ~64 images_unavailable, 0 errors). Two fixes + found: (a) backfill worker needed the runner's per-frame fsspec cache + cleanup (8 workers leaked ~53 GB and got the run killed); (b) resumability + via commit-every-100 worked exactly as designed. +4. Saturation pass: 11,778 finite margins; season fits incl. 2018 DC8 img1 + ceiling −24.0 dB, cross-cap step −4.1 dB. +5. Verification: 28/28 original arrays byte-identical; only the 4 calibration + arrays added. Uploaded to s3://opr-radar-metrics/icechunk/ase (old store + backed up locally at outputs/ase_s3_backup_20260831 and partially at + icechunk/ase-backup-20260831); S3 copy verified identical, opens anonymous. + +Operational notes for the big runs: +- This machine's wifi drops every ~2 min (systemd wifi-rescue service); + background tasks doing long S3 transfers get killed. Fix: detached + (setsid nohup) retry-loop wrappers (scratchpad resilient_sync.sh) that + re-run `aws s3 sync` until two consecutive clean zero-remaining passes. +- `aws s3 ls/rm s3://…/ase` without trailing slash prefix-matches + `ase-backup-*` too — always use the trailing slash for rm. + +## Greenland + Antarctica local calibration runs (2026-09-01, complete) + +Both pipelines: VERIFY_PASS — 28/28 original arrays byte-identical, only the 4 +calibration fields added; zero load_error across 11,845 frames. + +Statuses: GL 3861 ok / 1762 insufficient_overlap / 55 no_combine / 20 +images_unavailable; Ant 5178 ok / 633 insufficient_overlap / 287 +images_unavailable / 49 weight_recovery_failed (to inspect). +Note: 2013_Greenland_P3 has 663 ok frames — the 2-frame availability probe +(404s) badly underestimated that season's published img files. + +Seam-offset season means (ok frames): mostly |mean| < 2 dB, but real problem +seasons surfaced: 2016_GL +7.4±9.9 (160/272 frames >3 dB), 2017_Ant_Basler +−6.8±7.5 (210/337), 2013_Ant_Basler +4.8±6.1, 2017_Ant_P3 +3.0. Coverage: +per-trace offsets finite GL 118k/214k, Ant 201k/262k. + +Saturation (by_source_image fits): GL margins 97k/214k finite; 2014_GL both +regimes fit (img1 −31.6 / img2 −53.6, cross-cap 16.2 dB — matches predicted +bias); 2013_GL −50.6. Ant margins only 39k/262k: the source split fragments +several season populations below support — notably **2012_Antarctica_DC8** +(img1 span 0.195 dec < 0.25, 8 bins) loses its former flagship −33.1 dB / +17%-pileup detection. Decision needed: relax span again, or fall back to +all_traces fit (flagged as mixed) when both subpopulations fail support. +2018_GL img1 fit_ok(−58.2) is a weak flat envelope (slope −1.5, pileup 1.7%) — +review. Piecewise model: still no clear real-data winner (fits driven by the +flat rule); user-requested keep/drop decision now due. diff --git a/claude_notes/img2_surface_fraction.py b/claude_notes/img2_surface_fraction.py new file mode 100644 index 0000000..5f58bdc --- /dev/null +++ b/claude_notes/img2_surface_fraction.py @@ -0,0 +1,109 @@ +"""Estimate the fraction of store traces whose surface would be img2-sourced +(surface_twtt beyond the img1 gate cap T_end(img1) - T_guard). + +Rough per-season caps from 2 sampled frames' img1 gates (gates can vary by +segment/altitude, so min/max caps bound the estimate). Greenland + Antarctica +stores only (user scope). Read-only. +""" +import numpy as np +import icechunk +import zarr +from xopr import OPRConnection + + +def get_img_comb(attrs): + paths = [ + ("param_array", "array", "img_comb"), + ("param_records", "array", "img_comb"), + ("param_combine", "array_param", "img_comb"), + ("param_combine", "combine", "img_comb"), + ("param_array", "qlook", "img_comb"), + ("param_combine", "get_heights", "qlook", "img_comb"), + ] + for p in paths: + d = attrs + try: + for k in p: + d = d[k] + v = np.atleast_1d(np.asarray(d, dtype=float)) + if v.size >= 3: + return v + except (KeyError, TypeError, ValueError): + continue + return None + + +def get_tpd(attrs): + for pk in ("param_records", "param_array", "param_combine"): + try: + return np.atleast_1d(np.asarray(attrs[pk]["radar"]["wfs"]["Tpd"], dtype=float)) + except (KeyError, TypeError): + continue + return None + + +opr = OPRConnection() +print(f"{'season':32s} {'caps (us)':>16s} {'n_qc':>9s} {'frac>capmed':>11s} {'lo/hi':>13s}") + +for store_name in ["greenland", "antarctica"]: + storage = icechunk.s3_storage(bucket="opr-radar-metrics", + prefix=f"icechunk/{store_name}", + region="us-west-2", anonymous=True) + repo = icechunk.Repository.open(storage=storage) + root = zarr.open_group(repo.readonly_session(branch="main").store, mode="r") + surf_twtt = root["surface_twtt"][:] + qc = root["qc_surface_pass"][:].astype(bool) + frame_index = root["frame_index"][:] + colls = np.array(root.attrs.asdict()["frame_collections"]) + seasons = colls[frame_index] + good = qc & np.isfinite(surf_twtt) + + print(f"--- store: {store_name} (n_qc={good.sum()})") + store_affected_med = store_affected_lo = store_affected_hi = 0 + store_total = 0 + for season in np.unique(seasons): + m = good & (seasons == season) + n = int(m.sum()) + store_total += n + caps = [] + try: + frames = opr.query_frames(collections=[season], max_items=12) + # prefer two frames from different segments + segs = {} + for fid in frames.index: + segs.setdefault(fid.rsplit("_", 1)[0], fid) + picks = list(segs.values())[:2] or list(frames.index[:2]) + for fid in picks: + item = frames.loc[fid] + try: + img1 = opr.load_frame(item, data_product="CSARP_standard", + image=1, allow_unlisted_products=True) + except Exception: + continue + ic = get_img_comb(img1.attrs) + if ic is not None: + t_guard = ic[2] + else: + tpd = get_tpd(img1.attrs) + t_guard = tpd[0] if tpd is not None else 0.0 + caps.append(float(img1.twtt.values[-1]) - float(t_guard)) + except Exception as e: + print(f"{season:32s} query failed: {type(e).__name__}") + continue + if not caps: + print(f"{season:32s} {'no img1 files':>16s} {n:9d} {'unknown':>11s}") + continue + cap_med, cap_lo, cap_hi = np.median(caps), min(caps), max(caps) + st = surf_twtt[m] + f_med = float((st > cap_med).mean()) + f_lo = float((st > cap_hi).mean()) # highest cap -> fewest affected + f_hi = float((st > cap_lo).mean()) + store_affected_med += f_med * n + store_affected_lo += f_lo * n + store_affected_hi += f_hi * n + print(f"{season:32s} {cap_med*1e6:7.2f} [{cap_lo*1e6:.2f},{cap_hi*1e6:.2f}]" + f" {n:9d} {f_med:11.1%} {f_lo:6.1%}/{f_hi:6.1%}") + if store_total: + print(f"{'STORE TOTAL':32s} {'':16s} {store_total:9d} " + f"{store_affected_med/store_total:11.1%} " + f"{store_affected_lo/store_total:6.1%}/{store_affected_hi/store_total:6.1%}") diff --git a/claude_notes/plot_2012_dc8_saturation.py b/claude_notes/plot_2012_dc8_saturation.py new file mode 100644 index 0000000..10365a5 --- /dev/null +++ b/claude_notes/plot_2012_dc8_saturation.py @@ -0,0 +1,86 @@ +"""Visualize the 2012_Antarctica_DC8 saturation situation after the +by-source-image split (from the calibrated local antarctica store).""" +import numpy as np +import icechunk +import zarr +import scipy.constants +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +storage = icechunk.local_filesystem_storage("outputs/icechunk_store_antarctica_local") +root = zarr.open_group(icechunk.Repository.open(storage=storage) + .readonly_session(branch="main").store, mode="r") + +colls = np.array(root.attrs.asdict()["frame_collections"]) +seasons = colls[root["frame_index"][:]] +m = (seasons == "2012_Antarctica_DC8") +qc = root["qc_surface_pass"][:].astype(bool) +pw = root["surface_power_dB"][:] +r = scipy.constants.c * root["surface_twtt"][:] / 2.0 +src = root["surface_source_image_index"][:] +good = m & qc & np.isfinite(pw) & np.isfinite(r) & (r > 100) + +pops = { + "img1-sourced": good & (src == 1), + "img2-sourced": good & (src >= 2), + "unknown (-1)": good & (src < 0), +} +print({k: int(v.sum()) for k, v in pops.items()}) + +fig, axes = plt.subplots(1, 3, figsize=(20, 6), sharey=True) + +def pct_overlay(ax, rr, p, color): + bins = np.logspace(np.log10(rr.min()), np.log10(rr.max()), 20) + bidx = np.digitize(rr, bins) + for b in range(1, len(bins)): + sel = bidx == b + if sel.sum() > 50: + ax.plot(np.sqrt(bins[b - 1] * bins[b]), np.percentile(p[sel], 99), + ".", color=color, ms=10, mec="k", mew=0.5) + +# Panel 1: all traces (the pre-split view that gave fit_ok -33.1) +ax = axes[0] +rr, p = r[good], pw[good] +ax.hexbin(rr, p, gridsize=70, bins="log", cmap="viridis", xscale="log") +pct_overlay(ax, rr, p, "w") +ax.axhline(-33.1, color="r", ls="--", lw=1.5, label="old all-traces ceiling −33.1 dB") +ax.set_title(f"All QC traces (n={good.sum()})\npre-split fit: fit_ok(−33.1), pileup 17%") +ax.legend(loc="lower left", fontsize=9) + +# Panels 2-3: split populations +for ax, (label, mask), color in zip(axes[1:], list(pops.items())[:2], ["c", "orange"]): + rr, p = r[mask], pw[mask] + if mask.sum() < 10: + ax.set_title(f"{label}: n={mask.sum()}") + continue + ax.hexbin(rr, p, gridsize=70, bins="log", cmap="viridis", xscale="log") + pct_overlay(ax, rr, p, color) + span = np.log10(np.percentile(rr, 99)) - np.log10(np.percentile(rr, 1)) + ax.axhline(-33.1, color="r", ls="--", lw=1) + ax.set_title(f"{label} (n={mask.sum()})\nrange span {span:.2f} decades " + f"({'FAILS' if span < 0.25 else 'passes'} 0.25 min)") + +for ax in axes: + ax.set_xlabel("range to surface (m)") +axes[0].set_ylabel("surface power (dB)") +fig.suptitle("2012_Antarctica_DC8: saturation populations after surface-source split", y=1.02) +fig.tight_layout() +out = "outputs/figures/sat_2012dc8_split.png" +fig.savefig(out, dpi=110, bbox_inches="tight") +print("saved", out) + +# Range histogram per population (why the span collapsed) +fig2, ax2 = plt.subplots(figsize=(10, 4)) +for (label, mask), color in zip(list(pops.items())[:2], ["tab:blue", "tab:orange"]): + if mask.sum() > 10: + ax2.hist(np.log10(r[mask]), bins=60, alpha=0.6, label=f"{label} (n={mask.sum()})", color=color) +ax2.axvline(np.log10(1200), color="k", ls=":", lw=1) +ax2.set_xlabel("log10(range to surface [m])") +ax2.set_ylabel("traces") +ax2.legend() +ax2.set_title("2012_Antarctica_DC8 range distribution by surface source image") +fig2.tight_layout() +out2 = "outputs/figures/sat_2012dc8_range_hist.png" +fig2.savefig(out2, dpi=110) +print("saved", out2) diff --git a/claude_notes/plot_binning_flips.py b/claude_notes/plot_binning_flips.py new file mode 100644 index 0000000..ee5b74b --- /dev/null +++ b/claude_notes/plot_binning_flips.py @@ -0,0 +1,59 @@ +"""Inspect season fits that changed under quantile binning (v0.3.0).""" +import sys +import numpy as np +import icechunk +import zarr +import scipy.constants +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +sys.path.insert(0, "src") +from radar_return_statistics.calibration import _binned_upper_quantile, SAT_DEFAULTS, fit_ceiling + +CASES = [ + # (store, season, src_sel, label: old -> new) + ("greenland", "2014_Greenland_P3", 1, "img1: fit_ok(-31.6) -> insufficient_support"), + ("greenland", "2018_Greenland_P3", 2, "img2: insufficient_support -> fit_ok(-73.4)"), + ("antarctica", "2012_Antarctica_DC8", 1, "img1: insufficient_support -> fit_ok(-32.1)"), + ("antarctica", "2016_Antarctica_DC8", 1, "img1: no_plateau -> fit_ok(-26.8)"), + ("antarctica", "2014_Antarctica_DC8", 1, "img1: no_plateau -> fit_ok(-25.3)"), + ("antarctica", "2022_Antarctica_BaslerMKB", 1, "img1: no_plateau -> fit_ok(-42.7)"), +] + +roots = {} +for store in {c[0] for c in CASES}: + storage = icechunk.local_filesystem_storage(f"outputs/icechunk_store_{store}_local") + roots[store] = zarr.open_group(icechunk.Repository.open(storage=storage) + .readonly_session(branch="main").store, mode="r") + +fig, axes = plt.subplots(2, 3, figsize=(21, 11)) +for ax, (store, season, src_sel, label) in zip(axes.flat, CASES): + root = roots[store] + a = root.attrs.asdict() + seasons = np.array(a["frame_collections"])[root["frame_index"][:]] + qc = root["qc_surface_pass"][:].astype(bool) + pw = root["surface_power_dB"][:] + r = scipy.constants.c * root["surface_twtt"][:] / 2.0 + src = root["surface_source_image_index"][:] + m = (seasons == season) & qc & np.isfinite(pw) & np.isfinite(r) & (r > 100) + m &= (src == src_sel) if src_sel == 1 else (src >= 2) + rr, p = r[m], pw[m] + ax.hexbin(rr, p, gridsize=60, bins="log", cmap="viridis", xscale="log") + xs, qs, ns = _binned_upper_quantile(p, rr, SAT_DEFAULTS) + ax.plot(10 ** xs, qs, "wo-", ms=6, mec="k", lw=1, label="equal-count bin q99") + fit = fit_ceiling(p, rr) + if np.isfinite(fit["level"]): + ax.axhline(fit["level"], color="r", ls="--", + label=f"level {fit['level']:.1f} (pileup {fit['pileup_fraction']:.1%})") + ax.set_title(f"{season} src={src_sel} (n={m.sum()})\n{label}\n" + f"status={fit['status']} span={fit['span_decades']:.2f} " + f"slope={fit['single_slope'] if fit['single_slope'] is not None else float('nan'):.1f}", + fontsize=9) + ax.legend(fontsize=8, loc="lower left") + ax.set_xlabel("range (m)") + ax.set_ylabel("surface power (dB)") +fig.tight_layout() +out = "outputs/figures/binning_flips.png" +fig.savefig(out, dpi=100, bbox_inches="tight") +print("saved", out) diff --git a/claude_notes/radiometric_calibration_experiments.md b/claude_notes/radiometric_calibration_experiments.md new file mode 100644 index 0000000..00ef07a --- /dev/null +++ b/claude_notes/radiometric_calibration_experiments.md @@ -0,0 +1,124 @@ +# Radiometric calibration check experiments (2026-08-31) + +Exploration supporting `claude_plans/20260831-radiometric-calibration-checks.md`. +Scripts: `calib_explore1.py` … `calib_explore5.py`, `calib_saturation_explore.py`. + +## Loading individual images via xopr + +`opr.load_frame(item, data_product="CSARP_standard", image=N, allow_unlisted_products=True)` +works (xopr 0.5.1.dev27). Image files must exist on the public OPR site +(`Data_img_NN_...mat` next to the combined file). + +Availability spot check (first frame per collection): + +| Collection | img files | img_comb vector | notes | +|---|---|---|---| +| 2016_Antarctica_DC8 | yes (3) | `[3e-6, -inf, 1e-6, 1e-5, -inf, 3e-6]` in `param_combine.{array_param,combine}.img_comb` | Tpd = [1, 3, 10] µs | +| 2014_Greenland_P3 | yes (≥2) | **empty** (`param_combine` present, `img_comb=[]`, 7 imgs listed) | older param format | +| 2013_Greenland_P3 | **no** (404) | empty | neither check input available | +| 2018_Greenland_P3 | yes | `param_array.array.img_comb` (also `param_records.array/qlook`) | Tpd = [1,1,3,3,10,10] µs (wf-adc pairs) | +| 2019_Antarctica_GV | yes | same as 2018 layout | | +| 2022_Antarctica_BaslerMKB | yes | same layout | img1/img2 span identical twtt — full-range images | + +So the param attr path varies by season: `param_combine.*` (2016 DC8 era), +`param_array.array.img_comb` / `param_records.array.img_comb` (newer). +`radar.wfs.Tpd` under `param_records` gives per-waveform pulse durations. + +## Combine-point semantics (OPR guide §img_comb) + +Per adjacent image pair, three numbers `[T_comb, T_blank, T_guard]`: +transition at `max(T_blank, T_comb + td_surface)` (surface NaN → 0), limited by +`T_end(img_earlier) - T_guard`. T_comb ≈ Tpd of the later (higher-gain) image to +skip its saturated surface response; T_guard ≈ Tpd of earlier image to skip +pulse-compression rolloff. `img_comb_weights_mode='auto'` (if set) estimates and +applies a per-transition weight — would *hide* gain mismatch in the combined +product; not seen in the frames checked. + +## Inter-image overlap offset (healthy frame, Data_20161014_03_001) + +Per-trace median dB offset in valid overlap window +(`[max(T_blank, T_comb+td_surf), T_end_a - T_guard]`, linear interp onto common grid): + +- img1 vs img2: median +1.76 dB, std 0.22 dB +- img2 vs img3: median +1.76 dB, std 0.31 dB + +Very precise per trace (~0.2–0.3 dB). Note nonzero systematic offset even on a +healthy frame — thresholds must be relative to season baselines, not zero. +A-scope figure: `outputs/figures/img_combine_ascope.png` — combined follows +img1 → img2 → img3 exactly at predicted transition times; img3 shows the +saturated-surface bump above img2 at 8–10 µs (why T_comb exists). + +## Combined-only step detector + +Naive above/below-window median difference is dominated by the power decay +trend (−3.5 dB "step" on a clean frame). Detrended version (robust linear fits +±[0.3, 1.5] µs on each side, extrapolated to boundary): + +- 2016 DC8 healthy frames: frame medians +0.5 to +2.1 dB, per-trace MAD 1.5–3.4 dB + +Usable only as a frame-level aggregate and as fallback where img files are +missing; the overlap measure is ~10x more precise. + +## Surface saturation signature (Greenland store, existing data) + +`outputs/figures/saturation_power_vs_range_greenland.png`: surface_power_dB vs +range-to-surface (c·surface_twtt/2), per season, with per-bin 99th percentile +overlay. + +- **2014_Greenland_P3**: unmistakable hard ceiling ~−27 dB; upper envelope flat + vs range over 200 m – 5 km (expected −20 dB/decade absent). +- **2013_Greenland_P3**: ceiling ~−50 dB *plus* a strongly depressed + low-altitude population (power drops to −140 dB at short range — possibly + receiver blanking or saturation-distorted pulse compression). +- 2016–2019: upper envelope roughly follows −20 dB/decade; look unsaturated. + +Saturation check can run entirely from existing store variables +(surface_power_dB, surface_twtt, qc_surface_pass, frame_index) — no echogram +loading needed. + +## Image availability across all store seasons (2026-08-31 probe) + +`calib_img_availability_probe.py`: HEAD requests on img_01 URLs for 2 frames of +each of the 19 unique seasons across the antarctica/greenland/ase/utig stores +(crosssystem store predates the `frame_collections` attr): + +- All 200 (images published): 15 seasons. +- Season-wide 404: **2013_Greenland_P3** only. +- Mixed 200/404 (one of two frames): 2013_Antarctica_Basler, + 2014_Greenland_P3, 2019_Greenland_P3. Some 404s may be single-image frames + (per the OPR guide, a lone image is written straight to the combined file + with no img_01) — check params to distinguish `no_combine` from + `images_unavailable`. + +## Combine weights: declared vs effective (2026-08-31) + +`calib_weights_survey.py` (one frame per season, all param_* attrs searched): +most seasons declare no `img_comb_weights*` fields; several declare +`img_comb_weights=[0 0]` with empty mode; **`img_comb_weights_mode='auto'`** +declared on the array/combine path only for **2018_Antarctica_DC8** and +**2017_Antarctica_Basler** (2019_Antarctica_GV: qlook paths only). + +`calib_effective_weights.py`: effective per-image scaling is recoverable by +differencing the combined product against each image over the section the +combined sourced from it (`w_i = median(combined − img_i)`, per trace, then +median over traces). It is a per-frame constant with **zero** measured spread: + +- 2016_Antarctica_DC8: w1=w2=w3=+0.000 dB exactly — combined is the raw + stitch; the +1.76 dB raw inter-image offset is a *real* seam step in the + combined product. +- 2018_Antarctica_DC8 (auto): w1=0, **w2=+107.607 dB, w3=+114.840 dB** (MAD + 0.000). The raw img files there are not on a common gain reference at all + (auto weights and/or undeclared unit differences), so raw img-vs-img + offsets are meaningless without weight recovery. + +`calib_residual_step.py`: residual seam step = overlap offset of +weight-corrected images `(dB_a + w_a) − (dB_b + w_b)` = the step actually in +the combined product: + +- 2016 DC8: +1.764 dB (MAD 0.22) and +1.757 dB (MAD 0.31) — equals raw, as + weights are zero. +- 2018 DC8: +0.971 dB (MAD 0.80) and +1.364 dB (MAD 2.85) — the >100 dB + corrections bring the seams to ~1 dB residual. + +This makes the check weight-mode-independent: recover w_i empirically, report +the residual as the headline metric, keep raw offsets + w_i as diagnostics. diff --git a/claude_notes/review_fixes_diff.md b/claude_notes/review_fixes_diff.md new file mode 100644 index 0000000..56bb77a --- /dev/null +++ b/claude_notes/review_fixes_diff.md @@ -0,0 +1,72 @@ +# PR-review fixes: local re-run diff (2026-09-02, method 0.4.1) + +Seven review fixes implemented (see PR discussion), then the calibration +machinery re-run against the three local stores: +`calibration_backfill --statuses weight_recovery_failed` on antarctica +(49 frames), `saturation_pass` on all three, calibration-variable attrs +stamped. Snapshots via `claude_notes/calib_state_snapshot.py` +(before/after in the session scratchpad). + +## Headline + +**The published *measurements* are unchanged.** Across all three stores: +zero frame-status transitions, zero per-trace `img_comb_offset_dB` changes, +zero `surface_source_image_index` changes. The only data change is the +intended fix-5 effect: + +| store | change | +|---|---| +| greenland | none (margins 97,877 unchanged; method attr 0.4.0→0.4.1; array attrs added) | +| antarctica | 2018_Antarctica_DC8/img2: fit_ok(−18.8) → **no_plateau** → 1,428 margins NaN'd (91,867→90,439) | +| ase | 2018_Antarctica_DC8/img2: fit_ok(−21.2) → **no_plateau** → 651 margins NaN'd (16,199→15,548) | + +Science verification: PASS on all three stores (original arrays +byte-identical; snapshot verify). + +## Per-fix outcomes + +1. **Staleness ordering** — backfill now marks `saturation.stale` when it + rewrites `surface_source_image_index`; `saturation_pass` clears it. + Verified live: antarctica read stale=True after the 49-frame backfill, + fresh after the pass. Docs order swapped (backfill → saturation). +2. **All-or-nothing backfill writes** — any array-length mismatch now aborts + the frame update with nothing written (status untouched, frame stays + retryable). No real-data effect (no mismatches exist). +3. **Pair coverage vs expected images** — the pair loop now spans + max(expected, loaded) images (bounded by the declared img_comb length), so + a trailing missing image yields `partial_images`, never `ok`. Zero + real-data relabels in this re-run (only the 49 wrf frames were + re-measured; the sample audit predicted ≈0 store-wide). +4. **NaN-trace weight poisoning** — per-trace and aggregate medians now + filter non-finite values. **All 49 antarctica `weight_recovery_failed` + frames re-ran and all 49 remain weight_recovery_failed with identical + values** — they are genuine large-spread cases (frame-varying effective + weights, mostly 2017_Antarctica_Basler-era), not NaN artifacts. The fix is + defensive only. +5. **Two-sided flat band** (−12 ≤ slope ≤ +12 dB/decade) — rising envelopes + are no longer ceilings. Exactly one real fit affected: 2018_Antarctica_DC8 + img2 (slope +13.6), which flips to no_plateau in both stores that carry it. + All other fits identical. +6. **Single-session saturation pass** — fit and write now share one writable + session (no fit-from-one-snapshot / stamp-another race). No data effect. +7. **Metadata/docs** — CALIBRATION_VAR_ATTRS templates (units, sign + convention, provenance framing) now applied at every array-creation path + and stamped onto all three local stores (0→4 arrays with attrs each); + changelog gains `partial_images` + corrected NaN-offset claim + the + "img1 vs pooled higher-gain (index ≥ 2)" clarification; stale piecewise + module docstring fixed; METHOD_VERSION 0.4.1. + +## Tests + +81 passed (76 existing + 5 new: staleness marking/clearing, mismatch +all-or-nothing, trailing-404 → partial_images, NaN-trace weights, rising +envelope rejection). One existing test edited: the mismatch test's no-op +session cannot be committed (icechunk refuses empty commits), asserted +explicitly. + +## Not yet done (post-review) + +- Commits + push of the fixes; S3 stores still carry method 0.4.0 and + attr-less arrays until the small syncs are re-run after approval. +- The changelog rollout table will need its method reference refreshed when + the S3 sync happens. diff --git a/claude_notes/two_regime_investigation.md b/claude_notes/two_regime_investigation.md new file mode 100644 index 0000000..c99df60 --- /dev/null +++ b/claude_notes/two_regime_investigation.md @@ -0,0 +1,82 @@ +# Two-regime surface-power investigation (2026-08-31) + +Question: why do several Greenland P3 seasons show two flat-topped +surface-power bands with a ~15–20 dB step near ~800 m–1 km range +(baseline finding, `calibration_baseline_results.md`)? Candidates: +(A) image combining, (B) receiver/params changes. +Scripts: `two_regime_step{1..4}_*.py`; figures in +`outputs/calibration/baseline/figures/` (`two_regime_gl_seasons.png`, +`two_regime_ascope_2014.png`). + +## Verdict: Hypothesis A — the seam cap forces the combined product's surface +## sample to come from img2, whose surface response is corrupted + +**Evidence (2014_Greenland_P3, 5 frames across 5 segments, incl. 3 mixed- +altitude frames):** + +- Per-trace `surface_source_image_index` (from `calibration.check_img_combine`) + splits every mixed frame exactly at 780 m — precisely + `T_end(img1) − T_guard = 7.85 − 2.64 ≈ 5.21 µs`. Regime membership is + predicted perfectly by source image *within single frames*: + + | frame | src=1 (≤780 m) p50 | src=2 (>780 m) p50 | + |---|---|---| + | Data_20140502_01_031 | −40 to −47 dB range | −58 to −65 dB range | + | Data_20140505_01_002 | −42.6 | −60.0 | + | Data_20140412_03_001 | −40.7 | −64.9 | + | Data_20140407_01_012 | −46.7 | −64.8 | + | Data_20140415_01_002 (high) | — | −58.4 | + +- **Hypothesis B eliminated**: all gain-affecting params are identical across + all 5 frames / 5 segments spanning both regimes (adc_gains = 22.387 on every + channel, presums, rx_paths, Vpp_scale all constant). Same settings on both + sides of the step within a single frame → settings cannot explain it. + +**Sign explanation** (`two_regime_ascope_2014.png`): at 700 m the surface +(4.65 µs) precedes the seam cap (5.21 µs), so the combined follows img1 — +sharp compressed surface peak (−53 dB in the example trace). At 900 m the +surface (6.0 µs) lands *past* the cap, so the combined takes img2's surface — +a broadened, suppressed hump ~15–20 dB lower. This is the corrupted +surface response `T_comb` exists to avoid ("avoid using the saturated surface +response in image 2", OPR guide): the seam cap `T_end(img1) − T_guard` +overrides that protection whenever td_surf > cap. The lower band is therefore +**not a second calibration regime — it is invalid surface-power data**: for +`surface_source_image_index ≥ 2` traces the surface sample is smeared/ +suppressed, biasing surface power low by ~15–20 dB. + +## 2017 / 2019 Greenland + +- **2017_Greenland_P3: same phenomenon confirmed.** Mixed frame + Data_20170424_01_061 splits at 1050 m (= its era's cap + `T_end(img1) − T_guard ≈ 7 µs`): src=1 p50 −21.3 dB vs src=2 p50 −40.5 dB + (19 dB step at the boundary). +- **2019_Greenland_P3: inconclusive from one frame.** The sampled mixed frame + (Data_20190417_02_047) has an img1 gate long enough to cover all its ranges + (all src=1 to 1756 m), so its population step at ~1 km must come from + segments with shorter img1 gates (gates vary per segment) or another cause. + The per-trace index at integration scale will resolve it season-wide. + +## Implications for the saturation second pass (prototype results) + +Store-level prototype on the full 2014 population (img1/img2 approximated by +the 780 m cap): + +| population | n | fit | +|---|---|---| +| full season (current) | 49 251 | `no_plateau` (single slope −24 averages the regimes) | +| img1-sourced (r < 731 m) | 44 904 | `insufficient_support` — span 0.31 < 0.5 decades; with `min_span_decades=0.25` → **fit_ok, level −31.5 dB**, flat (−4.8 dB/dec) | +| img2-sourced (r > 831 m) | 3 953 | fit_ok at −53.7 dB — but this is a ceiling on *corrupted* data, not a clip level | + +1. The already-planned restriction to `surface_source_image_index == 1` + cleanly isolates the valid population and recovers a credible plateau for + 2014_Greenland_P3 — **but only if `min_span_decades` is relaxed (0.5 → + ~0.25) for that subpopulation**, because low-altitude flying inherently + spans little range. (Interacts with the user's decision to leave the span + minimum as-is for 2019_GV — needs a user call.) +2. Traces with index ≥ 2 should never enter a ceiling fit, and their stored + `surface_power_dB` is biased low by ~15–20 dB — this affects *all* + surface-power uses (e.g. RSSNR), not just saturation. The planned per-trace + index variable is exactly the needed filter; whether to also flag these in + docs/QC guidance is a user decision. +3. No bug found in `calibration.py` — `surface_source_image_index` behaved + exactly as designed throughout. diff --git a/claude_plans/20260831-radiometric-calibration-checks.md b/claude_plans/20260831-radiometric-calibration-checks.md new file mode 100644 index 0000000..055a47b --- /dev/null +++ b/claude_plans/20260831-radiometric-calibration-checks.md @@ -0,0 +1,410 @@ +# Radiometric calibration checks: image-combine errors and surface saturation + +Status: **COMPLETE 2026-09-01.** All checks implemented, integrated +(method 0.4.0: flat rule only, sparse-bin merging, piecewise dropped by user +decision), backfilled across ase/greenland/antarctica, verified (science +arrays byte-identical; only calibration fields added), and pushed to S3. +PR: https://github.com/englacial/radar-return-statistics/pull/4. +Dataset-facing summary: `docs/dataset_changelog.md`. Working notes: +`claude_notes/calibration_baseline_results.md`, +`claude_notes/two_regime_investigation.md`, +`claude_notes/calibration_integration_notes.md`. + +## User decisions from baseline review (2026-08-31) + +1. **Two-regime seasons — RESOLVED** (`claude_notes/two_regime_investigation.md`): + the second "regime" is surface-from-img2 traces forced in by the img1 gate + cap (`T_end(img1) − T_guard`); params ruled out (identical settings both + sides of the step within single frames). User decisions on the follow-ups + (2026-08-31): + - `surface_source_image_index` is a **provenance flag, not a validity + verdict**: img2-sourced surfaces are *more likely saturated* and may + carry a season-dependent low bias (~15–20 dB measured on 2014/2017 GL + P3's img2 surface response), but for high-altitude DC8 seasons the + surface landing in img2 is the normal operating geometry — blanket + exclusion would discard ~19–27% of the 2012/2014/2016 Antarctica DC8 + seasons (`claude_notes/img2_surface_fraction.py`: greenland store ~2% + affected overall, antarctica ~9%, concentrated in DC8 seasons). Document + the likelihood/bias in docs + variable attrs; downstream filters. + - The saturation second pass fits the **img1-sourced and img2-sourced + populations separately** per season (each with its own level/status/ + support in the season dict), and reports the cross-cap step (img1 vs + img2 population offset at the cap) as the per-season empirical bias + estimate — measured, not assumed. + - Per-trace `surface_ceiling_margin_dB` is computed against the ceiling of + the trace's own source population when that population's fit is + `fit_ok`, else NaN. + - `min_span_decades` = **0.25 everywhere** (user decision at integration + review: one uniform minimum, no special-casing by population; supersedes + the earlier 0.5 whole-season minimum — 2019_GV's 0.42-decade span now + passes support as a side effect). +2. **insufficient_overlap**: keep honest NaN (no relaxed window). +3. **Piecewise model**: keep for now; **revisit after it has run on more + data** — remind the user at the next review checkpoint (post-integration + test runs / production baseline). +4. **2019_GV span minimum**: leave as is. +5. **No cross-store pooling** — no guarantee the same instrument flew both + ice sheets. **Only the antarctica and greenland stores matter for now**; + ignore ase/utig/crosssystem in calibration work. +Branch: `calibration-checks`. **No pushes to the online S3 stores** — all runs +against local stores / `outputs/` until explicitly approved. +Experiment details: `claude_notes/radiometric_calibration_experiments.md`. + +## Goal + +Detect two classes of radiometric calibration problems in OPR echograms and +**expose them as per-trace/per-frame values in the store** so downstream users +can choose their own filtering thresholds (we report values and status flags, +not pass/fail decisions): + +1. **Image-combine errors** — gain mismatch between the "waveform playlist" + images stitched together by `img_combine.m`, visible as a power step at the + per-trace combine point `max(T_blank, T_comb + td_surface)`. Reported as a + mean dB offset per frame; downstream rejection threshold likely ~3 dB. +2. **Surface saturation** — receiver/ADC clipping of the surface return, + visible as a hard maximum in the surface-power-vs-altitude relationship + instead of the expected ~−20 dB/decade range dependence. + +The img-combine check runs **on by default** in the pipeline (config flag to +disable): the dataset is processed once and used many times, so the extra +download cost per frame is acceptable. Cost is benchmarked before any +production run (see step 5). + +## Key facts from experiments + +- Individual images load via + `opr.load_frame(item, data_product=..., image=N, allow_unlisted_products=True)`. +- Availability probe over all 19 store seasons (HEAD on img_01 URLs, + `claude_notes/calib_img_availability_probe.py`): images published for every + season except **2013_Greenland_P3** (season-wide 404); scattered per-frame + 404s in 2013_Antarctica_Basler, 2014_Greenland_P3, 2019_Greenland_P3 (some + may be single-image frames where no combining occurred — distinguishable + from params). Frames without images get status `images_unavailable` and NaN + offsets rather than an estimate from the combined product alone. +- `img_comb` vector location varies by season: + `param_array.array.img_comb`, `param_records.array.img_comb` (newer), + `param_combine.{array_param,combine}.img_comb` (2016-era), qlook variants; + sometimes present but empty (2014_Greenland_P3). `param_records.radar.wfs.Tpd` + gives per-waveform pulse durations. +- Direct img-vs-img overlap offset is precise (~0.2–0.3 dB per-trace std). A + healthy 2016 DC8 frame still shows a systematic ~+1.8 dB offset — not + surprising per user (combine thresholds are hard to get exactly right); + reinforces reporting values rather than baking in a threshold. +- Saturation is already visible in the existing stores: 2014_Greenland_P3 has a + flat hard ceiling (~−27 dB) across all ranges; 2013 has a ceiling plus a + depressed low-altitude population. 2016–2019 look clean. +- **Effective combine weights are recoverable by differencing** the combined + product against each image over the section it sourced + (`claude_notes/calib_effective_weights.py` / `calib_residual_step.py`): + `w_i = median(combined − img_i)` over section i is a per-frame constant with + ~zero spread. 2016 DC8: all w=0 (combined = raw stitch, so its +1.76 dB raw + offset is a *real* seam step). 2018_Antarctica_DC8 (weights_mode='auto' on + the array path): w2=+107.6 dB, w3=+114.8 dB — the raw img files aren't even + on a common gain reference there, so raw img-vs-img offsets are meaningless + without weight recovery; the weight-corrected residual seam step is + +0.97 / +1.36 dB. +- `img_comb_weights_mode='auto'` declared on the relevant (array/combine) path + in 2 of 19 seasons (2018_Antarctica_DC8, 2017_Antarctica_Basler; 2019_GV + qlook-only); most seasons declare no weights fields at all + (`claude_notes/calib_weights_survey.py`). Because weights are recovered + empirically, the declared mode is recorded for provenance only. + +## Design + +New module `src/radar_return_statistics/calibration.py`. Two integration +stages: (1) a standalone CLI for baseline runs writing to +`outputs/calibration/`, (2) hooks in `processing.py` adding store variables +(local test stores only until go-ahead). + +### Check 1: image-combine consistency (per frame, during processing) + +1. **Param extraction** — `find_img_comb(ds)` / `find_tpd(ds)` searching the + known attr paths in priority order; return the vector plus a provenance + string. +2. **Image count and waveform mapping** — the image count comes from the + params `imgs` list (primary), each entry being a wf-adc matrix; the + waveform index for image i is the wf entries of `imgs[i]` (e.g. 2016 DC8 + `imgs[0] = [[1,...],[1..6]]` → wf 1, adcs 1–6). `Tpd_img_i = + radar.wfs.Tpd[wf−1]` — **never index Tpd by image number**: Tpd arrays can + carry per-wf-adc duplicates (2018 P3: `[1,1,3,3,10,10] µs` for images with + wf pairs). Warn if the wfs of one image disagree on Tpd. Image-count + precedence: `imgs` list (primary); else `len(img_comb)/3 + 1` when a valid + `img_comb` vector exists; else a bounded probe capped at **3 images** (per + user: 3 is the normal maximum, some seasons use 2, none known to use + more — the "7 imgs" seen in 2014-era params is the old per-wf-adc format, + not 7 real images, and is treated as suspect → verify at baseline). Typed + error handling throughout: HTTP 404 = image permanently absent; + timeouts/5xx = transient → bounded retry, then status `load_error` + (retryable later), so a transient failure or a missing intermediate file + is never recorded as "no more images". +3. **Effective weight recovery** — per image i, compute the per-trace + `median(combined − img_i)` over the section of the combined product + sourced from img_i (between the adjacent per-trace combine boundaries, + trimmed 0.3 µs clear of the blend zone); the frame's effective weight + `w_i` is the median over traces. Empirically this is a per-frame constant + with ~zero spread (see key facts), and it captures *whatever* scaling the + combine applied — declared static weights, 'auto' weights, or undeclared + unit/gain differences between the img files and the combined product + (2018 DC8: >100 dB). Sanity check: if the within-section spread + (MAD over traces) exceeds ~0.1 dB, weight recovery is unreliable → status + `weight_recovery_failed`, raw offsets reported to parquet only. + `img_comb_weights` / `img_comb_weights_mode` from params are recorded for + provenance but are not needed for the measurement. +4. **Metric — residual seam offset**, per adjacent image pair (a, b): the + overlap offset of the *weight-corrected* images, + `(dB_a + w_a) − (dB_b + w_b)`, which equals the step actually present in + the combined product at the seam — the quantity downstream users care + about, valid regardless of weights mode. The raw img-vs-img offset is + derivable (`residual − (w_a − w_b)`) and goes to the parquet as a + diagnostic of the underlying pre-combine calibration. Details: + - Alignment: images are subset to the same decimated slow_time selection + as the combined frame, matched by index (verify identical trace counts; + else nearest-time with the pipeline's 1 s tolerance, unmatched → NaN). + Data orientation normalized to (twtt, trace) as elsewhere in + `processing.py`; the two images keep their own twtt grids and image a is + interpolated onto image b's grid (linear power). + - Valid overlap per trace: `[max(T_blank, T_comb + td_surf), T_end_a − T_guard]`. + If `img_comb` is empty/missing but multiple images exist, use + `[td_surf + Tpd_b, T_end_a − Tpd_a]` (what T_comb/T_guard default to). + `td_surf` is the frame's `Surface` variable (the value `img_combine` + used), NaN → 0 per the OPR convention. + - Per-trace `offset = median((dB_a + w_a) − (dB_b + w_b))` over bins where + both images are >6 dB above their own noise floor. The per-image noise + floor uses the pipeline's record-tail convention — window + `[end − 12 µs, end − 7 µs]` per image, matching `record_tail` defaults — + further clipped to end before `T_end − T_guard`, since images show the + same end-of-record rolloff the tail-window study found in combined + products (a last-2µs window would be biased low). If the image is too + short for that window, fall back to its global 10th percentile; + validate the fallback during the baseline run. The weight shifts signal + and floor equally, so the SNR gate is weight-independent. Require ≥10 + valid bins per trace (else NaN) and ≥30 valid traces per pair (else + pair status `insufficient_overlap`). + - Sign convention (documented in variable attrs): positive = the earlier + (lower-index, shallower) image is brighter than the later one. +5. **Computed metrics**, per frame and pair: + - *per-trace offset* — as above (median over bins within the trace). + - *mean offset* — mean of the per-trace offsets along the frame (the + headline value; downstream threshold ~3 dB). Mean is a user decision; + each per-trace value is already a median over bins, which absorbs + within-trace outliers. The baseline run reports mean vs median per frame + as a sensitivity check — if heavy-tailed frames make them diverge enough + to move worst-pair selection or cross a ~3 dB threshold, revisit (MAD + flags such frames regardless). + - *offset MAD* — median absolute deviation (×1.4826) of the per-trace + offsets: a robust spread / measurement-quality indicator, distinct from + the mean (mean +1 dB, MAD 0.3 dB = clean and stable; mean +1 dB, MAD + 4 dB = untrustworthy measurement). + - *drift* — difference between the mean offsets of the two frame halves. + - *worst pair* — the pair with the largest |mean offset|, chosen **once + per frame**; the store carries that single pair's per-trace series plus + its pair index, so frame-level MAD/drift for the reported pair are + exactly recomputable from the store (no per-trace pair switching, no + max-selection bias). All pairs' full tables go to the baseline parquet. +6. **Status** (per frame; per-pair statuses in parquet): + `ok` / `no_combine` (single image or never combined — check not + applicable) / `images_unavailable` (404 on required images) / + `partial_images` (some pairs measurable, some not) / `load_error` + (transient failure, retryable) / `invalid_params` (params present but + inconsistent) / `params_missing` / `insufficient_overlap` / + `weight_recovery_failed` (residuals unreliable; raw offsets in parquet + only) / `disabled`. + **Calibration never fails a frame**: every calibration step is wrapped so + an error degrades to a status, and the frame's science metrics are written + regardless. Because `processed_frames` blocks automatic retries, a + calibration-only backfill command re-runs the check for frames with + retryable statuses (`load_error`, `disabled`) and updates the calibration + variables in place without touching science metrics. + +### Check 2: surface saturation + +**Which image supplied the combined surface sample matters.** The combined +product takes the surface from img1 (low-gain) unless the surface falls +outside img1's time gate — the transition is capped at `T_end(img1) − +T_guard`, so at high AGL the surface sample comes from a higher-gain image +and is *expected* to saturate; that's precisely why the low-gain image +exists. These traces must not drive the ceiling fit. + +Per-trace, during processing (images/params already loaded for check 1): +- `surface_source_image_index` — the image the combined product's surface + sample came from: smallest i with `T_start(img_i) ≤ td_surf ≤ T_end(img_i) + − T_guard_i`, evaluated with the frame's original `Surface` value (what + `img_combine` used — not the later peak-refined `surface_twtt`), NaN → 0 + per the OPR convention. 1 for single-image frames (2022/2023 BaslerMKB + images share a full-range gate → always 1). −1 when it cannot be + determined: surface outside every gate, or images unavailable — the gate + ends are properties of the image files, not the params, so frames without + image files (2013_Greenland_P3) get −1 unless single-waveform. +- Optional (evaluate during implementation): `surface_peak_width_us` — width + of the img1 surface peak (e.g. −6 dB width around the pick). Clipping + broadens the compressed pulse and raises sidelobes, giving direct per-trace + saturation evidence from the radargram. Cheap since images are already + loaded; validate against the 2014_Greenland_P3 ceiling population before + adding to the schema. + +**Ceiling detection stays on the output products, not the radargrams**: a hard +max is a property of the surface-power-vs-range *population* across a season — +no single frame spans enough range diversity to reveal it — so it is a second +pass over a completed store (cheap, no echogram loads; runs from +`surface_power_dB`, `surface_twtt`, `qc_surface_pass`, `frame_index`). + +1. Range to surface `r = c·surface_twtt/2`; where + `surface_source_image_index` is populated, fit the img1-sourced and + img2-sourced populations **separately** (per the two-regime resolution + above), each with its own level/status/support recorded in the season + dict, plus the cross-cap step (img1-vs-img2 population offset at the gate + cap) as the season's empirical img2 bias estimate. Where the index is + unavailable (pre-schema stores, 2013_Greenland_P3's −1s), fall back to a + single all-traces fit; `ceiling_fit_population` records which mode ran + (`by_source_image` / `all_traces`). +2. **Ceiling estimate per season**: binned upper-quantile (99th) power vs + `log10(r)`. + - *Plateau detection first*: a clean season has no physical ceiling, so a + season-wide max would mislabel ordinary close-range returns. The + physical model is a constant clip level: partial saturation shows as a + flat upper envelope at short range transitioning to ~−20 dB/decade + where the unsaturated envelope drops below the clip level — a single + whole-season slope would average the two regimes and miss it. Use a + **contiguous low-range plateau detector / piecewise fit** (plateau + level + breakpoint + free slope beyond), compared against a pure-slope + model; declare a plateau only when the piecewise model wins with the + plateau segment ≫ shallower than −20 dB/decade over sufficient support: + minimum traces per bin, minimum occupied bins, and minimum log-range + span (defaults to fix during implementation via the baseline runs; all + recorded as metadata, with bootstrap CI on the plateau level). The + plateau level is the season ceiling; margins against it remain valid at + all ranges (distance below the clip level). + - Season fit status: `fit_ok` / `no_plateau` / `insufficient_support`. + - *Pile-up test*: fraction of traces within 1 dB of the ceiling + (only when `fit_ok`). + - *Regime check*: during the baseline run, repeat the fit per segment and + inspect ceiling multimodality within a season (mixed receiver/waveform + settings); if a season is multimodal, fall back to per-segment ceilings + before writing per-trace margins. In that case the second pass writes a + `segment_saturation` dict (keyed by segment ID, same fields as the + season dict: level, status, support, CI) so per-trace margins remain + explainable downstream; season attrs then record + `ceiling_scope='segment'`. +3. **Per-trace value** `surface_ceiling_margin_dB` = ceiling − + surface_power_dB, written **only when the season (or segment) fit is + `fit_ok`**; NaN otherwise — no credible plateau means no margin, never a + misleading number. +4. Diagnostics: per-season hexbin figures with −20 dB/decade reference and + binned upper-quantile overlay. +5. Later/optional: investigate the 2013-style depressed low-altitude + population (blanking?). + +### Store schema (local stores only until approved) + +Per-trace variables on `slow_time` (all **always emitted**, with NaN/−1 +placeholders when disabled or unavailable, so appends stay schema-uniform): +- `img_comb_offset_dB` — float32: the frame's worst pair's per-trace signed + **residual seam offset** (weight-corrected — the step actually present in + the combined product; NaN where undefined). Attrs document the sign + convention and that the pair is fixed per frame. +- `img_comb_pair` — int8: which image pair the offset series refers to + (1 = img1/img2, 2 = img2/img3, …; −1 undefined). +- `surface_source_image_index` — int8, −1 unknown. +- `surface_ceiling_margin_dB` — float32, created at schema time (all-NaN) and + filled by the second pass. +- (optional, pending validation) `surface_peak_width_us` — float32. + +Per-frame attrs arrays (parallel to `frame_names`): +- `frame_img_comb_offset_dB` — signed mean residual offset of the worst pair + (float). +- `frame_img_comb_status` — string (statuses above). +- `frame_img_comb_weights_mode` — string ('' / 'fixed' / 'auto'), declared + params mode, provenance only (measurement is weight-mode-independent). +Recovered effective weights `w_i` and raw (uncorrected) offsets per pair go +to the baseline parquet, not the store. + +Season-level attrs from the second pass, as a dict keyed by season name (no +bare parallel lists): per season — `ceiling_dB`, `slope_dB_per_decade`, +`slope_ci`, `pileup_fraction`, `fit_status`, `ceiling_fit_population`, +`n_traces`, `n_bins`, `log_range_span`, plus global `calibration_method_version` +and every fixed threshold used. Staleness is tracked by a **science-data +fingerprint**, not a commit ID (the second pass's own write would otherwise +make its result look stale immediately): the second pass records +`(len(slow_time), sha256(sorted processed_frames))` at fit time; any science +append/removal/reprocess changes the fingerprint, while calibration-only +writes do not. The runner warns when the current fingerprint differs from the +recorded one (second pass needs re-running). + +None of these are QC-masking; they are values downstream users filter on. + +**Migration/backfill**: `_zarr_append` silently skips variables absent from an +existing store (`store.py:88`) — that silent skip is removed, not worked +around: **append auto-migrates**. When an incoming dataset carries a variable +the store lacks, `_zarr_append` creates the array sized to the store's current +`slow_time` filled with the variable's fill value (NaN / −1), logs it loudly, +then appends — so an ordinary processing run against a pre-calibration store +migrates itself and can never silently drop calibration columns. (The +standalone backfill command remains for retryable-status recomputation, not as +a migration prerequisite.) `runner.py`'s frame-attr plumbing +(`batch_frame_attrs`, currently `dict[str, dict[str, float]]` with a `float()` +cast) is extended to typed per-frame attrs so string statuses ride along. +Entry-point defaults are unified by a shared config-normalization helper used +by both `load_config` and `process_frame`: calibration defaults **on** +everywhere, so behavior cannot depend on the entry point; callers that need it +off (unit tests, quick runs) set `processing.calibration.img_combine: false` +explicitly. + +Config: `processing.calibration.img_combine: true` (default true; disable for +quick test runs) and `processing.calibration.image_load_retries`. + +## Implementation steps + +1. `calibration.py`: param extraction + wf→image Tpd mapping + image loading + (typed 404/transient handling) + effective-weight recovery + residual + overlap offset + worst-pair selection + `surface_source_image_index`, + pure and testable. Unit tests (`tests/unit/test_calibration.py`): + synthetic combined+image frames with known injected seam step and known + per-image scalings (weights must be recovered and the residual must equal + the injected step, independent of the scalings); differing twtt grids and + dimension orders; `weight_recovery_failed` on inconsistent sections; + short-image noise-floor fallback; insufficient-overlap and + transient-vs-404 statuses. +2. Saturation second-pass functions (piecewise plateau fit / pile-up / + margin) on plain numpy arrays; unit tests with synthetic clean, fully + clipped, **partially clipped** (plateau at short range, −20 dB/decade + beyond the breakpoint), mixed/insufficient-support, and multimodal + populations (clean season must yield `no_plateau` and all-NaN margins; + partial clipping must recover the plateau level). +3. Standalone CLI (`--check {img-combine,saturation}`, + `--sample-per-segment N`) writing parquet + figures to + `outputs/calibration//`. Sampling: frames grouped by OPR segment + (`YYYYMMDD_SS`), N chosen per segment with a fixed seed; a run manifest + (JSON) records method version, store snapshot, parameters, seed, selected + frames, and per-frame failures. Add `pyarrow` to `pyproject.toml` (parquet + currently has no engine dependency). +4. Baseline runs: saturation on all five stores (cheap); img-combine on a + sample of Greenland P3 + ASE DC8 + one UTIG season, including at least one + real frame per known param layout (2016-era `param_combine`, newer + `param_array`, empty-`img_comb` 2014 era). Review distributions with user; + set plateau-slope threshold; evaluate `surface_peak_width_us` against the + 2014 ceiling population; run the per-segment regime check. +5. Integrate: `processing.py` (behind config flag, default on), + `runner.py` (typed frame attrs), `store.py` (schema creation, migration + backfill, second-pass in-place updates). Benchmark download volume, + runtime, and worker temp-disk growth on a sample before default-on + production use (images multiply per-frame downloads ~×N; load pairs + sequentially and release each image after use; mind the recently fixed + /tmp cache growth). Tests: append to a pre-calibration store auto-migrates + (arrays created, backfilled with fill values, loud log), enabled→ + disabled→enabled transitions, backfill of retryable statuses, fingerprint + staleness warning (calibration-only write stays fresh; science append goes + stale). Test against a local store (`test_config.yaml`). + **Stop before any production store update.** +6. Docs & viewer (with integration): `docs/architecture.md`, + `docs/data_access.md`, `docs/processing.md`; viewer handles the int8 + categorical variables (`surface_source_image_index`, `img_comb_pair`) + as integers, not via the float typed-array path. + +## Open questions + +- Is the ~1.8 dB img-to-img offset on a healthy DC8 frame instrumental + (Tsys/deconv per waveform) or partly an interpolation artifact? Baseline run + will show whether it is season-stable. +- 2014_Greenland_P3 lists 7 images with empty `img_comb` — confirm what the + combined file actually contains for that era before trusting the + Tpd-default overlap window there. diff --git a/docs/architecture.md b/docs/architecture.md index ffdb7ec..144c9cf 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -39,6 +39,29 @@ Per-trace (resampled) values stored with `slow_time` dimension: from missingness analyses as segment-edge quirks * `bed_pick_quality` - OPR layer quality flag (1 good / 2 moderate / 3 derived), -1 where no pick or flag unavailable +* `img_comb_offset_dB` - residual seam offset of the frame's worst image pair: + the power step actually present in the combined product at the image-combine + transition, after recovering the effective per-image weights the combine + applied. Positive = the earlier (shallower) image is brighter. The pair is + fixed per frame (`img_comb_pair`); NaN where unmeasured (see + `frame_img_comb_status`). A downstream rejection threshold around 3 dB on + the frame mean is typical; no threshold is applied in the store. +* `img_comb_pair` - which image pair `img_comb_offset_dB` refers to + (1 = img1/img2, 2 = img2/img3), -1 undefined +* `surface_source_image_index` - which image the combined product's surface + sample came from (the combine transition is capped at the earlier image's + gate end minus its guard, so at high AGL the surface comes from a + higher-gain image). **Provenance flag, not a validity verdict**: index >= 2 + surfaces are *more likely to be saturated* and may carry a season-dependent + low bias (~15-20 dB measured on 2014/2017 Greenland P3's img2 surface + response), but for high-altitude DC8 seasons surface-in-img2 is the normal + operating geometry. Downstream users decide filtering; the per-season + `cross_cap_step_db` attr is the measured bias estimate. -1 unknown + (images unavailable, e.g. 2013_Greenland_P3). +* `surface_ceiling_margin_dB` - season ceiling (clip level) minus surface + power, computed against the trace's own source-image population's ceiling. + Written by the saturation second pass; NaN until it runs and wherever no + credible ceiling fit exists (a clean season has no ceiling). * `record_end_twtt` - twtt (s) of the last sample in the record; lets users reconstruct the record-end-relative windows (post-bed, record tail). Never QC-masked. @@ -50,9 +73,21 @@ Per-trace (resampled) values stored with `slow_time` dimension: Coordinates: `latitude`, `longitude`, `elevation` Root attributes `frame_names`, `frame_collections`, `frame_bed_pick_fraction`, -and `segment_bed_pick_fraction` are parallel lists (one entry per unique frame) -giving each frame's season and bed-picking effort per frame / per segment. The -per-trace `frame_index` array indexes into them. +`segment_bed_pick_fraction`, `frame_img_comb_offset_dB` (worst-pair mean +residual, None where unmeasured), `frame_img_comb_status` (`ok` / +`no_combine` / `images_unavailable` / `partial_images` / `load_error` / +`invalid_params` / `params_missing` / `insufficient_overlap` / +`weight_recovery_failed` / `disabled`), and `frame_img_comb_weights_mode` +are parallel lists (one entry per unique frame). The per-trace `frame_index` +array indexes into them. + +The saturation second pass adds a root `saturation` attr: per-season fits +keyed by season name (per population `img1` / `img2` — or `all` where the +source index is unavailable — level, status, support, CI; plus +`cross_cap_step_db`, the measured img1-vs-img2 population offset at the gate +cap), the fit parameters, a method version, and a science-data fingerprint +(trace count + hash of `processed_frames`). Science appends make the second +pass stale (the runner warns); calibration-only writes do not. ### Missing bed picks as censored observations @@ -83,11 +118,20 @@ handles versioning and incremental tracking (processed frame IDs stored in the z ### Modules -* `config.py` - loads YAML config +* `config.py` - loads YAML config (`normalize_config` applies shared defaults) * `processing.py` - per-frame metric extraction (ports `extract_layer_peak_power` algorithm) -* `store.py` - icechunk read/write, frame tracking, commits +* `calibration.py` - radiometric calibration checks: image-combine residual + seam offsets (effective-weight recovery + weight-corrected overlap) and + saturation ceiling fits (per-source-image populations) +* `store.py` - icechunk read/write, frame tracking, commits; appends + auto-migrate (missing variables are created backfilled with fill values) * `runner.py` - orchestration: query, diff, process (parallel), write, commit * `__main__.py` - CLI entry point via click +* `saturation_pass.py` - second pass: fit season ceilings over a completed + store, write `surface_ceiling_margin_dB` + season attrs +* `calibration_backfill.py` - re-run the img-combine check for frames with + retryable statuses, updating calibration variables in place +* `calibration_cli.py` - standalone read-only baseline runs (parquet + figures) ### How to run @@ -117,6 +161,21 @@ Set `AWS_PROFILE` for local development. 6. Commit with summary message; `processing.checkpoint_every: N` additionally commits every N frames so long runs are resumable +During step 4 each frame also runs the image-combine calibration check +(`processing.calibration.img_combine`, default on; ~3x wall time and ~2.5x +download volume per frame from the individual image loads). Calibration +failures degrade to a per-frame status and never fail a frame. After the +science data is complete, first backfill any retryable calibration statuses, +then run the saturation second pass — the pass fits per surface-source-image +population (img1 vs the pooled higher-gain images, index >= 2), so backfills +that rewrite `surface_source_image_index` mark existing saturation results +stale until the pass is re-run: + +```bash +uv run python -m radar_return_statistics.calibration_backfill config/config_antarctica.yaml # retry load_error/disabled frames +uv run python -m radar_return_statistics.saturation_pass config/config_antarctica.yaml +``` + ### Testing `tests/unit` covers processing, config, and store logic against synthetic diff --git a/docs/data_access.md b/docs/data_access.md index 2f924e3..708502f 100644 --- a/docs/data_access.md +++ b/docs/data_access.md @@ -28,6 +28,17 @@ represented for downstream missingness analyses. The `processed_frames` array tracks which frame IDs have been ingested, and `frame_index` maps each trace to the `frame_names` root attribute. +Stores written with the calibration schema also carry radiometric calibration +values to filter on: per-trace `img_comb_offset_dB` (residual image-combine +seam step; a ~3 dB frame-mean cut is a typical downstream rejection +threshold), `surface_source_image_index` (which image supplied the surface +sample — a provenance flag: index >= 2 surfaces are more likely saturated, +not inherently invalid), and `surface_ceiling_margin_dB` (distance below the +season's fitted clip level; NaN where no credible ceiling exists), plus +per-frame `frame_img_comb_status` / `frame_img_comb_offset_dB` attrs and the +season-level `saturation` attr with the measured img2 bias +(`cross_cap_step_db`). See architecture.md for full semantics. + ## Python access ### Using icechunk + zarr directly diff --git a/docs/data_citations.md b/docs/data_citations.md new file mode 100644 index 0000000..11b5810 --- /dev/null +++ b/docs/data_citations.md @@ -0,0 +1,153 @@ +# Radar Data Citations + +Requested citations and acknowledgments for the underlying radar data, covering +the union of all seasons in the icechunk stores (antarctica, greenland, ase, +utig, crosssystem). Compiled 2026-08-07 from the data providers' citation +pages; verify against the linked sources before publication, as providers +occasionally update grant lists and citation formats. + +Seasons covered: + +| Season(s) | Program / operator | Instrument | See section | +|---|---|---|---| +| 2013, 2014, 2016, 2017, 2018, 2019 Greenland_P3 | NASA Operation IceBridge (CReSIS) | MCoRDS | 1, 2 | +| 2012, 2014, 2016, 2018 Antarctica_DC8 | NASA Operation IceBridge (CReSIS) | MCoRDS | 1, 2 | +| 2013_Antarctica_P3, 2017_Antarctica_P3 | NASA Operation IceBridge (CReSIS) | MCoRDS | 1, 2 | +| 2019_Antarctica_GV | NASA Operation IceBridge (CReSIS) | MCoRDS | 1, 2 | +| 2013_Antarctica_Basler, 2017_Antarctica_Basler | CReSIS Basler surveys | CReSIS RDS | 1 | +| 2008_Antarctica_BaslerJKB | UTIG ICECAP | HiCARS 1 | 3 | +| 2017_Antarctica_BaslerJKB | UTIG ICECAP-2 era | HiCARS 2 | 3 | +| 2022, 2023 Antarctica_BaslerMKB | NSF COLDEX | MARFA | 4 | + +## 1. Open Polar Radar / CReSIS (all seasons) + +All data were accessed through [Open Polar Radar](https://data.cresis.ku.edu/) +(via the xopr library). OPR/CReSIS request the following. + +Data citation (replace radar name as appropriate, e.g. "RDS"): + +> CReSIS. 2024. RDS Data, Lawrence, Kansas, USA. Digital Media. +> http://data.cresis.ku.edu/. + +Requested acknowledgment for CReSIS data: + +> We acknowledge the use of data and/or data products from CReSIS generated +> with support from the University of Kansas, NASA Operation IceBridge grant +> NNX16AH54G, NSF grants ACI-1443054, OPP-1739003, and IIS-1838230, Lilly +> Endowment Incorporated, and Indiana METACyt Initiative. + +OPR toolbox (the software that produced the CSARP products): + +> Open Polar Radar. (2024). opr (Version 3.0.1) [Computer software]. +> https://gitlab.com/openpolarradar/opr/. https://doi.org/10.5281/zenodo.5683959 + +(Note: the Zenodo record behind that DOI is titled "CReSIS Toolbox version +3.01" — the opr toolbox's predecessor name; the citation format above is what +OPR's own site requests alongside this DOI.) + +Requested acknowledgment for the OPR toolbox: + +> We acknowledge the use of software from Open Polar Radar generated with +> support from the University of Kansas, NASA grants 80NSSC20K1242 and +> 80NSSC21K0753, and NSF grants OPP-2027615, OPP-2019719, OPP-1739003, +> IIS-1838230, RISE-2126503, RISE-2127606, and RISE-2126468. + +Contact for citation questions: opr@openpolarradar.org. + +## 2. NASA Operation IceBridge MCoRDS (P3 / DC8 / GV seasons) + +The Greenland P3 (2013–2019), Antarctica DC8 (2012–2018), Antarctica P3 +(2013, 2017), and Antarctica GV (2019) seasons were flown as NASA Operation +IceBridge campaigns with the CReSIS MCoRDS radar. The archival NSIDC datasets: + +L1B radar echograms ([IRMCR1B v2](https://nsidc.org/data/irmcr1b/versions/2)): + +> Paden, J., Li, J., Leuschen, C., Rodriguez-Morales, F. & Hale, R. (2014). +> IceBridge MCoRDS L1B Geolocated Radar Echo Strength Profiles, Version 2 +> (IRMCR1B). Boulder, Colorado USA. NASA National Snow and Ice Data Center +> Distributed Active Archive Center. https://doi.org/10.5067/90S1XZRBAX5N + +L2 surface/bottom/thickness ([IRMCR2 v1](https://nsidc.org/data/irmcr2/versions/1)) — +relevant since the layer picks used here derive from the same picks: + +> Paden, J., Li, J., Leuschen, C., Rodriguez-Morales, F. & Hale, R. (2010). +> IceBridge MCoRDS L2 Ice Thickness, Version 1 (IRMCR2). Boulder, Colorado +> USA. NASA National Snow and Ice Data Center Distributed Active Archive +> Center. https://doi.org/10.5067/GDQ0CUCVTE2Q + +Also acknowledge NASA Operation IceBridge for these campaigns. The +2013/2017 Antarctica_Basler seasons are CReSIS-operated surveys distributed +through OPR; they are covered by the CReSIS citation/acknowledgment in +section 1. + +## 3. UTIG ICECAP / HiCARS (BaslerJKB seasons) + +2008_Antarctica_BaslerJKB (HiCARS 1 era) and 2017_Antarctica_BaslerJKB +(HiCARS 2 era; 3 frames from January 2018) were collected by the University +of Texas Institute for Geophysics ICECAP program. + +HiCARS 1 ([IR1HI1B v1](https://nsidc.org/data/ir1hi1b/versions/1)): + +> Blankenship, D. D., Kempf, S. D., Young, D. A., Richter, T. G., Schroeder, +> D. M., Greenbaum, J. S., van Ommen, T., Warner, R. C., Roberts, J. L., +> Young, N. W., Lemeur, E., Siegert, M. J. & Holt, J. W. (2017). IceBridge +> HiCARS 1 L1B Time-Tagged Echo Strength Profiles, Version 1 (IR1HI1B). +> https://doi.org/10.5067/W2KXX0MYNJ9G + +HiCARS 2 ([IR2HI1B v1](https://nsidc.org/data/ir2hi1b/versions/1)); note the +January 2018 flights postdate the NSIDC archive's coverage and are +distributed via OPR: + +> Blankenship, D. D., Kempf, S. D., Young, D. A., Richter, T. G., Schroeder, +> D. M., Ng, G., Greenbaum, J. S., van Ommen, T., Warner, R. C., Roberts, +> J. L., Young, N. W., Lemeur, E. & Siegert, M. J. (2017). IceBridge HiCARS 2 +> L1B Time-Tagged Echo Strength Profiles, Version 1 (IR2HI1B). +> https://doi.org/10.5067/0I7PFBVQOGO5 + +ICECAP acknowledgment (per the NSIDC dataset documentation): data were +collected by the International Collaborative Exploration of the Cryosphere +through Airborne Profiling (ICECAP) project, funded by NSF, the Antarctic +Climate and Ecosystems Collaborative Research Centre, and NERC, with +additional support from NASA Operation IceBridge. + +## 4. NSF COLDEX / MARFA (BaslerMKB seasons) + +2022_Antarctica_BaslerMKB and 2023_Antarctica_BaslerMKB are the NSF Center +for Oldest Ice Exploration (COLDEX) airborne seasons (CXA1 2022/23, CXA2 +2023/24) with the UTIG MARFA radar. Citations from the +[COLDEX published datasets page](https://coldex.org/published-datasets): + +Raw radar data (USAP Data Center): + +> Young, D. A., Blankenship, D. D., Buhl, D., Chan, K., Greenbaum, J., +> Kempf, S. D., et al. (2024). NSF COLDEX Raw MARFA Ice Penetrating Radar +> data. U.S. Antarctic Program (USAP) Data Center. +> https://doi.org/10.15784/601768 + +OPR-processed radargrams (Texas Data Repository): + +> Young, D. A., Paden, J., Greenbaum, J., Blankenship, D. D., Kerr, M., +> Singh, S., Kaundinya, S., Chan, K., Chan, R., Ng, G., & Kempf, S. D. +> (2024). COLDEX Open Polar Radar MARFA Airborne Radar Data. Texas Data +> Repository. https://doi.org/10.18738/T8/J38CO5 + +(Related: COLDEX VHF MARFA Open Polar Radar radargrams, +https://doi.org/10.18738/T8/NEF2XM.) + +Acknowledge the NSF Science and Technology Center award supporting COLDEX +(NSF OPP-2019719, which also appears in the OPR grant list in section 1). + +## 5. Optional instrument references + +Commonly cited instrument papers for the radars involved: + +- MCoRDS: Rodriguez-Morales, F., et al. (2014). Advanced Multifrequency Radar + Instrumentation for Polar Research. IEEE Transactions on Geoscience and + Remote Sensing, 52(5), 2824–2842. https://doi.org/10.1109/TGRS.2013.2266415 +- HiCARS: Peters, M. E., Blankenship, D. D., & Morse, D. L. (2005). Analysis + techniques for coherent airborne radar sounding: Application to West + Antarctic ice streams. Journal of Geophysical Research, 110, B06303. + https://doi.org/10.1029/2004JB003222 + +All DOIs and links in this document were verified against their registrars' +metadata (DataCite/Crossref content negotiation) and live pages on 2026-08-07. diff --git a/docs/dataset_changelog.md b/docs/dataset_changelog.md new file mode 100644 index 0000000..62b34e4 --- /dev/null +++ b/docs/dataset_changelog.md @@ -0,0 +1,130 @@ +# Dataset changelog + +## 2026-09 — Radiometric calibration fields (calibration method 0.4.1) + +New per-trace variables and per-frame/season metadata for two radiometric +calibration checks, added to the **ase**, **greenland**, and **antarctica** +stores. All values are diagnostics for downstream filtering — nothing is +QC-masked by them, and **no pre-existing values changed** (verified per store: +every original array byte-identical before/after the update). + +### New per-trace variables (dimension `slow_time`) + +| Variable | Type | Meaning | +|---|---|---| +| `img_comb_offset_dB` | float32 | Residual seam offset (dB) of the frame's worst image pair — the power step actually present in the combined product at the image-combine boundary, measured by comparing the individual images (weight-corrected) in their overlap. Positive = the shallower image is brighter. NaN where unmeasurable (see statuses). A downstream rejection threshold of ~3 dB is suggested; season-level means vary (see below). | +| `img_comb_pair` | int8 | Which image pair the offset series refers to (1 = img1/img2, 2 = img2/img3), fixed per frame (largest \|mean offset\|); −1 undefined. | +| `surface_source_image_index` | int8 | Which image supplied the combined product's surface sample (1 = low-gain img1; ≥2 = a higher-gain image; −1 unknown). **Provenance flag, not a validity verdict**: img2-sourced surfaces are more likely saturated and may carry a season-dependent low bias (~15–20 dB measured on 2014/2017 Greenland P3), but for high-altitude DC8 operations the surface landing in img2 is normal geometry. Filter as your application requires. | +| `surface_ceiling_margin_dB` | float32 | Season ceiling minus surface power (dB), computed against the ceiling of the trace's own source-image population. Small/zero ⇒ at the ceiling (likely saturated). NaN wherever no credible ceiling was fitted — a clean season has no ceiling and therefore no margins. | + +### New per-frame attributes (parallel to `frame_names`) + +`frame_img_comb_offset_dB` (worst-pair mean residual, dB), +`frame_img_comb_status`, `frame_img_comb_weights_mode` (declared OPR combine +weights mode, provenance only). + +Statuses: `ok` (every expected image pair was assessed — measured, or +honestly `insufficient_overlap`), `partial_images` (an expected pair was +never assessed because its image file is missing; offsets from the measured +pairs are still reported), `no_combine` (single image — check not +applicable), `images_unavailable` (img files not published), +`insufficient_overlap` (no artifact-free overlap window at the frame's +geometry — e.g. 2014-era Greenland P3 at typical AGL), +`weight_recovery_failed`, `load_error` (transient, retryable), +`invalid_params`, `params_missing`, `disabled`. Except for `ok` and +`partial_images`, offsets are NaN; no estimate is ever substituted. + +### New root attribute: `saturation` + +Season-keyed dict from the saturation second pass. Populations are split as +**img1-sourced vs higher-gain-sourced** — every `surface_source_image_index +>= 2` trace is pooled into one population reported as `img2` (individual +img3+ surfaces are too rare to fit separately). Per population — `status` +(`fit_ok` / `no_plateau` / `insufficient_support`), +`level` (dB), `level_ci` (bootstrap 95%), `single_slope` (dB/decade), +`pileup_fraction`, `n_traces`, `n_bins_occupied`, `span_decades` — plus +`cross_cap_step_db` (measured img1-vs-img2 envelope offset: the season's +empirical img2 surface-power bias), `method_version`, every fixed threshold, +and a science-data fingerprint for staleness tracking. Consumers wanting +stricter acceptance than `fit_ok` can filter on `span_decades`, +`pileup_fraction`, and `level_ci`. + +### Suggested QC filtering + +Starting points — tighten or loosen per application: + +```python +seam_ok = np.abs(img_comb_offset_dB) < 3.0 # or np.isnan(...) to keep unmeasured +surf_ok = surface_source_image_index == 1 # drop likely-saturated img2 surfaces +not_sat = ~(surface_ceiling_margin_dB < 2.0) # NaN margin (no ceiling) passes +``` + +- **Seam step**: reject traces (or frames, via `frame_img_comb_offset_dB`) + with |offset| > ~3 dB. Decide explicitly how to treat NaN offsets + (statuses like `images_unavailable` / `insufficient_overlap`): the check + didn't run there, which is not evidence the frame is clean. +- **Surface saturation**: traces with `surface_ceiling_margin_dB` below + ~1–2 dB sit at their season's clip level — exclude them from surface-power + and RSSNR analyses. NaN margins mean no credible ceiling exists for that + population, not that the trace is clean or dirty. +- **Surface provenance**: for absolute surface-power work, prefer + `surface_source_image_index == 1`; where img2-sourced surfaces are kept, + the season's `cross_cap_step_db` estimates their low bias. +- **Trusting a ceiling**: before leaning on a season's `level`, check its + `span_decades` (≥0.4 is comfortable), `pileup_fraction` (well above zero), + and `level_ci` width — weak fits (e.g. 2018_Greenland_P3 img1: pileup + 1.7%) are reported, not suppressed. + +### Method notes + +- **Seam offsets** are measured on the individual OPR images after + empirically recovering each image's effective combine scaling + (`w_i = median(combined − img_i)` over the section the combined product + sourced from img_i), so the metric reflects the shipped combined product + regardless of `img_comb_weights_mode` (including 'auto' seasons, where + per-frame scalings of >100 dB were recovered on 2018_Antarctica_DC8). +- **Ceiling detection** uses a binned 99th-percentile envelope vs log range + with equal-width bins and sparse-adjacent-bin merging (0.3.0), fitted per + source-image population; a ceiling is declared when the Theil–Sen slope is + within a near-flat band (−12 to +12 dB/decade, vs the unsaturated ~−20; + a strongly rising envelope is not a ceiling either — 0.4.1). A piecewise + partial-saturation model was removed (0.4.0) after deciding zero fits on + real data across all stores. + +### Headline findings shipped with this release + +- Seam steps are typically ≤2 dB, but several seasons carry systematic + image-combine errors: 2016_Greenland_P3 (+7.4±9.9 dB), 2017_Antarctica_Basler + (−6.8±7.5 dB), 2013_Antarctica_Basler (+4.8±6.1 dB). +- Clear surface-saturation ceilings (img1 populations): 2012_Antarctica_DC8 + (−30 dB), 2014_Greenland_P3 (−31.6 dB, with a measured 16.2 dB img2 bias + step), 2013_Greenland_P3 (−50.6 dB), 2017_Antarctica_Basler (−25.9 dB), + among others; 2016–2019 Greenland P3 and most Basler seasons show no + ceiling. + +### Store rollout + +| Store | Status | +|---|---| +| ase | updated 2026-08-31; method 0.4.1 refresh 2026-09-03 | +| greenland | updated 2026-09-01; method 0.4.1 refresh 2026-09-03 | +| antarctica | updated 2026-09-01; method 0.4.1 refresh 2026-09-03 | +| utig / crosssystem | not updated (out of current scope) | + +The 0.4.1 refresh (PR-review fixes) changed one fit — 2018_Antarctica_DC8's +img2 population, a rising envelope, is no longer reported as a ceiling +(margins now NaN there) — and added descriptive attrs to the calibration +arrays; no other values moved. + +Old ASE store preserved at `icechunk/ase-backup-20260831`; the greenland and +antarctica updates were additive (full icechunk history retained in place, so +the pre-calibration snapshots remain in each store's history). + +### Code + +Implemented on branch `calibration-checks` +([PR #4](https://github.com/englacial/radar-return-statistics/pull/4)): +`93da0ec` (checks, pipeline integration, docs, viewer), `e80d852` +(sparse-adjacent-bin merging, method 0.3.0), `d5cd19e` (piecewise model +removed, method 0.4.0), `cc16499` (this changelog), `c398ff7` (PR-review +fixes, method 0.4.1). diff --git a/docs/processing.md b/docs/processing.md index a3ec786..8842f4a 100644 --- a/docs/processing.md +++ b/docs/processing.md @@ -86,6 +86,40 @@ processing: Set `opr.cache_dir` to a path to cache downloaded frames locally (fast reruns, large disk footprint) or `null` to stream (used by the full-continent configs). +### Calibration settings + +```yaml +processing: + calibration: + img_combine: true # default; set false for quick test runs + image_load_retries: 2 # transient-failure retries per image +``` + +With `img_combine: true` (the default) each frame additionally loads its +individual images (`Data_img_NN_*`) and measures the residual seam offset at +each image-combine transition plus `surface_source_image_index` — roughly 3x +wall time and 2.5x download volume per frame. Failures degrade to a per-frame +`frame_img_comb_status`; they never skip a frame. Frames left with a +retryable status (`load_error`, `disabled`) can be filled in later: + +```bash +uv run python -m radar_return_statistics.calibration_backfill config/config_antarctica.yaml +``` + +The saturation ceiling fit is a **second pass** over the completed store +(population-level; no echogram loads): + +```bash +uv run python -m radar_return_statistics.saturation_pass config/config_antarctica.yaml +uv run python -m radar_return_statistics.saturation_pass config/config_antarctica.yaml --dry-run +``` + +It writes `surface_ceiling_margin_dB` and the season-level `saturation` root +attr, stamped with a science-data fingerprint; the runner warns when later +science appends make it stale. Appending to a store written before the +calibration schema auto-migrates it (missing variables are created backfilled +with NaN / -1). + ### QC settings ```yaml diff --git a/pyproject.toml b/pyproject.toml index a9b8984..13ecb49 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,6 +19,7 @@ dependencies = [ "click", "cftime>=1.6.5", "seaborn>=0.13.2", + "pyarrow>=23.0.1", ] [build-system] diff --git a/scripts/reprocess_regression_check.py b/scripts/reprocess_regression_check.py new file mode 100644 index 0000000..dd2ec96 --- /dev/null +++ b/scripts/reprocess_regression_check.py @@ -0,0 +1,137 @@ +"""Verify reprocessing with the calibration-enabled code does NOT change any +existing per-trace values in the production stores. + +Picks frames from the greenland and antarctica stores, re-runs process_frame +with the branch code (calibration ON, production configs), and compares every +pre-existing variable trace-by-trace against the stored values. Read-only. +""" +import sys +import numpy as np +import icechunk +import zarr +from xopr import OPRConnection + +sys.path.insert(0, "src") +from radar_return_statistics.config import load_config +from radar_return_statistics.processing import process_frame + +CASES = [ + ("antarctica", "config/config_antarctica.yaml", + ["Data_20161104_05_037", "Data_20121025_03_022"]), +] + +CALIB_VARS = {"img_comb_offset_dB", "img_comb_pair", + "surface_source_image_index", "surface_ceiling_margin_dB"} + + +def open_store(name): + storage = icechunk.s3_storage(bucket="opr-radar-metrics", + prefix=f"icechunk/{name}", + region="us-west-2", anonymous=True) + repo = icechunk.Repository.open(storage=storage) + return zarr.open_group(repo.readonly_session(branch="main").store, mode="r") + + +def compare(name, stored, computed): + stored = np.asarray(stored) + computed = np.asarray(computed) + if stored.shape != computed.shape: + return f"SHAPE {stored.shape} vs {computed.shape}" + if stored.dtype.kind in "fc" or computed.dtype.kind in "fc": + s = stored.astype(float) + c = computed.astype(float) + if np.array_equal(np.isnan(s), np.isnan(c)) and np.allclose( + s[~np.isnan(s)], c[~np.isnan(c)], rtol=1e-5, atol=1e-8): + return None + both = ~(np.isnan(s) | np.isnan(c)) + max_diff = np.abs(s[both] - c[both]).max() if both.any() else 0.0 + nan_mismatch = int((np.isnan(s) != np.isnan(c)).sum()) + return f"DIFF max={max_diff:.3e}, nan_mismatch={nan_mismatch}" + if stored.dtype.kind in "US" or computed.dtype.kind in "US": + return None if np.array_equal(stored.astype(str), computed.astype(str)) else "STRING DIFF" + return None if np.array_equal(stored, computed) else f"DIFF ({int((stored != computed).sum())} elems)" + + +any_fail = False +for store_name, config_path, frame_ids in CASES: + root = open_store(store_name) + config = load_config(config_path) + store_fids = root["frame_id"] + frame_names = list(root.attrs.asdict().get("frame_names", [])) + + # auto-pick a second frame from a different season than the first + resolved = [] + for fid in frame_ids: + if fid is not None: + resolved.append(fid) + else: + first_season = resolved[0].split("_")[1][:4] if resolved else "" + pick = next((f for f in frame_names + if f not in resolved and f.split("_")[1][:4] != first_season), None) + if pick: + resolved.append(pick) + + fid_arr = store_fids[:].astype(str) + opr = OPRConnection(cache_dir=config.get("opr", {}).get("cache_dir")) + for fid in resolved: + mask = fid_arr == fid + n_store = int(mask.sum()) + if n_store == 0: + print(f"[{store_name}] {fid}: NOT IN STORE, skipping") + continue + collection = None + fi = root["frame_index"][:][mask] + colls = root.attrs.asdict().get("frame_collections") + if colls is not None and fi.size: + collection = colls[int(fi[0])] + try: + frames = opr.query_frames(collections=[collection] if collection else None, + max_items=None) + item = frames.loc[fid] + except Exception as e: + print(f"[{store_name}] {fid}: STAC lookup failed ({type(e).__name__}: {e})") + any_fail = True + continue + ds = process_frame(opr, item, config) + if ds is None: + print(f"[{store_name}] {fid}: process_frame returned None — MISMATCH vs stored frame") + any_fail = True + continue + print(f"[{store_name}] {fid} ({collection}): store n={n_store}, new n={len(ds.slow_time)}, " + f"calib status={ds.attrs.get('frame_img_comb_status', '?')}") + if n_store != len(ds.slow_time): + print(f" TRACE COUNT MISMATCH") + any_fail = True + continue + checked = failed = 0 + for var in sorted(set(root.array_keys())): + if var in CALIB_VARS or var in ("slow_time", "processed_frames", "frame_index"): + continue + arr = root[var] + if arr.ndim != 1 or arr.shape[0] != fid_arr.shape[0]: + continue + if var in ds: + computed = ds[var].values + elif var in ds.coords: + computed = ds.coords[var].values + elif var in ("latitude", "longitude", "elevation"): + computed = ds.coords[var].values if var in ds.coords else None + else: + continue + res = compare(var, arr[:][mask], computed) + checked += 1 + if res: + print(f" {var}: {res}") + failed += 1 + any_fail = True + # frame-level attr + try: + fb = root.attrs.asdict()["frame_bed_pick_fraction"][int(fi[0])] + if not np.isclose(fb, ds.attrs["frame_bed_pick_fraction"], rtol=1e-6): + print(f" frame_bed_pick_fraction: {fb} vs {ds.attrs['frame_bed_pick_fraction']}") + any_fail = True + except (KeyError, IndexError, TypeError): + pass + print(f" {checked} variables compared, {failed} mismatched") + +print("\nRESULT:", "FAIL — differences found" if any_fail else "PASS — reprocessing reproduces all stored values") diff --git a/scripts/store_science_snapshot.py b/scripts/store_science_snapshot.py new file mode 100644 index 0000000..df87c94 --- /dev/null +++ b/scripts/store_science_snapshot.py @@ -0,0 +1,85 @@ +"""Snapshot / verify a store's science arrays. + + uv run python scripts/store_science_snapshot.py snapshot + uv run python scripts/store_science_snapshot.py verify + +snapshot: records sha256 of every array's bytes (NaNs included) + attrs. +verify: recomputes and reports, per array: unchanged / MODIFIED / new. +Calibration additions are expected as 'new' or attr additions; anything +MODIFIED is a failure. Used to prove backfill + second pass touched nothing +but the new fields. +""" +import hashlib +import json +import sys + +import icechunk +import numpy as np +import zarr + +CALIB_VARS = {"img_comb_offset_dB", "img_comb_pair", + "surface_source_image_index", "surface_ceiling_margin_dB", + "surface_peak_width_us"} +CALIB_ATTRS_PREFIXES = ("frame_img_comb", "saturation", "calibration_") + + +def open_root(path): + storage = icechunk.local_filesystem_storage(path) + repo = icechunk.Repository.open(storage=storage) + return zarr.open_group(repo.readonly_session(branch="main").store, mode="r") + + +def array_hash(arr): + a = arr[:] + return hashlib.sha256(np.ascontiguousarray(a).tobytes()).hexdigest() + + +def main(): + mode, path, out = sys.argv[1], sys.argv[2], sys.argv[3] + root = open_root(path) + current = {k: array_hash(root[k]) for k in sorted(root.array_keys())} + attrs = root.attrs.asdict() + non_calib_attrs = { + k: v for k, v in attrs.items() + if not any(k.startswith(p) for p in CALIB_ATTRS_PREFIXES) + } + attrs_hash = hashlib.sha256( + json.dumps(non_calib_attrs, sort_keys=True, default=str).encode() + ).hexdigest() + + if mode == "snapshot": + with open(out, "w") as f: + json.dump({"arrays": current, "attrs_hash": attrs_hash}, f, indent=1) + print(f"snapshot: {len(current)} arrays, attrs hash {attrs_hash[:12]}") + return + + with open(out) as f: + base = json.load(f) + fail = False + for k, h in base["arrays"].items(): + if k not in current: + print(f"MISSING array: {k}") + fail = True + elif current[k] != h: + if k in CALIB_VARS: + print(f"calibration array changed (expected): {k}") + else: + print(f"MODIFIED science array: {k}") + fail = True + new = [k for k in current if k not in base["arrays"]] + print(f"new arrays: {sorted(new)}") + unexpected_new = [k for k in new if k not in CALIB_VARS] + if unexpected_new: + print(f"UNEXPECTED new arrays: {unexpected_new}") + fail = True + if attrs_hash != base["attrs_hash"]: + print("MODIFIED non-calibration attrs") + fail = True + n_unchanged = sum(1 for k, h in base["arrays"].items() if current.get(k) == h) + print(f"{n_unchanged}/{len(base['arrays'])} original arrays unchanged") + print("RESULT:", "FAIL" if fail else "PASS — science data untouched; only calibration fields added") + sys.exit(1 if fail else 0) + + +if __name__ == "__main__": + main() diff --git a/src/radar_return_statistics/calibration.py b/src/radar_return_statistics/calibration.py new file mode 100644 index 0000000..c3db171 --- /dev/null +++ b/src/radar_return_statistics/calibration.py @@ -0,0 +1,862 @@ +"""Radiometric calibration checks. + +Check 1 — image-combine consistency: recover the effective per-image weights +the OPR combine applied (by differencing the combined product against each +image over the section it sourced), then measure the residual seam offset +actually present in the combined product at each image transition. + +Check 2 — surface saturation: population-level detection of a hard ceiling in +surface power vs range-to-surface (flat-envelope rule on a binned +upper-quantile curve, fitted per surface-source-image population), run as a +second pass over store outputs. + +Sign convention: positive offset = the earlier (lower-index, shallower) image +is brighter than the later one in the combined product. +""" +import logging +import time + +import numpy as np + +logger = logging.getLogger(__name__) + +# Frame-level statuses for the image-combine check +STATUS_OK = "ok" +STATUS_NO_COMBINE = "no_combine" +STATUS_IMAGES_UNAVAILABLE = "images_unavailable" +STATUS_PARTIAL_IMAGES = "partial_images" +STATUS_LOAD_ERROR = "load_error" +STATUS_INVALID_PARAMS = "invalid_params" +STATUS_PARAMS_MISSING = "params_missing" +STATUS_INSUFFICIENT_OVERLAP = "insufficient_overlap" +STATUS_WEIGHT_RECOVERY_FAILED = "weight_recovery_failed" +STATUS_DISABLED = "disabled" + +MAX_IMAGES = 3 # per user: 3 is the normal max; no known season uses more + +DEFAULTS = { + "trim_s": 0.3e-6, # blend-zone trim around combine boundaries + "snr_gate_db": 6.0, # bins must exceed noise floor by this much + "min_bins_per_trace": 10, + "min_traces_per_pair": 30, + "min_weight_bins": 5, + "weight_spread_max_db": 0.1, # MAD*1.4826 of per-trace weights above this fails recovery + "noise_start_us": 12.0, # record-tail-convention noise window offsets + "noise_end_us": 7.0, + "noise_min_bins": 5, + "noise_fallback_percentile": 10.0, +} + + +# --------------------------------------------------------------------------- +# Param extraction +# --------------------------------------------------------------------------- + +def _get_path(attrs, path): + d = attrs + for k in path: + if not isinstance(d, dict) or k not in d: + return None + d = d[k] + return d + + +IMG_COMB_PATHS = [ + ("param_array", "array", "img_comb"), + ("param_records", "array", "img_comb"), + ("param_combine", "array_param", "img_comb"), + ("param_combine", "combine", "img_comb"), + ("param_array", "qlook", "img_comb"), + ("param_records", "qlook", "img_comb"), + ("param_combine", "get_heights", "qlook", "img_comb"), + # 2012-era attr; carries the OLD 2-per-pair [T_comb, T_guard] format, + # normalized to the modern 3-per-pair layout in find_img_comb. + ("param_combine_wf_chan", "img_comb"), +] + +TPD_PATHS = [ + ("param_records", "radar", "wfs", "Tpd"), + ("param_array", "radar", "wfs", "Tpd"), + ("param_combine", "radar", "wfs", "Tpd"), + ("param_sar", "radar", "wfs", "Tpd"), +] + +# 2012-era wfs is a list of per-waveform dicts (each with a scalar Tpd) rather +# than a dict of arrays. +WFS_LIST_PATHS = [ + ("param_csarp", "radar", "wfs"), + ("param_records", "radar", "wfs"), +] + +IMGS_PATHS = [ + ("param_array", "array", "imgs"), + ("param_records", "array", "imgs"), + ("param_combine", "combine", "imgs"), + ("param_combine", "csarp", "imgs"), + ("param_combine_wf_chan", "array", "imgs"), +] + +WEIGHTS_MODE_PATHS = [ + ("param_array", "array", "img_comb_weights_mode"), + ("param_combine", "combine", "img_comb_weights_mode"), + ("param_records", "array", "img_comb_weights_mode"), + ("param_array", "qlook", "img_comb_weights_mode"), +] + + +def find_img_comb(attrs): + """Return (img_comb vector, provenance string) or (None, None). + + A present-but-empty vector returns (empty array, provenance) so callers can + distinguish 'declared empty' (2014-era) from 'missing'. + """ + for path in IMG_COMB_PATHS: + v = _get_path(attrs, path) + if v is None: + continue + try: + arr = np.atleast_1d(np.asarray(v, dtype=float)) + except (TypeError, ValueError): + continue + # 2012-era param_combine_wf_chan uses the old 2-per-pair + # [T_comb, T_guard] format (no T_blank; verified empirically on + # 2012_Antarctica_DC8: transition at surf + T_comb, weights ~0). + # Normalize to the modern [T_comb, T_blank, T_guard] layout. + if path[0] == "param_combine_wf_chan" and arr.size and arr.size % 2 == 0: + pairs = arr.reshape(-1, 2) + arr = np.column_stack( + [pairs[:, 0], np.full(len(pairs), -np.inf), pairs[:, 1]] + ).ravel() + return arr, ".".join(path) + return None, None + + +def find_tpd(attrs): + """Per-waveform pulse durations (seconds), or None.""" + for path in TPD_PATHS: + v = _get_path(attrs, path) + if v is None: + continue + try: + return np.atleast_1d(np.asarray(v, dtype=float)) + except (TypeError, ValueError): + continue + # 2012-era: wfs is a list of per-waveform dicts with scalar Tpd + for path in WFS_LIST_PATHS: + v = _get_path(attrs, path) + if isinstance(v, (list, tuple)) and v and all(isinstance(w, dict) for w in v): + try: + return np.array([float(w["Tpd"]) for w in v]) + except (KeyError, TypeError, ValueError): + continue + return None + + +def find_imgs(attrs): + """The params imgs list (wf-adc matrices), or None.""" + for path in IMGS_PATHS: + v = _get_path(attrs, path) + if isinstance(v, (list, tuple)) and len(v) > 0: + return list(v) + return None + + +def find_weights_mode(attrs): + """Declared img_comb_weights_mode ('' if absent/empty).""" + for path in WEIGHTS_MODE_PATHS: + v = _get_path(attrs, path) + if v is not None and str(v).strip(): + return str(v).strip() + return "" + + +def tpd_for_images(imgs, tpd): + """Map each image's wf entries to a pulse duration. + + imgs entries are wf-adc matrices with waveform indices in the first row + (1-based into the radar.wfs arrays). Returns (list of Tpd seconds, + consistent: bool) — consistent False if any image's waveforms disagree. + """ + out, consistent = [], True + for entry in imgs: + m = np.atleast_2d(np.asarray(entry, dtype=float)) + # Modern layout: (2, N) with wf indices in row 0. 2012-era layout: + # (N, 2) rows of [wf, adc] pairs — wf indices in column 0. The wf + # axis is the one whose values are constant per image; prefer it. + wf_vals = m[0] + if m.shape[0] > 2 and m.shape[1] == 2: + wf_vals = m[:, 0] + elif m.shape[0] == 2 and m.shape[1] == 2 and len(set(m[0])) > 1 and len(set(m[:, 0])) == 1: + wf_vals = m[:, 0] + wfs = np.unique(wf_vals.astype(int)) + vals = [tpd[w - 1] for w in wfs if 1 <= w <= len(tpd)] + if not vals: + out.append(np.nan) + consistent = False + continue + if len(set(np.round(vals, 12))) > 1: + consistent = False + out.append(float(vals[0])) + return out, consistent + + +def expected_image_count(attrs): + """Expected image count: imgs list, else img_comb-derived, else None. + + Counts above MAX_IMAGES (seen in 2014-era per-wf-adc imgs lists) are + treated as suspect and capped. + """ + imgs = find_imgs(attrs) + if imgs is not None: + return min(len(imgs), MAX_IMAGES), "imgs" + img_comb, _ = find_img_comb(attrs) + if img_comb is not None and img_comb.size >= 3 and img_comb.size % 3 == 0: + return min(img_comb.size // 3 + 1, MAX_IMAGES), "img_comb" + return None, None + + +# --------------------------------------------------------------------------- +# Image loading +# --------------------------------------------------------------------------- + +def _load_image(opr, stac_item, data_product, i, retries=2, retry_sleep_s=2.0): + """Load image i with typed error handling. + + Returns (Dataset|None, error|None) where error is 'missing' (HTTP 404, + permanent) or 'transient' (retried, still failing). + """ + for attempt in range(retries + 1): + try: + return opr.load_frame( + stac_item, data_product=data_product, image=i, + allow_unlisted_products=True, + ), None + except FileNotFoundError: + return None, "missing" + except Exception: + if attempt == retries: + logger.warning("image %d: transient load failure", i, exc_info=True) + return None, "transient" + time.sleep(retry_sleep_s) + + +def _align_traces(ds_i, combined_ds): + """Align an image's traces to the combined frame's slow_time. + + The pipeline decimates the combined frame, so images (always full + resolution) are subset by nearest slow_time within the pipeline's 1 s + pick-alignment tolerance; unmatched traces become NaN rows. Identical + trace axes pass through untouched. + """ + if ds_i.sizes.get("slow_time") == combined_ds.sizes["slow_time"] and np.array_equal( + ds_i.slow_time.values, combined_ds.slow_time.values + ): + return ds_i + import pandas as pd + + return ds_i.sortby("slow_time").reindex( + slow_time=combined_ds.slow_time, + method="nearest", + tolerance=pd.Timedelta(seconds=1), + ) + + +# --------------------------------------------------------------------------- +# Geometry helpers (pure numpy) +# --------------------------------------------------------------------------- + +def normalize_data_db(data, twtt): + """abs->dB, oriented (n_traces, n_twtt).""" + d = np.abs(np.asarray(data)) + if d.ndim != 2: + raise ValueError("expected 2-D data") + if d.shape[1] != twtt.size: + if d.shape[0] == twtt.size: + d = d.T + else: + raise ValueError("data shape does not match twtt") + return 10.0 * np.log10(np.maximum(d, 1e-30)) + + +def combine_boundaries(img_comb, td_surf, img_t_ends): + """Per-pair, per-trace combine boundary times. + + img_comb: (n_pairs*3,) [T_comb, T_blank, T_guard] per pair. + img_t_ends: {image index: last twtt} for the earlier image of each pair. + Returns list of per-trace arrays (one per pair). + """ + td = np.where(np.isfinite(td_surf), td_surf, 0.0) + bounds = [] + for p in range(len(img_comb) // 3): + T_comb, T_blank, T_guard = img_comb[3 * p: 3 * p + 3] + b = np.maximum(T_blank, T_comb + td) + t_end = img_t_ends.get(p + 1) + if t_end is not None: + b = np.minimum(b, t_end - T_guard) + bounds.append(b) + return bounds + + +def default_img_comb(tpds, n_images): + """Tpd-default img_comb vector for when the params vector is empty/missing: + T_comb = Tpd of the later image, T_blank = -inf, T_guard = Tpd of the + earlier image.""" + vec = [] + for p in range(n_images - 1): + vec.extend([tpds[p + 1], -np.inf, tpds[p]]) + return np.asarray(vec, dtype=float) + + +def recover_weights(tc, Dc_db, images_db, bounds, trim_s=DEFAULTS["trim_s"], + min_bins=DEFAULTS["min_weight_bins"]): + """Effective weight per image: median(combined - img_i) over the section + the combined product sourced from img_i. + + images_db: {idx: (twtt, dB array (n_traces, n_twtt))}. + Returns ({idx: w}, {idx: robust spread}, {idx: n_traces_used}). + """ + n_traces = Dc_db.shape[0] + w, spread, n_used = {}, {}, {} + for i, (ti, Di) in images_db.items(): + lo_b = bounds[i - 2] if i >= 2 else None + hi_b = bounds[i - 1] if i - 1 < len(bounds) else None + vals = [] + for r in range(n_traces): + lo = (lo_b[r] + trim_s) if lo_b is not None else max(tc[0], ti[0]) + hi = (hi_b[r] - trim_s) if hi_b is not None else min(tc[-1], ti[-1]) + lo, hi = max(lo, ti[0], tc[0]), min(hi, ti[-1], tc[-1]) + if hi <= lo: + continue + m = (tc >= lo) & (tc <= hi) + if m.sum() < min_bins: + continue + diff = Dc_db[r, m] - np.interp(tc[m], ti, Di[r]) + diff = diff[np.isfinite(diff)] + if diff.size: + vals.append(np.median(diff)) + # one NaN trace (unmatched alignment, NaN samples in the section) must + # not poison the whole image's weight + vals = np.asarray(vals) + vals = vals[np.isfinite(vals)] + if vals.size == 0: + w[i], spread[i], n_used[i] = np.nan, np.nan, 0 + else: + med = np.median(vals) + w[i] = float(med) + spread[i] = float(1.4826 * np.median(np.abs(vals - med))) + n_used[i] = int(vals.size) + return w, spread, n_used + + +def image_noise_floor_db(ti, Di_db, t_guard=0.0, + start_us=DEFAULTS["noise_start_us"], + end_us=DEFAULTS["noise_end_us"], + min_bins=DEFAULTS["noise_min_bins"], + fallback_percentile=DEFAULTS["noise_fallback_percentile"]): + """Per-trace noise floor (dB): the record-tail convention (window + [end-12us, end-7us], clipped to end before T_end - T_guard for the + pulse-compression rolloff), floored by the per-trace 10th percentile of + the whole record. + + The tail convention assumes the record end lies past the bed — true for + the combined product but not for short-gate images whose record ends + mid-ice (2012 DC8 img1 ends at 18.8 us: its tail window sits in strong + englacial signal, inflating the floor ~45 dB and starving the SNR gate). + Taking the elementwise min of the two estimates keeps whichever is + uncontaminated; a slightly-low floor only admits extra bins into a + median, which is safe.""" + fallback = np.percentile(Di_db, fallback_percentile, axis=1) + t_end = ti[-1] + hi = min(t_end - end_us * 1e-6, t_end - t_guard) + lo = t_end - start_us * 1e-6 + if lo >= ti[0]: # a truncated window would slide into top-of-record signal + m = (ti >= lo) & (ti <= hi) + if m.sum() >= min_bins: + return np.minimum(np.median(Di_db[:, m], axis=1), fallback) + return fallback + + +def residual_offsets(pair_a, pair_b, w_a, w_b, boundary, t_guard, + snr_gate_db=DEFAULTS["snr_gate_db"], + min_bins=DEFAULTS["min_bins_per_trace"]): + """Per-trace residual seam offset for one image pair. + + pair_a/pair_b: (twtt, dB (n_traces, n_twtt), noise floor (n_traces,)). + boundary: per-trace combine boundary (start of the valid overlap); + the overlap ends at T_end(a) - T_guard. Offset = median over gated bins + of (dB_a + w_a) - (dB_b + w_b); positive = earlier image brighter. + """ + ta, Da, fa = pair_a + tb, Db, fb = pair_b + n_traces = Da.shape[0] + out = np.full(n_traces, np.nan) + hi = ta[-1] - t_guard + for r in range(n_traces): + lo = boundary[r] + if not np.isfinite(lo) or hi <= lo: + continue + m = (tb >= lo) & (tb <= hi) & (tb >= ta[0]) + if not m.any(): + continue + a_interp = np.interp(tb[m], ta, Da[r]) + b_vals = Db[r, m] + gate = (a_interp > fa[r] + snr_gate_db) & (b_vals > fb[r] + snr_gate_db) + if gate.sum() < min_bins: + continue + out[r] = np.median((a_interp[gate] + w_a) - (b_vals[gate] + w_b)) + return out + + +def pair_summary(offsets): + """Mean / median / MAD / drift / n_valid of a per-trace offset series.""" + ok = np.isfinite(offsets) + n = int(ok.sum()) + if n == 0: + return {"mean": np.nan, "median": np.nan, "mad": np.nan, + "drift": np.nan, "n_valid": 0} + v = offsets[ok] + med = float(np.median(v)) + half = len(offsets) // 2 + first, second = offsets[:half], offsets[half:] + drift = np.nan + if np.isfinite(first).sum() >= 5 and np.isfinite(second).sum() >= 5: + drift = float(np.nanmean(second) - np.nanmean(first)) + return { + "mean": float(np.mean(v)), + "median": med, + "mad": float(1.4826 * np.median(np.abs(v - med))), + "drift": drift, + "n_valid": n, + } + + +def surface_source_image_index(td_surf, img_gates, guards): + """Which image supplied the combined product's surface sample. + + img_gates: {idx: (t_start, t_end)}; guards: {idx: T_guard for the pair in + which idx is the earlier image} (0 / absent for the last image). + Smallest index whose gate contains td (NaN -> 0); -1 if none. + """ + td = np.where(np.isfinite(td_surf), td_surf, 0.0) + out = np.full(td.shape, -1, dtype=np.int8) + for i in sorted(img_gates): + t0, t1 = img_gates[i] + g = guards.get(i, 0.0) + hit = (out == -1) & (td >= t0) & (td <= t1 - g) + out[hit] = i + return out + + +# --------------------------------------------------------------------------- +# Frame-level orchestration +# --------------------------------------------------------------------------- + +def check_img_combine(opr, stac_item, combined_ds, data_product, + params=None): + """Run the image-combine consistency check for one frame. + + combined_ds: the already-loaded combined product. Returns a dict with + frame status, per-pair summaries, recovered weights, worst-pair per-trace + offsets, and surface_source_image_index (aligned to combined_ds traces). + """ + p = dict(DEFAULTS) + if params: + p.update(params) + attrs = combined_ds.attrs + n_traces = combined_ds.sizes["slow_time"] + result = { + "status": STATUS_OK, "pairs": [], "weights": {}, "weight_spread": {}, + "worst_pair": -1, "offsets": np.full(n_traces, np.nan), + "surface_source_image_index": np.full(n_traces, -1, dtype=np.int8), + "weights_mode": find_weights_mode(attrs), + "img_comb": None, "img_comb_provenance": None, + "n_images_expected": None, "n_images_loaded": 0, "image_errors": {}, + } + + img_comb, provenance = find_img_comb(attrs) + result["img_comb_provenance"] = provenance + expected, count_source = expected_image_count(attrs) + result["n_images_expected"] = expected + if expected == 1: + result["status"] = STATUS_NO_COMBINE + result["surface_source_image_index"][:] = 1 + return result + + tc = combined_ds.twtt.values + + # Load images sequentially; convert each to numpy immediately and release + # the dataset so at most one image's file/dataset is held at a time. + images_db, errors = {}, {} + for i in range(1, (expected or MAX_IMAGES) + 1): + ds_i, err = _load_image(opr, stac_item, data_product, i, + retries=p.get("image_load_retries", 2)) + if ds_i is None: + errors[i] = err + continue + try: + aligned = _align_traces(ds_i, combined_ds) + ti = np.asarray(aligned.twtt.values, dtype=float) + images_db[i] = (ti, normalize_data_db(aligned.Data.values, ti)) + except Exception: + logger.warning("image %d: could not align/convert", i, exc_info=True) + errors[i] = "transient" + finally: + try: + ds_i.close() + except Exception: + pass + result["n_images_loaded"] = len(images_db) + result["image_errors"] = errors + if len(images_db) <= 1: + # A true single-image frame has no img_NN files at all AND params + # declaring one image (handled above); anything else with <2 loadable + # images means combining occurred but the files aren't usable. + if any(e == "transient" for e in errors.values()): + result["status"] = STATUS_LOAD_ERROR + else: + result["status"] = STATUS_IMAGES_UNAVAILABLE + return result + + # Pair coverage must reflect EXPECTED images, not merely loaded ones: a + # trailing missing image (img3 404 when params declare 3) must surface as + # an unassessed pair (-> partial_images), never silently shorten the pair + # list and leave the frame "ok". The declared img_comb length still bounds + # the pair count below, which protects against inflated 2014-era + # per-wf-adc imgs lists. + n_images = max(max(images_db), min(expected or 0, MAX_IMAGES)) + # Resolve the img_comb vector (Tpd defaults when declared empty/missing) + if img_comb is None or img_comb.size < 3: + tpd = find_tpd(attrs) + imgs = find_imgs(attrs) + if tpd is None: + result["status"] = STATUS_PARAMS_MISSING + return result + if imgs is not None and len(imgs) >= n_images: + tpds, consistent = tpd_for_images(imgs[:n_images], tpd) + if not consistent: + logger.info("imgs entries disagree on Tpd; using first wf per image") + elif len(tpd) >= n_images: + tpds = [float(t) for t in tpd[:n_images]] + else: + result["status"] = STATUS_INVALID_PARAMS + return result + if not np.all(np.isfinite(tpds)): + result["status"] = STATUS_INVALID_PARAMS + return result + img_comb = default_img_comb(tpds, n_images) + result["img_comb_provenance"] = (provenance or "") + "+tpd_default" + result["img_comb"] = img_comb + + n_pairs_possible = min(len(img_comb) // 3, n_images - 1) + if n_pairs_possible < 1: + result["status"] = STATUS_INVALID_PARAMS + return result + + Dc = normalize_data_db(combined_ds.Data.values, tc) + td_surf = combined_ds.Surface.values if "Surface" in combined_ds else np.full(n_traces, np.nan) + + img_t_ends = {i: images_db[i][0][-1] for i in images_db} + bounds = combine_boundaries(img_comb, td_surf, img_t_ends) + + # surface source image index + guards = {p_ + 1: img_comb[3 * p_ + 2] for p_ in range(len(img_comb) // 3)} + gates = {i: (images_db[i][0][0], images_db[i][0][-1]) for i in images_db} + result["surface_source_image_index"] = surface_source_image_index( + td_surf, gates, guards) + + # effective weights + w, spread, n_used = recover_weights(tc, Dc, images_db, bounds, + trim_s=p["trim_s"], + min_bins=p["min_weight_bins"]) + result["weights"], result["weight_spread"] = w, spread + usable = {i for i in w if np.isfinite(w[i])} + recovery_ok = all( + np.isfinite(spread[i]) and spread[i] <= p["weight_spread_max_db"] + for i in usable + ) and len(usable) >= 2 + if not recovery_ok: + result["status"] = STATUS_WEIGHT_RECOVERY_FAILED + w = {i: 0.0 for i in images_db} # raw offsets only (parquet diagnostic) + + # residual offsets per adjacent pair + floors = {i: image_noise_floor_db( + images_db[i][0], images_db[i][1], + t_guard=guards.get(i, 0.0), + start_us=p["noise_start_us"], end_us=p["noise_end_us"], + min_bins=p["noise_min_bins"], + fallback_percentile=p["noise_fallback_percentile"], + ) for i in images_db} + + pair_series = {} + any_ok, any_missing = False, False + # min_traces_per_pair is tuned against full-resolution frames (thousands of + # traces); cap it proportionally so decimated pipeline frames (~tens of + # traces) are not held to a far stricter relative bar. + min_traces = min(p["min_traces_per_pair"], max(5, int(0.25 * n_traces))) + for pr in range(n_pairs_possible): + a, b = pr + 1, pr + 2 + if a not in images_db or b not in images_db: + any_missing = True + result["pairs"].append({"pair": pr + 1, "a": a, "b": b, + "status": STATUS_IMAGES_UNAVAILABLE, + **pair_summary(np.full(n_traces, np.nan))}) + continue + offs = residual_offsets( + (images_db[a][0], images_db[a][1], floors[a]), + (images_db[b][0], images_db[b][1], floors[b]), + w.get(a, 0.0), w.get(b, 0.0), + bounds[pr], img_comb[3 * pr + 2], + snr_gate_db=p["snr_gate_db"], min_bins=p["min_bins_per_trace"], + ) + summ = pair_summary(offs) + status = STATUS_OK if summ["n_valid"] >= min_traces else STATUS_INSUFFICIENT_OVERLAP + result["pairs"].append({"pair": pr + 1, "a": a, "b": b, "status": status, **summ}) + if status == STATUS_OK: + any_ok = True + pair_series[pr + 1] = offs + + if result["status"] == STATUS_OK: + if not any_ok: + result["status"] = STATUS_INSUFFICIENT_OVERLAP + elif any_missing: + result["status"] = STATUS_PARTIAL_IMAGES + + if pair_series: + worst = max(pair_series, key=lambda k: abs(result["pairs"][k - 1]["mean"])) + result["worst_pair"] = worst + result["offsets"] = pair_series[worst] + return result + + +# --------------------------------------------------------------------------- +# Check 2: surface saturation (second pass, plain numpy) +# --------------------------------------------------------------------------- + +SAT_DEFAULTS = { + "quantile": 99.0, + "n_bins": 24, + "min_traces_per_bin": 50, + "min_occupied_bins": 6, + # 0.25 decades for every fit — source-restricted subpopulations inherently + # span less range, and the user chose one uniform minimum over + # special-casing (2026-08-31). + "min_span_decades": 0.25, + # Flat rule: the envelope slope must sit in a near-flat band to count as a + # ceiling. The unsaturated envelope falls ~-20 dB/decade; a strongly + # RISING envelope is not a clip ceiling either (gain/geometry artifact), + # so the band is two-sided. + "flat_slope_db_per_decade": -12.0, + "flat_slope_max_db_per_decade": 12.0, + "pileup_delta_db": 1.0, + "n_boot": 200, +} + +FIT_OK = "fit_ok" +NO_PLATEAU = "no_plateau" +INSUFFICIENT_SUPPORT = "insufficient_support" + + +def _theil_sen(x, y): + """Theil-Sen slope + intercept (median of pairwise slopes).""" + n = len(x) + slopes = [(y[j] - y[i]) / (x[j] - x[i]) + for i in range(n) for j in range(i + 1, n) if x[j] != x[i]] + if not slopes: + return np.nan, np.nan + slope = np.median(slopes) + return float(slope), float(np.median(y - slope * x)) + + +def _binned_upper_quantile(power_db, range_m, params): + """Equal-width bins in log10(range), with adjacent sparse bins merged. + + Equal-width base bins keep bin spacing (and hence slope estimation) + faithful to the range geometry rather than the altitude histogram. Bins + under min_traces_per_bin are greedily merged with the next bin(s) until + they reach the floor, so a thin ceiling-following tail (2012 DC8 img1: + 95% of traces at 400-700 m, real tail to 3 km) contributes a few wide + bins that extend the fitted span instead of being discarded — without + concentrating extra bins in the dense band, which would dilute the + Theil-Sen slope and admit false flats. Bin x = median log10(r) of the + (merged) bin's traces; a trailing merged bin still under the floor is + dropped. + """ + logr = np.log10(range_m) + lo, hi = np.percentile(logr, [1, 99]) + edges = np.linspace(lo, hi, params["n_bins"] + 1) + idx = np.digitize(logr, edges) + xs, qs, ns = [], [], [] + pend_sel = np.zeros(logr.size, dtype=bool) + for b in range(1, params["n_bins"] + 1): + pend_sel |= idx == b + n = int(pend_sel.sum()) + if n >= params["min_traces_per_bin"]: + xs.append(float(np.median(logr[pend_sel]))) + qs.append(float(np.percentile(power_db[pend_sel], params["quantile"]))) + ns.append(n) + pend_sel = np.zeros(logr.size, dtype=bool) + return np.asarray(xs), np.asarray(qs), np.asarray(ns) + + +def fit_ceiling(power_db, range_m, params=None, rng=None): + """Season/segment-level ceiling (clip level) detection via the + flat-envelope rule: a binned upper-quantile envelope whose Theil-Sen slope + is far shallower than the unsaturated ~-20 dB/decade is a ceiling. + + (A piecewise plateau+decline model for partial saturation was removed in + method 0.4.0: it never decided a fit on any real season across all + stores — every detection came from the flat rule — and it complicated the + two-regime hardening. See docs/dataset_changelog.md.) + + Returns a dict: status, level, level_ci, single_slope, pileup_fraction, + n_traces, n_bins_occupied, span_decades. + """ + p = dict(SAT_DEFAULTS) + if params: + p.update(params) + rng = rng or np.random.default_rng(0) + ok = np.isfinite(power_db) & np.isfinite(range_m) & (range_m > 1) + power_db, range_m = np.asarray(power_db)[ok], np.asarray(range_m)[ok] + out = {"status": INSUFFICIENT_SUPPORT, "level": np.nan, "level_ci": (np.nan, np.nan), + "single_slope": np.nan, + "pileup_fraction": np.nan, "n_traces": int(ok.sum()), + "n_bins_occupied": 0, "span_decades": np.nan} + if power_db.size < p["min_traces_per_bin"] * p["min_occupied_bins"]: + return out + xs, qs, ns = _binned_upper_quantile(power_db, range_m, p) + out["n_bins_occupied"] = len(xs) + if len(xs) < p["min_occupied_bins"]: + return out + span = xs[-1] - xs[0] + out["span_decades"] = float(span) + if span < p["min_span_decades"]: + return out + + single_slope, _ = _theil_sen(xs, qs) + out["single_slope"] = single_slope + + def decide(xs_, qs_): + s, _ = _theil_sen(xs_, qs_) + if p["flat_slope_db_per_decade"] <= s <= p["flat_slope_max_db_per_decade"]: + return float(np.median(qs_)) # near-flat envelope = ceiling + return None + + level = decide(xs, qs) + if level is None: + out["status"] = NO_PLATEAU + return out + out["level"] = level + out["status"] = FIT_OK + + # bootstrap CI on the level (resample occupied bins) + levels = [] + for _ in range(p["n_boot"]): + sel = np.sort(rng.integers(0, len(xs), len(xs))) + lvl = decide(xs[sel], qs[sel]) + if lvl is not None: + levels.append(lvl) + if len(levels) >= 20: + out["level_ci"] = (float(np.percentile(levels, 2.5)), + float(np.percentile(levels, 97.5))) + out["pileup_fraction"] = pileup_fraction(power_db, out["level"], p["pileup_delta_db"]) + return out + + +def pileup_fraction(power_db, level, delta_db=SAT_DEFAULTS["pileup_delta_db"]): + """Fraction of finite traces within delta of the ceiling level.""" + ok = np.isfinite(power_db) + if not ok.any() or not np.isfinite(level): + return np.nan + return float((power_db[ok] >= level - delta_db).mean()) + + +def ceiling_margin(power_db, fit): + """Per-trace ceiling margin (dB): level - power. All-NaN unless fit_ok.""" + power_db = np.asarray(power_db, dtype=float) + if fit.get("status") != FIT_OK: + return np.full(power_db.shape, np.nan) + return fit["level"] - power_db + + +def cross_cap_step(power_db, range_m, m1, m2, params=None): + """Empirical img2-population bias: img1 minus img2 upper-quantile curves on + a common log-range binning, over overlapping bins; if the populations are + range-disjoint (the usual case — the gate cap splits by range), difference + of the two adjacent boundary bins instead. NaN when either side is empty. + """ + p = dict(SAT_DEFAULTS) + if params: + p.update(params) + ok = np.isfinite(power_db) & np.isfinite(range_m) & (range_m > 1) + m1, m2 = m1 & ok, m2 & ok + if m1.sum() < p["min_traces_per_bin"] or m2.sum() < p["min_traces_per_bin"]: + return np.nan + logr = np.log10(np.asarray(range_m, dtype=float)) + lo, hi = np.percentile(logr[m1 | m2], [1, 99]) + edges = np.linspace(lo, hi, p["n_bins"] + 1) + + def q_per_bin(mask): + idx = np.digitize(logr[mask], edges) + pw = np.asarray(power_db)[mask] + out = {} + for b in range(1, p["n_bins"] + 1): + sel = idx == b + if sel.sum() >= p["min_traces_per_bin"]: + out[b] = float(np.percentile(pw[sel], p["quantile"])) + return out + + q1, q2 = q_per_bin(m1), q_per_bin(m2) + if not q1 or not q2: + return np.nan + common = sorted(set(q1) & set(q2)) + if common: + return float(np.median([q1[b] - q2[b] for b in common])) + return q1[max(q1)] - q2[min(q2)] # adjacent boundary bins + + +def fit_season_by_source(power_db, range_m, source_index, params=None): + """Season-level saturation fits split by surface source image. + + Where the source index is populated, fit the img1-sourced (index==1) and + img2-sourced (index>=2) populations separately (relaxed span minimum) and + report the cross-cap step as the empirical img2 bias. Where it is not + (pre-schema stores, images unavailable), fall back to one all-traces fit. + """ + p = dict(SAT_DEFAULTS) + if params: + p.update(params) + src = np.asarray(source_index) + power_db = np.asarray(power_db, dtype=float) + range_m = np.asarray(range_m, dtype=float) + if not (src > 0).any(): + return { + "ceiling_fit_population": "all_traces", + "populations": {"all": fit_ceiling(power_db, range_m, p)}, + "cross_cap_step_db": np.nan, + } + m1, m2 = src == 1, src >= 2 + return { + "ceiling_fit_population": "by_source_image", + "populations": { + "img1": fit_ceiling(power_db[m1], range_m[m1], p), + "img2": fit_ceiling(power_db[m2], range_m[m2], p), + }, + "cross_cap_step_db": cross_cap_step(power_db, range_m, m1, m2, p), + } + + +def margins_by_source(power_db, source_index, season_fit): + """Per-trace ceiling margin against the trace's own source population's + fit when that fit is fit_ok; NaN otherwise (including unknown source in + by_source_image mode).""" + power_db = np.asarray(power_db, dtype=float) + src = np.asarray(source_index) + out = np.full(power_db.shape, np.nan) + pops = season_fit["populations"] + if season_fit["ceiling_fit_population"] == "all_traces": + return ceiling_margin(power_db, pops["all"]) + for key, mask in (("img1", src == 1), ("img2", src >= 2)): + fit = pops.get(key, {}) + if fit.get("status") == FIT_OK: + out[mask] = fit["level"] - power_db[mask] + return out diff --git a/src/radar_return_statistics/calibration_backfill.py b/src/radar_return_statistics/calibration_backfill.py new file mode 100644 index 0000000..bdb0c99 --- /dev/null +++ b/src/radar_return_statistics/calibration_backfill.py @@ -0,0 +1,196 @@ +"""Calibration-only backfill: run the image-combine check for frames whose +frame_img_comb_status is retryable — or absent entirely (pseudo-status +``missing``, for stores written before the calibration schema) — and update +their calibration variables in place. Science metrics and processed_frames are +never touched. + + uv run python -m radar_return_statistics.calibration_backfill + uv run python -m radar_return_statistics.calibration_backfill \\ + --statuses missing --workers 8 # populate a pre-calibration store + +Commits every --commit-every updated frames, so an interrupted run resumes by +re-running the same command (already-backfilled frames have a status and drop +out of the ``missing`` selection). +""" +import logging +import multiprocessing +from concurrent.futures import ProcessPoolExecutor, as_completed + +import click +import numpy as np +import zarr +from xopr import OPRConnection + +from . import store as store_mod +from .config import load_config +from .processing import decimate_frame, run_frame_calibration + +logger = logging.getLogger(__name__) + +RETRYABLE = ("load_error", "disabled") +MISSING = "missing" + + +def _calibrate_frame_worker(stac_item_row, frame_id, config): + """Load one frame and run calibration; returns (per_trace, frame_attrs). + + Streaming mode (no cache_dir) spools downloads into a per-frame temp dir + removed when the frame is done — same disk-leak guard as + runner._process_frame_worker; without it each pool worker's cache grows + unbounded for the life of the run. + """ + cache_dir = config["opr"].get("cache_dir") + if cache_dir: + return _calibrate_frame(OPRConnection(cache_dir=cache_dir), + stac_item_row, frame_id, config) + + import shutil + import tempfile + + import fsspec + + workdir = tempfile.mkdtemp(prefix="rrs_calib_cache_") + fsspec.config.conf["simplecache"] = {"cache_storage": workdir} + try: + return _calibrate_frame(OPRConnection(cache_dir=None), + stac_item_row, frame_id, config) + finally: + fsspec.config.conf.pop("simplecache", None) + shutil.rmtree(workdir, ignore_errors=True) + + +def _calibrate_frame(opr, stac_item_row, frame_id, config): + frame = opr.load_frame(stac_item_row, data_product=config["processing"]["data_product"]) + frame = frame.sortby("slow_time") + frame = decimate_frame(frame, config["processing"].get("decimate_interval")) + calib = run_frame_calibration(opr, stac_item_row, frame, + config["processing"], frame_id=frame_id) + frame.close() + per_trace = { + "img_comb_offset_dB": calib["img_comb_offset_dB"], + "img_comb_pair": calib["img_comb_pair"], + "surface_source_image_index": calib["surface_source_image_index"], + } + frame_attrs = { + "frame_img_comb_status": calib["status"], + "frame_img_comb_weights_mode": calib["weights_mode"], + "frame_img_comb_offset_dB": calib["frame_offset"], + } + return per_trace, frame_attrs + + +@click.command() +@click.argument("config_path", type=click.Path(exists=True)) +@click.option("--statuses", default=",".join(RETRYABLE), show_default=True, + help=f"Comma-separated statuses to retry ('{MISSING}' selects " + "frames with no calibration status at all)") +@click.option("--limit", default=0, show_default=True, help="Max frames (0 = all)") +@click.option("--workers", default=0, show_default=True, + help="Parallel frame workers (0 = processing.max_workers)") +@click.option("--commit-every", default=100, show_default=True, + help="Commit after this many updated frames") +@click.option("-v", "--verbose", is_flag=True) +def main(config_path: str, statuses: str, limit: int, workers: int, + commit_every: int, verbose: bool) -> None: + """Run the img-combine check for frames with retryable/missing status.""" + logging.basicConfig(level=logging.DEBUG if verbose else logging.INFO, + format="%(asctime)s %(levelname)s %(message)s") + retry_set = {s.strip() for s in statuses.split(",") if s.strip()} + config = load_config(config_path) + repo = store_mod.open_or_create_repo(config["store"]) + root = zarr.open_group(repo.readonly_session(branch="main").store, mode="r") + + names = list(root.attrs.get("frame_names", []) or []) + status_list = list(root.attrs.get("frame_img_comb_status", []) or []) + colls = list(root.attrs.get("frame_collections", []) or []) + + def effective_status(i): + if i < len(status_list) and status_list[i]: + return status_list[i] + return MISSING + + targets = [ + (name, colls[i] if i < len(colls) else "") + for i, name in enumerate(names) + if effective_status(i) in retry_set + ] + if limit: + targets = targets[:limit] + if not targets: + click.echo("No frames with retryable calibration status") + return + n_workers = workers or config["processing"].get("max_workers", 4) + logger.info("Backfilling calibration for %d frames (%d workers)", + len(targets), n_workers) + + opr = OPRConnection(cache_dir=config["opr"].get("cache_dir")) + by_coll: dict[str, list[str]] = {} + for name, coll in targets: + by_coll.setdefault(coll, []).append(name) + + # Resolve STAC rows up front (queries are per-collection and cheap). + jobs = [] # (frame_id, stac_row) + for coll, frame_ids in by_coll.items(): + try: + frames_gdf = opr.query_frames(collections=[coll] if coll else None, + max_items=None) + except Exception: + logger.exception("Query failed for collection %r", coll) + continue + for fid in frame_ids: + if fid not in frames_gdf.index: + logger.warning("Frame %s not found in STAC query; skipping", fid) + continue + jobs.append((fid, frames_gdf.loc[fid])) + + session = repo.writable_session("main") + n_done = n_failed = n_since_commit = 0 + + def write_result(fid, per_trace, frame_attrs): + nonlocal n_done, n_since_commit, session + n = store_mod.update_frame_calibration(session, fid, per_trace, frame_attrs) + if n: + n_done += 1 + n_since_commit += 1 + logger.info("Frame %s: backfilled %d traces (status %s) [%d/%d]", + fid, n, frame_attrs["frame_img_comb_status"], n_done, len(jobs)) + else: + logger.warning("Frame %s: no traces updated (missing or trace-count " + "mismatch)", fid) + if n_since_commit >= commit_every: + store_mod.commit_session(session, f"[calibration] backfill {n_since_commit} frames") + n_since_commit = 0 + session = repo.writable_session("main") + + if n_workers <= 1: + for fid, row in jobs: + try: + per_trace, frame_attrs = _calibrate_frame_worker(row, fid, config) + except Exception: + logger.exception("Frame %s: backfill failed; leaving status", fid) + n_failed += 1 + continue + write_result(fid, per_trace, frame_attrs) + else: + # spawn: workers must not inherit thread state (icechunk tokio runtime) + mp_ctx = multiprocessing.get_context("spawn") + with ProcessPoolExecutor(max_workers=n_workers, mp_context=mp_ctx) as ex: + futures = {ex.submit(_calibrate_frame_worker, row, fid, config): fid + for fid, row in jobs} + for fut in as_completed(futures): + fid = futures[fut] + try: + per_trace, frame_attrs = fut.result() + except Exception: + logger.exception("Frame %s: backfill failed; leaving status", fid) + n_failed += 1 + continue + write_result(fid, per_trace, frame_attrs) + + if n_since_commit: + store_mod.commit_session(session, f"[calibration] backfill {n_since_commit} frames") + click.echo(f"backfilled {n_done}/{len(jobs)} frames ({n_failed} failed)") + + +if __name__ == "__main__": + main() diff --git a/src/radar_return_statistics/calibration_cli.py b/src/radar_return_statistics/calibration_cli.py new file mode 100644 index 0000000..baf9930 --- /dev/null +++ b/src/radar_return_statistics/calibration_cli.py @@ -0,0 +1,281 @@ +"""Standalone CLI for baseline calibration runs (no store writes). + + uv run python -m radar_return_statistics.calibration_cli --check img-combine \ + --collection 2018_Greenland_P3 --sample-per-segment 2 --max-segments 5 --name baseline + uv run python -m radar_return_statistics.calibration_cli --check saturation \ + --store greenland --store ase --name baseline + +Writes parquet tables, figures, and a JSON run manifest to +outputs/calibration//. +""" +import json +import logging +import shutil +import subprocess +import tempfile +import time +from datetime import datetime, timezone +from pathlib import Path + +import click +import numpy as np +import pandas as pd + +from . import calibration as cal + +logger = logging.getLogger(__name__) + +METHOD_VERSION = "0.1.0" +S3_BUCKET = "opr-radar-metrics" +S3_REGION = "us-west-2" +STORES = ["antarctica", "greenland", "ase", "utig", "crosssystem"] + + +def _git_rev(): + try: + return subprocess.run(["git", "rev-parse", "HEAD"], capture_output=True, + text=True, timeout=10).stdout.strip() + except Exception: + return None + + +def _segment_of(frame_id): + # Data_20180426_03_001 -> 20180426_03 + parts = frame_id.replace("Data_", "").rsplit("_", 1) + return parts[0] + + +def sample_frames(frame_ids, sample_per_segment, max_segments, seed): + """Deterministic per-segment sampling.""" + rng = np.random.default_rng(seed) + by_seg = {} + for fid in sorted(frame_ids): + by_seg.setdefault(_segment_of(fid), []).append(fid) + segs = sorted(by_seg) + if max_segments and len(segs) > max_segments: + segs = sorted(rng.choice(segs, max_segments, replace=False)) + out = [] + for seg in segs: + frames = sorted(by_seg[seg]) + if len(frames) > sample_per_segment: + frames = sorted(rng.choice(frames, sample_per_segment, replace=False)) + out.extend(frames) + return out + + +def run_img_combine(collections, sample_per_segment, max_segments, seed, + data_product, out_dir, manifest): + from xopr import OPRConnection + + cache_dir = Path(tempfile.mkdtemp(prefix="opr_calib_cache_")) + frame_rows, pair_rows, failures = [], [], [] + try: + for coll in collections: + opr = OPRConnection(cache_dir=str(cache_dir)) + frames = opr.query_frames(collections=[coll], max_items=2000) + selected = sample_frames(frames.index.tolist(), sample_per_segment, + max_segments, seed) + manifest["selected_frames"][coll] = selected + logger.info("%s: %d frames selected", coll, len(selected)) + for fid in selected: + t0 = time.time() + try: + item = frames.loc[fid] + ds = opr.load_frame(item, data_product=data_product) + res = cal.check_img_combine(opr, item, ds, data_product) + ds.close() + except Exception as e: + logger.exception("%s: failed", fid) + failures.append({"frame": fid, "collection": coll, + "error": f"{type(e).__name__}: {e}"}) + continue + w = res["weights"] + worst = res["worst_pair"] + worst_summ = next((p for p in res["pairs"] if p["pair"] == worst), None) + frame_rows.append({ + "collection": coll, "frame": fid, "status": res["status"], + "worst_pair": worst, + "mean_offset_db": worst_summ["mean"] if worst_summ else np.nan, + "median_offset_db": worst_summ["median"] if worst_summ else np.nan, + "mad_db": worst_summ["mad"] if worst_summ else np.nan, + "drift_db": worst_summ["drift"] if worst_summ else np.nan, + "n_valid": worst_summ["n_valid"] if worst_summ else 0, + "w1": w.get(1, np.nan), "w2": w.get(2, np.nan), "w3": w.get(3, np.nan), + "w_spread_max": max([s for s in res["weight_spread"].values() + if np.isfinite(s)], default=np.nan), + "weights_mode": res["weights_mode"], + "img_comb_provenance": res["img_comb_provenance"], + "n_images_loaded": res["n_images_loaded"], + "n_images_expected": res["n_images_expected"], + "surface_img1_frac": float((res["surface_source_image_index"] == 1).mean()), + "seconds": round(time.time() - t0, 1), + }) + for pr in res["pairs"]: + raw = pr["mean"] - (w.get(pr["a"], 0.0) - w.get(pr["b"], 0.0)) \ + if np.isfinite(pr["mean"]) else np.nan + pair_rows.append({"collection": coll, "frame": fid, + "raw_mean_offset_db": raw, **pr}) + # bound cache growth: clear per frame + for f in cache_dir.glob("**/*"): + if f.is_file(): + f.unlink(missing_ok=True) + finally: + shutil.rmtree(cache_dir, ignore_errors=True) + + manifest["failures"] = failures + frames_df = pd.DataFrame(frame_rows) + pairs_df = pd.DataFrame(pair_rows) + frames_df.to_parquet(out_dir / "img_combine_frames.parquet") + pairs_df.to_parquet(out_dir / "img_combine_pairs.parquet") + + if len(frames_df): + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + fig, ax = plt.subplots(figsize=(10, 5)) + for coll, g in pairs_df[pairs_df.status == cal.STATUS_OK].groupby("collection"): + ax.hist(g["mean"], bins=40, alpha=0.6, label=f"{coll} (n={len(g)})") + ax.set_xlabel("residual seam offset, pair mean (dB)") + ax.set_ylabel("frame-pairs") + ax.axvline(0, color="k", lw=0.5) + ax.legend(fontsize=8) + ax.grid(True, alpha=0.3) + fig.savefig(out_dir / "figures" / "img_combine_offsets.png", dpi=120, + bbox_inches="tight") + plt.close(fig) + return frames_df, pairs_df + + +def _open_store(prefix): + import icechunk + import zarr + storage = icechunk.s3_storage(bucket=S3_BUCKET, prefix=f"icechunk/{prefix}", + region=S3_REGION, anonymous=True) + repo = icechunk.Repository.open(storage=storage) + session = repo.readonly_session(branch="main") + snapshot = getattr(session, "snapshot_id", None) + return zarr.open_group(session.store, mode="r"), str(snapshot) + + +def run_saturation(stores, out_dir, manifest): + import scipy.constants + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + rows, seg_rows = [], [] + for prefix in stores: + try: + root, snapshot = _open_store(prefix) + except Exception as e: + manifest["failures"].append({"store": prefix, "error": str(e)}) + continue + manifest["store_snapshots"][prefix] = snapshot + attrs = root.attrs.asdict() + power = root["surface_power_dB"][:] + twtt = root["surface_twtt"][:] + qc_key = "qc_surface_pass" if "qc_surface_pass" in root else "qc_pass" + qc = root[qc_key][:].astype(bool) + r = scipy.constants.c * twtt / 2.0 + colls = attrs.get("frame_collections") + names = attrs.get("frame_names") + if colls is None: + manifest["failures"].append({"store": prefix, + "error": "no frame_collections attr"}) + continue + fidx = root["frame_index"][:] + seasons = np.asarray(colls)[fidx] + segs = np.asarray([_segment_of(n) for n in names])[fidx] if names else None + + for season in np.unique(seasons): + m = (seasons == season) & qc + fit = cal.fit_ceiling(power[m], r[m]) + rows.append({"store": prefix, "season": season, + **{k: (str(v) if isinstance(v, tuple) else v) + for k, v in fit.items()}}) + # per-segment regime check (segments are altitude-uniform, so + # two-regime seasons resolve into single-regime segments) + if segs is not None: + for seg in np.unique(segs[m]): + ms = m & (segs == seg) + if ms.sum() < 2000: + continue + sf = cal.fit_ceiling(power[ms], r[ms]) + seg_rows.append({"store": prefix, "season": season, "segment": seg, + "status": sf["status"], "level": sf["level"], + "n_traces": sf["n_traces"]}) + # figure + mm = m & np.isfinite(power) & np.isfinite(r) & (r > 100) + if mm.sum() > 500: + fig, ax = plt.subplots(figsize=(7, 5)) + rr, pp = r[mm], power[mm] + if mm.sum() > 60000: + sel = np.random.default_rng(0).choice(mm.sum(), 60000, replace=False) + rr, pp = rr[sel], pp[sel] + ax.hexbin(rr, pp, gridsize=70, bins="log", cmap="viridis", xscale="log") + rline = np.logspace(np.log10(rr.min()), np.log10(rr.max()), 50) + ax.plot(rline, np.median(pp) - 20 * (np.log10(rline) - np.log10(np.median(rr))), + "r--", lw=1, label="-20 dB/decade") + if fit["status"] == cal.FIT_OK: + ax.axhline(fit["level"], color="w", ls="-", lw=1.5, + label=f"ceiling {fit['level']:.1f} dB") + ax.set_title(f"{season} [{fit['status']}] slope={fit['single_slope']:.1f} dB/dec", + fontsize=9) + ax.set_xlabel("range to surface (m)") + ax.set_ylabel("surface power (dB)") + ax.legend(fontsize=8) + fig.savefig(out_dir / "figures" / f"saturation_{prefix}_{season}.png", + dpi=110, bbox_inches="tight") + plt.close(fig) + + df = pd.DataFrame(rows) + df.to_parquet(out_dir / "saturation_seasons.parquet") + if seg_rows: + pd.DataFrame(seg_rows).to_parquet(out_dir / "saturation_segments.parquet") + return df + + +@click.command() +@click.option("--check", "checks", type=click.Choice(["img-combine", "saturation"]), + multiple=True, required=True) +@click.option("--collection", "collections", multiple=True, + help="Collections for img-combine") +@click.option("--store", "stores", multiple=True, + help=f"Store prefixes for saturation ({', '.join(STORES)}); default all") +@click.option("--sample-per-segment", default=2, show_default=True) +@click.option("--max-segments", default=5, show_default=True) +@click.option("--seed", default=0, show_default=True) +@click.option("--data-product", default="CSARP_standard", show_default=True) +@click.option("--name", default="baseline", show_default=True) +@click.option("-v", "--verbose", is_flag=True) +def main(checks, collections, stores, sample_per_segment, max_segments, seed, + data_product, name, verbose): + """Run standalone calibration baseline checks (read-only; no store writes).""" + logging.basicConfig(level=logging.DEBUG if verbose else logging.INFO, + format="%(asctime)s %(levelname)s %(message)s") + out_dir = Path("outputs/calibration") / name + (out_dir / "figures").mkdir(parents=True, exist_ok=True) + manifest = { + "method_version": METHOD_VERSION, + "git_rev": _git_rev(), + "started": datetime.now(timezone.utc).isoformat(), + "args": {"checks": list(checks), "collections": list(collections), + "stores": list(stores) or STORES, "sample_per_segment": sample_per_segment, + "max_segments": max_segments, "seed": seed, + "data_product": data_product}, + "params": {"img_combine": cal.DEFAULTS, "saturation": cal.SAT_DEFAULTS}, + "selected_frames": {}, "store_snapshots": {}, "failures": [], + } + if "img-combine" in checks: + run_img_combine(list(collections), sample_per_segment, max_segments, + seed, data_product, out_dir, manifest) + if "saturation" in checks: + run_saturation(list(stores) or STORES, out_dir, manifest) + manifest["finished"] = datetime.now(timezone.utc).isoformat() + with open(out_dir / "manifest.json", "w") as f: + json.dump(manifest, f, indent=2, default=str) + click.echo(f"wrote {out_dir}") + + +if __name__ == "__main__": + main() diff --git a/src/radar_return_statistics/config.py b/src/radar_return_statistics/config.py index d16aae2..f007674 100644 --- a/src/radar_return_statistics/config.py +++ b/src/radar_return_statistics/config.py @@ -3,16 +3,12 @@ import yaml -def load_config(config_path: str | Path) -> dict: - """Load and return config from YAML file.""" - config_path = Path(config_path) - if not config_path.exists(): - raise FileNotFoundError(f"Config file not found: {config_path}") +def normalize_config(config: dict) -> dict: + """Apply defaults in place and return the config. - with open(config_path) as f: - config = yaml.safe_load(f) - - # Apply defaults + Shared by load_config and process_frame so behavior (notably the + calibration default) cannot depend on the entry point. + """ config.setdefault("opr", {}) config.setdefault("region", {}) config.setdefault("query", {}) @@ -36,6 +32,14 @@ def load_config(config_path: str | Path) -> dict: post = noise.setdefault("post_bed", {}) post.setdefault("start_offset_us", 5.0) post.setdefault("end_offset_us", 5.0) + + # Radiometric calibration checks. img_combine defaults ON: the dataset is + # processed once and used many times, so the extra image downloads are + # accepted. Disable explicitly for quick test runs. + calibration = config["processing"].setdefault("calibration", {}) + calibration.setdefault("img_combine", True) + calibration.setdefault("image_load_retries", 2) + config.setdefault("qc", {}) config["qc"].setdefault("max_heading_change_deg_per_km", None) config["qc"].setdefault("min_ice_thickness_m", None) @@ -47,3 +51,15 @@ def load_config(config_path: str | Path) -> dict: config["store"].setdefault("remove_out_of_scope", False) return config + + +def load_config(config_path: str | Path) -> dict: + """Load and return config from YAML file.""" + config_path = Path(config_path) + if not config_path.exists(): + raise FileNotFoundError(f"Config file not found: {config_path}") + + with open(config_path) as f: + config = yaml.safe_load(f) + + return normalize_config(config) diff --git a/src/radar_return_statistics/processing.py b/src/radar_return_statistics/processing.py index 2fcf3f2..33b9970 100644 --- a/src/radar_return_statistics/processing.py +++ b/src/radar_return_statistics/processing.py @@ -1,4 +1,5 @@ import logging +import time import warnings import numpy as np @@ -8,6 +9,9 @@ from xopr import OPRConnection from xopr import qc as xopr_qc +from . import calibration +from .config import normalize_config + logger = logging.getLogger(__name__) SURFACE_KEY = "standard:surface" @@ -293,8 +297,71 @@ def _build_qc_checks(qc_config: dict) -> dict: return checks +def decimate_frame(frame: xr.Dataset, decimate_interval) -> xr.Dataset: + """Keep the first trace of each decimation interval (deterministic, so a + later calibration backfill reproduces the pipeline's trace selection).""" + if not decimate_interval: + return frame + interval = pd.Timedelta(decimate_interval) + times = frame.slow_time.values + selected = [0] + last = times[0] + for idx in range(1, len(times)): + if times[idx] - last >= interval: + selected.append(idx) + last = times[idx] + return frame.isel(slow_time=selected) + + +def run_frame_calibration(opr, stac_item, frame, proc_config, frame_id="?"): + """Image-combine calibration check for one frame, degraded to placeholders + plus a status on any failure — calibration never fails a frame. + + Returns per-trace arrays (offsets/pair/source index) and frame-level attrs. + When weight recovery failed, the residuals are unreliable, so the store + gets NaN offsets (raw offsets live only in baseline-CLI parquet). + """ + n = frame.sizes["slow_time"] + out = { + "img_comb_offset_dB": np.full(n, np.nan, dtype=np.float32), + "img_comb_pair": np.full(n, -1, dtype=np.int8), + "surface_source_image_index": np.full(n, -1, dtype=np.int8), + "status": calibration.STATUS_DISABLED, + "frame_offset": None, + "weights_mode": "", + } + cfg = proc_config.get("calibration", {}) + if not cfg.get("img_combine", True): + return out + t0 = time.time() + try: + res = calibration.check_img_combine( + opr, stac_item, frame, proc_config["data_product"], + params={"image_load_retries": cfg.get("image_load_retries", 2)}, + ) + out["status"] = res["status"] + out["weights_mode"] = res["weights_mode"] + out["surface_source_image_index"] = np.asarray( + res["surface_source_image_index"], dtype=np.int8) + if res["worst_pair"] > 0 and res["status"] in ( + calibration.STATUS_OK, calibration.STATUS_PARTIAL_IMAGES, + ): + out["img_comb_offset_dB"] = np.asarray(res["offsets"], dtype=np.float32) + out["img_comb_pair"][:] = res["worst_pair"] + summ = next(pp for pp in res["pairs"] if pp["pair"] == res["worst_pair"]) + if np.isfinite(summ["mean"]): + out["frame_offset"] = float(summ["mean"]) + except Exception: + logger.exception("Frame %s: calibration check failed; recording load_error", frame_id) + out["status"] = calibration.STATUS_LOAD_ERROR + logger.debug("Frame %s: calibration status=%s (%.1fs)", + frame_id, out["status"], time.time() - t0) + return out + + def process_frame(opr: OPRConnection, stac_item, config: dict) -> xr.Dataset | None: """Process a single radar frame and return a Dataset of metrics, or None on failure.""" + config = normalize_config(config) proc = config["processing"] qc_config = config.get("qc", {}) frame_id = stac_item.name if hasattr(stac_item, "name") else stac_item.get("id", "unknown") @@ -303,17 +370,7 @@ def process_frame(opr: OPRConnection, stac_item, config: dict) -> xr.Dataset | N frame = opr.load_frame(stac_item, data_product=proc["data_product"]) frame = frame.sortby("slow_time") - decimate_interval = proc.get("decimate_interval") - if decimate_interval: - interval = pd.Timedelta(decimate_interval) - times = frame.slow_time.values - selected = [0] - last = times[0] - for idx in range(1, len(times)): - if times[idx] - last >= interval: - selected.append(idx) - last = times[idx] - frame = frame.isel(slow_time=selected) + frame = decimate_frame(frame, proc.get("decimate_interval")) try: layers = opr.get_layers(frame, include_geometry=False) @@ -506,6 +563,23 @@ def process_frame(opr: OPRConnection, stac_item, config: dict) -> xr.Dataset | N qc_pass = qc_mask if qc_mask is not None else all_true + # Radiometric calibration (image-combine residual seam offsets + + # surface source image). Not QC-masking and never QC-masked: these are + # provenance/quality values downstream users filter on. + calib = run_frame_calibration(opr, stac_item, frame, proc, frame_id=frame_id) + # Attribute templates live in store.CALIBRATION_VAR_ATTRS (single + # source of truth with append auto-migration and backfill creation). + from .store import CALIBRATION_VAR_ATTRS + img_comb_offset = as_da(calib["img_comb_offset_dB"]).assign_attrs( + **CALIBRATION_VAR_ATTRS["img_comb_offset_dB"]) + img_comb_pair = as_da(calib["img_comb_pair"]).assign_attrs( + **CALIBRATION_VAR_ATTRS["img_comb_pair"]) + surface_source_image_index = as_da(calib["surface_source_image_index"]).assign_attrs( + **CALIBRATION_VAR_ATTRS["surface_source_image_index"]) + surface_ceiling_margin_dB = as_da( + np.full(len(frame.slow_time), np.nan, dtype=np.float32)).assign_attrs( + **CALIBRATION_VAR_ATTRS["surface_ceiling_margin_dB"]) + # Surface-side metrics are masked only by pick-independent QC so that # traces missing a bed pick (or failing thin-ice / bed-SNR checks) keep # the surface power and noise floor needed to treat the missing bed as @@ -543,6 +617,10 @@ def process_frame(opr: OPRConnection, stac_item, config: dict) -> xr.Dataset | N "bed_pick_attempted": bed_pick_attempted, "bed_pick_quality": bed_pick_quality, "record_end_twtt": record_end_twtt, + "img_comb_offset_dB": img_comb_offset, + "img_comb_pair": img_comb_pair, + "surface_source_image_index": surface_source_image_index, + "surface_ceiling_margin_dB": surface_ceiling_margin_dB, "frame_id": ("slow_time", [str(frame_id)] * len(frame.slow_time)), }, coords={ @@ -552,6 +630,10 @@ def process_frame(opr: OPRConnection, stac_item, config: dict) -> xr.Dataset | N ) ds.attrs["frame_bed_pick_fraction"] = float(bed_pick_available.values.mean()) ds.attrs["segment_bed_pick_fraction"] = segment_bed_pick_fraction + ds.attrs["frame_img_comb_status"] = calib["status"] + ds.attrs["frame_img_comb_weights_mode"] = calib["weights_mode"] + if calib["frame_offset"] is not None: + ds.attrs["frame_img_comb_offset_dB"] = calib["frame_offset"] if "Elevation" in frame: ds.coords["elevation"] = frame.Elevation diff --git a/src/radar_return_statistics/runner.py b/src/radar_return_statistics/runner.py index 117c8de..8d38138 100644 --- a/src/radar_return_statistics/runner.py +++ b/src/radar_return_statistics/runner.py @@ -1,10 +1,12 @@ import logging +import math import multiprocessing from concurrent.futures import ProcessPoolExecutor, as_completed import aiohttp import fsspec import xarray as xr +import zarr from xopr import OPRConnection from xopr import geometry as xopr_geometry @@ -15,6 +17,23 @@ logger = logging.getLogger(__name__) +def _finite_or_none(v): + f = float(v) + return f if math.isfinite(f) else None + + +# Per-frame attrs carried from process_frame output into frame_names-parallel +# root attribute lists, with a cast per attr (strings ride along uncast; NaN +# floats become None so the JSON attrs stay valid). +FRAME_ATTR_CASTS = { + "frame_bed_pick_fraction": float, + "segment_bed_pick_fraction": float, + "frame_img_comb_offset_dB": _finite_or_none, + "frame_img_comb_status": str, + "frame_img_comb_weights_mode": str, +} + + # fsspec's HTTPFileSystem inherits aiohttp's default ~5min total timeout, but # the default sock_read timeout (5 min) can fire on slow CReSIS downloads when # many workers compete for bandwidth on 200 MB+ frames. Set generous totals. @@ -205,10 +224,7 @@ def run(config_path: str | None = None, *, config: dict | None = None, reprocess # Per-checkpoint state — flushed at each checkpoint commit. batch_count = 0 batch_frame_collections: dict[str, str] = {} - batch_frame_attrs: dict[str, dict[str, float]] = { - "frame_bed_pick_fraction": {}, - "segment_bed_pick_fraction": {}, - } + batch_frame_attrs: dict[str, dict[str, object]] = {k: {} for k in FRAME_ATTR_CASTS} # Use 'spawn' so workers don't inherit thread state (icechunk's tokio # runtime, aiohttp clients, etc.) from the parent. Fork+threads can @@ -242,9 +258,9 @@ def run(config_path: str | None = None, *, config: dict | None = None, reprocess col = str(ds.attrs["collection"]) batch_frame_collections[fid] = col all_collections.add(col) - for attr_name in batch_frame_attrs: + for attr_name, cast in FRAME_ATTR_CASTS.items(): if attr_name in ds.attrs: - batch_frame_attrs[attr_name][fid] = float(ds.attrs[attr_name]) + batch_frame_attrs[attr_name][fid] = cast(ds.attrs[attr_name]) logger.info("Completed frame %s (%d/%d)", fid, total_frames_written, total_to_process) if checkpoint_every and batch_count >= checkpoint_every: @@ -304,3 +320,14 @@ def run(config_path: str | None = None, *, config: dict | None = None, reprocess message = f"[run] {summary}" store.commit_session(session, message) logger.info("Done: %s", message) + + # The saturation second pass fits season populations over the whole store; + # this run just changed the science data, so any prior results are stale. + try: + root = zarr.open_group(repo.readonly_session(branch="main").store, mode="r") + if store.saturation_stale(root): + logger.warning( + "Saturation second-pass results are stale after this run — re-run " + "`uv run python -m radar_return_statistics.saturation_pass `") + except Exception: + pass diff --git a/src/radar_return_statistics/saturation_pass.py b/src/radar_return_statistics/saturation_pass.py new file mode 100644 index 0000000..9b4d23a --- /dev/null +++ b/src/radar_return_statistics/saturation_pass.py @@ -0,0 +1,125 @@ +"""Saturation second pass: fit season-level ceilings over a completed store and +write per-trace surface_ceiling_margin_dB + season attrs back to it. + + uv run python -m radar_return_statistics.saturation_pass + +Fits the img1-sourced and higher-gain-sourced (index >= 2, one pooled +population reported as "img2") populations separately where +surface_source_image_index is populated, falling back to a single all-traces +fit otherwise; the cross-cap step is the season's empirical higher-gain bias +estimate. Margins are computed against the trace's own source population's +ceiling when that fit is fit_ok, else NaN. +""" +import logging +import math +from datetime import datetime, timezone + +import click +import numpy as np +import scipy.constants +import zarr + +from . import calibration as cal +from . import store as store_mod +from .config import load_config + +logger = logging.getLogger(__name__) + +# 0.4.1: two-sided flat-slope band (a rising envelope is no longer a ceiling) +# 0.4.0: piecewise partial-saturation model removed (flat rule only) +# 0.3.0: sparse-adjacent-bin merging in range binning +METHOD_VERSION = "0.4.1" + + +def _jsonable(obj): + """np scalars -> python, NaN -> None, tuples -> lists, recursively.""" + if isinstance(obj, dict): + return {k: _jsonable(v) for k, v in obj.items()} + if isinstance(obj, (list, tuple)): + return [_jsonable(v) for v in obj] + if isinstance(obj, (np.floating, float)): + f = float(obj) + return f if math.isfinite(f) else None + if isinstance(obj, np.integer): + return int(obj) + return obj + + +def run_saturation_pass(root: zarr.Group) -> tuple[np.ndarray, dict]: + """Compute per-trace margins and season attrs from store arrays.""" + power = root["surface_power_dB"][:] + twtt = root["surface_twtt"][:] + r = scipy.constants.c * twtt / 2.0 + qc_key = "qc_surface_pass" if "qc_surface_pass" in root else "qc_pass" + qc = root[qc_key][:].astype(bool) if qc_key in root else np.ones(power.shape, bool) + if "surface_source_image_index" in root: + src = root["surface_source_image_index"][:] + else: + src = np.full(power.shape, -1, dtype=np.int8) + + attrs = root.attrs.asdict() + colls = attrs.get("frame_collections") + if colls and "frame_index" in root: + seasons = np.asarray(colls)[root["frame_index"][:]] + else: + seasons = np.full(power.shape, "", dtype=object) + + margins = np.full(power.shape, np.nan, dtype=np.float32) + season_results = {} + for season in np.unique(seasons): + m = seasons == season + fit_mask = m & qc # fits use QC-passing traces only + season_fit = cal.fit_season_by_source(power[fit_mask], r[fit_mask], src[fit_mask]) + # margins for every trace of the season with finite power + margins[m] = cal.margins_by_source(power[m], src[m], season_fit).astype(np.float32) + season_results[str(season)] = season_fit + pops = season_fit["populations"] + logger.info("%s [%s]: %s; cross-cap step %s dB", season or "(no season)", + season_fit["ceiling_fit_population"], + {k: f"{v['status']}({v['level']:.1f})" if math.isfinite(v["level"]) + else v["status"] for k, v in pops.items()}, + f"{season_fit['cross_cap_step_db']:.1f}" + if math.isfinite(season_fit["cross_cap_step_db"]) else "n/a") + + saturation_attrs = _jsonable({ + "method_version": METHOD_VERSION, + "fitted": datetime.now(timezone.utc).isoformat(), + "params": cal.SAT_DEFAULTS, + "seasons": season_results, + }) + return margins, saturation_attrs + + +@click.command() +@click.argument("config_path", type=click.Path(exists=True)) +@click.option("--dry-run", is_flag=True, help="Fit and report, but write nothing") +@click.option("-v", "--verbose", is_flag=True) +def main(config_path: str, dry_run: bool, verbose: bool) -> None: + """Fit season saturation ceilings and write margins back to the store.""" + logging.basicConfig(level=logging.DEBUG if verbose else logging.INFO, + format="%(asctime)s %(levelname)s %(message)s") + config = load_config(config_path) + repo = store_mod.open_or_create_repo(config["store"]) + + if dry_run: + root = zarr.open_group(repo.readonly_session(branch="main").store, mode="r") + if "surface_power_dB" not in root: + raise click.ClickException("Store has no surface_power_dB — run the pipeline first") + run_saturation_pass(root) + click.echo("dry run: no store writes") + return + + # Fit and write through ONE session so the margins cannot be computed from + # one snapshot and stamped fresh against a concurrently updated one. + session = repo.writable_session("main") + root = zarr.open_group(session.store, mode="a") + if "surface_power_dB" not in root: + raise click.ClickException("Store has no surface_power_dB — run the pipeline first") + margins, saturation_attrs = run_saturation_pass(root) + store_mod.write_saturation_results(session, margins, saturation_attrs) + store_mod.commit_session(session, "[calibration] saturation second pass") + click.echo(f"wrote surface_ceiling_margin_dB ({np.isfinite(margins).sum()} finite margins)") + + +if __name__ == "__main__": + main() diff --git a/src/radar_return_statistics/store.py b/src/radar_return_statistics/store.py index ae241c5..13b098e 100644 --- a/src/radar_return_statistics/store.py +++ b/src/radar_return_statistics/store.py @@ -1,3 +1,4 @@ +import hashlib import logging import numpy as np @@ -8,6 +9,71 @@ logger = logging.getLogger(__name__) + +# Attribute templates for the calibration variables, applied wherever the +# arrays are created (fresh store writes, append auto-migration, backfill). +CALIBRATION_VAR_ATTRS = { + "img_comb_offset_dB": { + "units": "dB", + "description": ( + "Residual image-combine seam offset of the frame's worst image " + "pair — the power step present in the combined product at the " + "combine boundary, measured on the weight-corrected individual " + "images in their overlap. Positive = the earlier (lower-index, " + "shallower) image is brighter than the later one. NaN = " + "unmeasured; see frame_img_comb_status. Suggested downstream " + "rejection threshold ~3 dB."), + }, + "img_comb_pair": { + "description": ( + "Image pair img_comb_offset_dB refers to (1 = img1/img2, " + "2 = img2/img3), chosen once per frame as the pair with the " + "largest |mean offset|; -1 = undefined."), + }, + "surface_source_image_index": { + "description": ( + "Image that supplied the combined product's surface sample " + "(1 = low-gain img1; >=2 = a higher-gain image; -1 = unknown). " + "Provenance flag, not a validity verdict: higher-gain-sourced " + "surfaces are more likely saturated and may carry a " + "season-dependent low bias (see the season's cross_cap_step_db), " + "but for high-altitude operations they are the normal geometry."), + }, + "surface_ceiling_margin_dB": { + "units": "dB", + "description": ( + "Season ceiling of the trace's surface-source population minus " + "surface_power_dB. Small/zero => at the clip ceiling (likely " + "saturated). NaN => no credible ceiling was fitted for that " + "population (a clean season has no ceiling), or power undefined."), + }, +} + + +def apply_calibration_var_attrs(root: zarr.Group) -> int: + """Stamp the calibration attribute templates onto existing arrays that + lack them. Returns the number of arrays updated.""" + n = 0 + for name, attrs in CALIBRATION_VAR_ATTRS.items(): + if name in root and not dict(root[name].attrs): + root[name].attrs.update(attrs) + n += 1 + return n + + +def _fill_value_for(dtype: np.dtype): + """Fill value used when backfilling a variable that predates the store.""" + dtype = np.dtype(dtype) + if np.issubdtype(dtype, np.floating): + return np.nan + if np.issubdtype(dtype, np.bool_): + return False + if np.issubdtype(dtype, np.signedinteger): + return -1 + if dtype.kind in ("U", "S"): + return "" + return 0 + # Per-trace zarr arrays are appended to over the lifetime of the store. xarray's # default chunking on first write is the length of the first frame (~38 traces), # which makes the viewer fire thousands of HTTP requests per variable. Force a @@ -79,14 +145,24 @@ def _zarr_append(root: zarr.Group, ds: xr.Dataset) -> None: encoded_times, _, _ = encode_cf_datetime(ds.slow_time.values, units=units, calendar=calendar) arrays["slow_time"] = np.asarray(encoded_times) + existing_size = st_arr.shape[0] for name, data in arrays.items(): - if name in root: - arr = root[name] - old_size = arr.shape[0] - arr.resize(old_size + n_new) - arr[old_size:] = data - # Variables present in new frames but not yet in store (e.g. required_surface_snr_dB) - # are silently skipped — they won't appear until a --reprocess run. + if name not in root: + # Auto-migrate: a variable the store predates is created sized to + # the current store, backfilled with its fill value, then appended + # to like any other — never silently dropped. + fill = _fill_value_for(data.dtype) + logger.warning( + "Store migration: creating missing variable %r (%s), backfilling " + "%d existing traces with %r", name, data.dtype, existing_size, fill) + root.create_array(name, shape=(existing_size,), dtype=data.dtype, + chunks=(PER_TRACE_CHUNK_SIZE,), fill_value=fill) + if name in CALIBRATION_VAR_ATTRS: + root[name].attrs.update(CALIBRATION_VAR_ATTRS[name]) + arr = root[name] + old_size = arr.shape[0] + arr.resize(old_size + n_new) + arr[old_size:] = data def write_frame_results( @@ -169,7 +245,7 @@ def clear_store(session: icechunk.Session) -> None: def update_frame_index( session: icechunk.Session, frame_collections: dict[str, str] | None = None, - frame_scalar_attrs: dict[str, dict[str, float]] | None = None, + frame_scalar_attrs: dict[str, dict[str, object]] | None = None, ) -> None: """Rebuild frame_index (uint16 per trace) and frame_names root attribute from frame_id. @@ -237,3 +313,113 @@ def commit_session(session: icechunk.Session, message: str) -> str: snapshot_id = session.commit(message) logger.info("Committed: %s (snapshot: %s)", message, snapshot_id) return snapshot_id + + +def science_fingerprint(root: zarr.Group) -> dict: + """Fingerprint of the science data: trace count + hash of the processed + frame set. Science appends/removals/reprocessing change it; calibration-only + writes (saturation second pass, backfill) do not — so second-pass staleness + can be judged without being invalidated by its own commit.""" + n = int(root["slow_time"].shape[0]) if "slow_time" in root else 0 + frames = sorted(root["processed_frames"][:].tolist()) if "processed_frames" in root else [] + digest = hashlib.sha256("\n".join(frames).encode()).hexdigest() + return {"n_traces": n, "frames_sha256": digest} + + +def saturation_stale(root: zarr.Group) -> bool | None: + """True if the stored saturation second-pass results predate the current + science data OR were explicitly marked stale (e.g. a calibration backfill + rewrote surface_source_image_index, a fit input); None if no second pass + has been run.""" + meta = root.attrs.get("saturation") + if not meta: + return None + if meta.get("stale"): + return True + return meta.get("fingerprint") != science_fingerprint(root) + + +def mark_saturation_stale(root: zarr.Group) -> None: + """Flag existing saturation results as stale (cleared by the next + write_saturation_results). No-op when no second pass has run yet.""" + meta = root.attrs.get("saturation") + if meta and not meta.get("stale"): + root.attrs["saturation"] = {**meta, "stale": True} + + +def write_saturation_results( + session: icechunk.Session, + margins: np.ndarray, + saturation_attrs: dict, +) -> None: + """Write second-pass results: per-trace surface_ceiling_margin_dB plus the + root `saturation` attr (season fits, params, method version), stamped with + the current science fingerprint.""" + root = zarr.open_group(session.store, mode="a") + n = root["slow_time"].shape[0] + margins = np.asarray(margins, dtype=np.float32) + if margins.shape != (n,): + raise ValueError(f"margins shape {margins.shape} != store traces ({n},)") + attrs = dict(CALIBRATION_VAR_ATTRS["surface_ceiling_margin_dB"]) + if "surface_ceiling_margin_dB" in root: + attrs.update(dict(root["surface_ceiling_margin_dB"].attrs)) + root.create_array("surface_ceiling_margin_dB", data=margins, + chunks=(PER_TRACE_CHUNK_SIZE,), overwrite=True) + root["surface_ceiling_margin_dB"].attrs.update(attrs) + root.attrs["saturation"] = {**saturation_attrs, + "fingerprint": science_fingerprint(root)} + + +def update_frame_calibration( + session: icechunk.Session, + frame_id: str, + per_trace: dict[str, np.ndarray], + frame_attrs: dict[str, object], +) -> int: + """Rewrite one frame's calibration values in place (backfill of retryable + statuses) without touching science metrics. Returns traces updated (0 if + the frame is absent or ANY array's trace count mismatches — nothing is + written in that case, so the frame stays retryable). Marks existing + saturation results stale when surface_source_image_index (a saturation + fit input) is rewritten.""" + root = zarr.open_group(session.store, mode="a") + if "frame_id" not in root: + return 0 + idx = np.where(root["frame_id"][:] == frame_id)[0] + if idx.size == 0: + return 0 + # All-or-nothing: validate every array before any mutation. A partial + # write recorded as success would drop the frame from retry selection. + per_trace = {name: np.asarray(vals) for name, vals in per_trace.items()} + for name, vals in per_trace.items(): + if vals.shape != idx.shape: + logger.warning("Backfill %s: %s has %d values for %d traces; frame " + "left unwritten and retryable", + frame_id, name, vals.size, idx.size) + return 0 + n_total = root["frame_id"].shape[0] + for name, vals in per_trace.items(): + if name not in root: + fill = _fill_value_for(vals.dtype) + root.create_array(name, shape=(n_total,), dtype=vals.dtype, + chunks=(PER_TRACE_CHUNK_SIZE,), fill_value=fill) + if name in CALIBRATION_VAR_ATTRS: + root[name].attrs.update(CALIBRATION_VAR_ATTRS[name]) + arr = root[name] + if idx.size and np.all(np.diff(idx) == 1): + arr[idx[0]: idx[-1] + 1] = vals + else: + arr.oindex[idx] = vals + + frame_names = list(root.attrs.get("frame_names", []) or []) + if frame_id in frame_names: + pos = frame_names.index(frame_id) + for attr_name, value in frame_attrs.items(): + prev = root.attrs.get(attr_name) + lst = list(prev) if isinstance(prev, (list, tuple)) else [None] * len(frame_names) + lst += [None] * (len(frame_names) - len(lst)) + lst[pos] = value + root.attrs[attr_name] = lst + if "surface_source_image_index" in per_trace: + mark_saturation_stale(root) + return int(idx.size) diff --git a/tests/unit/test_calibration.py b/tests/unit/test_calibration.py new file mode 100644 index 0000000..9e9dcb6 --- /dev/null +++ b/tests/unit/test_calibration.py @@ -0,0 +1,366 @@ +import numpy as np +import pytest +import xarray as xr + +from radar_return_statistics import calibration as cal + + +# --------------------------------------------------------------------------- +# Synthetic frame builder: combined + 2 images with known per-image scalings +# (s1, s2) and a known seam step injected into the combined product. +# --------------------------------------------------------------------------- + +def _profile_db(t): + return -40.0 - 3.5e6 * t # smooth decay, dB + + +def _to_linear(db): + return 10.0 ** (db / 10.0) + + +def _dataset(data_db, twtt, surface, attrs, transpose=False): + n_traces = data_db.shape[0] + data = _to_linear(data_db) + if transpose: + data = data.T + dims = ("twtt", "slow_time") + else: + dims = ("slow_time", "twtt") + return xr.Dataset( + { + "Data": (dims, data), + "Surface": (("slow_time",), surface), + }, + coords={ + "slow_time": np.arange(n_traces).astype("datetime64[s]"), + "twtt": twtt, + }, + attrs=attrs, + ) + + +def make_synthetic(s1=0.0, s2=0.0, step=2.0, n_traces=100, img1_end=8e-6, + img_comb=(2e-6, -np.inf, 1e-6), combined_extra=None, + transpose_img2=False, imgs_entries=2): + """Combined product sections: [start, b) from img1+s1, [b, end] from + img2+s2-step (so residual seam step = +step, earlier image brighter).""" + tc = np.arange(0, 20e-6, 20e-9) + t1 = np.arange(0, img1_end, 20e-9) + t2 = np.arange(0, 20e-6, 25e-9) # deliberately different grid + surface = np.full(n_traces, 1.0e-6) + + A = _profile_db(tc) + # noise floors ride each image's own gain reference + img1_db = np.tile(np.maximum(_profile_db(t1) - s1, -130.0 - s1), (n_traces, 1)) + img2_db = np.tile(np.maximum(_profile_db(t2) - step - s2, -135.0 - s2), (n_traces, 1)) + + T_comb, T_blank, T_guard = img_comb + b = np.minimum(np.maximum(T_blank, T_comb + surface), t1[-1] - T_guard) + combined_db = np.empty((n_traces, tc.size)) + for r in range(n_traces): + sec2 = tc >= b[r] + combined_db[r] = A + combined_db[r, sec2] = A[sec2] - step + if combined_extra is not None: + combined_db[r, sec2] += combined_extra[r] + + attrs = { + "param_array": {"array": { + "img_comb": np.asarray(img_comb, dtype=float), + "imgs": [np.array([[i + 1.0], [1.0]]) for i in range(imgs_entries)], + }}, + "param_records": {"radar": {"wfs": {"Tpd": np.array([1e-6, 2e-6])}}}, + } + combined = _dataset(combined_db, tc, surface, attrs) + images = { + 1: _dataset(img1_db, t1, surface, {}), + 2: _dataset(img2_db, t2, surface, {}, transpose=transpose_img2), + } + return combined, images + + +class StubOPR: + def __init__(self, images, raises=None): + self.images = images + self.raises = raises or {} + + def load_frame(self, item, data_product=None, image=None, + allow_unlisted_products=False): + if image in self.raises: + raise self.raises[image] + if image in self.images: + return self.images[image] + raise FileNotFoundError(f"no image {image}") + + +def run_check(combined, images, raises=None): + return cal.check_img_combine(StubOPR(images, raises), None, combined, + "CSARP_standard") + + +# --------------------------------------------------------------------------- +# Check 1 tests +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("s2", [0.0, 107.6]) +def test_seam_step_recovered_independent_of_scalings(s2): + combined, images = make_synthetic(s1=1.5, s2=s2, step=2.0) + res = run_check(combined, images) + assert res["status"] == cal.STATUS_OK + assert res["weights"][1] == pytest.approx(1.5, abs=0.05) + assert res["weights"][2] == pytest.approx(s2, abs=0.05) + assert res["pairs"][0]["mean"] == pytest.approx(2.0, abs=0.15) + assert res["worst_pair"] == 1 + ok = np.isfinite(res["offsets"]) + assert ok.sum() >= 30 + assert np.nanmean(res["offsets"]) == pytest.approx(2.0, abs=0.15) + + +def test_differing_grids_and_orientation(): + combined, images = make_synthetic(step=3.0, transpose_img2=True) + res = run_check(combined, images) + assert res["status"] == cal.STATUS_OK + assert res["pairs"][0]["mean"] == pytest.approx(3.0, abs=0.15) + + +def test_weight_recovery_failed_on_inconsistent_sections(): + n = 100 + drift = np.linspace(0.0, 1.0, n) # w2 varies along the frame + combined, images = make_synthetic(step=2.0, combined_extra=drift) + res = run_check(combined, images) + assert res["status"] == cal.STATUS_WEIGHT_RECOVERY_FAILED + # raw offsets still computed for the parquet diagnostic + assert np.isfinite(res["offsets"]).sum() > 0 + + +def test_short_image_noise_floor_fallback(): + t_short = np.arange(0, 6e-6, 20e-9) + Di = np.tile(np.maximum(_profile_db(t_short), -130.0), (10, 1)) + floor = cal.image_noise_floor_db(t_short, Di, t_guard=1e-6) + # record shorter than the 12us window start -> percentile fallback, + # never the top-of-record signal median (~ -41 dB here) + assert np.all(floor < -55) + combined, images = make_synthetic(step=2.0, img1_end=6e-6) + res = run_check(combined, images) + assert res["status"] == cal.STATUS_OK + assert res["pairs"][0]["mean"] == pytest.approx(2.0, abs=0.2) + + +def test_signal_contaminated_tail_noise_floor(): + """A record long enough for the tail window whose end lies mid-ice + (2012 DC8 img1): the tail median sits in signal, so the floor must be + capped by the global low percentile, not the contaminated tail.""" + t = np.arange(0, 20e-6, 20e-9) + profile = np.full(t.size, -120.0) # true noise + profile[(t > 3e-6) & (t < 19e-6)] = -80.0 # englacial signal through the end + Di = np.tile(profile, (10, 1)) + floor = cal.image_noise_floor_db(t, Di, t_guard=2e-6) + # tail window [8, 13] us medians -80 (signal); pre-surface noise keeps the + # 10th percentile near -120 + assert np.all(floor < -110) + + +def test_insufficient_overlap(): + combined, images = make_synthetic(step=2.0, + img_comb=(2e-6, -np.inf, 4.95e-6)) + res = run_check(combined, images) + assert res["status"] == cal.STATUS_INSUFFICIENT_OVERLAP + assert res["pairs"][0]["status"] == cal.STATUS_INSUFFICIENT_OVERLAP + + +def test_transient_vs_missing_images(): + combined, images = make_synthetic() + res = run_check(combined, {1: images[1]}, + raises={2: ConnectionError("timeout")}) + assert res["status"] == cal.STATUS_LOAD_ERROR + res = run_check(combined, {1: images[1]}) # image 2 -> FileNotFoundError + assert res["status"] == cal.STATUS_IMAGES_UNAVAILABLE + assert res["image_errors"][2] == "missing" + + +def test_no_combine_single_image(): + combined, images = make_synthetic(imgs_entries=1) + res = run_check(combined, images) + assert res["status"] == cal.STATUS_NO_COMBINE + assert np.all(res["surface_source_image_index"] == 1) + + +def test_empty_img_comb_uses_tpd_defaults(): + combined, images = make_synthetic(step=2.0, img_comb=(2e-6, -np.inf, 1e-6)) + combined.attrs["param_array"]["array"]["img_comb"] = np.array([]) + res = run_check(combined, images) + # Tpd = [1us, 2us] -> default [T_comb=2us, -inf, T_guard=1us], same window + assert res["status"] == cal.STATUS_OK + assert "tpd_default" in res["img_comb_provenance"] + assert res["pairs"][0]["mean"] == pytest.approx(2.0, abs=0.15) + + +def test_surface_source_image_index(): + n = 60 + surface = np.full(n, 1e-6) + surface[40:] = 7.5e-6 # beyond img1 end (8us) - guard (1us) + surface[50:] = np.nan # NaN -> 0 -> img1 + gates = {1: (0.0, 8e-6), 2: (0.0, 20e-6)} + guards = {1: 1e-6} + idx = cal.surface_source_image_index(surface, gates, guards) + assert np.all(idx[:40] == 1) + assert np.all(idx[40:50] == 2) + assert np.all(idx[50:] == 1) + + +def test_tpd_for_images_wf_mapping(): + tpd = np.array([1e-6, 1e-6, 3e-6, 3e-6, 10e-6, 10e-6]) + imgs = [np.array([[1.0, 2.0], [1.0, 2.0]]), + np.array([[3.0, 4.0], [1.0, 2.0]]), + np.array([[5.0, 6.0], [1.0, 2.0]])] + tpds, consistent = cal.tpd_for_images(imgs, tpd) + assert tpds == [1e-6, 3e-6, 10e-6] + assert consistent + + +def test_2012_era_param_layout(): + """param_combine_wf_chan attr, old 2-per-pair img_comb, wfs as list of + dicts, imgs entries as (N, 2) wf-adc rows (2012_Antarctica_DC8 layout).""" + attrs = { + "param_combine_wf_chan": { + "img_comb": [1e-5, 2e-6], + "array": {"imgs": [ + np.array([[1, 1], [1, 2], [1, 3], [1, 4], [1, 5]]), + np.array([[2, 1], [2, 2], [2, 3], [2, 4], [2, 5]]), + ]}, + }, + "param_csarp": {"radar": {"wfs": [{"Tpd": 1e-6}, {"Tpd": 1e-5}]}}, + } + img_comb, prov = cal.find_img_comb(attrs) + np.testing.assert_allclose(img_comb, [1e-5, -np.inf, 2e-6]) + assert prov == "param_combine_wf_chan.img_comb" + tpd = cal.find_tpd(attrs) + np.testing.assert_allclose(tpd, [1e-6, 1e-5]) + imgs = cal.find_imgs(attrs) + assert len(imgs) == 2 + tpds, consistent = cal.tpd_for_images(imgs, tpd) + assert tpds == [1e-6, 1e-5] + assert consistent + + +# --------------------------------------------------------------------------- +# Check 2 (saturation) tests +# --------------------------------------------------------------------------- + +def _population(n=20000, clip=None, seed=0): + rng = np.random.default_rng(seed) + r = 10 ** rng.uniform(np.log10(300), np.log10(4000), n) + p = 20.0 - 20.0 * np.log10(r) + rng.normal(0, 2.0, n) + if clip is not None: + p = np.minimum(p, clip) + return p, r + + +def test_trailing_missing_image_partial_not_ok(): + """3 images expected, img3's file 404s: pair 2 was never assessed, so the + frame must be partial_images (with pair 1's offsets intact), never ok.""" + combined, images = make_synthetic(step=2.0, imgs_entries=3) + # declare a second combine transition + third waveform in the params + combined.attrs["param_array"]["array"]["img_comb"] = np.array( + [2e-6, -np.inf, 1e-6, 4e-6, -np.inf, 2e-6]) + combined.attrs["param_records"]["radar"]["wfs"]["Tpd"] = np.array( + [1e-6, 2e-6, 4e-6]) + res = run_check(combined, images) # StubOPR raises FileNotFoundError for img3 + assert res["status"] == cal.STATUS_PARTIAL_IMAGES + pair_statuses = {pr["pair"]: pr["status"] for pr in res["pairs"]} + assert pair_statuses[1] == cal.STATUS_OK + assert pair_statuses[2] == cal.STATUS_IMAGES_UNAVAILABLE + assert np.isfinite(res["offsets"]).sum() > 0 # measured pair still reported + + +def test_nan_trace_does_not_poison_weights(): + """One all-NaN trace in the combined product must not turn the image + weights NaN (weight_recovery_failed); the frame stays measurable.""" + combined, images = make_synthetic(s1=1.5, step=2.0) + data = combined["Data"].values + data[0, :] = np.nan + res = run_check(combined, images) + assert res["status"] == cal.STATUS_OK + assert np.isfinite(res["weights"][1]) and res["weights"][1] == pytest.approx(1.5, abs=0.05) + assert res["pairs"][0]["mean"] == pytest.approx(2.0, abs=0.2) + + +def test_clean_population_no_plateau(): + p, r = _population() + fit = cal.fit_ceiling(p, r) + assert fit["status"] == cal.NO_PLATEAU + margins = cal.ceiling_margin(p, fit) + assert np.all(np.isnan(margins)) + + +def test_dense_blob_thin_tail_clipped(): + """2012 DC8 img1 shape: ~98% of traces packed into a narrow range band, + a real ceiling-following tail so thin its equal-width bins individually + fail the per-bin floor. Sparse-bin merging must keep the tail's support + so the clipped population fits, with the full span reported.""" + rng = np.random.default_rng(1) + n_blob, n_tail = 19000, 350 + r = np.concatenate([ + 10 ** rng.uniform(np.log10(400), np.log10(700), n_blob), + 10 ** rng.uniform(np.log10(700), np.log10(3000), n_tail), + ]) + p = np.minimum(20.0 - 20.0 * np.log10(r) + rng.normal(0, 2.0, r.size), -33.0) + fit = cal.fit_ceiling(p, r) + assert fit["status"] == cal.FIT_OK + assert fit["level"] == pytest.approx(-33.0, abs=1.0) + # tail support retained (blob alone spans 0.24; merged tail bins extend it) + assert fit["span_decades"] > 0.3 + # same blob without the clip must not produce a ceiling + p_clean = 20.0 - 20.0 * np.log10(r) + rng.normal(0, 2.0, r.size) + assert cal.fit_ceiling(p_clean, r)["status"] != cal.FIT_OK + + +def test_rising_envelope_is_not_a_ceiling(): + """A steadily rising upper envelope is a gain/geometry artifact, not a + clip ceiling: the two-sided flat band must reject it.""" + rng = np.random.default_rng(3) + r = 10 ** rng.uniform(np.log10(300), np.log10(4000), 20000) + p = -80.0 + 20.0 * np.log10(r) + rng.normal(0, 2.0, r.size) # +20 dB/decade + fit = cal.fit_ceiling(p, r) + assert fit["status"] == cal.NO_PLATEAU + assert np.all(np.isnan(cal.ceiling_margin(p, fit))) + + +def test_fully_clipped_population(): + p, r = _population(clip=-45.0) + # clip at -45: unsaturated envelope is ~-30 at 300m, ~-52 at 4km -> flat + # over most of the range + fit = cal.fit_ceiling(p, r) + assert fit["status"] == cal.FIT_OK + assert fit["level"] == pytest.approx(-45.0, abs=1.0) + assert fit["pileup_fraction"] > 0.05 + margins = cal.ceiling_margin(p, fit) + assert np.isfinite(margins).all() + assert np.nanmin(margins) >= -1.5 + + +def test_partially_clipped_population(): + p, r = _population(clip=-52.0) + # envelope crosses -52 dB around r ~ 2km: plateau below, slope beyond + fit = cal.fit_ceiling(p, r) + assert fit["status"] == cal.FIT_OK + assert fit["level"] == pytest.approx(-52.0, abs=1.0) + + +def test_insufficient_support(): + p, r = _population(n=200) + fit = cal.fit_ceiling(p, r) + assert fit["status"] == cal.INSUFFICIENT_SUPPORT + assert np.all(np.isnan(cal.ceiling_margin(p, fit))) + + +def test_two_regime_step_not_misfit(): + # Two flat levels with a cliff between (two gain/altitude regimes, as in + # 2014/2018 Greenland P3) must not be reported as a single-season plateau: + # the second segment does not fall like an unsaturated envelope. + rng = np.random.default_rng(1) + n = 20000 + r = 10 ** rng.uniform(np.log10(300), np.log10(4000), n) + p = np.where(r < 1000, -31.0, -53.0) + rng.normal(0, 1.5, n) + fit = cal.fit_ceiling(p, r) + assert fit["status"] != cal.FIT_OK diff --git a/tests/unit/test_calibration_integration.py b/tests/unit/test_calibration_integration.py new file mode 100644 index 0000000..c92df0e --- /dev/null +++ b/tests/unit/test_calibration_integration.py @@ -0,0 +1,340 @@ +"""Step-5 integration tests: schema auto-migration, calibration placeholders and +enable/disable transitions, backfill-in-place, saturation second pass with +split populations, and science-fingerprint staleness.""" +import logging +import types + +import numpy as np +import pandas as pd +import pytest +import xarray as xr +import zarr + +from radar_return_statistics import calibration as cal +from radar_return_statistics import store as store_mod +from radar_return_statistics.processing import process_frame +from radar_return_statistics.saturation_pass import run_saturation_pass + + +def _make_ds(n_traces=5, frame_id="F", hour_offset=0, with_calibration=False, + offset_value=1.5): + slow_times = pd.date_range(f"2012-01-01 {hour_offset:02d}:00:00", + periods=n_traces, freq="10s") + data = { + "surface_twtt": ("slow_time", np.full(n_traces, 5e-6)), + "surface_power_dB": ("slow_time", np.full(n_traces, 10.0)), + "qc_pass": ("slow_time", np.ones(n_traces, dtype=bool)), + "frame_id": ("slow_time", [frame_id] * n_traces), + } + if with_calibration: + data.update({ + "img_comb_offset_dB": ("slow_time", + np.full(n_traces, offset_value, dtype=np.float32)), + "img_comb_pair": ("slow_time", np.ones(n_traces, dtype=np.int8)), + "surface_source_image_index": ("slow_time", np.ones(n_traces, dtype=np.int8)), + "surface_ceiling_margin_dB": ("slow_time", + np.full(n_traces, np.nan, dtype=np.float32)), + }) + return xr.Dataset(data, coords={ + "slow_time": slow_times, + "latitude": ("slow_time", np.full(n_traces, -75.0)), + "longitude": ("slow_time", np.full(n_traces, 160.0)), + }) + + +@pytest.fixture +def local_repo(tmp_path): + return store_mod.open_or_create_repo({"backend": "local", "path": str(tmp_path / "s")}) + + +def _root(repo): + return zarr.open_group(repo.readonly_session(branch="main").store, mode="r") + + +# --------------------------------------------------------------------------- +# Auto-migration on append +# --------------------------------------------------------------------------- + +def test_append_to_precalibration_store_automigrates(local_repo, caplog): + session = local_repo.writable_session("main") + store_mod.write_frame_results(session, "F1", _make_ds(5, "F1", 0)) + with caplog.at_level(logging.WARNING): + store_mod.write_frame_results( + session, "F2", _make_ds(3, "F2", 1, with_calibration=True)) + store_mod.commit_session(session, "test") + + assert any("Store migration" in r.message for r in caplog.records) + root = _root(local_repo) + off = root["img_comb_offset_dB"][:] + assert off.shape == (8,) + assert np.all(np.isnan(off[:5])) # backfilled fill value + assert np.allclose(off[5:], 1.5) + pair = root["img_comb_pair"][:] + assert np.all(pair[:5] == -1) and np.all(pair[5:] == 1) + assert np.all(root["surface_source_image_index"][:5] == -1) + + +def test_enabled_disabled_enabled_transitions(local_repo): + # Disabled frames still emit placeholder columns (schema stays uniform), + # so appends in any enabled/disabled order line up. + session = local_repo.writable_session("main") + store_mod.write_frame_results(session, "F1", _make_ds(2, "F1", 0, True, 2.0)) + disabled = _make_ds(2, "F2", 1, True) + disabled["img_comb_offset_dB"][:] = np.nan + disabled["img_comb_pair"][:] = -1 + disabled["surface_source_image_index"][:] = -1 + store_mod.write_frame_results(session, "F2", disabled) + store_mod.write_frame_results(session, "F3", _make_ds(2, "F3", 2, True, 3.0)) + store_mod.commit_session(session, "test") + + root = _root(local_repo) + off = root["img_comb_offset_dB"][:] + assert np.allclose(off[[0, 1]], 2.0) and np.allclose(off[[4, 5]], 3.0) + assert np.all(np.isnan(off[[2, 3]])) + + +# --------------------------------------------------------------------------- +# Processing emits placeholders / real columns +# --------------------------------------------------------------------------- + +def _minimal_config(img_combine): + return { + "processing": { + "data_product": "CSARP_standard", "decimate_interval": None, + "layer_margin_m": 50, "ice_permittivity": 3.17, "max_workers": 1, + "calibration": {"img_combine": img_combine}, + }, + "qc": {"max_heading_change_deg_per_km": None, "min_ice_thickness_m": None, + "min_agl_m": None, "min_bed_snr_db": None, "min_traces_after_qc": 1}, + } + + +def test_process_frame_calibration_disabled_placeholders( + mocker, synthetic_frame, synthetic_layers): + opr = mocker.MagicMock() + opr.load_frame.return_value = synthetic_frame + opr.get_layers.return_value = synthetic_layers + ds = process_frame(opr, types.SimpleNamespace(name="F"), _minimal_config(False)) + assert np.all(np.isnan(ds["img_comb_offset_dB"].values)) + assert np.all(ds["img_comb_pair"].values == -1) + assert np.all(ds["surface_source_image_index"].values == -1) + assert ds.attrs["frame_img_comb_status"] == cal.STATUS_DISABLED + assert "frame_img_comb_offset_dB" not in ds.attrs + # disabled: the (mock) opr must not be asked for image products + assert not any(kw.get("image") for _, kw in opr.load_frame.call_args_list) + + +def test_process_frame_calibration_enabled_degrades_to_status( + mocker, synthetic_frame, synthetic_layers): + # Synthetic frames carry no OPR params: the check must degrade to a status + # (params_missing), never fail the frame. + opr = mocker.MagicMock() + opr.load_frame.return_value = synthetic_frame + opr.get_layers.return_value = synthetic_layers + ds = process_frame(opr, types.SimpleNamespace(name="F"), _minimal_config(True)) + assert ds is not None + assert ds.attrs["frame_img_comb_status"] == cal.STATUS_PARAMS_MISSING + assert np.all(np.isnan(ds["img_comb_offset_dB"].values)) + + +# --------------------------------------------------------------------------- +# Backfill in place +# --------------------------------------------------------------------------- + +def test_update_frame_calibration_in_place(local_repo): + session = local_repo.writable_session("main") + store_mod.write_frame_results(session, "F1", _make_ds(3, "F1", 0, True, np.nan)) + store_mod.write_frame_results(session, "F2", _make_ds(4, "F2", 1, True, np.nan)) + store_mod.update_frame_index(session) + root_w = zarr.open_group(session.store, mode="a") + root_w.attrs["frame_img_comb_status"] = ["load_error", "load_error"] + store_mod.commit_session(session, "initial") + + session = local_repo.writable_session("main") + n = store_mod.update_frame_calibration( + session, "F2", + per_trace={ + "img_comb_offset_dB": np.full(4, 0.7, dtype=np.float32), + "img_comb_pair": np.full(4, 1, dtype=np.int8), + "surface_source_image_index": np.full(4, 2, dtype=np.int8), + }, + frame_attrs={"frame_img_comb_status": "ok", + "frame_img_comb_offset_dB": 0.7}, + ) + store_mod.commit_session(session, "backfill") + assert n == 4 + + root = _root(local_repo) + off = root["img_comb_offset_dB"][:] + assert np.all(np.isnan(off[:3])) and np.allclose(off[3:], 0.7) + names = list(root.attrs["frame_names"]) + statuses = list(root.attrs["frame_img_comb_status"]) + assert dict(zip(names, statuses)) == {"F1": "load_error", "F2": "ok"} + offsets_attr = list(root.attrs["frame_img_comb_offset_dB"]) + assert offsets_attr[names.index("F2")] == pytest.approx(0.7) + + +# --------------------------------------------------------------------------- +# Split-population second pass + cross-cap step +# --------------------------------------------------------------------------- + +def _two_population_arrays(n=30000, seed=0, level1=-31.0, level2=-53.0): + """img1 clipped at level1 below the 1 km cap, img2 clipped at level2 beyond + — the resolved two-regime structure with a known cross-cap step.""" + rng = np.random.default_rng(seed) + r = 10 ** rng.uniform(np.log10(300), np.log10(4000), n) + src = np.where(r < 1000, 1, 2).astype(np.int8) + # a clip truncates from above: noise spreads the population downward only + p = np.where(src == 1, level1, level2) - np.abs(rng.normal(0, 1.0, n)) + return p, r, src + + +def test_fit_season_by_source_split_populations(): + p, r, src = _two_population_arrays() + fit = cal.fit_season_by_source(p, r, src) + assert fit["ceiling_fit_population"] == "by_source_image" + assert fit["populations"]["img1"]["status"] == cal.FIT_OK + assert fit["populations"]["img1"]["level"] == pytest.approx(-31.0, abs=1.5) + assert fit["populations"]["img2"]["status"] == cal.FIT_OK + assert fit["populations"]["img2"]["level"] == pytest.approx(-53.0, abs=1.5) + assert fit["cross_cap_step_db"] == pytest.approx(22.0, abs=3.0) + + margins = cal.margins_by_source(p, src, fit) + assert np.isfinite(margins).all() + # each trace's margin is against its own population's ceiling + assert np.median(margins[src == 1]) == pytest.approx( + fit["populations"]["img1"]["level"] - np.median(p[src == 1]), abs=0.3) + + +def test_fit_season_by_source_falls_back_without_index(): + p, r, _ = _two_population_arrays() + fit = cal.fit_season_by_source(p, r, np.full(p.shape, -1, dtype=np.int8)) + assert fit["ceiling_fit_population"] == "all_traces" + assert "all" in fit["populations"] + + +def test_saturation_pass_writes_margins_and_attrs(local_repo): + n = 12000 + p, r, src = _two_population_arrays(n=n, seed=1) + twtt = 2.0 * r / 299792458.0 + slow_times = pd.date_range("2012-01-01", periods=n, freq="s") + ds = xr.Dataset( + { + "surface_twtt": ("slow_time", twtt), + "surface_power_dB": ("slow_time", p), + "qc_pass": ("slow_time", np.ones(n, dtype=bool)), + "surface_source_image_index": ("slow_time", src), + "frame_id": ("slow_time", ["F1"] * n), + }, + coords={"slow_time": slow_times, + "latitude": ("slow_time", np.full(n, -75.0)), + "longitude": ("slow_time", np.full(n, 160.0))}, + ) + session = local_repo.writable_session("main") + store_mod.write_frame_results(session, "F1", ds) + store_mod.update_frame_index(session, frame_collections={"F1": "test_season"}) + store_mod.commit_session(session, "data") + + root = _root(local_repo) + margins, sat_attrs = run_saturation_pass(root) + season = sat_attrs["seasons"]["test_season"] + assert season["ceiling_fit_population"] == "by_source_image" + assert season["cross_cap_step_db"] == pytest.approx(22.0, abs=3.0) + assert np.isfinite(margins).sum() > 0.9 * n + + session = local_repo.writable_session("main") + store_mod.write_saturation_results(session, margins, sat_attrs) + store_mod.commit_session(session, "second pass") + root = _root(local_repo) + assert np.isfinite(root["surface_ceiling_margin_dB"][:]).sum() > 0.9 * n + assert root.attrs["saturation"]["seasons"]["test_season"]["cross_cap_step_db"] is not None + + +# --------------------------------------------------------------------------- +# Fingerprint staleness +# --------------------------------------------------------------------------- + +def test_fingerprint_staleness(local_repo): + session = local_repo.writable_session("main") + store_mod.write_frame_results(session, "F1", _make_ds(5, "F1", 0, True)) + store_mod.commit_session(session, "science") + + root = _root(local_repo) + assert store_mod.saturation_stale(root) is None # no second pass yet + + # Second pass -> fresh (its own write must not look stale) + session = local_repo.writable_session("main") + store_mod.write_saturation_results( + session, np.full(5, np.nan, dtype=np.float32), {"seasons": {}}) + store_mod.commit_session(session, "second pass") + assert store_mod.saturation_stale(_root(local_repo)) is False + + # Calibration-only backfill write -> still fresh + session = local_repo.writable_session("main") + store_mod.update_frame_calibration( + session, "F1", + per_trace={"img_comb_offset_dB": np.full(5, 0.1, dtype=np.float32)}, + frame_attrs={}) + store_mod.commit_session(session, "backfill") + assert store_mod.saturation_stale(_root(local_repo)) is False + + # Science append -> stale + session = local_repo.writable_session("main") + store_mod.write_frame_results(session, "F2", _make_ds(3, "F2", 1, True)) + store_mod.commit_session(session, "more science") + assert store_mod.saturation_stale(_root(local_repo)) is True + + +def test_backfill_source_index_marks_saturation_stale(local_repo): + """A backfill that rewrites surface_source_image_index (a saturation fit + input) must invalidate existing saturation results; the next second pass + clears the flag.""" + session = local_repo.writable_session("main") + store_mod.write_frame_results(session, "F1", _make_ds(5, "F1", 0, True)) + store_mod.write_saturation_results( + session, np.full(5, np.nan, dtype=np.float32), {"seasons": {}}) + store_mod.commit_session(session, "science + second pass") + assert store_mod.saturation_stale(_root(local_repo)) is False + + session = local_repo.writable_session("main") + n = store_mod.update_frame_calibration( + session, "F1", + per_trace={"surface_source_image_index": np.full(5, 2, dtype=np.int8)}, + frame_attrs={}) + store_mod.commit_session(session, "backfill with source index") + assert n == 5 + assert store_mod.saturation_stale(_root(local_repo)) is True + + session = local_repo.writable_session("main") + store_mod.write_saturation_results( + session, np.full(5, np.nan, dtype=np.float32), {"seasons": {}}) + store_mod.commit_session(session, "re-run second pass") + assert store_mod.saturation_stale(_root(local_repo)) is False + + +def test_update_frame_calibration_mismatch_writes_nothing(local_repo): + """Any array-length mismatch aborts the whole frame update (no partial + write recorded as success), leaving the frame retryable.""" + session = local_repo.writable_session("main") + store_mod.write_frame_results(session, "F1", _make_ds(4, "F1", 0, True, np.nan)) + store_mod.update_frame_index(session) + root_w = zarr.open_group(session.store, mode="a") + root_w.attrs["frame_img_comb_status"] = ["load_error"] + store_mod.commit_session(session, "initial") + + session = local_repo.writable_session("main") + n = store_mod.update_frame_calibration( + session, "F1", + per_trace={ + "img_comb_offset_dB": np.full(4, 0.7, dtype=np.float32), # correct + "surface_source_image_index": np.full(3, 1, dtype=np.int8), # WRONG length + }, + frame_attrs={"frame_img_comb_status": "ok"}, + ) + assert n == 0 + # nothing to commit — the session must hold zero changes + with pytest.raises(Exception, match="no changes"): + store_mod.commit_session(session, "attempted backfill") + root = _root(local_repo) + assert np.all(np.isnan(root["img_comb_offset_dB"][:])) # nothing written + assert list(root.attrs["frame_img_comb_status"]) == ["load_error"] # retryable diff --git a/tests/unit/test_processing.py b/tests/unit/test_processing.py index 22f38af..15a7a18 100644 --- a/tests/unit/test_processing.py +++ b/tests/unit/test_processing.py @@ -113,6 +113,8 @@ def test_process_frame_output_variables(mocker, synthetic_frame, synthetic_layer "record_end_twtt", "qc_pass", "qc_surface_pass", "qc_heading_pass", "qc_agl_pass", "bed_pick_available", "bed_pick_attempted", "bed_pick_quality", "frame_id", + "img_comb_offset_dB", "img_comb_pair", "surface_source_image_index", + "surface_ceiling_margin_dB", } assert ds.attrs["frame_bed_pick_fraction"] == 1.0 assert ds.attrs["segment_bed_pick_fraction"] == 1.0 diff --git a/uv.lock b/uv.lock index 2dd3558..86911d2 100644 --- a/uv.lock +++ b/uv.lock @@ -2366,6 +2366,7 @@ dependencies = [ { name = "matplotlib" }, { name = "numpy" }, { name = "pandas" }, + { name = "pyarrow" }, { name = "pyproj" }, { name = "pyyaml" }, { name = "scipy" }, @@ -2391,6 +2392,7 @@ requires-dist = [ { name = "matplotlib" }, { name = "numpy" }, { name = "pandas" }, + { name = "pyarrow", specifier = ">=23.0.1" }, { name = "pyproj" }, { name = "pytest", marker = "extra == 'test'", specifier = ">=8" }, { name = "pytest-mock", marker = "extra == 'test'", specifier = ">=3" }, diff --git a/web/index.html b/web/index.html index a14ddc2..a0a4ea0 100644 --- a/web/index.html +++ b/web/index.html @@ -399,6 +399,9 @@

Variable

+ + + diff --git a/web/src/config.ts b/web/src/config.ts index a105243..09f4952 100644 --- a/web/src/config.ts +++ b/web/src/config.ts @@ -50,6 +50,14 @@ export const BED_SIDE_VARIABLES = new Set([ "post_bed_noise_dB", ]); +// int8 categorical variables using -1 as the "unknown" sentinel; loaded via +// element-wise integer conversion and shown as integers (never through the +// float byte-reinterpret path). +export const INT_SENTINEL_VARIABLES = new Set([ + "surface_source_image_index", + "img_comb_pair", +]); + export interface VariableInfo { label: string; cmap: string; @@ -57,6 +65,8 @@ export interface VariableInfo { // Multiplier applied when formatting values for display (legend + tooltip). // Stored data is unchanged; e.g. TWTT is stored in seconds but shown as µs. displayScale?: number; + // Categorical/integer variable: format values without decimals. + integer?: boolean; } export const VARIABLES: Record = { @@ -115,6 +125,22 @@ export const VARIABLES: Record = { cmap: "turbo", unit: "dB", }, + img_comb_offset_dB: { + label: "Img-Combine Seam Offset", + cmap: "turbo", + unit: "dB", + }, + surface_source_image_index: { + label: "Surface Source Image", + cmap: "turbo", + unit: "", + integer: true, + }, + surface_ceiling_margin_dB: { + label: "Surface Ceiling Margin", + cmap: "turbo", + unit: "dB", + }, surface_twtt: { label: "Surface TWTT", cmap: "turbo", diff --git a/web/src/map.ts b/web/src/map.ts index b622abe..fdea3f5 100644 --- a/web/src/map.ts +++ b/web/src/map.ts @@ -257,6 +257,7 @@ export function tracesInPolygon(data: StoreData): number[] { export function formatScaledValue(value: number, info: VariableInfo): string { const scaled = (info.displayScale ?? 1) * value; + if (info.integer) return Math.round(scaled).toString(); return scaled.toPrecision(4); } diff --git a/web/src/store.ts b/web/src/store.ts index 11a00f1..edf6da6 100644 --- a/web/src/store.ts +++ b/web/src/store.ts @@ -1,6 +1,6 @@ import { IcechunkStore } from "@carbonplan/icechunk-js"; import * as zarr from "zarrita"; -import { VARIABLE_SOURCE } from "./config"; +import { INT_SENTINEL_VARIABLES, VARIABLE_SOURCE } from "./config"; export interface StoreData { latitude: Float64Array; @@ -59,9 +59,10 @@ async function loadArray( function toFloat64Array(chunk: zarr.Chunk): Float64Array { const data = chunk.data; if (data instanceof Float64Array) return data; - if (data instanceof Float32Array) return new Float64Array(data); + // Element-wise conversion: integer typed arrays (e.g. the int8 categorical + // variables) must be converted by value, never byte-reinterpreted. if (ArrayBuffer.isView(data)) - return new Float64Array(data.buffer, data.byteOffset, data.byteLength / 8); + return Float64Array.from(data as unknown as ArrayLike, Number); return new Float64Array(data as unknown as ArrayLike); } @@ -161,7 +162,14 @@ export async function loadVariables( const source = VARIABLE_SOURCE[name] ?? name; try { const chunk = await loadArray(store, source); - return [name, toFloat64Array(chunk)] as const; + const values = toFloat64Array(chunk); + // Integer variables use -1 as "unknown"; map it to NaN so the + // existing missing-value handling applies. + if (INT_SENTINEL_VARIABLES.has(name)) { + for (let i = 0; i < values.length; i++) + if (values[i] === -1) values[i] = NaN; + } + return [name, values] as const; } catch { return null; }