diff --git a/lib/microwaveprop/weather/grib2/complex_packing.ex b/lib/microwaveprop/weather/grib2/complex_packing.ex index 978aadfc..85eb386e 100644 --- a/lib/microwaveprop/weather/grib2/complex_packing.ex +++ b/lib/microwaveprop/weather/grib2/complex_packing.ex @@ -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. diff --git a/lib/microwaveprop/weather/grib2/extractor.ex b/lib/microwaveprop/weather/grib2/extractor.ex index 71872925..cffc9512 100644 --- a/lib/microwaveprop/weather/grib2/extractor.ex +++ b/lib/microwaveprop/weather/grib2/extractor.ex @@ -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 diff --git a/lib/microwaveprop/weather/grib2/simple_packing.ex b/lib/microwaveprop/weather/grib2/simple_packing.ex index 7c2b97e4..fa197d2d 100644 --- a/lib/microwaveprop/weather/grib2/simple_packing.ex +++ b/lib/microwaveprop/weather/grib2/simple_packing.ex @@ -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) diff --git a/test/microwaveprop/weather/grib2/complex_packing_test.exs b/test/microwaveprop/weather/grib2/complex_packing_test.exs index e0b40081..19bedcd6 100644 --- a/test/microwaveprop/weather/grib2/complex_packing_test.exs +++ b/test/microwaveprop/weather/grib2/complex_packing_test.exs @@ -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 diff --git a/test/microwaveprop/weather/grib2/extractor_test.exs b/test/microwaveprop/weather/grib2/extractor_test.exs index 6acf016e..79d08d9c 100644 --- a/test/microwaveprop/weather/grib2/extractor_test.exs +++ b/test/microwaveprop/weather/grib2/extractor_test.exs @@ -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" diff --git a/test/microwaveprop/weather/grib2/simple_packing_test.exs b/test/microwaveprop/weather/grib2/simple_packing_test.exs index f2d4a863..9915838d 100644 --- a/test/microwaveprop/weather/grib2/simple_packing_test.exs +++ b/test/microwaveprop/weather/grib2/simple_packing_test.exs @@ -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