Graham McIntire
5a9359191f
Eliminates the external wgrib2 C tool dependency that blocked HRRR processing. Implements Lambert Conformal projection, simple packing (Template 5.0), complex packing with spatial differencing (Template 5.3), and GRIB2 section parsing — enough to extract point values from HRRR grid data using only Elixir.
243 lines
6.4 KiB
Elixir
243 lines
6.4 KiB
Elixir
defmodule Microwaveprop.Weather.Grib2.ExtractorTest do
|
|
use ExUnit.Case, async: true
|
|
|
|
alias Microwaveprop.Weather.Grib2.Extractor
|
|
|
|
# HRRR grid parameters for building synthetic messages
|
|
@hrrr_nx 1799
|
|
@hrrr_ny 1059
|
|
|
|
describe "split_messages/1" do
|
|
test "splits concatenated GRIB2 messages" do
|
|
msg1 = build_synthetic_grib2_message(0, 0, 0, 103, 2)
|
|
msg2 = build_synthetic_grib2_message(0, 0, 6, 103, 2)
|
|
blob = msg1 <> msg2
|
|
|
|
messages = Extractor.split_messages(blob)
|
|
assert length(messages) == 2
|
|
assert Enum.at(messages, 0) == msg1
|
|
assert Enum.at(messages, 1) == msg2
|
|
end
|
|
|
|
test "handles single message" do
|
|
msg = build_synthetic_grib2_message(0, 0, 0, 103, 2)
|
|
messages = Extractor.split_messages(msg)
|
|
assert length(messages) == 1
|
|
end
|
|
|
|
test "handles empty binary" do
|
|
assert Extractor.split_messages(<<>>) == []
|
|
end
|
|
end
|
|
|
|
describe "extract_points/3" do
|
|
test "extracts single message at first grid point" do
|
|
msg = build_synthetic_grib2_message(0, 0, 0, 103, 2)
|
|
|
|
# la1/lo1 = first grid point, so index (0,0) → linear index 0
|
|
assert {:ok, result} = Extractor.extract_points(msg, 21.138123, 237.280472)
|
|
|
|
assert Map.has_key?(result, "TMP:2 m above ground")
|
|
assert is_float(result["TMP:2 m above ground"])
|
|
end
|
|
|
|
test "extracts multiple concatenated messages" do
|
|
msg1 = build_synthetic_grib2_message(0, 0, 0, 103, 2)
|
|
msg2 = build_synthetic_grib2_message(0, 0, 6, 103, 2)
|
|
msg3 = build_synthetic_grib2_message(0, 3, 0, 1, 0)
|
|
blob = msg1 <> msg2 <> msg3
|
|
|
|
assert {:ok, result} = Extractor.extract_points(blob, 21.138123, 237.280472)
|
|
|
|
assert Map.has_key?(result, "TMP:2 m above ground")
|
|
assert Map.has_key?(result, "DPT:2 m above ground")
|
|
assert Map.has_key?(result, "PRES:surface")
|
|
assert map_size(result) == 3
|
|
end
|
|
|
|
test "handles negative longitude" do
|
|
msg = build_synthetic_grib2_message(0, 0, 0, 103, 2)
|
|
|
|
# -122.719528 == 237.280472
|
|
assert {:ok, result} = Extractor.extract_points(msg, 21.138123, -122.719528)
|
|
assert Map.has_key?(result, "TMP:2 m above ground")
|
|
end
|
|
|
|
test "returns error for point outside grid" do
|
|
msg = build_synthetic_grib2_message(0, 0, 0, 103, 2)
|
|
|
|
# London — way outside HRRR domain
|
|
assert {:error, _} = Extractor.extract_points(msg, 51.5, -0.12)
|
|
end
|
|
end
|
|
|
|
describe "linear_index/3" do
|
|
test "scan mode 64 (j positive, i positive, i consecutive)" do
|
|
nx = 100
|
|
scan_mode = 64
|
|
assert Extractor.linear_index({5, 10}, nx, scan_mode) == 10 * 100 + 5
|
|
end
|
|
|
|
test "scan mode 0 (j negative, i positive, i consecutive)" do
|
|
# j decreasing means row 0 is top, so index = (ny - 1 - j) * nx + i
|
|
# But we just return j * nx + i for now since we don't flip
|
|
nx = 100
|
|
scan_mode = 0
|
|
assert Extractor.linear_index({5, 10}, nx, scan_mode) == 10 * 100 + 5
|
|
end
|
|
end
|
|
|
|
describe "extract_points/3 with real HRRR fixture" do
|
|
@fixture_path "test/fixtures/grib2/hrrr_tmp_2m.grib2"
|
|
@multi_fixture_path "test/fixtures/grib2/hrrr_multi.grib2"
|
|
|
|
@tag :fixture
|
|
@tag timeout: 120_000
|
|
test "extracts TMP:2m from real single-message GRIB2" do
|
|
if File.exists?(@fixture_path) do
|
|
fixture = File.read!(@fixture_path)
|
|
|
|
assert {:ok, result} = Extractor.extract_points(fixture, 32.78, -96.8)
|
|
assert Map.has_key?(result, "TMP:2 m above ground")
|
|
|
|
temp_k = result["TMP:2 m above ground"]
|
|
temp_c = temp_k - 273.15
|
|
assert temp_c > -10.0 and temp_c < 45.0, "Dallas temp #{temp_c}C unreasonable"
|
|
end
|
|
end
|
|
|
|
@tag :fixture
|
|
@tag timeout: 300_000
|
|
test "extracts from multi-message GRIB2" do
|
|
if File.exists?(@multi_fixture_path) do
|
|
fixture = File.read!(@multi_fixture_path)
|
|
|
|
assert {:ok, result} = Extractor.extract_points(fixture, 32.78, -96.8)
|
|
assert Map.has_key?(result, "TMP:2 m above ground")
|
|
assert Map.has_key?(result, "DPT:2 m above ground")
|
|
assert map_size(result) == 2
|
|
end
|
|
end
|
|
end
|
|
|
|
# --- Synthetic GRIB2 builder ---
|
|
|
|
defp build_synthetic_grib2_message(discipline, param_category, param_number, surface_type, surface_value) do
|
|
num_points = @hrrr_nx * @hrrr_ny
|
|
|
|
sec1 = section_wrap(1, build_section1())
|
|
|
|
sec3 =
|
|
section_wrap(3, build_section3(@hrrr_nx, @hrrr_ny))
|
|
|
|
sec4 =
|
|
section_wrap(
|
|
4,
|
|
build_section4(param_category, param_number, surface_type, surface_value)
|
|
)
|
|
|
|
sec5 = section_wrap(5, build_section5(num_points, 16))
|
|
|
|
sec6 = section_wrap(6, <<255::8>>)
|
|
|
|
# Data: num_points * 2 bytes (16-bit values), all set to 42
|
|
data_bytes = :binary.copy(<<42::16-big>>, num_points)
|
|
sec7 = section_wrap(7, data_bytes)
|
|
|
|
sec8 = "7777"
|
|
|
|
body = sec1 <> sec3 <> sec4 <> sec5 <> sec6 <> sec7 <> sec8
|
|
total_length = 16 + byte_size(body)
|
|
|
|
sec0 = <<"GRIB", 0::16, discipline::8, 2::8, total_length::64-big>>
|
|
sec0 <> body
|
|
end
|
|
|
|
defp build_section1 do
|
|
<<
|
|
0::16-big,
|
|
0::16-big,
|
|
2::8,
|
|
1::8,
|
|
1::8,
|
|
2026::16-big,
|
|
3::8,
|
|
28::8,
|
|
18::8,
|
|
0::8,
|
|
0::8,
|
|
0::8,
|
|
1::8
|
|
>>
|
|
end
|
|
|
|
defp build_section3(nx, ny) do
|
|
<<
|
|
0::8,
|
|
nx * ny::32-big,
|
|
0::8,
|
|
0::8,
|
|
30::16-big,
|
|
6::8,
|
|
0::8,
|
|
0::32-big,
|
|
0::8,
|
|
0::32-big,
|
|
0::8,
|
|
0::32-big,
|
|
nx::32-big,
|
|
ny::32-big,
|
|
21_138_123::32-big,
|
|
237_280_472::32-big,
|
|
0b00111000::8,
|
|
38_500_000::32-big,
|
|
262_500_000::32-big,
|
|
3_000_000::32-big,
|
|
3_000_000::32-big,
|
|
0::8,
|
|
64::8,
|
|
38_500_000::32-big,
|
|
38_500_000::32-big,
|
|
0::32-big,
|
|
0::32-big
|
|
>>
|
|
end
|
|
|
|
defp build_section4(param_category, param_number, surface_type, surface_value) do
|
|
<<
|
|
0::16-big,
|
|
0::16-big,
|
|
param_category::8,
|
|
param_number::8,
|
|
2::8,
|
|
0::8,
|
|
0::8,
|
|
0::16-big,
|
|
0::8,
|
|
1::8,
|
|
0::32-big,
|
|
surface_type::8,
|
|
0::8,
|
|
surface_value::32-big,
|
|
255::8,
|
|
0::8,
|
|
0::32-big
|
|
>>
|
|
end
|
|
|
|
defp build_section5(num_points, bits_per_value) do
|
|
<<
|
|
num_points::32-big,
|
|
0::16-big,
|
|
0.0::float-32-big,
|
|
0::16-big,
|
|
0::16-big,
|
|
bits_per_value::8
|
|
>>
|
|
end
|
|
|
|
defp section_wrap(section_number, body) do
|
|
length = 5 + byte_size(body)
|
|
<<length::32-big, section_number::8, body::binary>>
|
|
end
|
|
end
|