fix(wgrib2): account for 8-byte Fortran record overhead between messages

build_messages_per_message/3 was slicing the wgrib2 -lola bin output as
a packed float32 stream, ignoring the 4-byte record-length header + 4-byte
trailer that Fortran unformatted records prepend/append around each
message's data block. First value of message 0 came back as ~2.24e-44
(little-endian float32 of the int-16 record length), and every subsequent
message was shifted by 8 bytes so values spilled into neighbouring cells.

Fix: propagate the correct offset math already used by parse_lola_binary/3
(record_overhead = 8, stride = bytes_per_message + 8, data_offset =
msg_idx * stride + 4).

Flipped the characterization test from metadata-only to per-cell physical
asserts (200 K < TMP, DPT < 340 K; DPT <= TMP) and a 0.01 K cross-check
against extract_grid/3 output on the same fixture + bbox. Latent bug —
no production callers of extract_grid_messages*.
This commit is contained in:
Graham McIntire 2026-04-16 14:52:32 -05:00
parent 8c3c097009
commit 221d5516bc
No known key found for this signature in database
GPG key ID: F4ABF488E6029E59
2 changed files with 35 additions and 17 deletions

View file

@ -438,13 +438,20 @@ 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. Same math as
# parse_lola_binary/3 — see that function for context.
record_overhead = 8
stride = bytes_per_message + record_overhead
messages
|> Enum.with_index()
|> Enum.flat_map(fn {msg, msg_idx} ->
offset = msg_idx * bytes_per_message
data_offset = msg_idx * stride + 4
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)
values = extract_message_values(chunk, grid_spec)
[Map.put(msg, :values, values)]

View file

@ -185,20 +185,11 @@ defmodule Microwaveprop.Weather.Grib2.Wgrib2Test do
@describetag :slow
@describetag timeout: 60_000
# NOTE: This test characterizes CURRENT (BUGGY) behavior. See moduledoc
# above for the details. The metadata (var/level/datetime) is correct,
# but the per-cell values are corrupted because
# `build_messages_per_message` does not account for the 8-byte Fortran
# record overhead used by wgrib2's `-lola ... bin` output format.
# The sibling function `parse_lola_binary` does account for it
# (`record_overhead = 8` and `data_offset = msg_idx * stride + 4`),
# and the matching `extract_grid/3` assertions above confirm correct
# values for the same fixture + grid. Fixing this is out of scope for
# characterization — flagged in the report to the caller.
test "returns one entry per matched message with metadata (values currently corrupted by record-overhead bug)" do
test "returns one entry per matched message with real physical values" do
fixture = File.read!(@multi_fixture)
{:ok, messages} = Wgrib2.extract_grid_messages(fixture, ":(TMP|DPT):", @dallas_grid)
{:ok, grid} = Wgrib2.extract_grid(fixture, ":(TMP|DPT):", @dallas_grid)
# hrrr_multi.grib2 contains exactly TMP:2m and DPT:2m.
assert length(messages) == 2
@ -206,21 +197,42 @@ defmodule Microwaveprop.Weather.Grib2.Wgrib2Test do
vars = messages |> Enum.map(& &1.var) |> Enum.sort()
assert vars == ["DPT", "TMP"]
# Physical-value assertions per message — verifies the Fortran record
# overhead (8 bytes/message) is accounted for in the parser.
tmp_msg = Enum.find(messages, &(&1.var == "TMP"))
dpt_msg = Enum.find(messages, &(&1.var == "DPT"))
for msg <- messages do
assert msg.level == "2 m above ground"
assert %DateTime{} = msg.datetime
# parse_wgrib2_date always truncates to :second.
assert msg.datetime == DateTime.truncate(msg.datetime, :second)
assert msg.datetime.time_zone == "Etc/UTC"
assert is_map(msg.values)
assert map_size(msg.values) == 4
for {{lat, lon}, v} <- msg.values do
assert lat in [32.78, 32.88]
assert lon in [-96.8, -96.7]
assert is_float(v)
assert v > 200.0 and v < 340.0, "#{msg.var} value #{v} outside physical range"
end
end
# Dew point must be <= dry-bulb temp per point.
for {point, tmp} <- tmp_msg.values do
dpt = Map.fetch!(dpt_msg.values, point)
assert dpt <= tmp, "DPT #{dpt} > TMP #{tmp} at #{inspect(point)}"
end
# Cross-check: same inputs as extract_grid/3 — values must agree within 0.01 K.
for {point, tmp} <- tmp_msg.values do
assert_in_delta grid[point]["TMP:2 m above ground"], tmp, 0.01
end
for {point, dpt} <- dpt_msg.values do
assert_in_delta grid[point]["DPT:2 m above ground"], dpt, 0.01
end
end
end
@ -229,8 +241,7 @@ defmodule Microwaveprop.Weather.Grib2.Wgrib2Test do
@describetag timeout: 60_000
# Both the binary and file variants funnel into `build_messages_per_message`,
# so they produce the same (currently buggy — see note above) output.
# This test just asserts that the two code paths agree with each other.
# so they must produce identical output for identical inputs.
test "mirrors extract_grid_messages/3 on the same data" do
fixture = File.read!(@multi_fixture)
{:ok, from_bin} = Wgrib2.extract_grid_messages(fixture, ":(TMP|DPT):", @dallas_grid)