prop/docs/plans/2026-04-15-merra2-historical-backfill.md

# NARR Historical Backfill Implementation Plan

> **For Claude:** REQUIRED SUB-SKILL: Use superpowers:executing-plans to implement this plan task-by-task.

> **History:** This file originally specified a NASA MERRA-2 path. While executing
> Phase 1 it turned out the MERRA-2 AWS bucket (`gesdisc-cumulus-prod-protected`)
> requires Earthdata credentials, contradicting the initial research. While
> investigating alternatives I discovered NCEP NARR is anonymously available at
> NCEI with native GRIB1 + wgrib1-format `.inv` files. Plan was rewritten around
> NARR. A second correction landed during the spike: wgrib2 cannot decode GRIB1
> messages (`rd_grib2_seq_msg, grib1 message ignored, use wgrib`), so we *do*
> need a separate decoder. cdo handles GRIB1 + Lambert-conformal projection
> + nearest-neighbor remap natively, so the plan now uses cdo for decoding.
> Filename kept for git continuity.

**Goal:** Replace the broken ERA5/CDS historical backfill with NCEP NARR fetched anonymously from NCEI, parsed via the existing `wgrib2` pipeline, and inserted into the existing `era5_profiles` table.

**Architecture:** New `Microwaveprop.Weather.NarrClient` mirrors the existing `HrrrClient` byte-range fetch pattern: download the small `.inv` file, parse it into a map of `{var, level} => {byte_offset, length}`, issue HTTP `Range: bytes=START-END` GETs for the ~10 records we need, write each record to a per-record temp file, then run `cdo -merge` to combine them into a single GRIB1 and `cdo -outputtab,name,lev,value -remapnn,lon=X_lat=Y` to extract values at the QSO point. New `NarrFetchWorker` is a single-stage Oban worker (no submit/poll split — there's no upstream queue) that takes a `contact_id`, picks the nearest 3-hourly analysis hour to the QSO timestamp, calls `NarrClient`, and inserts an `era5_profiles` row. `ContactWeatherEnqueueWorker` dispatches `NarrFetchWorker` for pre-2014 contacts in place of `Era5FetchWorker`. The `era5_profiles` table is reused 1:1 — schema rename can be a follow-up if we feel like it.

**Tech Stack:** Elixir, Ecto, Oban, `Req` for HTTPS byte-range fetches, **cdo** (Climate Data Operators) CLI for GRIB1 record handling and Lambert-conformal nearest-neighbor extraction. The existing `wgrib2` CLI does *not* decode GRIB1 messages so we add `cdo` to the Dockerfile (one new system dep, ~30–60 MB image growth via libnetcdf/libhdf5/libudunits2 transitive deps).

**cdo parameter-number → semantic-name lookup (verified empirically against NARR 2010-06-15 12Z):**

| cdo `name` | NCEP GRIB1 param | Meaning | Where in NARR |
|---|---|---|---|
| `var1` | 1 | PRES (Pressure) | sfc / 10 m / hybrid |
| `var7` | 7 | HGT (Geopotential height) | pressure levels |
| `var11` | 11 | TMP (Temperature) | sfc / 2 m / 10 m / pressure levels |
| `var17` | 17 | DPT (Dew point) | 2 m only |
| `var33` | 33 | UGRD (u wind) | 10 m / pressure levels |
| `var34` | 34 | VGRD (v wind) | 10 m / pressure levels |
| `var51` | 51 | SPFH (Specific humidity) | 2 m / pressure levels |
| `var54` | 54 | PWAT (Precipitable water) | atmos col |
| `var221` | 221 | HPBL (PBL height) | sfc |

cdo's `lev` column for pressure-level records is in **Pa** (85000 = 850 mb). For 2 m above ground the `lev` is `2`, for 10 m the `lev` is `10`, for surface/atmos col the `lev` is `0`. The parser must use *both* the var ID and the lev value to disambiguate (e.g. `var11 lev=2` is 2 m TMP, `var11 lev=85000` is 850 mb TMP).

**Why NARR over MERRA-2 / ERA5:**
- Native GRIB1 with wgrib1-format `.inv` (already verified — `narr-a_221_20100615_1200_000.inv` parses cleanly)
- `Accept-Ranges: bytes` confirmed on `https://www.ncei.noaa.gov/data/north-american-regional-reanalysis/access/3-hourly/...` — same byte-range fetch pattern `HrrrClient` already uses
- Anonymous HTTPS — no auth, no quota, no job queue
- 32 km Lambert conformal grid covering CONUS (Grid 221), 3-hourly, 1979-01-01 → 2014-10-02 — exactly the pre-HRRR window
- DPT is a direct surface dewpoint record, no specific-humidity → Td derivation
- 56 MB per analysis hour total, but a per-point fetch only pulls ~10 records ≈ 1.5 MB max
- Same atmospheric variables we need on the same surface + pressure-level layout as the existing `era5_profiles` schema

**Spike findings (2026-04-15, verified by hand):**
- For 2010-06-15 12Z at DFW (32.9°N, -97.0°W): TMP_2m=296.94 K (23.8°C), DPT_2m=294.49 K (21.3°C), HPBL=703.5 m, PRES_sfc=99383.6 Pa (993.8 mb), PWAT=45.1 mm, TMP_850mb=291.52 K (18.4°C), HGT_850mb=1547.8 m, SPFH_850mb=10.4 g/kg — all physically sane Texas summer values
- cdo `-merge` is the right operator for combining GRIB1 records of different variables. cdo `-cat` is for time concatenation only and silently drops records with mismatched levels. Shell `cat` produces an invalid composite (`GRIB message error`). Always use `cdo -merge`.

**Verified facts (gathered live during plan rewrite):**
- File URL pattern: `https://www.ncei.noaa.gov/data/north-american-regional-reanalysis/access/3-hourly/{YYYYMM}/{YYYYMMDD}/narr-a_221_{YYYYMMDD}_{HHMM}_000.grb`
- Inventory file: same path with `.inv` extension
- HHMM ∈ {0000, 0300, 0600, 0900, 1200, 1500, 1800, 2100} (3-hourly analyses)
- Inventory format: `record:offset:D=YYYYMMDDHH:VAR:level:kpds=...:anl:...:"description [units]`
- Records we need (verified present in real inventory for 2010-06-15 12Z):
  - 15: `HPBL:sfc`
  - 265: `PRES:sfc`
  - 288: `TMP:2 m above gnd`
  - 290: `DPT:2 m above gnd`
  - 293: `UGRD:10 m above gnd`
  - 294: `VGRD:10 m above gnd`
  - 317: `PWAT:atmos col`
  - Pressure levels 100, 125, 150, 175, 200, 225, 250, 275, 300, 350, 400, 450, 500, 550, 600, 650, 700, 725, 750, 775, 800, 825, 850, 875, 900, 925, 950, 975, 1000 mb each have `TMP`, `HGT`, `SPFH` records

---

## Phase 1: Spike — DONE 2026-04-15

Verified by hand. Findings recorded in the spike-findings block above. Key deltas from the original assumption:

- wgrib2 cannot read GRIB1 messages — it errors out with `rd_grib2_seq_msg, grib1 message ignored, use wgrib`.
- cdo handles GRIB1 + Lambert-conformal projection + nearest-neighbor remap natively. `cdo -merge` for combining records, `cdo -outputtab,name,lev,value -remapnn,lon=X_lat=Y` for point extraction.
- cdo emits variable names as `varNNN` (NCEP GRIB1 parameter numbers). Mapping table in the architecture section.
- All 5 surface vars + 3 pressure-level vars at 850 mb extract to sane physical values for DFW June 12Z.
- Adding cdo to the Dockerfile is now a Phase 2 prerequisite (see Task 2.5).

---

## Phase 2: Build NarrClient (TDD)

### Task 2.5 (PREREQUISITE): Add cdo to Dockerfile runtime stage

**Files:**
- Modify: `Dockerfile:137-141` (final runtime apt-get install)

Add `cdo` to the `apt-get install` line. cdo on Debian trixie pulls in libnetcdf, libhdf5, libudunits2 transitively. Verify with `docker build --target final` then `docker run --rm <image> cdo --version`.

```dockerfile
RUN apt-get update \
  && apt-get install -y --no-install-recommends \
     libstdc++6 openssl libncurses6 locales ca-certificates snmp \
     libgfortran5 libaec0 zlib1g libpng16-16t64 libopenjp2-7 \
     cdo \
  && rm -rf /var/lib/apt/lists/*
```

This is a prerequisite because every other Phase 2/3 task depends on `cdo` being available at runtime. Local dev already has cdo from `brew install cdo` during the spike — no local change needed there. Commit this Dockerfile change before the NarrClient tasks so prod CI builds with the new image.

### Task 2: NarrClient.url_for/2 + url_for_inventory/2

**Files:**
- Create: `lib/microwaveprop/weather/narr_client.ex`
- Create: `test/microwaveprop/weather/narr_client_test.exs`

**Step 1: Failing test**

```elixir
defmodule Microwaveprop.Weather.NarrClientTest do
  use ExUnit.Case, async: true
  alias Microwaveprop.Weather.NarrClient

  describe "url_for/1" do
    test "builds the NCEI NARR HTTPS URL for a given analysis time" do
      vt = ~U[2010-06-15 12:00:00Z]
      assert NarrClient.url_for(vt) ==
        "https://www.ncei.noaa.gov/data/north-american-regional-reanalysis/access/3-hourly/201006/20100615/narr-a_221_20100615_1200_000.grb"
    end

    test "raises on non-3-hourly analysis times (NARR is 00/03/06/.../21Z only)" do
      vt = ~U[2010-06-15 13:00:00Z]
      assert_raise ArgumentError, fn -> NarrClient.url_for(vt) end
    end
  end

  describe "url_for_inventory/1" do
    test "swaps .grb for .inv" do
      vt = ~U[2010-06-15 12:00:00Z]
      assert NarrClient.url_for_inventory(vt) ==
        "https://www.ncei.noaa.gov/data/north-american-regional-reanalysis/access/3-hourly/201006/20100615/narr-a_221_20100615_1200_000.inv"
    end
  end

  describe "snap_to_analysis_hour/1" do
    test "rounds an arbitrary timestamp down to the nearest 3-hourly NARR slot" do
      assert NarrClient.snap_to_analysis_hour(~U[2010-06-15 13:42:00Z]) == ~U[2010-06-15 12:00:00Z]
      assert NarrClient.snap_to_analysis_hour(~U[2010-06-15 14:30:00Z]) == ~U[2010-06-15 12:00:00Z]
      assert NarrClient.snap_to_analysis_hour(~U[2010-06-15 15:00:00Z]) == ~U[2010-06-15 15:00:00Z]
      assert NarrClient.snap_to_analysis_hour(~U[2010-06-15 23:59:59Z]) == ~U[2010-06-15 21:00:00Z]
    end
  end
end
```

**Step 2:** Run, see fail.

**Step 3:** Implement minimally.

**Step 4:** Run, green.

**Step 5: Commit**

### Task 3: NarrClient.fetch_inventory/1

**Files:** modify `lib/microwaveprop/weather/narr_client.ex` + test file.

**Step 1: Failing test** using a fixture inventory file.

Capture a real NARR `.inv` file from Phase 1 and commit it as a fixture:

```bash
cp /tmp/narr.inv test/fixtures/narr/narr-a_221_20100615_1200_000.inv
```

Then test the parser:

```elixir
describe "parse_inventory/2" do
  test "parses the inv text into a map of {var, level} => {offset, length}" do
    inv = File.read!("test/fixtures/narr/narr-a_221_20100615_1200_000.inv")
    file_size = 56_221_430

    {:ok, index} = NarrClient.parse_inventory(inv, file_size)

    # Spot-check a surface record
    assert {37_836_396, _len_1} = index[{"TMP", "2 m above gnd"}]
    # Spot-check a pressure-level record
    assert {_off, _len} = index[{"TMP", "850 mb"}]
    # File-end record's length is computed against file_size
    assert {_off, len_last} = index[{"PWAT", "atmos col"}]
    assert len_last > 0
  end
end
```

**Step 2: Run, see fail.**

**Step 3: Implement.**

```elixir
@spec parse_inventory(String.t(), pos_integer()) :: {:ok, %{{String.t(), String.t()} => {non_neg_integer(), pos_integer()}}}
def parse_inventory(inv_text, file_size) do
  lines =
    inv_text
    |> String.split("\n", trim: true)
    |> Enum.map(&parse_line/1)
    |> Enum.reject(&is_nil/1)

  index =
    lines
    |> Enum.zip(Enum.drop(lines, 1) ++ [{nil, file_size, nil, nil}])
    |> Map.new(fn {{_n, off, var, level}, {_, next_off, _, _}} ->
      {{var, level}, {off, next_off - off}}
    end)

  {:ok, index}
end

defp parse_line(line) do
  case String.split(line, ":") do
    [num, off, _date, var, level | _] ->
      {String.to_integer(num), String.to_integer(off), var, level}
    _ ->
      nil
  end
end
```

**Step 4: Run, green.**

**Step 5: Add `fetch_inventory/1` that wraps `parse_inventory` with an `Req.get` of the `.inv` URL** — this part needs a `Req.Test.stub` test to avoid hitting NCEI in the test suite.

**Step 6: Commit.**

### Task 4: NarrClient.fetch_profile_at/2

**Files:** modify `narr_client.ex` + test.

This is the big one. Given a `valid_time` and a `{lat, lon}` point, fetch all the records we need via byte-range, write each record to its own per-record temp file, run `cdo -merge` to combine them, run `cdo -outputtab,name,lev,value -remapnn,lon=X_lat=Y` to extract values at the QSO point, and return a `%{surface_temp_c: ..., profile: [...], ...}` map shaped exactly like what `Era5BatchClient.build_profile_attrs/5` produces.

**Step 1: Failing test using a real captured composite GRIB1 fixture.**

Capture a tiny composite GRIB1 (~1 MB total) containing all needed records for one analysis hour from the spike, commit as `test/fixtures/narr/narr_dfw_2010-06-15_12z.grb`. The test reads the fixture, calls `NarrClient.extract_profile_from_file/3` (a public-but-internal helper exposed for testability), and asserts the returned map has all expected keys with sane values.

```elixir
test "extract_profile_from_file/3 builds a profile from a captured fixture" do
  profile = NarrClient.extract_profile_from_file(
    "test/fixtures/narr/narr_dfw_2010-06-15_12z.grb",
    32.9,
    -97.0
  )

  # Surface — values come from the spike (DFW 2010-06-15 12Z)
  assert_in_delta profile.surface_temp_c, 23.8, 1.0
  assert_in_delta profile.surface_dewpoint_c, 21.3, 1.0
  assert profile.surface_pressure_mb > 950 and profile.surface_pressure_mb < 1050
  assert profile.hpbl_m >= 0 and profile.hpbl_m < 5000
  assert profile.pwat_mm > 30 and profile.pwat_mm < 60

  # Profile array
  assert length(profile.profile) >= 5
  Enum.each(profile.profile, fn level ->
    assert is_float(level["pres"])
    assert is_float(level["tmpc"])
    assert is_float(level["hght"])
  end)
end
```

**Step 2: Run, see fail.**

**Step 3: Implement.**

Two layers:
- `extract_profile_from_file/3` — testable, takes a composite GRIB1 path + lat + lon, runs `cdo -outputtab,name,lev,value -remapnn,lon=#{lon}_lat=#{lat} <path>`, parses the text output, maps `varNNN`+`lev` to semantic field names via the lookup table in the Architecture section, builds the era5_profiles-shaped map. Reuses `ThetaE.dewpoint_from_spfh/2` for pressure-level dewpoint derivation (NARR has DPT only at 2 m).
- `fetch_profile_at/2` — non-testable (does network), calls `fetch_inventory/1`, picks the records to fetch, byte-range-fetches each into its own `/tmp/narr_<tag>_<i>.grb` file, calls `cdo -merge` to combine them into `/tmp/narr_<tag>_merged.grb`, calls `extract_profile_from_file/3`, then deletes all temp files.

cdo invocation pattern (verified in spike):

```bash
cdo -merge /tmp/r1.grb /tmp/r2.grb /tmp/r3.grb /tmp/merged.grb
cdo -outputtab,name,lev,value -remapnn,lon=-97.0_lat=32.9 /tmp/merged.grb
```

Output format (one line per record):

```
#    name    lev    value 
   var11      2  296.939 
   var17      2 294.4932 
  var221      0  703.519 
    var1      0 99383.64 
   var54      0 45.10547 
   var11  85000 291.5186 
    var7  85000 1547.793 
   var51  85000 0.01040225
```

Parse: skip the header line (starts with `#`), split each data line on whitespace, ignore empty fields, take `[name, lev, value]`, look up `name` in the parameter table, build the row.

**Step 4: Run, green.**

**Step 5: Commit.**

---

## Phase 3: NarrFetchWorker (TDD)

### Task 5: Build NarrFetchWorker

**Files:**
- Create: `lib/microwaveprop/workers/narr_fetch_worker.ex`
- Create: `test/microwaveprop/workers/narr_fetch_worker_test.exs`

Per-contact Oban worker. Args: `%{contact_id: id}`. Looks up the contact, picks the nearest 3-hourly analysis hour to `qso_timestamp`, calls `NarrClient.fetch_profile_at(snapped_vt, {contact.lat, contact.lon})`, inserts an `era5_profiles` row via `Era5BatchClient.bulk_insert_profiles/1`, marks the contact as enriched.

**Step 1: Failing test** with `Req.Test.stub` for the network calls. Asserts the era5_profiles row is created and the contact's status is updated.

**Step 2: Implement.**

**Step 3: Commit.**

### Task 6: Add :narr queue to dev.exs and runtime.exs

**Files:**
- Modify: `config/dev.exs` (queue list)
- Modify: `config/runtime.exs` (queue list)

`narr: 6` in both. Network-bound, so high concurrency is fine.

### Task 7: Wire ContactWeatherEnqueueWorker to NarrFetchWorker

**Files:**
- Modify: `lib/microwaveprop/workers/contact_weather_enqueue_worker.ex` (`Era5FetchWorker.new(%{...})` call site, currently around line 196)
- Modify: matching test file

**Step 1: Update the failing test** to assert `NarrFetchWorker` is enqueued instead of `Era5FetchWorker`.

**Step 2: Implement the swap. The pre-2014 cutoff stays the same — anything with `qso_timestamp < 2014-10-02 00:00 UTC` gets NarrFetchWorker, after that gets HrrrFetchWorker.**

**Step 3: Run all tests.**

**Step 4: Commit.**

### Task 8: End-to-end smoke test in dev

**Files:** none — manual verification.

**Step 1:** Pick one pre-2014 contact in dev with `hrrr_status == :unavailable`, manually enqueue `NarrFetchWorker` for it, and verify a row lands in `era5_profiles` and the propagation map shows weather data for that historical timestamp.

---

## Phase 4: Decommission ERA5

### Task 9: Permanently document Era5* deprecation in runtime.exs

**Files:**
- Modify: `config/runtime.exs` (era5* paused-queue comment block)

Update the existing "TEMPORARILY PAUSED while we diagnose" comment to "DEPRECATED in favor of NARR via NarrFetchWorker — see plan 2026-04-15. Slated for deletion in a follow-up PR."

### Task 10 (follow-up PR, do not execute in this pass): Drop Era5* code paths

This is a separate PR after NARR has been running cleanly for a few days. Scope:
- Drop the `era5_cds_jobs` table (Ecto migration)
- Delete `Era5SubmitWorker`, `Era5PollWorker`, `Era5FetchWorker`, `Era5MonthBatchWorker`, `Era5BatchClient`, `Era5Client`, `Era5CdsJob` schema, and the `era5_*` queue config blocks
- Maybe rename `era5_profiles` → `historical_profiles` (or leave it — the misnomer is a footgun but harmless)

---

## Open Questions / Risks

1. ~~**`wgrib2 -lola` on the Lambert conformal grid.**~~ Resolved by spike: wgrib2 doesn't read GRIB1 at all, so we use cdo's `-remapnn,lon=X_lat=Y` instead. cdo handles the Lambert conformal projection natively.

2. **NARR records sometimes have different metadata strings across years.** The 1979 records may have slightly different inventory strings than 2014 records (kpds tables changed once or twice). The Phase 1 spike used 2010 — also test 1995 and 2013 in Phase 3 Task 8 to confirm the parser handles the full date range.

3. **NARR is 3-hourly, ERA5 was hourly.** That's a 3× temporal coarsening, but the propagation algorithm doesn't care — it picks the nearest sample to a QSO timestamp anyway. Acceptable.

4. **No tile-month batching.** ERA5 had it because CDS jobs were month-tile shaped. NARR fetches are per-point, much smaller, much faster — no batching needed. The flow is simpler: enqueue worker → worker calls client → client byte-range-fetches → done.

5. **NARR file 56 MB total per analysis hour, but byte-range fetches only pull ~1.5 MB max** for our 10 records. NCEI's HTTP server supports ranges (verified). If for some reason a range request fails, fall back to the full file download — but expect this to be rare.

6. **`narr-a` vs `narr-b`.** `narr-a` is the analysis (what we want). `narr-b` is a 3-hour forecast interpolation (filling in the off-hour times). Use `narr-a` exclusively — analyses are higher quality than forecast interpolations.