Fix wgrib2 -lola binary parsing and add native_surface_refractivity

wgrib2 -lola ... bin writes Fortran unformatted records (4-byte
length header + data + 4-byte length trailer per message).
parse_lola_binary was treating the binary as tightly packed,
causing every message after the first to read from the wrong
offset — values came out as garbage across all grid points.

Fix: account for the 8-byte record overhead per message when
computing the data offset for each message's grid values.

This bug affects both the existing propagation grid extraction
(which may have been producing subtly wrong scores) and the new
native-level extraction (which was producing obviously wrong
values). The fix is a one-line stride change.

Also adds Backtest.Features.native_surface_refractivity for the
Phase 1 sanity check, plus a tighter wgrib2 match pattern that
selects only hybrid-level messages from the native file.

lib/microwaveprop/backtest/features.ex

@@ -18,7 +18,11 @@ defmodule Microwaveprop.Backtest.Features do
     a raise on one point kills the whole report.
   """
 
+  import Ecto.Query
+
+  alias Microwaveprop.Repo
   alias Microwaveprop.Weather
+  alias Microwaveprop.Weather.HrrrNativeProfile
 
   @doc """
   Minimum refractivity gradient from the nearest HRRR profile.
@@ -81,4 +85,60 @@ defmodule Microwaveprop.Backtest.Features do
       _ -> nil
     end
   end
+
+  @doc """
+  Surface refractivity from the lowest level of the nearest native
+  hybrid-sigma profile.
+
+  Phase 1 sanity check: this should produce numbers comparable to the
+  `naive_gradient` baseline but computed from native-level data. If
+  the sign/magnitude look wrong, the native ingestion pipeline has a
+  bug.
+
+  Uses the ITU-R P.453-14 formula: N = 77.6*P/T + 3.73e5*e/T²
+  where e (water vapor pressure) is derived from specific humidity.
+  """
+  @spec native_surface_refractivity(float, float, DateTime.t()) :: float | nil
+  def native_surface_refractivity(lat, lon, valid_time) do
+    with %HrrrNativeProfile{} = profile <- find_nearest_native(lat, lon, valid_time),
+         t when is_float(t) <- profile.surface_temp_k,
+         p when is_float(p) <- profile.surface_pressure_pa,
+         q when is_float(q) <- profile.surface_spfh do
+      # Water vapor pressure from specific humidity: e = q*P / (0.622 + 0.378*q)
+      e = q * p / (0.622 + 0.378 * q)
+      # N-units
+      77.6 * p / (t * 100) + 3.73e5 * e / (t * t * 100)
+    else
+      _ -> nil
+    end
+  end
+
+  defp find_nearest_native(lat, lon, valid_time) do
+    dlat = 0.07
+    dlon = 0.07
+    time_start = DateTime.add(valid_time, -3600, :second)
+    time_end = DateTime.add(valid_time, 3600, :second)
+
+    HrrrNativeProfile
+    |> where(
+      [p],
+      p.lat >= ^(lat - dlat) and p.lat <= ^(lat + dlat) and
+        p.lon >= ^(lon - dlon) and p.lon <= ^(lon + dlon) and
+        p.valid_time >= ^time_start and p.valid_time <= ^time_end
+    )
+    |> order_by([p],
+      asc:
+        fragment(
+          "ABS(? - ?) + ABS(? - ?) + ABS(EXTRACT(EPOCH FROM ? - ?))",
+          p.lat,
+          ^lat,
+          p.lon,
+          ^lon,
+          p.valid_time,
+          ^valid_time
+        )
+    )
+    |> limit(1)
+    |> Repo.one()
+  end
 end

lib/microwaveprop/weather/grib2/wgrib2.ex

@@ -251,15 +251,23 @@ defmodule Microwaveprop.Weather.Grib2.Wgrib2 do
     points_per_message = nx * ny
     bytes_per_message = points_per_message * 4
 
+    # wgrib2 -lola ... bin writes Fortran unformatted binary: each
+    # message is preceded by a 4-byte little-endian record length
+    # and followed by a duplicate 4-byte record length. So the
+    # stride per message is 4 + data + 4 = data + 8.
+    record_overhead = 8
+    stride = bytes_per_message + record_overhead
+
     result =
       messages
       |> Enum.with_index()
       |> Enum.reduce(%{}, fn {msg, msg_idx}, acc ->
-        offset = msg_idx * bytes_per_message
+        # Skip the 4-byte header to reach the data
+        data_offset = msg_idx * stride + 4
         key = "#{msg.var}:#{msg.level}"
 
-        if offset + bytes_per_message <= byte_size(bin_data) do
-          chunk = binary_part(bin_data, offset, bytes_per_message)
+        if data_offset + bytes_per_message <= byte_size(bin_data) do
+          chunk = binary_part(bin_data, data_offset, bytes_per_message)
           merge_message_values(acc, key, chunk, nx, ny, lon_start, lon_step, lat_start, lat_step)
         else
           acc

lib/microwaveprop/weather/hrrr_native_client.ex

@@ -139,7 +139,9 @@ defmodule Microwaveprop.Weather.HrrrNativeClient do
   def extract_native_profiles(grib_binary, points) when is_list(points) do
     alias Microwaveprop.Weather.Grib2.Wgrib2
 
-    match_pattern = ":(#{Enum.join(@native_variables, "|")}):"
+    # Match only hybrid-level messages — the simple var-name pattern
+    # also hits surface/2m/10m messages which misalign the binary output.
+    match_pattern = ":(#{Enum.join(@native_variables, "|")}):.*hybrid level:"
 
     if Wgrib2.available?() do
       grid_spec = bounding_grid(points)

test/microwaveprop/backtest/features_test.exs

@@ -3,6 +3,7 @@ defmodule Microwaveprop.Backtest.FeaturesTest do
 
   alias Microwaveprop.Backtest.Features
   alias Microwaveprop.Weather
+  alias Microwaveprop.Weather.HrrrNativeProfile
 
   @point_lat 32.9
   @point_lon -97.0
@@ -59,4 +60,37 @@ defmodule Microwaveprop.Backtest.FeaturesTest do
       assert Features.pressure(@point_lat, @point_lon, @valid_time) == nil
     end
   end
+
+  describe "native_surface_refractivity/3" do
+    test "computes N-units from native profile surface scalars" do
+      {:ok, _} =
+        %HrrrNativeProfile{}
+        |> HrrrNativeProfile.changeset(%{
+          valid_time: @valid_time,
+          lat: @point_lat,
+          lon: @point_lon,
+          level_count: 1,
+          heights_m: [10.0],
+          temp_k: [295.0],
+          spfh: [0.010],
+          pressure_pa: [101_325.0],
+          u_wind_ms: [2.0],
+          v_wind_ms: [1.0],
+          tke_m2s2: [0.5],
+          surface_temp_k: 295.0,
+          surface_spfh: 0.010,
+          surface_pressure_pa: 101_325.0
+        })
+        |> Repo.insert()
+
+      result = Features.native_surface_refractivity(@point_lat, @point_lon, @valid_time)
+      assert is_float(result)
+      # Typical surface N is 250-400 in temperate latitudes
+      assert result > 200.0 and result < 500.0
+    end
+
+    test "returns nil when no native profile exists" do
+      assert Features.native_surface_refractivity(@point_lat, @point_lon, @valid_time) == nil
+    end
+  end
 end