Merge branch 'worktree-agent-a791a51b' into feature/propagation-map

lib/microwaveprop/weather/grib2/extractor.ex

@@ -26,6 +26,34 @@ defmodule Microwaveprop.Weather.Grib2.Extractor do
     results
   end
 
+  @doc """
+  Extract weather values from a GRIB2 binary blob for multiple lat/lon points.
+
+  Takes a binary and a list of `{lat, lon}` tuples.
+  Returns `{:ok, %{{lat, lon} => %{"VAR:LEVEL" => float}}}` or `{:error, term}`.
+  """
+  def extract_grid(binary, points) do
+    messages = split_messages(binary)
+
+    Enum.reduce_while(messages, {:ok, init_grid(points)}, fn msg, {:ok, acc} ->
+      case extract_single_grid(msg, points) do
+        {:ok, point_values} ->
+          merged =
+            Enum.reduce(point_values, acc, fn {point, key, value}, grid ->
+              Map.update!(grid, point, &Map.put(&1, key, value))
+            end)
+
+          {:cont, {:ok, merged}}
+
+        {:error, :outside_grid} ->
+          {:halt, {:error, :outside_grid}}
+
+        {:error, reason} ->
+          {:halt, {:error, reason}}
+      end
+    end)
+  end
+
   @doc """
   Split a binary blob into individual GRIB2 messages by scanning for "GRIB" magic bytes
   and reading the total length from the indicator section.
@@ -82,6 +110,40 @@ defmodule Microwaveprop.Weather.Grib2.Extractor do
     split_messages(rest, acc)
   end
 
+  defp init_grid(points) do
+    Map.new(points, fn point -> {point, %{}} end)
+  end
+
+  defp extract_single_grid(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}"
+
+      results =
+        Enum.reduce_while(points, {:ok, []}, fn {lat, lon} = point, {:ok, acc} ->
+          case LambertConformal.to_grid_index(grid, lat, lon) do
+            {:ok, {i, j}} ->
+              index = linear_index({i, j}, grid.nx, grid.scan_mode)
+
+              case unpack_value(packing, data, index) do
+                {:ok, value} -> {:cont, {:ok, [{point, key, value} | acc]}}
+                {:error, reason} -> {:halt, {:error, reason}}
+              end
+
+            {:error, :outside_grid} ->
+              {:halt, {:error, :outside_grid}}
+          end
+        end)
+
+      case results do
+        {:ok, point_values} -> {:ok, point_values}
+        {:error, reason} -> {:error, reason}
+      end
+    end
+  rescue
+    e -> {:error, "GRIB2 grid extraction failed: #{inspect(e)}"}
+  end
+
   defp extract_single(msg, lat, lon) do
     with {:ok, parsed} <- Section.parse_message(msg),
          %{grid_params: grid, product: prod, packing_params: packing, data: data} = parsed,

lib/microwaveprop/weather/grib2/section.ex

@@ -28,6 +28,10 @@ defmodule Microwaveprop.Weather.Grib2.Section do
   def identify_var(0, 3, 5), do: "HGT"
   def identify_var(0, 3, 18), do: "HPBL"
   def identify_var(0, 1, 3), do: "PWAT"
+  def identify_var(0, 1, 8), do: "APCP"
+  def identify_var(0, 2, 2), do: "UGRD"
+  def identify_var(0, 2, 3), do: "VGRD"
+  def identify_var(0, 6, 1), do: "TCDC"
   def identify_var(_discipline, cat, num), do: "#{cat}:#{num}"
 
   @doc """
@@ -37,6 +41,7 @@ defmodule Microwaveprop.Weather.Grib2.Section do
   def identify_level(100, value_pa), do: "#{div(value_pa, 100)} mb"
   def identify_level(103, value_m), do: "#{value_m} m above ground"
 
+  def identify_level(10, _value), do: "entire atmosphere"
   def identify_level(200, _value), do: "entire atmosphere (considered as a single layer)"
 
   def identify_level(type, value), do: "unknown:#{type}:#{value}"

lib/microwaveprop/weather/hrrr_client.ex

@@ -14,11 +14,36 @@ defmodule Microwaveprop.Weather.HrrrClient do
     %{var: "DPT", level: "2 m above ground"},
     %{var: "PRES", level: "surface"},
     %{var: "HPBL", level: "surface"},
-    %{var: "PWAT", level: "entire atmosphere (considered as a single layer)"}
+    %{var: "PWAT", level: "entire atmosphere (considered as a single layer)"},
+    %{var: "UGRD", level: "10 m above ground"},
+    %{var: "VGRD", level: "10 m above ground"},
+    %{var: "TCDC", level: "entire atmosphere"},
+    %{var: "APCP", level: "surface"}
   ]
 
   # --- Public API ---
 
+  def surface_messages, do: @surface_messages
+
+  def fetch_grid(points, valid_time, opts \\ []) do
+    hour_dt = nearest_hrrr_hour(valid_time)
+    date = DateTime.to_date(hour_dt)
+    hour = hour_dt.hour
+
+    with {:ok, sfc_grid} <- fetch_product_grid(date, hour, :surface, points),
+         {:ok, prs_grid} <- maybe_fetch_pressure_grid(date, hour, points, opts) do
+      merged = merge_grid_data(sfc_grid, prs_grid)
+
+      profiles =
+        Map.new(merged, fn {point, data} ->
+          profile = build_profile(data)
+          {point, Map.put(profile, :run_time, hour_dt)}
+        end)
+
+      {:ok, profiles}
+    end
+  end
+
   def fetch_profile(lat, lon, valid_time) do
     hour_dt = nearest_hrrr_hour(valid_time)
     date = DateTime.to_date(hour_dt)
@@ -92,6 +117,12 @@ defmodule Microwaveprop.Weather.HrrrClient do
     end)
   end
 
+  def merge_grid_data(sfc_grid, prs_grid) do
+    Map.merge(sfc_grid, prs_grid, fn _point, sfc_data, prs_data ->
+      Map.merge(sfc_data, prs_data)
+    end)
+  end
+
   def build_profile(parsed) do
     sfc_temp_k = parsed["TMP:2 m above ground"]
     sfc_dpt_k = parsed["DPT:2 m above ground"]
@@ -124,6 +155,10 @@ defmodule Microwaveprop.Weather.HrrrClient do
       surface_pressure_mb: if(sfc_pres_pa, do: sfc_pres_pa / 100.0),
       hpbl_m: parsed["HPBL:surface"],
       pwat_mm: parsed["PWAT:entire atmosphere (considered as a single layer)"],
+      wind_u: parsed["UGRD:10 m above ground"],
+      wind_v: parsed["VGRD:10 m above ground"],
+      cloud_cover_pct: parsed["TCDC:entire atmosphere"],
+      precip_mm: parsed["APCP:surface"],
       profile: profile
     }
   end
@@ -152,6 +187,36 @@ defmodule Microwaveprop.Weather.HrrrClient do
     end
   end
 
+  defp fetch_product_grid(date, hour, product, points) do
+    url = hrrr_url(date, hour, product)
+    idx_url = url <> ".idx"
+
+    wanted =
+      case product do
+        :surface -> @surface_messages
+        :pressure -> pressure_messages()
+      end
+
+    Logger.info("HRRR grid fetching #{product} idx from #{idx_url}")
+
+    with {:ok, idx_text} <- fetch_idx(idx_url),
+         idx_entries = parse_idx(idx_text),
+         ranges = byte_ranges_for_messages(idx_entries, wanted),
+         _ = Logger.info("HRRR grid downloading #{length(ranges)} GRIB ranges for #{product}"),
+         {:ok, grib_binary} <- download_grib_ranges(url, ranges) do
+      Logger.info("HRRR grid #{product} downloaded #{byte_size(grib_binary)} bytes, extracting")
+      Extractor.extract_grid(grib_binary, points)
+    end
+  end
+
+  defp maybe_fetch_pressure_grid(date, hour, points, opts) do
+    if Keyword.get(opts, :pressure, true) do
+      fetch_product_grid(date, hour, :pressure, points)
+    else
+      {:ok, %{}}
+    end
+  end
+
   defp pressure_messages do
     for level <- @pressure_levels, var <- ["TMP", "DPT", "HGT"] do
       %{var: var, level: "#{level} mb"}

test/microwaveprop/weather/hrrr_client_test.exs

@@ -187,6 +187,44 @@ defmodule Microwaveprop.Weather.HrrrClientTest do
       assert first_level["hght"] == 110.0
     end
 
+    test "includes wind, cloud cover, and precip fields" do
+      parsed = %{
+        "TMP:2 m above ground" => 299.0,
+        "DPT:2 m above ground" => 292.0,
+        "PRES:surface" => 101_350.0,
+        "HPBL:surface" => 1500.0,
+        "PWAT:entire atmosphere (considered as a single layer)" => 25.0,
+        "UGRD:10 m above ground" => 3.5,
+        "VGRD:10 m above ground" => -2.1,
+        "TCDC:entire atmosphere" => 75.0,
+        "APCP:surface" => 1.2
+      }
+
+      result = HrrrClient.build_profile(parsed)
+
+      assert result.wind_u == 3.5
+      assert result.wind_v == -2.1
+      assert result.cloud_cover_pct == 75.0
+      assert result.precip_mm == 1.2
+    end
+
+    test "wind, cloud, and precip fields are nil when missing" do
+      parsed = %{
+        "TMP:2 m above ground" => 299.0,
+        "DPT:2 m above ground" => 292.0,
+        "PRES:surface" => 101_350.0,
+        "HPBL:surface" => 1500.0,
+        "PWAT:entire atmosphere (considered as a single layer)" => 25.0
+      }
+
+      result = HrrrClient.build_profile(parsed)
+
+      assert result.wind_u == nil
+      assert result.wind_v == nil
+      assert result.cloud_cover_pct == nil
+      assert result.precip_mm == nil
+    end
+
     test "skips pressure levels with missing data" do
       parsed = %{
         "TMP:1000 mb" => 298.0,
@@ -206,4 +244,112 @@ defmodule Microwaveprop.Weather.HrrrClientTest do
       assert length(result.profile) == 1
     end
   end
+
+  describe "surface_messages/0" do
+    test "returns list of surface message descriptors" do
+      messages = HrrrClient.surface_messages()
+
+      assert is_list(messages)
+      assert length(messages) == 9
+
+      vars = Enum.map(messages, & &1.var)
+      assert "TMP" in vars
+      assert "DPT" in vars
+      assert "PRES" in vars
+      assert "HPBL" in vars
+      assert "PWAT" in vars
+      assert "UGRD" in vars
+      assert "VGRD" in vars
+      assert "TCDC" in vars
+      assert "APCP" in vars
+    end
+
+    test "includes wind messages at 10 m above ground" do
+      messages = HrrrClient.surface_messages()
+
+      ugrd = Enum.find(messages, &(&1.var == "UGRD"))
+      assert ugrd.level == "10 m above ground"
+
+      vgrd = Enum.find(messages, &(&1.var == "VGRD"))
+      assert vgrd.level == "10 m above ground"
+    end
+
+    test "includes cloud cover for entire atmosphere" do
+      messages = HrrrClient.surface_messages()
+      tcdc = Enum.find(messages, &(&1.var == "TCDC"))
+      assert tcdc.level == "entire atmosphere"
+    end
+
+    test "includes precipitation at surface" do
+      messages = HrrrClient.surface_messages()
+      apcp = Enum.find(messages, &(&1.var == "APCP"))
+      assert apcp.level == "surface"
+    end
+  end
+
+  describe "merge_grid_data/2" do
+    test "merges surface and pressure grids by point key" do
+      point_a = {32.90, -97.04}
+      point_b = {33.10, -96.80}
+
+      sfc_grid = %{
+        point_a => %{
+          "TMP:2 m above ground" => 299.0,
+          "DPT:2 m above ground" => 292.0,
+          "PRES:surface" => 101_350.0,
+          "HPBL:surface" => 1500.0,
+          "PWAT:entire atmosphere (considered as a single layer)" => 25.0,
+          "UGRD:10 m above ground" => 3.5,
+          "VGRD:10 m above ground" => -2.1,
+          "TCDC:entire atmosphere" => 75.0,
+          "APCP:surface" => 1.2
+        },
+        point_b => %{
+          "TMP:2 m above ground" => 297.0,
+          "DPT:2 m above ground" => 290.0,
+          "PRES:surface" => 101_200.0,
+          "HPBL:surface" => 1200.0,
+          "PWAT:entire atmosphere (considered as a single layer)" => 22.0,
+          "UGRD:10 m above ground" => 2.0,
+          "VGRD:10 m above ground" => -1.5,
+          "TCDC:entire atmosphere" => 50.0,
+          "APCP:surface" => 0.0
+        }
+      }
+
+      prs_grid = %{
+        point_a => %{
+          "TMP:1000 mb" => 298.0,
+          "DPT:1000 mb" => 291.0,
+          "HGT:1000 mb" => 110.0
+        },
+        point_b => %{
+          "TMP:1000 mb" => 296.0,
+          "DPT:1000 mb" => 289.0,
+          "HGT:1000 mb" => 105.0
+        }
+      }
+
+      merged = HrrrClient.merge_grid_data(sfc_grid, prs_grid)
+
+      assert map_size(merged) == 2
+      assert merged[point_a]["TMP:2 m above ground"] == 299.0
+      assert merged[point_a]["TMP:1000 mb"] == 298.0
+      assert merged[point_b]["TMP:2 m above ground"] == 297.0
+      assert merged[point_b]["TMP:1000 mb"] == 296.0
+    end
+
+    test "handles empty pressure data" do
+      point = {32.90, -97.04}
+
+      sfc_grid = %{
+        point => %{"TMP:2 m above ground" => 299.0}
+      }
+
+      merged = HrrrClient.merge_grid_data(sfc_grid, %{})
+
+      assert map_size(merged) == 1
+      assert merged[point]["TMP:2 m above ground"] == 299.0
+    end
+  end
 end