Merge branch 'worktree-agent-abacc39e' into feature/propagation-map
This commit is contained in:
commit
5c344ccedb
6 changed files with 402 additions and 0 deletions
|
|
@ -19,6 +19,25 @@ defmodule Microwaveprop.Weather.Grib2.ComplexPacking do
|
|||
end
|
||||
end
|
||||
|
||||
@doc """
|
||||
Extract multiple values from GRIB2 complex-packed data at the given grid indices.
|
||||
|
||||
Decodes all values (since spatial differencing requires sequential access),
|
||||
then returns the values at the requested indices.
|
||||
|
||||
Returns `{:ok, %{index => value}}`.
|
||||
"""
|
||||
def extract_values(params, data, indices) do
|
||||
case decode_all(params, data) do
|
||||
{:ok, array} ->
|
||||
results = Map.new(indices, fn index -> {index, :array.get(index, array)} end)
|
||||
{:ok, results}
|
||||
|
||||
{:error, _} = err ->
|
||||
err
|
||||
end
|
||||
end
|
||||
|
||||
@doc """
|
||||
Decode all values from complex-packed data with spatial differencing.
|
||||
|
||||
|
|
|
|||
|
|
@ -41,6 +41,28 @@ defmodule Microwaveprop.Weather.Grib2.Extractor do
|
|||
j * nx + i
|
||||
end
|
||||
|
||||
@doc """
|
||||
Extract weather values from a GRIB2 binary blob at multiple lat/lon points.
|
||||
|
||||
For each message in the binary, converts all points to grid indices, then
|
||||
batch-extracts values. Points outside the grid are silently skipped.
|
||||
|
||||
Returns `{:ok, %{{lat, lon} => %{"VAR:LEVEL" => float}}}` or `{:error, term}`.
|
||||
"""
|
||||
def extract_grid(binary, points) when is_list(points) do
|
||||
messages = split_messages(binary)
|
||||
|
||||
Enum.reduce_while(messages, {:ok, %{}}, fn msg, {:ok, acc} ->
|
||||
case extract_grid_single(msg, points) do
|
||||
{:ok, key, point_values} ->
|
||||
{:cont, {:ok, merge_point_values(acc, key, point_values)}}
|
||||
|
||||
{:error, reason} ->
|
||||
{:halt, {:error, reason}}
|
||||
end
|
||||
end)
|
||||
end
|
||||
|
||||
# --- Private ---
|
||||
|
||||
defp split_messages(<<>>, acc), do: Enum.reverse(acc)
|
||||
|
|
@ -73,6 +95,54 @@ defmodule Microwaveprop.Weather.Grib2.Extractor do
|
|||
e -> {:error, "GRIB2 extraction failed: #{inspect(e)}"}
|
||||
end
|
||||
|
||||
defp extract_grid_single(msg, points) do
|
||||
with {:ok, parsed} <- Section.parse_message(msg) do
|
||||
%{grid_params: grid, product: prod, packing_params: packing, data: data} = parsed
|
||||
key = "#{prod.var}:#{prod.level}"
|
||||
point_indices = resolve_grid_indices(grid, points)
|
||||
batch_extract_points(packing, data, key, point_indices)
|
||||
end
|
||||
rescue
|
||||
e -> {:error, "GRIB2 grid extraction failed: #{inspect(e)}"}
|
||||
end
|
||||
|
||||
defp resolve_grid_indices(grid, points) do
|
||||
points
|
||||
|> Enum.reduce([], fn {lat, lon} = point, acc ->
|
||||
case LambertConformal.to_grid_index(grid, lat, lon) do
|
||||
{:ok, {i, j}} ->
|
||||
index = linear_index({i, j}, grid.nx, grid.scan_mode)
|
||||
[{point, index} | acc]
|
||||
|
||||
{:error, :outside_grid} ->
|
||||
acc
|
||||
end
|
||||
end)
|
||||
|> Enum.reverse()
|
||||
end
|
||||
|
||||
defp batch_extract_points(_packing, _data, key, []), do: {:ok, key, []}
|
||||
|
||||
defp batch_extract_points(packing, data, key, point_indices) do
|
||||
unique_indices = point_indices |> Enum.map(&elem(&1, 1)) |> Enum.uniq()
|
||||
|
||||
case unpack_values(packing, data, unique_indices) do
|
||||
{:ok, values_map} ->
|
||||
point_values = Enum.map(point_indices, fn {point, index} -> {point, Map.fetch!(values_map, index)} end)
|
||||
{:ok, key, point_values}
|
||||
|
||||
{:error, _} = err ->
|
||||
err
|
||||
end
|
||||
end
|
||||
|
||||
defp merge_point_values(acc, key, point_values) do
|
||||
Enum.reduce(point_values, acc, fn {point, value}, a ->
|
||||
existing = Map.get(a, point, %{})
|
||||
Map.put(a, point, Map.put(existing, key, value))
|
||||
end)
|
||||
end
|
||||
|
||||
defp unpack_value(%{template: 3} = params, data, index) do
|
||||
ComplexPacking.extract_value(params, data, index)
|
||||
end
|
||||
|
|
@ -80,4 +150,12 @@ defmodule Microwaveprop.Weather.Grib2.Extractor do
|
|||
defp unpack_value(params, data, index) do
|
||||
SimplePacking.extract_value(params, data, index)
|
||||
end
|
||||
|
||||
defp unpack_values(%{template: 3} = params, data, indices) do
|
||||
ComplexPacking.extract_values(params, data, indices)
|
||||
end
|
||||
|
||||
defp unpack_values(params, data, indices) do
|
||||
SimplePacking.extract_values(params, data, indices)
|
||||
end
|
||||
end
|
||||
|
|
|
|||
|
|
@ -31,6 +31,42 @@ defmodule Microwaveprop.Weather.Grib2.SimplePacking do
|
|||
end
|
||||
end
|
||||
|
||||
@doc """
|
||||
Extract multiple values from GRIB2 simple-packed data at the given grid indices.
|
||||
|
||||
Returns `{:ok, %{index => value}}`. Out-of-range indices are silently skipped.
|
||||
"""
|
||||
def extract_values(%{bits_per_value: 0} = params, _data, indices) do
|
||||
value = params.reference_value * :math.pow(10, -params.decimal_scale)
|
||||
{:ok, Map.new(indices, fn i -> {i, value} end)}
|
||||
end
|
||||
|
||||
def extract_values(params, data, indices) do
|
||||
%{
|
||||
reference_value: r,
|
||||
binary_scale: e,
|
||||
decimal_scale: d,
|
||||
bits_per_value: n,
|
||||
num_data_points: num
|
||||
} = params
|
||||
|
||||
factor_2e = :math.pow(2, e)
|
||||
factor_10d = :math.pow(10, -d)
|
||||
|
||||
results =
|
||||
Enum.reduce(indices, %{}, fn index, acc ->
|
||||
if index >= 0 and index < num do
|
||||
x = extract_bits(data, index, n)
|
||||
value = (r + x * factor_2e) * factor_10d
|
||||
Map.put(acc, index, value)
|
||||
else
|
||||
acc
|
||||
end
|
||||
end)
|
||||
|
||||
{:ok, results}
|
||||
end
|
||||
|
||||
defp extract_bits(data, index, bits_per_value) do
|
||||
bit_offset = index * bits_per_value
|
||||
byte_offset = div(bit_offset, 8)
|
||||
|
|
|
|||
|
|
@ -92,4 +92,52 @@ defmodule Microwaveprop.Weather.Grib2.ComplexPackingTest do
|
|||
end
|
||||
end
|
||||
end
|
||||
|
||||
describe "extract_values/3" do
|
||||
@tag :fixture
|
||||
test "batch extracts multiple values matching individual calls" do
|
||||
if File.exists?(@fixture_path) do
|
||||
fixture = File.read!(@fixture_path)
|
||||
{:ok, parsed} = Section.parse_message(fixture)
|
||||
|
||||
indices = [0, 1000, 100_000, 500_000, 1_000_000]
|
||||
|
||||
{:ok, batch_result} =
|
||||
ComplexPacking.extract_values(parsed.packing_params, parsed.data, indices)
|
||||
|
||||
assert map_size(batch_result) == length(indices)
|
||||
|
||||
# Each batch value should match individual extraction
|
||||
for idx <- indices do
|
||||
{:ok, individual} =
|
||||
ComplexPacking.extract_value(parsed.packing_params, parsed.data, idx)
|
||||
|
||||
assert_in_delta batch_result[idx],
|
||||
individual,
|
||||
0.001,
|
||||
"Index #{idx}: batch #{batch_result[idx]} != individual #{individual}"
|
||||
end
|
||||
else
|
||||
:ok
|
||||
end
|
||||
end
|
||||
|
||||
@tag :fixture
|
||||
test "returns values in meteorological range" do
|
||||
if File.exists?(@fixture_path) do
|
||||
fixture = File.read!(@fixture_path)
|
||||
{:ok, parsed} = Section.parse_message(fixture)
|
||||
|
||||
indices = [0, 500_000, 1_905_140]
|
||||
{:ok, result} = ComplexPacking.extract_values(parsed.packing_params, parsed.data, indices)
|
||||
|
||||
for idx <- indices do
|
||||
assert result[idx] > 200.0 and result[idx] < 340.0,
|
||||
"Index #{idx}: expected 200-340K, got #{result[idx]}"
|
||||
end
|
||||
else
|
||||
:ok
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
|
|
|||
|
|
@ -87,6 +87,121 @@ defmodule Microwaveprop.Weather.Grib2.ExtractorTest do
|
|||
end
|
||||
end
|
||||
|
||||
describe "extract_grid/2" do
|
||||
test "extracts values for multiple points from a single message" do
|
||||
msg = build_synthetic_grib2_message(0, 0, 0, 103, 2)
|
||||
|
||||
# la1/lo1 is the first grid point at (21.138123, 237.280472)
|
||||
# Use the same point twice to verify both get values
|
||||
points = [{21.138123, -122.719528}, {21.138123, -122.719528}]
|
||||
|
||||
assert {:ok, result} = Extractor.extract_grid(msg, points)
|
||||
|
||||
# Both points are the same, so we should get 1 unique key
|
||||
assert map_size(result) == 1
|
||||
assert Map.has_key?(result, {21.138123, -122.719528})
|
||||
assert Map.has_key?(result[{21.138123, -122.719528}], "TMP:2 m above ground")
|
||||
end
|
||||
|
||||
test "extracts from concatenated messages for multiple points" do
|
||||
msg1 = build_synthetic_grib2_message(0, 0, 0, 103, 2)
|
||||
msg2 = build_synthetic_grib2_message(0, 0, 6, 103, 2)
|
||||
blob = msg1 <> msg2
|
||||
|
||||
points = [{21.138123, -122.719528}]
|
||||
|
||||
assert {:ok, result} = Extractor.extract_grid(blob, points)
|
||||
point_data = result[{21.138123, -122.719528}]
|
||||
|
||||
assert Map.has_key?(point_data, "TMP:2 m above ground")
|
||||
assert Map.has_key?(point_data, "DPT:2 m above ground")
|
||||
end
|
||||
|
||||
test "skips points outside the grid" do
|
||||
msg = build_synthetic_grib2_message(0, 0, 0, 103, 2)
|
||||
|
||||
# Mix valid and invalid points
|
||||
points = [
|
||||
{21.138123, -122.719528},
|
||||
{51.5, -0.12}
|
||||
]
|
||||
|
||||
assert {:ok, result} = Extractor.extract_grid(msg, points)
|
||||
|
||||
# Only the valid point should be in the result
|
||||
assert Map.has_key?(result, {21.138123, -122.719528})
|
||||
refute Map.has_key?(result, {51.5, -0.12})
|
||||
end
|
||||
|
||||
test "returns empty map for empty points list" do
|
||||
msg = build_synthetic_grib2_message(0, 0, 0, 103, 2)
|
||||
|
||||
assert {:ok, result} = Extractor.extract_grid(msg, [])
|
||||
assert result == %{}
|
||||
end
|
||||
|
||||
test "handles all points outside grid" do
|
||||
msg = build_synthetic_grib2_message(0, 0, 0, 103, 2)
|
||||
|
||||
points = [{51.5, -0.12}, {48.8, 2.3}]
|
||||
|
||||
assert {:ok, result} = Extractor.extract_grid(msg, points)
|
||||
assert result == %{}
|
||||
end
|
||||
end
|
||||
|
||||
describe "extract_grid/2 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: 300_000
|
||||
test "extracts grid values matching extract_points for same location" do
|
||||
if File.exists?(@fixture_path) do
|
||||
fixture = File.read!(@fixture_path)
|
||||
|
||||
lat = 32.78
|
||||
lon = -96.8
|
||||
|
||||
{:ok, single_result} = Extractor.extract_points(fixture, lat, lon)
|
||||
{:ok, grid_result} = Extractor.extract_grid(fixture, [{lat, lon}])
|
||||
|
||||
point_data = grid_result[{lat, lon}]
|
||||
|
||||
for {key, value} <- single_result do
|
||||
assert_in_delta point_data[key], value, 0.001, "#{key}: grid #{point_data[key]} != single #{value}"
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
@tag :fixture
|
||||
@tag timeout: 300_000
|
||||
test "extracts grid values for multiple US cities" do
|
||||
if File.exists?(@fixture_path) do
|
||||
fixture = File.read!(@fixture_path)
|
||||
|
||||
points = [
|
||||
{32.78, -96.8},
|
||||
{40.71, -74.01},
|
||||
{34.05, -118.24}
|
||||
]
|
||||
|
||||
{:ok, result} = Extractor.extract_grid(fixture, points)
|
||||
|
||||
for point <- points do
|
||||
assert Map.has_key?(result, point), "Missing point #{inspect(point)}"
|
||||
assert Map.has_key?(result[point], "TMP:2 m above ground")
|
||||
|
||||
temp_k = result[point]["TMP:2 m above ground"]
|
||||
temp_c = temp_k - 273.15
|
||||
|
||||
assert temp_c > -30.0 and temp_c < 50.0,
|
||||
"#{inspect(point)} temp #{temp_c}C unreasonable"
|
||||
end
|
||||
end
|
||||
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"
|
||||
|
|
|
|||
|
|
@ -165,4 +165,110 @@ defmodule Microwaveprop.Weather.Grib2.SimplePackingTest do
|
|||
assert {:error, :index_out_of_range} = SimplePacking.extract_value(params, data, 5)
|
||||
end
|
||||
end
|
||||
|
||||
describe "extract_values/3" do
|
||||
test "extracts multiple values in a single batch call" do
|
||||
params = %{
|
||||
reference_value: 0.0,
|
||||
binary_scale: 0,
|
||||
decimal_scale: 0,
|
||||
bits_per_value: 16,
|
||||
num_data_points: 4
|
||||
}
|
||||
|
||||
data = <<100::16-big, 200::16-big, 300::16-big, 400::16-big>>
|
||||
|
||||
assert {:ok, result} = SimplePacking.extract_values(params, data, [0, 2, 3])
|
||||
assert map_size(result) == 3
|
||||
assert_in_delta result[0], 100.0, 0.001
|
||||
assert_in_delta result[2], 300.0, 0.001
|
||||
assert_in_delta result[3], 400.0, 0.001
|
||||
end
|
||||
|
||||
test "skips out-of-range indices" do
|
||||
params = %{
|
||||
reference_value: 0.0,
|
||||
binary_scale: 0,
|
||||
decimal_scale: 0,
|
||||
bits_per_value: 16,
|
||||
num_data_points: 2
|
||||
}
|
||||
|
||||
data = <<100::16-big, 200::16-big>>
|
||||
|
||||
assert {:ok, result} = SimplePacking.extract_values(params, data, [0, 5, -1, 1])
|
||||
assert map_size(result) == 2
|
||||
assert_in_delta result[0], 100.0, 0.001
|
||||
assert_in_delta result[1], 200.0, 0.001
|
||||
end
|
||||
|
||||
test "returns all same value when bits_per_value is 0" do
|
||||
params = %{
|
||||
reference_value: 293.5,
|
||||
binary_scale: 0,
|
||||
decimal_scale: 0,
|
||||
bits_per_value: 0,
|
||||
num_data_points: 10
|
||||
}
|
||||
|
||||
assert {:ok, result} = SimplePacking.extract_values(params, <<>>, [0, 3, 7])
|
||||
assert map_size(result) == 3
|
||||
assert_in_delta result[0], 293.5, 0.001
|
||||
assert_in_delta result[3], 293.5, 0.001
|
||||
assert_in_delta result[7], 293.5, 0.001
|
||||
end
|
||||
|
||||
test "applies scaling correctly to batch" do
|
||||
# value = (R + X * 2^E) * 10^(-D)
|
||||
# R=100.0, E=1, D=1 -> (100 + X*2) / 10
|
||||
params = %{
|
||||
reference_value: 100.0,
|
||||
binary_scale: 1,
|
||||
decimal_scale: 1,
|
||||
bits_per_value: 16,
|
||||
num_data_points: 2
|
||||
}
|
||||
|
||||
data = <<50::16-big, 75::16-big>>
|
||||
|
||||
assert {:ok, result} = SimplePacking.extract_values(params, data, [0, 1])
|
||||
# (100 + 50*2)/10 = 20.0
|
||||
assert_in_delta result[0], 20.0, 0.001
|
||||
# (100 + 75*2)/10 = 25.0
|
||||
assert_in_delta result[1], 25.0, 0.001
|
||||
end
|
||||
|
||||
test "handles empty indices list" do
|
||||
params = %{
|
||||
reference_value: 0.0,
|
||||
binary_scale: 0,
|
||||
decimal_scale: 0,
|
||||
bits_per_value: 16,
|
||||
num_data_points: 2
|
||||
}
|
||||
|
||||
data = <<100::16-big, 200::16-big>>
|
||||
|
||||
assert {:ok, result} = SimplePacking.extract_values(params, data, [])
|
||||
assert result == %{}
|
||||
end
|
||||
|
||||
test "handles non-byte-aligned 12-bit batch extraction" do
|
||||
params = %{
|
||||
reference_value: 0.0,
|
||||
binary_scale: 0,
|
||||
decimal_scale: 0,
|
||||
bits_per_value: 12,
|
||||
num_data_points: 3
|
||||
}
|
||||
|
||||
# Pack 3 12-bit values: 0xABC=2748, 0x123=291, 0xDEF=3567
|
||||
data = <<0xAB, 0xC1, 0x23, 0xDE, 0xF0>>
|
||||
|
||||
assert {:ok, result} = SimplePacking.extract_values(params, data, [0, 1, 2])
|
||||
assert_in_delta result[0], 2748.0, 0.001
|
||||
assert_in_delta result[1], 291.0, 0.001
|
||||
assert_in_delta result[2], 3567.0, 0.001
|
||||
end
|
||||
end
|
||||
end
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue