From e7f0f03bf0f8b2c22ed4f9b6f358a0afd27109ea Mon Sep 17 00:00:00 2001 From: Graham McIntire Date: Fri, 10 Apr 2026 08:13:32 -0500 Subject: [PATCH] Fix wgrib2 -lola binary parsing and add native_surface_refractivity MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 | 60 +++++++++++++++++++ lib/microwaveprop/weather/grib2/wgrib2.ex | 14 ++++- .../weather/hrrr_native_client.ex | 4 +- test/microwaveprop/backtest/features_test.exs | 34 +++++++++++ 4 files changed, 108 insertions(+), 4 deletions(-) diff --git a/lib/microwaveprop/backtest/features.ex b/lib/microwaveprop/backtest/features.ex index 2cc0d226..00c26c7a 100644 --- a/lib/microwaveprop/backtest/features.ex +++ b/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 diff --git a/lib/microwaveprop/weather/grib2/wgrib2.ex b/lib/microwaveprop/weather/grib2/wgrib2.ex index efa2aadd..cc56af07 100644 --- a/lib/microwaveprop/weather/grib2/wgrib2.ex +++ b/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 diff --git a/lib/microwaveprop/weather/hrrr_native_client.ex b/lib/microwaveprop/weather/hrrr_native_client.ex index b5fc31a4..faacb0c0 100644 --- a/lib/microwaveprop/weather/hrrr_native_client.ex +++ b/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) diff --git a/test/microwaveprop/backtest/features_test.exs b/test/microwaveprop/backtest/features_test.exs index a67d750f..b7430273 100644 --- a/test/microwaveprop/backtest/features_test.exs +++ b/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