diff --git a/lib/microwaveprop/weather/grib2/extractor.ex b/lib/microwaveprop/weather/grib2/extractor.ex index 71872925..00fc2d5b 100644 --- a/lib/microwaveprop/weather/grib2/extractor.ex +++ b/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. @@ -60,6 +88,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, diff --git a/lib/microwaveprop/weather/grib2/section.ex b/lib/microwaveprop/weather/grib2/section.ex index c3b1e0bc..67c5bc7d 100644 --- a/lib/microwaveprop/weather/grib2/section.ex +++ b/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}" diff --git a/lib/microwaveprop/weather/hrrr_client.ex b/lib/microwaveprop/weather/hrrr_client.ex index 21db10ef..41cad034 100644 --- a/lib/microwaveprop/weather/hrrr_client.ex +++ b/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"} diff --git a/test/microwaveprop/weather/hrrr_client_test.exs b/test/microwaveprop/weather/hrrr_client_test.exs index 49c4fa84..9598b622 100644 --- a/test/microwaveprop/weather/hrrr_client_test.exs +++ b/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