Add native HRRR duct detection to hourly propagation scoring

The PropagationGridWorker now fetches native hybrid-sigma levels
(TMP, SPFH, HGT, PRES × 50 levels) alongside the standard surface
and pressure products. Native data provides 10-50m vertical spacing
vs 250m from pressure levels, detecting thin surface ducts invisible
to the standard product.

Key design: cell-by-cell reducer in Wgrib2.extract_grid_from_file_mapped
processes each of the 95k CONUS cells through a duct analysis function
inline, keeping only scalar metrics per cell. Peak memory ~86 MB
instead of ~1.8 GB for the full grid map.

Per-cell output: native_min_gradient, best_duct_freq_ghz,
max_duct_thickness_m, duct_count. The scorer prefers the native
gradient over the pressure-level gradient when available.

Native fetch is optional — if it fails, scoring continues with
pressure-level data only.

lib/microwaveprop/propagation.ex

@@ -86,7 +86,7 @@ defmodule Microwaveprop.Propagation do
       pressure_mb: hrrr_profile.surface_pressure_mb,
       prev_pressure_mb: nil,
       rain_rate_mmhr: Scorer.precip_to_rate_mmhr(hrrr_profile[:precip_mm]),
-      min_refractivity_gradient: derived[:min_refractivity_gradient],
+      min_refractivity_gradient: hrrr_profile[:native_min_gradient] || derived[:min_refractivity_gradient],
       bl_depth_m: hrrr_profile[:hpbl_m],
       pwat_mm: hrrr_profile[:pwat_mm]
     }

lib/microwaveprop/weather/grib2/wgrib2.ex

@@ -40,6 +40,26 @@ defmodule Microwaveprop.Weather.Grib2.Wgrib2 do
     end
   end
 
+  @doc """
+  Like `extract_grid_from_file/3`, but instead of building the full
+  `%{{lat, lon} => %{"VAR:LEVEL" => float}}` map (which can be 1+ GB for
+  200+ messages × 95k cells), processes each grid cell through a reducer
+  function and keeps only the reduced output.
+
+  The reducer receives `%{"VAR:LEVEL" => float}` for one cell and returns
+  an arbitrary term. The result is `%{{lat, lon} => reducer_output}`.
+
+  Peak memory: binary output (~76 MB for 200 msgs × 95k cells) + one
+  cell's worth of data at a time + output map (~10 MB of scalars).
+  """
+  def extract_grid_from_file_mapped(grib_path, match_pattern, grid_spec, cell_reducer) do
+    if available?() do
+      extract_file_mapped_with_wgrib2(grib_path, match_pattern, grid_spec, cell_reducer)
+    else
+      {:error, :wgrib2_not_available}
+    end
+  end
+
   @doc "Check if wgrib2 is available on the system."
   def available?, do: wgrib2_path() != nil
 
@@ -82,6 +102,104 @@ defmodule Microwaveprop.Weather.Grib2.Wgrib2 do
     end
   end
 
+  defp extract_file_mapped_with_wgrib2(grib_path, match_pattern, grid_spec, cell_reducer) do
+    %{
+      lon_start: lon_start,
+      lon_count: lon_count,
+      lon_step: lon_step,
+      lat_start: lat_start,
+      lat_count: lat_count,
+      lat_step: lat_step
+    } = grid_spec
+
+    wgrib2_lon_start = normalize_lon(lon_start)
+    lon_spec = "#{wgrib2_lon_start}:#{lon_count}:#{lon_step}"
+    lat_spec = "#{lat_start}:#{lat_count}:#{lat_step}"
+
+    tmp_bin = grib_path <> ".lola.mapped.bin"
+
+    try do
+      args = [grib_path, "-match", match_pattern, "-lola", lon_spec, lat_spec, tmp_bin, "bin"]
+
+      case System.cmd(wgrib2_path(), args, stderr_to_stdout: true) do
+        {output, 0} ->
+          messages = parse_wgrib2_inventory(output)
+
+          case File.read(tmp_bin) do
+            {:ok, bin_data} ->
+              parse_lola_binary_mapped(bin_data, messages, grid_spec, cell_reducer)
+
+            {:error, :enoent} ->
+              {:ok, %{}}
+
+            {:error, reason} ->
+              {:error, "Failed to read wgrib2 output: #{inspect(reason)}"}
+          end
+
+        {output, exit_code} ->
+          {:error, "wgrib2 failed (exit #{exit_code}): #{String.slice(output, 0, 200)}"}
+      end
+    after
+      File.rm(tmp_bin)
+    end
+  end
+
+  # Cell-by-cell variant of parse_lola_binary. Instead of building the full
+  # %{{lat,lon} => %{"VAR:LEVEL" => val}} map, iterates grid cells and for
+  # each cell extracts its values from every message, passes the per-cell
+  # map to the reducer, and keeps only the reduced output.
+  defp parse_lola_binary_mapped(bin_data, messages, grid_spec, cell_reducer) do
+    %{lon_count: nx, lat_count: ny, lon_start: lon_start, lon_step: lon_step, lat_start: lat_start, lat_step: lat_step} =
+      grid_spec
+
+    points_per_message = nx * ny
+    bytes_per_message = points_per_message * 4
+    record_overhead = 8
+    stride = bytes_per_message + record_overhead
+
+    # Pre-compute message metadata: {key, data_offset} for each message
+    msg_meta =
+      messages
+      |> Enum.with_index()
+      |> Enum.map(fn {msg, idx} ->
+        {"#{msg.var}:#{msg.level}", idx * stride + 4}
+      end)
+      |> Enum.filter(fn {_key, offset} -> offset + bytes_per_message <= byte_size(bin_data) end)
+
+    # Iterate each grid cell, extract all message values, reduce
+    result =
+      Enum.reduce(0..(ny - 1), %{}, fn j, acc_outer ->
+        lat = Float.round(lat_start + j * lat_step, 3)
+
+        Enum.reduce(0..(nx - 1), acc_outer, fn i, acc_inner ->
+          cell_offset = (j * nx + i) * 4
+
+          # Extract all values for this cell across all messages
+          cell_data =
+            Enum.reduce(msg_meta, %{}, fn {key, data_offset}, cell_acc ->
+              offset = data_offset + cell_offset
+              <<_::binary-size(offset), value::float-little-32, _::binary>> = bin_data
+
+              if value > @undefined_value / 2 do
+                cell_acc
+              else
+                Map.put(cell_acc, key, value)
+              end
+            end)
+
+          if map_size(cell_data) > 0 do
+            lon = Float.round(denormalize_lon(lon_start + i * lon_step), 3)
+            reduced = cell_reducer.(cell_data)
+            Map.put(acc_inner, {lat, lon}, reduced)
+          else
+            acc_inner
+          end
+        end)
+      end)
+
+    {:ok, result}
+  end
+
   defp extract_with_wgrib2(grib_binary, match_pattern, grid_spec) do
     %{
       lon_start: lon_start,

lib/microwaveprop/weather/hrrr_native_client.ex

@@ -125,6 +125,116 @@ defmodule Microwaveprop.Weather.HrrrNativeClient do
     HrrrClient.byte_ranges_for_messages(idx_entries, native_messages())
   end
 
+  # The 4 variables needed for duct detection (N-profile + height).
+  # Skips UGRD, VGRD, TKE (shear/Richardson were dead features).
+  @duct_variables ~w(TMP SPFH HGT PRES)
+
+  @doc """
+  Messages needed for duct detection only (4 vars × 50 levels = 200 messages).
+  ~300 MB download vs ~530 MB for all 7 variables.
+  """
+  def duct_messages do
+    for level <- @native_levels, var <- @duct_variables do
+      %{var: var, level: "#{level} hybrid level"}
+    end
+  end
+
+  @doc """
+  Byte ranges for duct-detection variables only.
+  """
+  def duct_byte_ranges(idx_entries) do
+    HrrrClient.byte_ranges_for_messages(idx_entries, duct_messages())
+  end
+
+  @doc """
+  Fetch native duct metrics for the CONUS grid for one forecast hour.
+
+  Downloads TMP, SPFH, HGT, PRES on all 50 hybrid levels, extracts to
+  the CONUS grid via wgrib2, and computes duct metrics per cell using a
+  cell-by-cell reducer (peak memory ~86 MB instead of ~1.8 GB).
+
+  Returns `{:ok, %{{lat, lon} => duct_metrics}}` or `{:error, reason}`.
+  """
+  def fetch_native_duct_grid(date, hour, grid_spec, forecast_hour \\ 0) do
+    alias Microwaveprop.Weather.Grib2.Wgrib2
+
+    url = hrrr_native_url(date, hour, forecast_hour)
+    idx_url = url <> ".idx"
+
+    tmp_grib = Path.join(System.tmp_dir!(), "hrrr_native_duct_#{System.unique_integer([:positive])}.grib2")
+
+    try do
+      with {:ok, idx_text} <- fetch_idx(idx_url),
+           idx_entries = HrrrClient.parse_idx(idx_text),
+           ranges = duct_byte_ranges(idx_entries),
+           :ok <- HrrrClient.download_grib_ranges_to_file(url, ranges, tmp_grib) do
+        match_pattern = ":(#{Enum.join(@duct_variables, "|")}):.*hybrid level:"
+
+        Wgrib2.extract_grid_from_file_mapped(tmp_grib, match_pattern, grid_spec, &compute_duct_metrics/1)
+      end
+    after
+      File.rm(tmp_grib)
+    end
+  end
+
+  defp fetch_idx(url) do
+    case Req.get(url, receive_timeout: 120_000) do
+      {:ok, %{status: 200, body: body}} -> {:ok, body}
+      {:ok, %{status: status}} -> {:error, "HRRR native idx HTTP #{status}"}
+      {:error, reason} -> {:error, reason}
+    end
+  end
+
+  # Compute duct metrics from a single cell's native profile data.
+  # Input: %{"TMP:1 hybrid level" => 297.5, "SPFH:1 hybrid level" => 0.012, ...}
+  # Output: %{native_min_gradient: float, best_duct_freq_ghz: float, max_duct_thickness_m: float}
+  defp compute_duct_metrics(parsed) do
+    alias Microwaveprop.Propagation.Duct
+
+    profile = build_native_profile(parsed)
+
+    if profile.level_count >= 3 do
+      duct_result = Duct.analyze(profile)
+
+      %{
+        native_min_gradient: min_m_gradient(profile),
+        best_duct_freq_ghz: duct_result.best_duct_band_ghz,
+        max_duct_thickness_m: max_duct_thickness(duct_result.ducts),
+        duct_count: length(duct_result.ducts)
+      }
+    else
+      %{native_min_gradient: nil, best_duct_freq_ghz: nil, max_duct_thickness_m: nil, duct_count: 0}
+    end
+  end
+
+  defp min_m_gradient(%{heights_m: heights, temp_k: temps, spfh: spfhs, pressure_pa: pressures}) do
+    # Compute modified refractivity M at each level, find minimum dM/dh
+    ms =
+      Enum.zip([heights, temps, spfhs, pressures])
+      |> Enum.map(fn {h, t, q, p} ->
+        # N = 77.6 * P/T + 3.73e5 * e/T^2, where e = q*P/(0.622 + 0.378*q)
+        q = max(q || 0.0, 1.0e-8)
+        e = q * p / (0.622 + 0.378 * q) / 100.0
+        t_safe = max(t, 1.0)
+        n = 77.6 * (p / 100.0) / t_safe + 3.73e5 * e / (t_safe * t_safe)
+        m = n + 157.0 * h / 1000.0
+        {h, m}
+      end)
+
+    ms
+    |> Enum.chunk_every(2, 1, :discard)
+    |> Enum.map(fn [{h1, m1}, {h2, m2}] ->
+      dh = h2 - h1
+      if dh > 0, do: (m2 - m1) / dh * 1000.0, else: 0.0
+    end)
+    |> Enum.min(fn -> 0.0 end)
+  end
+
+  defp min_m_gradient(_), do: nil
+
+  defp max_duct_thickness([]), do: nil
+  defp max_duct_thickness(ducts), do: ducts |> Enum.map(& &1.thickness_m) |> Enum.max()
+
   @doc """
   Extract native profiles for a list of `{lat, lon}` points from a
   GRIB2 file on disk, using wgrib2. Avoids loading the entire file

lib/microwaveprop/workers/propagation_grid_worker.ex

@@ -14,6 +14,7 @@ defmodule Microwaveprop.Workers.PropagationGridWorker do
   alias Microwaveprop.Propagation.Grid
   alias Microwaveprop.Weather
   alias Microwaveprop.Weather.HrrrClient
+  alias Microwaveprop.Weather.HrrrNativeClient
   alias Microwaveprop.Weather.SoundingParams
 
   require Logger
@@ -78,6 +79,10 @@ defmodule Microwaveprop.Workers.PropagationGridWorker do
         store_hrrr_profiles(grid_data, valid_time)
         :erlang.garbage_collect()
 
+        # Fetch native duct metrics and merge into grid_data for scoring
+        grid_data = merge_native_duct_data(grid_data, run_time, forecast_hour)
+        :erlang.garbage_collect()
+
         Phoenix.PubSub.broadcast(
           Microwaveprop.PubSub,
           "weather:updated",
@@ -163,6 +168,31 @@ defmodule Microwaveprop.Workers.PropagationGridWorker do
     Logger.info("PropagationGrid: stored #{count} HRRR profiles")
   end
 
+  defp merge_native_duct_data(grid_data, run_time, forecast_hour) do
+    hour_dt = HrrrClient.nearest_hrrr_hour(run_time)
+    date = DateTime.to_date(hour_dt)
+    hour = hour_dt.hour
+    grid_spec = Grid.wgrib2_grid_spec()
+
+    case timed("native", fn ->
+           HrrrNativeClient.fetch_native_duct_grid(date, hour, grid_spec, forecast_hour)
+         end) do
+      {:ok, duct_grid} ->
+        Logger.info("PropagationGrid: merged #{map_size(duct_grid)} native duct cells")
+
+        Map.new(grid_data, fn {point, profile} ->
+          case Map.get(duct_grid, point) do
+            nil -> {point, profile}
+            duct -> {point, Map.merge(profile, duct)}
+          end
+        end)
+
+      {:error, reason} ->
+        Logger.warning("PropagationGrid: native duct fetch failed (continuing without): #{inspect(reason)}")
+        grid_data
+    end
+  end
+
   defp compute_scores(grid_data, valid_time) do
     # Algorithm is the primary scorer. ML score stored in factors for comparison.
     compute_scores_algorithm(grid_data, valid_time)