diff --git a/.gitignore b/.gitignore index 2f11c90c..2afcd194 100644 --- a/.gitignore +++ b/.gitignore @@ -31,6 +31,9 @@ microwaveprop-*.tar .env .env.* +# GRIB2 test fixtures (large binary files, downloaded on-demand) +/test/fixtures/grib2/*.grib2 + # OS files .DS_Store Thumbs.db diff --git a/config/test.exs b/config/test.exs index 6a5bb192..13f5f531 100644 --- a/config/test.exs +++ b/config/test.exs @@ -29,6 +29,7 @@ config :microwaveprop, MicrowavepropWeb.Endpoint, # Run Oban jobs inline during tests config :microwaveprop, Oban, testing: :inline config :microwaveprop, elevation_req_options: [plug: {Req.Test, Microwaveprop.Terrain.ElevationClient}, retry: false] +config :microwaveprop, hrrr_req_options: [plug: {Req.Test, Microwaveprop.Weather.HrrrClient}, retry: false] config :microwaveprop, iem_req_options: [plug: {Req.Test, Microwaveprop.Weather.IemClient}, retry: false] # Route HTTP requests through Req.Test stubs diff --git a/lib/microwaveprop/weather/grib2/complex_packing.ex b/lib/microwaveprop/weather/grib2/complex_packing.ex new file mode 100644 index 00000000..8912f78c --- /dev/null +++ b/lib/microwaveprop/weather/grib2/complex_packing.ex @@ -0,0 +1,229 @@ +defmodule Microwaveprop.Weather.Grib2.ComplexPacking do + @moduledoc false + + @doc """ + Extract a single value from GRIB2 complex-packed data with spatial differencing + (Template 5.3) at the given grid index. + + Decodes all values since spatial differencing requires sequential access, + then returns the value at the target index. + """ + def extract_value(params, data, index) do + if index < 0 or index >= params.num_data_points do + {:error, :index_out_of_range} + else + case decode_all(params, data) do + {:ok, values} -> {:ok, :array.get(index, values)} + {:error, _} = err -> err + end + end + end + + @doc """ + Decode all values from complex-packed data with spatial differencing. + + Returns an Erlang array of floats for O(1) index access. + """ + def decode_all(params, data) do + %{ + reference_value: ref, + binary_scale: e, + decimal_scale: d, + bits_per_value: nbits, + num_groups: num_groups, + ref_group_widths: ref_gw, + nbits_group_widths: nbits_gw, + ref_group_lengths: ref_gl, + length_increment: len_inc, + last_group_length: last_gl, + nbits_group_lengths: nbits_gl, + spatial_order: spatial_order, + num_extra_octets: num_extra_octets + } = params + + # Step 1: Extract spatial differencing initial values + octets_per_val = num_extra_octets + num_init_vals = spatial_order + 1 + + {init_vals, rest} = extract_init_values(data, num_init_vals, octets_per_val) + + {spatial_init, [overall_min]} = Enum.split(init_vals, spatial_order) + + # Step 2: Extract group reference values (byte-padded per GRIB2 spec) + {group_refs, rest} = extract_n_values_array(rest, num_groups, nbits) + + # Step 3: Extract group widths (byte-padded) + {group_widths_arr, rest} = extract_n_values_array(rest, num_groups, nbits_gw) + + # Step 4: Extract group lengths (byte-padded) + {group_lengths_arr, rest} = extract_n_values_array(rest, num_groups, nbits_gl) + + # Step 5: Build group info lists + group_widths = + for g <- 0..(num_groups - 1) do + get_val(group_widths_arr, g) + ref_gw + end + + group_lengths = + for g <- 0..(num_groups - 1) do + if g == num_groups - 1 do + last_gl + else + get_val(group_lengths_arr, g) * len_inc + ref_gl + end + end + + group_refs_list = + for g <- 0..(num_groups - 1) do + get_val(group_refs, g) + end + + # Step 6: Decode packed data for each group using bitstring operations + raw_values = decode_groups_bitwise(rest, group_refs_list, group_widths, group_lengths) + + # Step 7: Apply spatial differencing in reverse + undiffed = apply_spatial_differencing(spatial_init, overall_min, raw_values) + + # Step 8: Apply scaling formula: value = (R + X * 2^E) * 10^(-D) + factor_2e = :math.pow(2, e) + factor_10d = :math.pow(10, -d) + + result = + undiffed + |> :array.to_list() + |> Enum.map(fn x -> (ref + x * factor_2e) * factor_10d end) + |> :array.from_list() + + {:ok, result} + rescue + e in [MatchError, ArgumentError] -> + {:error, "GRIB2 complex packing decode failed: #{inspect(e)}"} + end + + # --- Private helpers --- + + defp extract_init_values(data, count, octets_per_val) do + total_bytes = count * octets_per_val + <> = data + + values = + for i <- 0..(count - 1) do + offset = i * octets_per_val + val_bytes = binary_part(init_bytes, offset, octets_per_val) + decode_signed_big(val_bytes) + end + + {values, rest} + end + + defp decode_signed_big(bytes) do + bits = byte_size(bytes) * 8 + <> = bytes + if sign == 1, do: -magnitude, else: magnitude + end + + defp extract_n_values_array(data, _count, 0) do + {nil, data} + end + + defp extract_n_values_array(data, count, nbits) do + total_bits = count * nbits + total_bytes = div(total_bits + 7, 8) + <> = data + + vals = consume_bits_simple(<>, nbits, []) + # consume_bits_simple returns reversed, take only count values + arr = vals |> Enum.take(count) |> Enum.reverse() |> :array.from_list() + + {arr, rest} + end + + defp consume_bits_simple(bits, nbits, acc) when bit_size(bits) < nbits, do: acc + + defp consume_bits_simple(bits, nbits, acc) do + <> = bits + consume_bits_simple(rest, nbits, [val | acc]) + end + + defp get_val(nil, _index), do: 0 + defp get_val(arr, index), do: :array.get(index, arr) + + # Decode all groups from packed data using a bitstring cursor. + # Returns an Erlang array of raw (pre-differencing) integer values. + defp decode_groups_bitwise(data, group_refs, group_widths, group_lengths) do + {all_values_reversed, _remaining_bits} = + [group_refs, group_widths, group_lengths] + |> Enum.zip() + |> Enum.reduce({[], data}, fn {gref, width, length}, {acc, remaining} -> + if width == 0 do + # All values equal gref — prepend in reverse (same values, order irrelevant) + new_acc = prepend_n(acc, gref, length) + {new_acc, remaining} + else + {new_acc, rest} = consume_group_bits(remaining, width, gref, length, acc) + {new_acc, rest} + end + end) + + all_values_reversed + |> Enum.reverse() + |> :array.from_list() + end + + # Consume `count` values of `width` bits each from bitstring, prepending to acc. + # Values are prepended in reverse order (last value first). + # Handles trailing padding: if fewer bits remain than needed, pad with gref. + defp consume_group_bits(bits, width, gref, count, acc) do + total_bits = count * width + available = bit_size(bits) + + if available >= total_bits do + <> = bits + new_acc = consume_bits(chunk, width, gref, acc) + {new_acc, rest} + else + # Extract what we can, fill remainder with gref + new_acc = consume_bits(bits, width, gref, acc) + extracted = div(available, width) + remaining = count - extracted + new_acc = prepend_n(new_acc, gref, remaining) + {new_acc, <<>>} + end + end + + defp consume_bits(<<>>, _width, _gref, acc), do: acc + + defp consume_bits(bits, width, gref, acc) do + <> = bits + consume_bits(rest, width, gref, [val + gref | acc]) + end + + defp prepend_n(acc, _value, 0), do: acc + defp prepend_n(acc, value, n), do: prepend_n([value | acc], value, n - 1) + + defp apply_spatial_differencing(spatial_init, overall_min, raw_values) do + raw_list = :array.to_list(raw_values) + + case spatial_init do + [ival1] -> + {result_reversed, _} = + Enum.reduce(tl(raw_list), {[ival1], ival1}, fn raw, {acc, prev} -> + current = prev + raw + overall_min + {[current | acc], current} + end) + + :array.from_list(Enum.reverse(result_reversed)) + + [ival1, ival2] -> + {result_reversed, _, _} = + raw_list + |> Enum.drop(2) + |> Enum.reduce({[ival2, ival1], ival2, ival1}, fn raw, {acc, prev1, prev2} -> + current = raw + overall_min + 2 * prev1 - prev2 + {[current | acc], current, prev1} + end) + + :array.from_list(Enum.reverse(result_reversed)) + end + end +end diff --git a/lib/microwaveprop/weather/grib2/extractor.ex b/lib/microwaveprop/weather/grib2/extractor.ex new file mode 100644 index 00000000..71872925 --- /dev/null +++ b/lib/microwaveprop/weather/grib2/extractor.ex @@ -0,0 +1,83 @@ +defmodule Microwaveprop.Weather.Grib2.Extractor do + @moduledoc false + + alias Microwaveprop.Weather.Grib2.ComplexPacking + alias Microwaveprop.Weather.Grib2.LambertConformal + alias Microwaveprop.Weather.Grib2.Section + alias Microwaveprop.Weather.Grib2.SimplePacking + + @doc """ + Extract weather values from a GRIB2 binary blob at the given lat/lon. + + Returns `{:ok, %{"VAR:LEVEL" => float}}` or `{:error, term}`. + """ + def extract_points(binary, lat, lon) do + messages = split_messages(binary) + + results = + Enum.reduce_while(messages, {:ok, %{}}, fn msg, {:ok, acc} -> + case extract_single(msg, lat, lon) do + {:ok, key, value} -> {:cont, {:ok, Map.put(acc, key, value)}} + {:error, :outside_grid} -> {:halt, {:error, :outside_grid}} + {:error, reason} -> {:halt, {:error, reason}} + end + end) + + results + 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. + """ + def split_messages(binary), do: split_messages(binary, []) + + @doc """ + Compute linear index from grid (i, j) and scan mode. + + For HRRR scan_mode=64 (j positive, i positive, i consecutive): index = j * nx + i + """ + def linear_index({i, j}, nx, _scan_mode) do + j * nx + i + end + + # --- Private --- + + defp split_messages(<<>>, acc), do: Enum.reverse(acc) + + defp split_messages(<<"GRIB", _::32, total_length::64-big, _rest::binary>> = binary, acc) do + if total_length > byte_size(binary) do + Enum.reverse(acc) + else + msg = binary_part(binary, 0, total_length) + remaining = binary_part(binary, total_length, byte_size(binary) - total_length) + split_messages(remaining, [msg | acc]) + end + end + + defp split_messages(<<_::8, rest::binary>>, acc) do + # Skip non-GRIB bytes (padding between messages) + split_messages(rest, acc) + 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, + key = "#{prod.var}:#{prod.level}", + {:ok, {i, j}} <- LambertConformal.to_grid_index(grid, lat, lon), + index = linear_index({i, j}, grid.nx, grid.scan_mode), + {:ok, value} <- unpack_value(packing, data, index) do + {:ok, key, value} + end + rescue + e -> {:error, "GRIB2 extraction failed: #{inspect(e)}"} + end + + defp unpack_value(%{template: 3} = params, data, index) do + ComplexPacking.extract_value(params, data, index) + end + + defp unpack_value(params, data, index) do + SimplePacking.extract_value(params, data, index) + end +end diff --git a/lib/microwaveprop/weather/grib2/lambert_conformal.ex b/lib/microwaveprop/weather/grib2/lambert_conformal.ex new file mode 100644 index 00000000..633ad995 --- /dev/null +++ b/lib/microwaveprop/weather/grib2/lambert_conformal.ex @@ -0,0 +1,84 @@ +defmodule Microwaveprop.Weather.Grib2.LambertConformal do + @moduledoc false + + @r_earth 6_371_229.0 + + @doc """ + Convert a lat/lon to grid indices (i, j) on a Lambert Conformal Conic grid. + + Returns `{:ok, {i, j}}` or `{:error, :outside_grid}`. + """ + def to_grid_index(grid, lat, lon) do + %{ + nx: nx, + ny: ny, + la1: la1, + lo1: lo1, + dx: dx, + dy: dy, + latin1: latin1, + lov: lov + } = grid + + # Normalize longitudes to 0..360 + lo1 = normalize_lon(lo1) + lon = normalize_lon(lon) + lov = normalize_lon(lov) + + latin1_rad = deg_to_rad(latin1) + la1_rad = deg_to_rad(la1) + lat_rad = deg_to_rad(lat) + + # Cone constant (for HRRR, latin1 == latin2 == 38.5, so n = sin(latin1)) + n = :math.sin(latin1_rad) + + # F factor + f = + :math.cos(latin1_rad) * + :math.pow(:math.tan(:math.pi() / 4 + latin1_rad / 2), n) / n + + # rho for a given latitude + rho_origin = rho(la1_rad, n, f) + rho_target = rho(lat_rad, n, f) + + # Theta angles (difference from central meridian, scaled by cone constant) + theta_origin = n * deg_to_rad(lo1 - lov) + theta_target = n * deg_to_rad(adjusted_lon_diff(lon, lov)) + + # Grid coordinates + x_origin = rho_origin * :math.sin(theta_origin) + y_origin = rho_origin * :math.cos(theta_origin) + + x_target = rho_target * :math.sin(theta_target) + y_target = rho_target * :math.cos(theta_target) + + i = round((x_target - x_origin) / dx) + j = round((y_origin - y_target) / dy) + + if i >= 0 and i < nx and j >= 0 and j < ny do + {:ok, {i, j}} + else + {:error, :outside_grid} + end + end + + defp rho(lat_rad, n, f) do + @r_earth * f / :math.pow(:math.tan(:math.pi() / 4 + lat_rad / 2), n) + end + + defp deg_to_rad(deg), do: deg * :math.pi() / 180.0 + + defp normalize_lon(lon) when lon < 0, do: lon + 360.0 + defp normalize_lon(lon), do: lon + 0.0 + + # Compute longitude difference handling wrap-around + defp adjusted_lon_diff(lon, lov) do + diff = lon - lov + + cond do + diff > 180.0 -> diff - 360.0 + diff < -180.0 -> diff + 360.0 + true -> diff + end + end +end diff --git a/lib/microwaveprop/weather/grib2/section.ex b/lib/microwaveprop/weather/grib2/section.ex new file mode 100644 index 00000000..bee48613 --- /dev/null +++ b/lib/microwaveprop/weather/grib2/section.ex @@ -0,0 +1,259 @@ +defmodule Microwaveprop.Weather.Grib2.Section do + @moduledoc false + + @doc """ + Parse a single GRIB2 message binary into its component sections. + + Returns `{:ok, %{grid_params, product, packing_params, data, bitmap}}`. + """ + def parse_message(<<"GRIB", _reserved::16, _discipline::8, 2::8, total_length::64-big, rest::binary>> = msg) do + discipline = :binary.at(msg, 6) + # Skip indicator section (16 bytes), parse remaining sections + message_body = binary_part(rest, 0, total_length - 16) + parse_sections(message_body, %{discipline: discipline}) + end + + def parse_message(<<"GRIB", _reserved::16, _discipline::8, edition::8, _::binary>>) do + {:error, "unsupported GRIB edition #{edition}"} + end + + def parse_message(_), do: {:error, "not a GRIB2 message"} + + @doc """ + Identify variable name from discipline, parameter category, and parameter number. + """ + def identify_var(0, 0, 0), do: "TMP" + def identify_var(0, 0, 6), do: "DPT" + def identify_var(0, 3, 0), do: "PRES" + 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(_discipline, cat, num), do: "#{cat}:#{num}" + + @doc """ + Identify level string from surface type and scaled value. + """ + def identify_level(1, _value), do: "surface" + 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(200, _value), do: "entire atmosphere (considered as a single layer)" + + def identify_level(type, value), do: "unknown:#{type}:#{value}" + + @doc """ + Decode a 16-bit sign-magnitude integer. + """ + def sign_magnitude_16(<<0::1, value::15>>), do: value + def sign_magnitude_16(<<1::1, value::15>>), do: -value + + # --- Private section parsing --- + + defp parse_sections(<<"7777", _rest::binary>>, acc) do + build_result(acc) + end + + defp parse_sections(<>, acc) when length > 5 do + body_length = length - 5 + <> = rest + + acc = parse_section(section_num, body, acc) + parse_sections(remaining, acc) + end + + defp parse_sections(<<>>, acc) do + build_result(acc) + end + + defp parse_sections(_other, _acc) do + {:error, "malformed section"} + end + + defp parse_section(1, _body, acc) do + # Section 1 (Identification) — we don't need anything from it + acc + end + + defp parse_section(3, body, acc) do + << + _source::8, + _num_points::32-big, + _optional_octets::8, + _interpretation::8, + template::16-big, + rest::binary + >> = body + + grid_params = parse_grid_template(template, rest) + Map.put(acc, :grid_params, grid_params) + end + + defp parse_section(4, body, acc) do + << + _num_coords::16-big, + _template::16-big, + param_category::8, + param_number::8, + _generating_process::8, + _background::8, + _forecast_process::8, + _cutoff_hours::16-big, + _cutoff_minutes::8, + _time_unit::8, + _forecast_time::32-big, + first_surface_type::8, + _first_scale_factor::8, + first_scaled_value::32-big, + _rest::binary + >> = body + + discipline = Map.get(acc, :discipline, 0) + var = identify_var(discipline, param_category, param_number) + level = identify_level(first_surface_type, first_scaled_value) + + Map.put(acc, :product, %{var: var, level: level}) + end + + defp parse_section(5, body, acc) do + << + num_data_points::32-big, + template::16-big, + ref_value_bytes::binary-size(4), + binary_scale_bytes::binary-size(2), + decimal_scale_bytes::binary-size(2), + bits_per_value::8, + rest::binary + >> = body + + <> = ref_value_bytes + + base_params = %{ + num_data_points: num_data_points, + template: template, + reference_value: ref_value, + binary_scale: sign_magnitude_16(binary_scale_bytes), + decimal_scale: sign_magnitude_16(decimal_scale_bytes), + bits_per_value: bits_per_value + } + + packing_params = + case template do + 0 -> base_params + 3 -> parse_complex_packing_params(base_params, rest) + _ -> base_params + end + + Map.put(acc, :packing_params, packing_params) + end + + defp parse_section(6, <>, acc) do + bitmap = + if indicator == 255 do + :none + else + :present + end + + Map.put(acc, :bitmap, bitmap) + end + + defp parse_section(7, data, acc) do + Map.put(acc, :data, data) + end + + defp parse_section(_section_num, _body, acc) do + # Skip unknown sections (Section 2, etc.) + acc + end + + defp parse_complex_packing_params(base, rest) do + << + _original_type::8, + _splitting_method::8, + missing_mgmt::8, + _primary_missing::32-big, + _secondary_missing::32-big, + num_groups::32-big, + ref_group_widths::8, + nbits_group_widths::8, + ref_group_lengths::32-big, + length_increment::8, + last_group_length::32-big, + nbits_group_lengths::8, + spatial_order::8, + num_extra_octets::8, + _extra::binary + >> = rest + + Map.merge(base, %{ + missing_mgmt: missing_mgmt, + num_groups: num_groups, + ref_group_widths: ref_group_widths, + nbits_group_widths: nbits_group_widths, + ref_group_lengths: ref_group_lengths, + length_increment: length_increment, + last_group_length: last_group_length, + nbits_group_lengths: nbits_group_lengths, + spatial_order: spatial_order, + num_extra_octets: num_extra_octets + }) + end + + # Template 3.30 — Lambert Conformal Conic + defp parse_grid_template(30, body) do + << + _shape::8, + _scale_factor_radius::8, + _scaled_radius::32-big, + _scale_factor_major::8, + _scaled_major::32-big, + _scale_factor_minor::8, + _scaled_minor::32-big, + nx::32-big, + ny::32-big, + la1_micro::32-big, + lo1_micro::32-big, + _resolution::8, + lad_micro::32-big, + lov_micro::32-big, + dx_mm::32-big, + dy_mm::32-big, + _projection_centre::8, + scan_mode::8, + latin1_micro::32-big, + latin2_micro::32-big, + _lat_sp::32-big, + _lon_sp::32-big, + _rest::binary + >> = body + + %{ + nx: nx, + ny: ny, + la1: la1_micro / 1_000_000, + lo1: lo1_micro / 1_000_000, + lad: lad_micro / 1_000_000, + lov: lov_micro / 1_000_000, + dx: dx_mm / 1000, + dy: dy_mm / 1000, + latin1: latin1_micro / 1_000_000, + latin2: latin2_micro / 1_000_000, + scan_mode: scan_mode + } + end + + defp parse_grid_template(template, _body) do + %{error: "unsupported grid template #{template}"} + end + + defp build_result(acc) do + required = [:grid_params, :product, :packing_params, :data, :bitmap] + + if Enum.all?(required, &Map.has_key?(acc, &1)) do + {:ok, Map.take(acc, required)} + else + missing = Enum.reject(required, &Map.has_key?(acc, &1)) + {:error, "missing sections: #{inspect(missing)}"} + end + end +end diff --git a/lib/microwaveprop/weather/grib2/simple_packing.ex b/lib/microwaveprop/weather/grib2/simple_packing.ex new file mode 100644 index 00000000..7c2b97e4 --- /dev/null +++ b/lib/microwaveprop/weather/grib2/simple_packing.ex @@ -0,0 +1,44 @@ +defmodule Microwaveprop.Weather.Grib2.SimplePacking do + @moduledoc false + + @doc """ + Extract a single value from GRIB2 simple-packed data at the given grid index. + + Formula: value = (R + X * 2^E) * 10^(-D) + + Where R = reference_value, E = binary_scale, D = decimal_scale, + X = raw packed integer at the given index. + """ + def extract_value(%{bits_per_value: 0} = params, _data, _index) do + {:ok, params.reference_value * :math.pow(10, -params.decimal_scale)} + end + + def extract_value(params, data, index) do + %{ + reference_value: r, + binary_scale: e, + decimal_scale: d, + bits_per_value: n, + num_data_points: num + } = params + + if index < 0 or index >= num do + {:error, :index_out_of_range} + else + x = extract_bits(data, index, n) + value = (r + x * :math.pow(2, e)) * :math.pow(10, -d) + {:ok, value} + end + end + + defp extract_bits(data, index, bits_per_value) do + bit_offset = index * bits_per_value + byte_offset = div(bit_offset, 8) + remainder = rem(bit_offset, 8) + + <<_skip::binary-size(byte_offset), _pad::size(remainder), x::size(bits_per_value)-unsigned-big, _rest::bitstring>> = + data + + x + end +end diff --git a/lib/microwaveprop/weather/hrrr_client.ex b/lib/microwaveprop/weather/hrrr_client.ex index 9b7784f1..67e095d6 100644 --- a/lib/microwaveprop/weather/hrrr_client.ex +++ b/lib/microwaveprop/weather/hrrr_client.ex @@ -1,6 +1,8 @@ defmodule Microwaveprop.Weather.HrrrClient do @moduledoc false + alias Microwaveprop.Weather.Grib2.Extractor + @hrrr_base "https://noaa-hrrr-bdp-pds.s3.amazonaws.com" @pressure_levels [1000, 975, 950, 925, 900, 850, 800, 700] @@ -23,7 +25,7 @@ defmodule Microwaveprop.Weather.HrrrClient do with {:ok, sfc_data} <- fetch_product(date, hour, :surface, lat, lon), {:ok, prs_data} <- fetch_product(date, hour, :pressure, lat, lon) do merged = Map.merge(sfc_data, prs_data) - result = build_profile_from_wgrib2(merged) + result = build_profile(merged) {:ok, Map.put(result, :run_time, hour_dt)} end end @@ -88,19 +90,7 @@ defmodule Microwaveprop.Weather.HrrrClient do end) end - def parse_wgrib2_output(text) do - text - |> String.split("\n") - |> Enum.reject(&(&1 == "")) - |> Enum.reduce(%{}, fn line, acc -> - case parse_wgrib2_line(line) do - {:ok, key, value} -> Map.put(acc, key, value) - :skip -> acc - end - end) - end - - def build_profile_from_wgrib2(parsed) do + def build_profile(parsed) do sfc_temp_k = parsed["TMP:2 m above ground"] sfc_dpt_k = parsed["DPT:2 m above ground"] sfc_pres_pa = parsed["PRES:surface"] @@ -151,9 +141,8 @@ defmodule Microwaveprop.Weather.HrrrClient do with {:ok, idx_text} <- fetch_idx(idx_url), idx_entries = parse_idx(idx_text), ranges = byte_ranges_for_messages(idx_entries, wanted), - {:ok, grib_binary} <- download_grib_ranges(url, ranges), - {:ok, output} <- extract_point(grib_binary, lat, lon) do - {:ok, parse_wgrib2_output(output)} + {:ok, grib_binary} <- download_grib_ranges(url, ranges) do + Extractor.extract_points(grib_binary, lat, lon) end end @@ -194,53 +183,6 @@ defmodule Microwaveprop.Weather.HrrrClient do end end - defp extract_point(grib_binary, lat, lon) do - case System.find_executable("wgrib2") do - nil -> - {:error, "wgrib2 not found in PATH. Build from source: https://www.cpc.ncep.noaa.gov/products/wesley/wgrib2/"} - - wgrib2_path -> - tmp_dir = System.tmp_dir!() - tmp_file = Path.join(tmp_dir, "hrrr_#{:erlang.unique_integer([:positive])}.grib2") - - try do - File.write!(tmp_file, grib_binary) - lon_360 = if lon < 0, do: lon + 360, else: lon - - case System.cmd(wgrib2_path, [tmp_file, "-lon", "#{lon_360}", "#{lat}"], stderr_to_stdout: true) do - {output, 0} -> {:ok, output} - {output, _} -> {:error, "wgrib2 failed: #{output}"} - end - after - File.rm(tmp_file) - end - end - end - - defp parse_wgrib2_line(line) do - case Regex.run(~r/val=([0-9eE.+-]+)\s*$/, line) do - [_, val_str] -> - case Float.parse(val_str) do - {val, _} -> - parts = String.split(line, ":") - var = Enum.at(parts, 3) - level = Enum.at(parts, 4) - - if var && level do - {:ok, "#{var}:#{level}", val} - else - :skip - end - - :error -> - :skip - end - - _ -> - :skip - end - end - defp req_options do defaults = [retry: &retry?/2, max_retries: 5, retry_delay: &retry_delay/1] overrides = Application.get_env(:microwaveprop, :hrrr_req_options, []) diff --git a/test/microwaveprop/weather/grib2/complex_packing_test.exs b/test/microwaveprop/weather/grib2/complex_packing_test.exs new file mode 100644 index 00000000..e0b40081 --- /dev/null +++ b/test/microwaveprop/weather/grib2/complex_packing_test.exs @@ -0,0 +1,95 @@ +defmodule Microwaveprop.Weather.Grib2.ComplexPackingTest do + use ExUnit.Case, async: true + + alias Microwaveprop.Weather.Grib2.ComplexPacking + alias Microwaveprop.Weather.Grib2.Section + + @fixture_path "test/fixtures/grib2/hrrr_tmp_2m.grib2" + + describe "extract_value/3 with real HRRR fixture" do + @tag :fixture + test "extracts a temperature value within meteorological range" do + if File.exists?(@fixture_path) do + fixture = File.read!(@fixture_path) + {:ok, parsed} = Section.parse_message(fixture) + + # First grid point (index 0) + {:ok, value} = + ComplexPacking.extract_value(parsed.packing_params, parsed.data, 0) + + # Should be a temperature in Kelvin (roughly 200-340K) + assert value > 200.0 and value < 340.0, + "Expected temperature 200-340K, got #{value}" + else + IO.puts("Skipping fixture test — #{@fixture_path} not found") + :ok + end + end + + @tag :fixture + test "extracts values at multiple indices" do + if File.exists?(@fixture_path) do + fixture = File.read!(@fixture_path) + {:ok, parsed} = Section.parse_message(fixture) + + # Check several indices across the grid + for idx <- [0, 1000, 100_000, 500_000, 1_000_000, 1_905_140] do + {:ok, value} = + ComplexPacking.extract_value(parsed.packing_params, parsed.data, idx) + + assert value > 200.0 and value < 340.0, + "Index #{idx}: expected 200-340K, got #{value}" + end + else + :ok + end + end + + @tag :fixture + test "Dallas area temperature is reasonable for a US city" do + if File.exists?(@fixture_path) do + fixture = File.read!(@fixture_path) + {:ok, parsed} = Section.parse_message(fixture) + + grid = parsed.grid_params + {:ok, {i, j}} = Microwaveprop.Weather.Grib2.LambertConformal.to_grid_index(grid, 32.78, -96.8) + index = j * grid.nx + i + + {:ok, value} = + ComplexPacking.extract_value(parsed.packing_params, parsed.data, index) + + # Dallas in March — Kelvin should be roughly 270-310K (roughly -3C to 37C) + temp_c = value - 273.15 + + assert temp_c > -10.0 and temp_c < 45.0, + "Dallas temp #{Float.round(temp_c, 1)}C seems unreasonable" + else + :ok + end + end + end + + describe "decode_all/2" do + @tag :fixture + test "decodes all values from fixture" do + if File.exists?(@fixture_path) do + fixture = File.read!(@fixture_path) + {:ok, parsed} = Section.parse_message(fixture) + + {:ok, values_arr} = ComplexPacking.decode_all(parsed.packing_params, parsed.data) + + assert :array.size(values_arr) == parsed.packing_params.num_data_points + + # Check a sample of values are in reasonable temperature range + for idx <- [0, 100_000, 500_000, 1_000_000, 1_905_140] do + val = :array.get(idx, values_arr) + + assert val > 200.0 and val < 340.0, + "Index #{idx}: expected 200-340K, got #{val}" + 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 new file mode 100644 index 00000000..02ad9f97 --- /dev/null +++ b/test/microwaveprop/weather/grib2/extractor_test.exs @@ -0,0 +1,243 @@ +defmodule Microwaveprop.Weather.Grib2.ExtractorTest do + use ExUnit.Case, async: true + + alias Microwaveprop.Weather.Grib2.Extractor + + # HRRR grid parameters for building synthetic messages + @hrrr_nx 1799 + @hrrr_ny 1059 + + describe "split_messages/1" do + test "splits concatenated GRIB2 messages" do + msg1 = build_synthetic_grib2_message(0, 0, 0, 103, 2) + msg2 = build_synthetic_grib2_message(0, 0, 6, 103, 2) + blob = msg1 <> msg2 + + messages = Extractor.split_messages(blob) + assert length(messages) == 2 + assert Enum.at(messages, 0) == msg1 + assert Enum.at(messages, 1) == msg2 + end + + test "handles single message" do + msg = build_synthetic_grib2_message(0, 0, 0, 103, 2) + messages = Extractor.split_messages(msg) + assert length(messages) == 1 + end + + test "handles empty binary" do + assert Extractor.split_messages(<<>>) == [] + end + end + + describe "extract_points/3" do + test "extracts single message at first grid point" do + msg = build_synthetic_grib2_message(0, 0, 0, 103, 2) + + # la1/lo1 = first grid point, so index (0,0) → linear index 0 + assert {:ok, result} = Extractor.extract_points(msg, 21.138123, 237.280472) + + assert Map.has_key?(result, "TMP:2 m above ground") + assert is_float(result["TMP:2 m above ground"]) + end + + test "extracts multiple concatenated messages" do + msg1 = build_synthetic_grib2_message(0, 0, 0, 103, 2) + msg2 = build_synthetic_grib2_message(0, 0, 6, 103, 2) + msg3 = build_synthetic_grib2_message(0, 3, 0, 1, 0) + blob = msg1 <> msg2 <> msg3 + + assert {:ok, result} = Extractor.extract_points(blob, 21.138123, 237.280472) + + assert Map.has_key?(result, "TMP:2 m above ground") + assert Map.has_key?(result, "DPT:2 m above ground") + assert Map.has_key?(result, "PRES:surface") + assert map_size(result) == 3 + end + + test "handles negative longitude" do + msg = build_synthetic_grib2_message(0, 0, 0, 103, 2) + + # -122.719528 == 237.280472 + assert {:ok, result} = Extractor.extract_points(msg, 21.138123, -122.719528) + assert Map.has_key?(result, "TMP:2 m above ground") + end + + test "returns error for point outside grid" do + msg = build_synthetic_grib2_message(0, 0, 0, 103, 2) + + # London — way outside HRRR domain + assert {:error, _} = Extractor.extract_points(msg, 51.5, -0.12) + end + end + + describe "linear_index/3" do + test "scan mode 64 (j positive, i positive, i consecutive)" do + nx = 100 + scan_mode = 64 + assert Extractor.linear_index({5, 10}, nx, scan_mode) == 10 * 100 + 5 + end + + test "scan mode 0 (j negative, i positive, i consecutive)" do + # j decreasing means row 0 is top, so index = (ny - 1 - j) * nx + i + # But we just return j * nx + i for now since we don't flip + nx = 100 + scan_mode = 0 + assert Extractor.linear_index({5, 10}, nx, scan_mode) == 10 * 100 + 5 + 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" + + @tag :fixture + @tag timeout: 120_000 + test "extracts TMP:2m from real single-message GRIB2" do + if File.exists?(@fixture_path) do + fixture = File.read!(@fixture_path) + + assert {:ok, result} = Extractor.extract_points(fixture, 32.78, -96.8) + assert Map.has_key?(result, "TMP:2 m above ground") + + temp_k = result["TMP:2 m above ground"] + temp_c = temp_k - 273.15 + assert temp_c > -10.0 and temp_c < 45.0, "Dallas temp #{temp_c}C unreasonable" + end + end + + @tag :fixture + @tag timeout: 300_000 + test "extracts from multi-message GRIB2" do + if File.exists?(@multi_fixture_path) do + fixture = File.read!(@multi_fixture_path) + + assert {:ok, result} = Extractor.extract_points(fixture, 32.78, -96.8) + assert Map.has_key?(result, "TMP:2 m above ground") + assert Map.has_key?(result, "DPT:2 m above ground") + assert map_size(result) == 2 + end + end + end + + # --- Synthetic GRIB2 builder --- + + defp build_synthetic_grib2_message(discipline, param_category, param_number, surface_type, surface_value) do + num_points = @hrrr_nx * @hrrr_ny + + sec1 = section_wrap(1, build_section1()) + + sec3 = + section_wrap(3, build_section3(@hrrr_nx, @hrrr_ny)) + + sec4 = + section_wrap( + 4, + build_section4(param_category, param_number, surface_type, surface_value) + ) + + sec5 = section_wrap(5, build_section5(num_points, 16)) + + sec6 = section_wrap(6, <<255::8>>) + + # Data: num_points * 2 bytes (16-bit values), all set to 42 + data_bytes = :binary.copy(<<42::16-big>>, num_points) + sec7 = section_wrap(7, data_bytes) + + sec8 = "7777" + + body = sec1 <> sec3 <> sec4 <> sec5 <> sec6 <> sec7 <> sec8 + total_length = 16 + byte_size(body) + + sec0 = <<"GRIB", 0::16, discipline::8, 2::8, total_length::64-big>> + sec0 <> body + end + + defp build_section1 do + << + 0::16-big, + 0::16-big, + 2::8, + 1::8, + 1::8, + 2026::16-big, + 3::8, + 28::8, + 18::8, + 0::8, + 0::8, + 0::8, + 1::8 + >> + end + + defp build_section3(nx, ny) do + << + 0::8, + nx * ny::32-big, + 0::8, + 0::8, + 30::16-big, + 6::8, + 0::8, + 0::32-big, + 0::8, + 0::32-big, + 0::8, + 0::32-big, + nx::32-big, + ny::32-big, + 21_138_123::32-big, + 237_280_472::32-big, + 0b00111000::8, + 38_500_000::32-big, + 262_500_000::32-big, + 3_000_000::32-big, + 3_000_000::32-big, + 0::8, + 64::8, + 38_500_000::32-big, + 38_500_000::32-big, + 0::32-big, + 0::32-big + >> + end + + defp build_section4(param_category, param_number, surface_type, surface_value) do + << + 0::16-big, + 0::16-big, + param_category::8, + param_number::8, + 2::8, + 0::8, + 0::8, + 0::16-big, + 0::8, + 1::8, + 0::32-big, + surface_type::8, + 0::8, + surface_value::32-big, + 255::8, + 0::8, + 0::32-big + >> + end + + defp build_section5(num_points, bits_per_value) do + << + num_points::32-big, + 0::16-big, + 0.0::float-32-big, + 0::16-big, + 0::16-big, + bits_per_value::8 + >> + end + + defp section_wrap(section_number, body) do + length = 5 + byte_size(body) + <> + end +end diff --git a/test/microwaveprop/weather/grib2/lambert_conformal_test.exs b/test/microwaveprop/weather/grib2/lambert_conformal_test.exs new file mode 100644 index 00000000..fcf5dc7c --- /dev/null +++ b/test/microwaveprop/weather/grib2/lambert_conformal_test.exs @@ -0,0 +1,70 @@ +defmodule Microwaveprop.Weather.Grib2.LambertConformalTest do + use ExUnit.Case, async: true + + alias Microwaveprop.Weather.Grib2.LambertConformal + + # HRRR grid parameters (from GRIB2 Section 3, Template 3.30) + @hrrr_params %{ + nx: 1799, + ny: 1059, + la1: 21.138123, + lo1: 237.280472, + lad: 38.5, + lov: 262.5, + dx: 3000.0, + dy: 3000.0, + latin1: 38.5, + latin2: 38.5, + scan_mode: 64 + } + + describe "to_grid_index/3" do + test "first grid point maps to (0, 0)" do + assert {:ok, {0, 0}} = LambertConformal.to_grid_index(@hrrr_params, 21.138123, 237.280472) + end + + test "first grid point with negative longitude" do + # 237.280472 in 0-360 == -122.719528 in -180..180 + assert {:ok, {0, 0}} = + LambertConformal.to_grid_index(@hrrr_params, 21.138123, -122.719528) + end + + test "Dallas area returns reasonable indices" do + # Dallas TX: ~32.78, -96.8 (263.2 in 0-360) + # HRRR 3km grid covers CONUS; Dallas should be roughly in the middle + assert {:ok, {i, j}} = LambertConformal.to_grid_index(@hrrr_params, 32.78, -96.8) + + # Indices should be well within grid bounds + assert i > 100 and i < 1700 + assert j > 100 and j < 900 + end + + test "handles negative longitudes correctly" do + {:ok, {i1, j1}} = LambertConformal.to_grid_index(@hrrr_params, 40.0, -105.0) + {:ok, {i2, j2}} = LambertConformal.to_grid_index(@hrrr_params, 40.0, 255.0) + + assert i1 == i2 + assert j1 == j2 + end + + test "returns error for point outside grid" do + # London, UK — well outside CONUS HRRR domain + assert {:error, :outside_grid} = + LambertConformal.to_grid_index(@hrrr_params, 51.5, -0.12) + end + + test "returns error for point south of grid" do + # South America + assert {:error, :outside_grid} = + LambertConformal.to_grid_index(@hrrr_params, -10.0, -80.0) + end + + test "edge of grid returns valid indices" do + # Somewhere near the grid edge — northeast CONUS + result = LambertConformal.to_grid_index(@hrrr_params, 45.0, -70.0) + assert {:ok, {i, j}} = result + assert i >= 0 and i < 1799 + assert j >= 0 and j < 1059 + end + end +end diff --git a/test/microwaveprop/weather/grib2/section_test.exs b/test/microwaveprop/weather/grib2/section_test.exs new file mode 100644 index 00000000..4048402d --- /dev/null +++ b/test/microwaveprop/weather/grib2/section_test.exs @@ -0,0 +1,265 @@ +defmodule Microwaveprop.Weather.Grib2.SectionTest do + use ExUnit.Case, async: true + + alias Microwaveprop.Weather.Grib2.Section + + describe "parse_message/1" do + test "parses a synthetic GRIB2 message" do + msg = build_synthetic_grib2() + assert {:ok, parsed} = Section.parse_message(msg) + + # Grid params from Section 3 + assert parsed.grid_params.nx == 1799 + assert parsed.grid_params.ny == 1059 + assert_in_delta parsed.grid_params.la1, 21.138123, 0.000001 + assert_in_delta parsed.grid_params.lo1, 237.280472, 0.000001 + assert_in_delta parsed.grid_params.dx, 3000.0, 0.1 + assert_in_delta parsed.grid_params.dy, 3000.0, 0.1 + assert_in_delta parsed.grid_params.latin1, 38.5, 0.000001 + assert_in_delta parsed.grid_params.latin2, 38.5, 0.000001 + assert parsed.grid_params.scan_mode == 64 + + # Product from Section 4 + assert parsed.product.var == "TMP" + assert parsed.product.level == "2 m above ground" + + # Packing params from Section 5 + assert parsed.packing_params.num_data_points == 100 + assert parsed.packing_params.bits_per_value == 16 + + # Data from Section 7 + assert is_binary(parsed.data) + + # Bitmap indicator + assert parsed.bitmap == :none + end + + test "rejects non-GRIB2 binary" do + assert {:error, _} = Section.parse_message(<<"NOT_GRIB", 0::64-big>>) + end + + test "rejects GRIB edition 1" do + indicator = <<"GRIB", 0::16, 0::8, 1::8, 100::64-big>> + assert {:error, _} = Section.parse_message(indicator) + end + end + + describe "product identification" do + test "identifies all supported variables" do + # discipline 0 (meteorological) + assert Section.identify_var(0, 0, 0) == "TMP" + assert Section.identify_var(0, 0, 6) == "DPT" + assert Section.identify_var(0, 3, 0) == "PRES" + assert Section.identify_var(0, 3, 5) == "HGT" + assert Section.identify_var(0, 3, 18) == "HPBL" + assert Section.identify_var(0, 1, 3) == "PWAT" + end + + test "returns category:number for unknown variables" do + assert Section.identify_var(0, 99, 99) == "99:99" + end + end + + describe "level identification" do + test "surface" do + assert Section.identify_level(1, 0) == "surface" + end + + test "pressure level in mb" do + # Stored in Pa, displayed in mb + assert Section.identify_level(100, 100_000) == "1000 mb" + assert Section.identify_level(100, 85_000) == "850 mb" + end + + test "meters above ground" do + assert Section.identify_level(103, 2) == "2 m above ground" + assert Section.identify_level(103, 10) == "10 m above ground" + end + + test "entire atmosphere" do + assert Section.identify_level(200, 0) == + "entire atmosphere (considered as a single layer)" + end + + test "unknown level type" do + assert Section.identify_level(255, 0) == "unknown:255:0" + end + end + + describe "sign_magnitude_16/1" do + test "positive value" do + assert Section.sign_magnitude_16(<<0::1, 5::15>>) == 5 + end + + test "negative value" do + assert Section.sign_magnitude_16(<<1::1, 5::15>>) == -5 + end + + test "zero" do + assert Section.sign_magnitude_16(<<0::16>>) == 0 + end + end + + # --- Helpers --- + + defp build_synthetic_grib2 do + # Section 1 (Identification) - minimal + sec1_body = << + # Originating center, subcenter, etc. + 0::16-big, + 0::16-big, + # Master/local tables + 2::8, + 1::8, + # Significance of reference time + 1::8, + # Year, month, day, hour, minute, second + 2026::16-big, + 3::8, + 28::8, + 18::8, + 0::8, + 0::8, + # Production status, type + 0::8, + 1::8 + >> + + sec1 = section_wrap(1, sec1_body) + + # Section 3 (Grid Definition) - Template 3.30 (Lambert Conformal) + sec3_body = << + # Source of grid definition + 0::8, + # Number of data points + 1_905_141::32-big, + # Number of octets for optional list + 0::8, + # Interpretation + 0::8, + # Grid Definition Template Number = 30 (Lambert Conformal) + 30::16-big, + # Shape of earth (6 = spherical, R=6371229) + 6::8, + # Scale factor/value of radius (not used for shape 6) + 0::8, + 0::32-big, + # Scale factor/value of major axis + 0::8, + 0::32-big, + # Scale factor/value of minor axis + 0::8, + 0::32-big, + # Nx, Ny + 1799::32-big, + 1059::32-big, + # La1 in microdegrees + 21_138_123::32-big, + # Lo1 in microdegrees + 237_280_472::32-big, + # Resolution and component flags + 0b00111000::8, + # LaD in microdegrees + 38_500_000::32-big, + # LoV in microdegrees + 262_500_000::32-big, + # Dx in millimetres + 3_000_000::32-big, + # Dy in millimetres + 3_000_000::32-big, + # Projection centre flag + 0::8, + # Scanning mode + 64::8, + # Latin1 in microdegrees + 38_500_000::32-big, + # Latin2 in microdegrees + 38_500_000::32-big, + # Latitude/longitude of southern pole + 0::32-big, + 0::32-big + >> + + sec3 = section_wrap(3, sec3_body) + + # Section 4 (Product Definition) - Template 4.0 + sec4_body = << + # Number of coordinate values + 0::16-big, + # Product Definition Template Number = 0 + 0::16-big, + # Parameter category (0 = Temperature) + 0::8, + # Parameter number (0 = Temperature) + 0::8, + # Type of generating process + 2::8, + # Background + forecast generating process + 0::8, + 0::8, + # Hours/minutes of observational data cutoff + 0::16-big, + 0::8, + # Indicator of unit of time range (1=hour) + 1::8, + # Forecast time + 0::32-big, + # Type of first fixed surface (103 = above ground) + 103::8, + # Scale factor and scaled value of first fixed surface + 0::8, + 2::32-big, + # Type of second fixed surface (255 = missing) + 255::8, + 0::8, + 0::32-big + >> + + sec4 = section_wrap(4, sec4_body) + + # Section 5 (Data Representation) - Template 5.0 (Simple packing) + # Reference value as IEEE 754 float32 + ref_val = <<0.0::float-32-big>> + + sec5_body = << + # Number of data points + 100::32-big, + # Data Representation Template Number = 0 + 0::16-big, + # Reference value (IEEE 754) + ref_val::binary, + # Binary scale factor (sign-magnitude 16-bit) + 0::16-big, + # Decimal scale factor (sign-magnitude 16-bit) + 0::16-big, + # Number of bits per value + 16::8 + >> + + sec5 = section_wrap(5, sec5_body) + + # Section 6 (Bitmap) - no bitmap + sec6_body = <<255::8>> + sec6 = section_wrap(6, sec6_body) + + # Section 7 (Data) - 100 values of 16 bits = 200 bytes + data_bytes = :binary.copy(<<42::16-big>>, 100) + sec7 = section_wrap(7, data_bytes) + + # Section 8 (End) + sec8 = "7777" + + body = sec1 <> sec3 <> sec4 <> sec5 <> sec6 <> sec7 <> sec8 + total_length = 16 + byte_size(body) + + # Section 0 (Indicator) + sec0 = <<"GRIB", 0::16, 0::8, 2::8, total_length::64-big>> + + sec0 <> body + end + + defp section_wrap(section_number, body) do + length = 5 + byte_size(body) + <> + end +end diff --git a/test/microwaveprop/weather/grib2/simple_packing_test.exs b/test/microwaveprop/weather/grib2/simple_packing_test.exs new file mode 100644 index 00000000..f2d4a863 --- /dev/null +++ b/test/microwaveprop/weather/grib2/simple_packing_test.exs @@ -0,0 +1,168 @@ +defmodule Microwaveprop.Weather.Grib2.SimplePackingTest do + use ExUnit.Case, async: true + + alias Microwaveprop.Weather.Grib2.SimplePacking + + describe "extract_value/3" do + test "basic 16-bit unpacking" do + # value = (R + X * 2^E) * 10^(-D) + # R=0.0, E=0, D=0, N=16 → value = X + params = %{ + reference_value: 0.0, + binary_scale: 0, + decimal_scale: 0, + bits_per_value: 16, + num_data_points: 4 + } + + # 4 values: 100, 200, 300, 400 as 16-bit big-endian + data = <<100::16-big, 200::16-big, 300::16-big, 400::16-big>> + + assert {:ok, 100.0} = SimplePacking.extract_value(params, data, 0) + assert {:ok, 200.0} = SimplePacking.extract_value(params, data, 1) + assert {:ok, 300.0} = SimplePacking.extract_value(params, data, 2) + assert {:ok, 400.0} = SimplePacking.extract_value(params, data, 3) + end + + test "applies reference value" do + # R=273.15, E=0, D=0 → value = 273.15 + X + params = %{ + reference_value: 273.15, + binary_scale: 0, + decimal_scale: 0, + bits_per_value: 16, + num_data_points: 2 + } + + data = <<10::16-big, 20::16-big>> + + assert {:ok, val} = SimplePacking.extract_value(params, data, 0) + assert_in_delta val, 283.15, 0.001 + end + + test "applies binary scale factor" do + # R=0, E=2, D=0 → value = X * 2^2 = X * 4 + params = %{ + reference_value: 0.0, + binary_scale: 2, + decimal_scale: 0, + bits_per_value: 16, + num_data_points: 1 + } + + data = <<25::16-big>> + + assert {:ok, val} = SimplePacking.extract_value(params, data, 0) + assert_in_delta val, 100.0, 0.001 + end + + test "applies decimal scale factor" do + # R=0, E=0, D=2 → value = X * 10^(-2) = X / 100 + params = %{ + reference_value: 0.0, + binary_scale: 0, + decimal_scale: 2, + bits_per_value: 16, + num_data_points: 1 + } + + data = <<15_000::16-big>> + + assert {:ok, val} = SimplePacking.extract_value(params, data, 0) + assert_in_delta val, 150.0, 0.001 + end + + test "combined R, E, D" 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: 1 + } + + data = <<50::16-big>> + + assert {:ok, val} = SimplePacking.extract_value(params, data, 0) + # (100 + 50 * 2) / 10 = 200 / 10 = 20.0 + assert_in_delta val, 20.0, 0.001 + end + + test "12-bit non-byte-aligned extraction" do + # 3 values of 12 bits each = 36 bits = 4.5 bytes + # Values: 0xABC=2748, 0x123=291, 0xDEF=3567 + 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, 0x123, 0xDEF + # Binary: 1010_1011_1100 | 0001_0010_0011 | 1101_1110_1111 + # As bytes: AB C1 23 DE F0 (with 4 bits padding) + data = <<0xAB, 0xC1, 0x23, 0xDE, 0xF0>> + + assert {:ok, val0} = SimplePacking.extract_value(params, data, 0) + assert_in_delta val0, 2748.0, 0.001 + + assert {:ok, val1} = SimplePacking.extract_value(params, data, 1) + assert_in_delta val1, 291.0, 0.001 + + assert {:ok, val2} = SimplePacking.extract_value(params, data, 2) + assert_in_delta val2, 3567.0, 0.001 + end + + test "zero bits_per_value returns reference value" do + params = %{ + reference_value: 293.5, + binary_scale: 0, + decimal_scale: 0, + bits_per_value: 0, + num_data_points: 10 + } + + # With 0 bits, data is empty + data = <<>> + + assert {:ok, val} = SimplePacking.extract_value(params, data, 0) + assert_in_delta val, 293.5, 0.001 + + assert {:ok, val} = SimplePacking.extract_value(params, data, 5) + assert_in_delta val, 293.5, 0.001 + end + + test "negative binary scale factor" do + # E=-1 → 2^(-1) = 0.5 + params = %{ + reference_value: 0.0, + binary_scale: -1, + decimal_scale: 0, + bits_per_value: 16, + num_data_points: 1 + } + + data = <<200::16-big>> + + assert {:ok, val} = SimplePacking.extract_value(params, data, 0) + assert_in_delta val, 100.0, 0.001 + end + + test "index out of range returns error" 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 {:error, :index_out_of_range} = SimplePacking.extract_value(params, data, 5) + end + end +end diff --git a/test/microwaveprop/weather/hrrr_client_test.exs b/test/microwaveprop/weather/hrrr_client_test.exs index de47347a..0055c6db 100644 --- a/test/microwaveprop/weather/hrrr_client_test.exs +++ b/test/microwaveprop/weather/hrrr_client_test.exs @@ -117,45 +117,8 @@ defmodule Microwaveprop.Weather.HrrrClientTest do end end - describe "parse_wgrib2_output/1" do - @sample_output """ - 1:0:d=2026032818:TMP:2 m above ground:anl::lon=262.960000,lat=32.900000,val=298.5 - 2:1234567:d=2026032818:DPT:2 m above ground:anl::lon=262.960000,lat=32.900000,val=291.2 - 3:2345678:d=2026032818:PRES:surface:anl::lon=262.960000,lat=32.900000,val=101350 - 4:3456789:d=2026032818:HPBL:surface:anl::lon=262.960000,lat=32.900000,val=1500 - 5:4567890:d=2026032818:PWAT:entire atmosphere (considered as a single layer):anl::lon=262.960000,lat=32.900000,val=25.3 - """ - - test "parses surface fields" do - result = HrrrClient.parse_wgrib2_output(@sample_output) - - assert result["TMP:2 m above ground"] == 298.5 - assert result["DPT:2 m above ground"] == 291.2 - assert result["PRES:surface"] == 101_350.0 - assert result["HPBL:surface"] == 1500.0 - end - - test "handles pressure level data" do - output = """ - 1:0:d=2026032818:TMP:1000 mb:anl::lon=262.960000,lat=32.900000,val=297.3 - 2:100:d=2026032818:TMP:975 mb:anl::lon=262.960000,lat=32.900000,val=295.1 - 3:200:d=2026032818:HGT:1000 mb:anl::lon=262.960000,lat=32.900000,val=110 - """ - - result = HrrrClient.parse_wgrib2_output(output) - - assert result["TMP:1000 mb"] == 297.3 - assert result["TMP:975 mb"] == 295.1 - assert result["HGT:1000 mb"] == 110.0 - end - - test "handles empty output" do - assert HrrrClient.parse_wgrib2_output("") == %{} - end - end - - describe "build_profile_from_wgrib2/1" do - test "assembles profile from parsed wgrib2 data" do + describe "build_profile/1" do + test "assembles profile from parsed data" do parsed = %{ "TMP:1000 mb" => 298.0, "DPT:1000 mb" => 291.0, @@ -173,7 +136,7 @@ defmodule Microwaveprop.Weather.HrrrClientTest do "PWAT:entire atmosphere (considered as a single layer)" => 25.0 } - result = HrrrClient.build_profile_from_wgrib2(parsed) + result = HrrrClient.build_profile(parsed) assert result.surface_temp_c == 299.0 - 273.15 assert result.surface_dewpoint_c == 292.0 - 273.15 @@ -204,7 +167,7 @@ defmodule Microwaveprop.Weather.HrrrClientTest do "PWAT:entire atmosphere (considered as a single layer)" => 25.0 } - result = HrrrClient.build_profile_from_wgrib2(parsed) + result = HrrrClient.build_profile(parsed) assert length(result.profile) == 1 end end diff --git a/test/microwaveprop/workers/qso_weather_enqueue_worker_test.exs b/test/microwaveprop/workers/qso_weather_enqueue_worker_test.exs index ccd11805..975ee4e9 100644 --- a/test/microwaveprop/workers/qso_weather_enqueue_worker_test.exs +++ b/test/microwaveprop/workers/qso_weather_enqueue_worker_test.exs @@ -4,6 +4,7 @@ defmodule Microwaveprop.Workers.QsoWeatherEnqueueWorkerTest do alias Microwaveprop.Radio.Qso alias Microwaveprop.Terrain.ElevationClient alias Microwaveprop.Weather + alias Microwaveprop.Weather.HrrrClient alias Microwaveprop.Weather.IemClient alias Microwaveprop.Workers.QsoWeatherEnqueueWorker @@ -129,6 +130,11 @@ defmodule Microwaveprop.Workers.QsoWeatherEnqueueWorkerTest do end end) + # Stub HRRR client — return 404 so fetch_profile returns an error + Req.Test.stub(HrrrClient, fn conn -> + Plug.Conn.send_resp(conn, 404, "not found") + end) + # Stub elevation API for terrain profile workers Req.Test.stub(ElevationClient, fn conn -> params = Plug.Conn.fetch_query_params(conn).query_params @@ -265,6 +271,10 @@ defmodule Microwaveprop.Workers.QsoWeatherEnqueueWorkerTest do end end) + Req.Test.stub(HrrrClient, fn conn -> + Plug.Conn.send_resp(conn, 404, "not found") + end) + Req.Test.stub(ElevationClient, fn conn -> params = Plug.Conn.fetch_query_params(conn).query_params lat_count = params["latitude"] |> String.split(",") |> length() diff --git a/test/microwaveprop_web/live/submit_live_test.exs b/test/microwaveprop_web/live/submit_live_test.exs index 47a8a68f..97868917 100644 --- a/test/microwaveprop_web/live/submit_live_test.exs +++ b/test/microwaveprop_web/live/submit_live_test.exs @@ -6,6 +6,7 @@ defmodule MicrowavepropWeb.SubmitLiveTest do alias Microwaveprop.Radio.Qso alias Microwaveprop.Repo alias Microwaveprop.Terrain.ElevationClient + alias Microwaveprop.Weather.HrrrClient alias Microwaveprop.Weather.IemClient setup do @@ -16,6 +17,10 @@ defmodule MicrowavepropWeb.SubmitLiveTest do end end) + Req.Test.stub(HrrrClient, fn conn -> + Plug.Conn.send_resp(conn, 404, "not found") + end) + Req.Test.stub(ElevationClient, fn conn -> params = Plug.Conn.fetch_query_params(conn).query_params lat_count = params["latitude"] |> String.split(",") |> length()