diff --git a/lib/microwaveprop/propagation.ex b/lib/microwaveprop/propagation.ex index dd2cbef7..f56964b2 100644 --- a/lib/microwaveprop/propagation.ex +++ b/lib/microwaveprop/propagation.ex @@ -86,7 +86,7 @@ defmodule Microwaveprop.Propagation do pressure_mb: hrrr_profile.surface_pressure_mb, prev_pressure_mb: nil, rain_rate_mmhr: Scorer.precip_to_rate_mmhr(hrrr_profile[:precip_mm]), - min_refractivity_gradient: derived[:min_refractivity_gradient], + min_refractivity_gradient: hrrr_profile[:native_min_gradient] || derived[:min_refractivity_gradient], bl_depth_m: hrrr_profile[:hpbl_m], pwat_mm: hrrr_profile[:pwat_mm] } diff --git a/lib/microwaveprop/weather/grib2/wgrib2.ex b/lib/microwaveprop/weather/grib2/wgrib2.ex index 8fe3a00a..e1d8af4d 100644 --- a/lib/microwaveprop/weather/grib2/wgrib2.ex +++ b/lib/microwaveprop/weather/grib2/wgrib2.ex @@ -40,6 +40,26 @@ defmodule Microwaveprop.Weather.Grib2.Wgrib2 do end end + @doc """ + Like `extract_grid_from_file/3`, but instead of building the full + `%{{lat, lon} => %{"VAR:LEVEL" => float}}` map (which can be 1+ GB for + 200+ messages × 95k cells), processes each grid cell through a reducer + function and keeps only the reduced output. + + The reducer receives `%{"VAR:LEVEL" => float}` for one cell and returns + an arbitrary term. The result is `%{{lat, lon} => reducer_output}`. + + Peak memory: binary output (~76 MB for 200 msgs × 95k cells) + one + cell's worth of data at a time + output map (~10 MB of scalars). + """ + def extract_grid_from_file_mapped(grib_path, match_pattern, grid_spec, cell_reducer) do + if available?() do + extract_file_mapped_with_wgrib2(grib_path, match_pattern, grid_spec, cell_reducer) + else + {:error, :wgrib2_not_available} + end + end + @doc "Check if wgrib2 is available on the system." def available?, do: wgrib2_path() != nil @@ -82,6 +102,104 @@ defmodule Microwaveprop.Weather.Grib2.Wgrib2 do end end + defp extract_file_mapped_with_wgrib2(grib_path, match_pattern, grid_spec, cell_reducer) do + %{ + lon_start: lon_start, + lon_count: lon_count, + lon_step: lon_step, + lat_start: lat_start, + lat_count: lat_count, + lat_step: lat_step + } = grid_spec + + wgrib2_lon_start = normalize_lon(lon_start) + lon_spec = "#{wgrib2_lon_start}:#{lon_count}:#{lon_step}" + lat_spec = "#{lat_start}:#{lat_count}:#{lat_step}" + + tmp_bin = grib_path <> ".lola.mapped.bin" + + try do + args = [grib_path, "-match", match_pattern, "-lola", lon_spec, lat_spec, tmp_bin, "bin"] + + case System.cmd(wgrib2_path(), args, stderr_to_stdout: true) do + {output, 0} -> + messages = parse_wgrib2_inventory(output) + + case File.read(tmp_bin) do + {:ok, bin_data} -> + parse_lola_binary_mapped(bin_data, messages, grid_spec, cell_reducer) + + {:error, :enoent} -> + {:ok, %{}} + + {:error, reason} -> + {:error, "Failed to read wgrib2 output: #{inspect(reason)}"} + end + + {output, exit_code} -> + {:error, "wgrib2 failed (exit #{exit_code}): #{String.slice(output, 0, 200)}"} + end + after + File.rm(tmp_bin) + end + end + + # Cell-by-cell variant of parse_lola_binary. Instead of building the full + # %{{lat,lon} => %{"VAR:LEVEL" => val}} map, iterates grid cells and for + # each cell extracts its values from every message, passes the per-cell + # map to the reducer, and keeps only the reduced output. + defp parse_lola_binary_mapped(bin_data, messages, grid_spec, cell_reducer) do + %{lon_count: nx, lat_count: ny, lon_start: lon_start, lon_step: lon_step, lat_start: lat_start, lat_step: lat_step} = + grid_spec + + points_per_message = nx * ny + bytes_per_message = points_per_message * 4 + record_overhead = 8 + stride = bytes_per_message + record_overhead + + # Pre-compute message metadata: {key, data_offset} for each message + msg_meta = + messages + |> Enum.with_index() + |> Enum.map(fn {msg, idx} -> + {"#{msg.var}:#{msg.level}", idx * stride + 4} + end) + |> Enum.filter(fn {_key, offset} -> offset + bytes_per_message <= byte_size(bin_data) end) + + # Iterate each grid cell, extract all message values, reduce + result = + Enum.reduce(0..(ny - 1), %{}, fn j, acc_outer -> + lat = Float.round(lat_start + j * lat_step, 3) + + Enum.reduce(0..(nx - 1), acc_outer, fn i, acc_inner -> + cell_offset = (j * nx + i) * 4 + + # Extract all values for this cell across all messages + cell_data = + Enum.reduce(msg_meta, %{}, fn {key, data_offset}, cell_acc -> + offset = data_offset + cell_offset + <<_::binary-size(offset), value::float-little-32, _::binary>> = bin_data + + if value > @undefined_value / 2 do + cell_acc + else + Map.put(cell_acc, key, value) + end + end) + + if map_size(cell_data) > 0 do + lon = Float.round(denormalize_lon(lon_start + i * lon_step), 3) + reduced = cell_reducer.(cell_data) + Map.put(acc_inner, {lat, lon}, reduced) + else + acc_inner + end + end) + end) + + {:ok, result} + end + defp extract_with_wgrib2(grib_binary, match_pattern, grid_spec) do %{ lon_start: lon_start, diff --git a/lib/microwaveprop/weather/hrrr_native_client.ex b/lib/microwaveprop/weather/hrrr_native_client.ex index 8be9682c..377ed6d8 100644 --- a/lib/microwaveprop/weather/hrrr_native_client.ex +++ b/lib/microwaveprop/weather/hrrr_native_client.ex @@ -125,6 +125,116 @@ defmodule Microwaveprop.Weather.HrrrNativeClient do HrrrClient.byte_ranges_for_messages(idx_entries, native_messages()) end + # The 4 variables needed for duct detection (N-profile + height). + # Skips UGRD, VGRD, TKE (shear/Richardson were dead features). + @duct_variables ~w(TMP SPFH HGT PRES) + + @doc """ + Messages needed for duct detection only (4 vars × 50 levels = 200 messages). + ~300 MB download vs ~530 MB for all 7 variables. + """ + def duct_messages do + for level <- @native_levels, var <- @duct_variables do + %{var: var, level: "#{level} hybrid level"} + end + end + + @doc """ + Byte ranges for duct-detection variables only. + """ + def duct_byte_ranges(idx_entries) do + HrrrClient.byte_ranges_for_messages(idx_entries, duct_messages()) + end + + @doc """ + Fetch native duct metrics for the CONUS grid for one forecast hour. + + Downloads TMP, SPFH, HGT, PRES on all 50 hybrid levels, extracts to + the CONUS grid via wgrib2, and computes duct metrics per cell using a + cell-by-cell reducer (peak memory ~86 MB instead of ~1.8 GB). + + Returns `{:ok, %{{lat, lon} => duct_metrics}}` or `{:error, reason}`. + """ + def fetch_native_duct_grid(date, hour, grid_spec, forecast_hour \\ 0) do + alias Microwaveprop.Weather.Grib2.Wgrib2 + + url = hrrr_native_url(date, hour, forecast_hour) + idx_url = url <> ".idx" + + tmp_grib = Path.join(System.tmp_dir!(), "hrrr_native_duct_#{System.unique_integer([:positive])}.grib2") + + try do + with {:ok, idx_text} <- fetch_idx(idx_url), + idx_entries = HrrrClient.parse_idx(idx_text), + ranges = duct_byte_ranges(idx_entries), + :ok <- HrrrClient.download_grib_ranges_to_file(url, ranges, tmp_grib) do + match_pattern = ":(#{Enum.join(@duct_variables, "|")}):.*hybrid level:" + + Wgrib2.extract_grid_from_file_mapped(tmp_grib, match_pattern, grid_spec, &compute_duct_metrics/1) + end + after + File.rm(tmp_grib) + end + end + + defp fetch_idx(url) do + case Req.get(url, receive_timeout: 120_000) do + {:ok, %{status: 200, body: body}} -> {:ok, body} + {:ok, %{status: status}} -> {:error, "HRRR native idx HTTP #{status}"} + {:error, reason} -> {:error, reason} + end + end + + # Compute duct metrics from a single cell's native profile data. + # Input: %{"TMP:1 hybrid level" => 297.5, "SPFH:1 hybrid level" => 0.012, ...} + # Output: %{native_min_gradient: float, best_duct_freq_ghz: float, max_duct_thickness_m: float} + defp compute_duct_metrics(parsed) do + alias Microwaveprop.Propagation.Duct + + profile = build_native_profile(parsed) + + if profile.level_count >= 3 do + duct_result = Duct.analyze(profile) + + %{ + native_min_gradient: min_m_gradient(profile), + best_duct_freq_ghz: duct_result.best_duct_band_ghz, + max_duct_thickness_m: max_duct_thickness(duct_result.ducts), + duct_count: length(duct_result.ducts) + } + else + %{native_min_gradient: nil, best_duct_freq_ghz: nil, max_duct_thickness_m: nil, duct_count: 0} + end + end + + defp min_m_gradient(%{heights_m: heights, temp_k: temps, spfh: spfhs, pressure_pa: pressures}) do + # Compute modified refractivity M at each level, find minimum dM/dh + ms = + Enum.zip([heights, temps, spfhs, pressures]) + |> Enum.map(fn {h, t, q, p} -> + # N = 77.6 * P/T + 3.73e5 * e/T^2, where e = q*P/(0.622 + 0.378*q) + q = max(q || 0.0, 1.0e-8) + e = q * p / (0.622 + 0.378 * q) / 100.0 + t_safe = max(t, 1.0) + n = 77.6 * (p / 100.0) / t_safe + 3.73e5 * e / (t_safe * t_safe) + m = n + 157.0 * h / 1000.0 + {h, m} + end) + + ms + |> Enum.chunk_every(2, 1, :discard) + |> Enum.map(fn [{h1, m1}, {h2, m2}] -> + dh = h2 - h1 + if dh > 0, do: (m2 - m1) / dh * 1000.0, else: 0.0 + end) + |> Enum.min(fn -> 0.0 end) + end + + defp min_m_gradient(_), do: nil + + defp max_duct_thickness([]), do: nil + defp max_duct_thickness(ducts), do: ducts |> Enum.map(& &1.thickness_m) |> Enum.max() + @doc """ Extract native profiles for a list of `{lat, lon}` points from a GRIB2 file on disk, using wgrib2. Avoids loading the entire file diff --git a/lib/microwaveprop/workers/propagation_grid_worker.ex b/lib/microwaveprop/workers/propagation_grid_worker.ex index 8e3b5cec..06f20663 100644 --- a/lib/microwaveprop/workers/propagation_grid_worker.ex +++ b/lib/microwaveprop/workers/propagation_grid_worker.ex @@ -14,6 +14,7 @@ defmodule Microwaveprop.Workers.PropagationGridWorker do alias Microwaveprop.Propagation.Grid alias Microwaveprop.Weather alias Microwaveprop.Weather.HrrrClient + alias Microwaveprop.Weather.HrrrNativeClient alias Microwaveprop.Weather.SoundingParams require Logger @@ -78,6 +79,10 @@ defmodule Microwaveprop.Workers.PropagationGridWorker do store_hrrr_profiles(grid_data, valid_time) :erlang.garbage_collect() + # Fetch native duct metrics and merge into grid_data for scoring + grid_data = merge_native_duct_data(grid_data, run_time, forecast_hour) + :erlang.garbage_collect() + Phoenix.PubSub.broadcast( Microwaveprop.PubSub, "weather:updated", @@ -163,6 +168,31 @@ defmodule Microwaveprop.Workers.PropagationGridWorker do Logger.info("PropagationGrid: stored #{count} HRRR profiles") end + defp merge_native_duct_data(grid_data, run_time, forecast_hour) do + hour_dt = HrrrClient.nearest_hrrr_hour(run_time) + date = DateTime.to_date(hour_dt) + hour = hour_dt.hour + grid_spec = Grid.wgrib2_grid_spec() + + case timed("native", fn -> + HrrrNativeClient.fetch_native_duct_grid(date, hour, grid_spec, forecast_hour) + end) do + {:ok, duct_grid} -> + Logger.info("PropagationGrid: merged #{map_size(duct_grid)} native duct cells") + + Map.new(grid_data, fn {point, profile} -> + case Map.get(duct_grid, point) do + nil -> {point, profile} + duct -> {point, Map.merge(profile, duct)} + end + end) + + {:error, reason} -> + Logger.warning("PropagationGrid: native duct fetch failed (continuing without): #{inspect(reason)}") + grid_data + end + end + defp compute_scores(grid_data, valid_time) do # Algorithm is the primary scorer. ML score stored in factors for comparison. compute_scores_algorithm(grid_data, valid_time)