Add batch GRIB2 extraction for multi-point grid queries

Add extract_values/3 to SimplePacking and ComplexPacking for batch
index extraction from a single GRIB2 message. Add extract_grid/2 to
Extractor which takes a list of {lat, lon} points and returns all
variable values for each point, skipping points outside the grid.

This enables extracting weather data for many grid points from a
single HRRR download instead of re-parsing per point.

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.
 

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

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)

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

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"

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